具体实施方式

为了使本发明的技术目的、技术方案和有益效果更加清楚，下面结合附图和具体实施例 对本发明的技术方案做出进一步的说明。

本发明提供了一种斜浪下船艏非对称砰击载荷的确定性方法，如图1所示，该方法包括 以下步骤：

步骤1):将船体假设为刚体模型，此时可忽略航行中波浪载荷对船体结构变形的影响。船 体重心处运动模型符合牛顿第二定律，即

式中，M为广义质量矩阵；F(t)为船体周围流体载荷；η(t)是刚体重心的运动。

步骤2):简化计算模型波浪载荷的求解。如图2所示，船体外围流域D包括四个边界：底 边界S_B，自由表面边界S_F，物面边界S_S和截断边界S_R。假定流场无旋、无粘不可压，流体 域满足Laplace方程，引入速度势，基本控制方程和边界条件表示为

为了得到速度势φ_j，引入格林函数G(P,Q)来建立边界积分方程，在物体表面分布的未知 源强度σ(q)满足

式中，P是流体域中的任意点；S为平均湿润面；Q是位于S的点。代入物面边界条件，可得 源强度的积分方程，为

采用面元法求解上述方程，其基本思想是利用三角/四边形面元离散湿表面，然后将方程 转化为未知源强度的线性方程组。

将作用在船体结构上的流体载荷F(t)可以分为两类：静水载荷F^S(t)和水动力载荷F^D(t)。 水动力载荷可以进一步分解为入射波力F_I(t)、绕射力F_D(t)和辐射力F_R(t)。因此，总流体荷载 可分解为

{F(t)}＝{F^S(t)}+{F_I(t)}+{F_D(t)}+{F_R(t)} (5)

对于静水载荷F^S(t)，假设船舶运动为线性，则简化为

{F^S(t)}＝-[C]{η}e^iωt (6) 式中，C为静水系数，可通过对沿船体平均湿表面的静水压力积分得到。入射波力F_I(t)和绕 射力F_D(t)共同构成波干涉力F_ID(t)，其与入射波的形式有关，记为

式中，ζ_a为波幅，和为单位振幅下的入射波势和绕射势。通过对湿表面的辐射压力进行 积分，得到辐射力，即

在船舶运动方程中，辐射力通常以三维水动力系数的形式表示，记为

式中，A为附加质量，B为阻尼系数。

步骤3):建立波浪下船体重心处的六自由度运动方程。将作用在船体上的各种流体载荷分 量(式((6)、(7)和(8))代入到船舶运动方程式(1)中，消去时间项，得到准静态方程系统为

([C]-ω²([M]+[A])+iω[B]){η}＝{f} (10)

上述运动方程为线性方程组，可采用高斯消去法进行求解，从而得到六自由度运动的实 部和虚部，此时船舶运动的幅值和相位可以表示为

步骤4):建立船波三维空间的相对运动预报模型。船艏的运动形式与重心处有所区别，假 定船体为刚体时，可根据刚体运动理论得到船体三维空间的运动形式，考虑到波浪的影响， 通常以相对运动形式求解三维空间的运动。假设船体遇到一个单位谐波，对于船体的任意点 (x,y,z)，则水平和垂直相对位移可表示为

式中，(x_b,y_b,z_b)为船体重心坐标；H_y、H_z为船舶与波浪的相对位移；β为波浪方向的角度；k 是波数。船波相对水平速度y_v和垂直速度z_v的线性分量表示为

式中，u_w和w_w分别为波浪的垂向和水平速度。

需要注意的是，横摇运动作为入水运动的单独量，由其引起的位移变化单独计及，在刚 体模型下，结构的转动在三维空间中与结构重心处的运动形式保持一致，即

w＝η₄ (16)

采用步骤1-4，即可得到船体特定剖面的横向运动、垂向运动和橫摇运动的时历，如图3所 示。

步骤5):建立非对称入水砰击模型。对船体进行切片，将三维砰击问题转化为二维砰击， 简化计算量。船体剖面以横向速度y_v(式(14))，垂向速度z_v(式(15))和横摇速度w入水，如图 4所示。整个流体域边界有上部S_T，底部S_B，远方S_C和物面S_S。自由面上部为空气，下部为 水。流域的上侧S_T设置为压力出口边界，底部S_B和两侧S_C为无滑移边界，物体表面S_S为无 滑移刚性边界。

步骤6):应用VOF方法求解入水砰击问题，预报砰击载荷。基于N-S方程建立入水冲击 问题的两相流模型。假设流体不可压，并忽略流体的粘性和重力效应，则连续性方程和动量 方程简化为

式中，p为压力，u为速度矢量，ρ为流体的密度。应用VOF方法求解上述方程，假设水为 基本相，引入相方程

式中，α为基本相的体积分数。则流体的密度可进一步表示为

ρ＝αρ₁+(1-α)ρ₂ (20)

式中，ρ₁和ρ₂分别为水和空气的密度，通常取为1.0×10³kg/m³和1.25kg/m³。方程(17)、(18) 和(19)构成了应用VOF法求解两相流模型的连续性方程、动量方程和相方程。入水速度时历 由步骤(4)给出的船波相对运动关系(式(14)-(16))确定。在数值模拟中，基于有限体积方法对动 量方程进行离散求解，时间项采用欧拉隐格式，式(18)各项离散为

式中，上标n和r分别表示当前时间步(值已知)和预测时间步(值未知)，下标f表示网格表面 的值，V_p为网格的体积，Δt为时间步长，F_f为通量，S表示网格单元的面矢量。网格表面的 速度可依据相邻体单元中心的速度获得，一般可采用中心差分方法获得，即 为表面速度的法向梯度，记为

式中，下标P和N分别表示当前网格和相邻网格的速度，|S_f|为面矢量的模，|d|为网格间的距 离。需要注意的是，在方程(18)对流项的离散中，其中通量F_f采用当前已知时间步的速度来 计算，并保留另一个预测速度为未知量。这种把其中一个速度已知化并保留另外一个未知速 度的过程即为线性化。将式(21)-(24)代入式(18)，并引入压力梯度和源项，有

其中

在求得预测速度之后，再次去掉公式(18)中的压力项和源项并进行移项，使用符号 HbyA代替预测时间步的速度，并按下式定义

可以看出，在某个时间步内A_P，A_N和E_P均保持不变，唯一的未知量可通过方程(29) 获得。值得注意的是此时并未用到连续性方程，这意味着预测的并不一定满足质量守恒， 可采用常用的PISO算法进行修正。

相方程(20)中时间项采用隐格式进行离散，体积分数采用标准有限差分方式离散，即

式中，上标n和n+1分别表示前一时间步和后一时间步，为表面通量。时间步进求解过程 中，当网格内为水或空气(α_a＝0或1)时，采用标准的中心差分模式获得面通量，当网格为空气 和水的交界面(0<α_a<1)，通过几何重构模式描述自由面。时间步进中，监测各个节点的压力， 即可得到各个部位的压力时历结果(图5)和整个剖面的压力分布(图6)。

最后所应说明的是：上述实施例仅用于说明而非限制本发明的技术方案，任何对本发明 进行的等同替换及不脱离本发明精神和范围的修改或局部替换，其均应涵盖在本发明权利要 求保护的范围之内。