一种任意二维剖面三自由度入水载荷确定的边界元方法与流程

1.本发明涉及流体动力学领域，尤其是一种任意二维剖面三自由度入水载荷确定的边界元方法。
技术背景
2.结构物多自由度入水具有广泛的工程背景，如海上救生艇投放、水上飞机降落及空投鱼雷入水等。结构入水过程中，会激起周围流体介质的运动；反过来，流体介质对结构又施加各种反作用力，是典型的流体与结构相互作用的动力学问题。
3.结构物多自由度入水通常从物体的底部与自由界面发生接触开始，随后能量迅速在流场中传递。随着时间的推移，运动体附近流场趋于稳定，在这个过程中往往伴随着毫秒量级的高峰值冲击载荷。一方面，巨大的冲击载荷会直接造成结构的变形，甚至引起结构损坏。另一方面，冲击载荷也会影响结构的运动，产生的加速度会改变物体的航行姿态。由此可见，无论是从结构安全性还是航行稳定性来看，结构物多自由度入水问题已经成为航海、航空等多个领域需要重点关注的问题，对其进行深入研究具有广阔的应用前景和重要的工程价值。
4.结构物多自由度入水是一个具有强瞬时性的复杂多相流动过程，其自由面演变和压力的时空分布特性与入水速度、几何形状、结构弹性、空气层、流体的压缩性等诸多因素有关。结构受冲击的同时也会激起流体的运动，撞击过程中结构的运动与变形又会影响自由面的演变，整个过程一直伴随着流体与结构的耦合作用。结构物多自由度入水中自由面的非线性效应占主导位置，包括细薄的射流、流体分离破碎、气泡等自由面特征，这些现象增加了问题的复杂性，使入水过程中冲击载荷的精确预测变得更加困难。
5.鉴于结构物多自由度入水问题的复杂性，载荷的预报极具挑战，本发明提出一种任意二维剖面三自由度入水载荷确定的边界元方法。本方法计及了垂向、横向和橫摇三种运动的影响，能够精确地预报斜向入水情况下入水载荷。采用边界元方法分析入水问题时，降低了计算的维度，与其他数值方法相比，有效的降低了入水问题的求解难度，显著提高了计算效率，为复杂几何结构多自由度入水载荷问题提供保障。


技术实现要素：

6.本发明目的是提供一种任意二维剖面三自由度入水载荷确定的边界元方法，弥补目前载荷预报中无法充分考虑非对称运动形式的缺陷，提供一种复杂几何结构多自由度入水载荷确定的高精度方法。
7.为实现发明目的，本发明采用如下技术方案：
8.步骤1：建立三自由度入水砰击的数学模型；
9.任意二维剖面以三自由度入水时包含横向速度u，垂向速度v，转动速度w。整个流体域d的基本边界有物面s
b
，自由面s
f
，远方边界s
c
和底部边界s
b
。
10.步骤2：建立计算模型的控制方程和边界条件；
11.假定流体无旋、无粘且不可压缩，则整个流体域d满足laplace方程
[0012][0013]
自由面s
f
上满足运动学方程和动力学方程
[0014][0015][0016]
在物面表面s
b
有
[0017][0018]
式中，(n
x
,n
y
)为物面单位法向量，指向流域外部，(x,y)为所示坐标系下物面的坐标。
[0019]
在远方控制面s
c
和底部s
b
上，流域未发生扰动，有
[0020][0021]
步骤3：建立入水初始时刻的理论解；
[0022]
初始入水阶段，自由面变化剧烈，直接利用边界元方法进行求解时，会导致计算量的急剧增大。通常的解决办法是引入一个近似的方法来模拟入水初始阶段的自由表面。对于二维入水砰击问题，引入的外域解来给出任意时刻自由面满足的方程，即
[0023][0024]
求解上述常微分方程，则自由表面给出以下形式
[0025][0026]
式中，c(t)为半浸湿宽度，当从t0时刻开始采用边界元方法时，此时自由面形状由上式给出。
[0027]
步骤4：边界离散建立未知量求解方程组；
[0028]
根据势流理论的推导结果，流域任意点p处的速度势可用边界上的速度势及其法向导数表示，即
[0029][0030]
式中，n
q
为q点处边界法向指向流域外部，r
pq
为p，q两点的距离。将整个流域边界离散为若干单元，边界条件转化为单元本身的边界条件，并采用常值单元，将上式离散为
[0031]
[0032]
式中，n为划分的单元总个数，s
j
为第j个单元边界，c
ij
和d
ij
为影响系数。对于每个单元的影响系数在该单元的局部坐标系下计算，局部坐标系原点在单元中心， x轴在单元上，y轴垂直于单元且指向流域内部，则
[0033][0034][0035][0036]
将上式写成矩阵形式，有
[0037][0038]
式中，c，d分别为元素c
ij
和d
ij
组成的矩阵，e为单位阵，维数均为n
×
n。根据边界条件分离已知量和未知量，上式改写为
[0039][0040]
采用gauss消去法求解上述方程组可得到物面上的速度势和自由面上的法向速度，之后对时间步进求解。
[0041]
步骤5：自由面更新、光顺与重构；
[0042]
对于自由面的更新，假设t时刻自由面速度势和速度均已知(初始时刻均为 0)，给定步长dt，则下个时刻的坐标及速度势可按下式进行更新
[0043][0044]
式中，自由面x向速度和y向速度可根据速度势的切向导数φ
s
和法向导数φ
n
求出，即
[0045]
[0046][0047]
式中，φ
s
为切向导数。本发明应用五阶三次平滑公式对自由面进行光顺。假设f
i
为更新后的y向坐标，则第一个节点处，f
i
更新为
[0048][0049]
第二个节点处，f
i
更新为
[0050][0051]
第三个至第n
‑
2个节点处，f
i
更新为
[0052][0053]
第n
‑
1个节点处，f
i
更新为
[0054][0055]
第n个节点处，f
i
更新为
[0056][0057]
步骤6：入水载荷的求解；
[0058]
对于速度势对时间的导数由时间差分得到，当第一个时间步时，记为
[0059][0060]
当第n个时间步时，记为
[0061][0062]
将所求压力沿物体表面积分，可得到剖面所受力和力矩，即
[0063]
f
x
＝∫
s
pn
x
ds
[0064]
f
y
＝∫
s
pn
y
ds
[0065]
m＝∫
s
p(r
×
n)ds
[0066]
式中，f
x
为垂向力，f
y
为水平力，m为所受力矩。
[0067]
有益效果
[0068]
目前入水冲击载荷的计算大多集中于楔形体的垂向入水问题，结构形式简单，运动形式单一。实际的结构多为不规则剖面形式且运动形式更为多样化。本方法将入水冲击载荷的预报进一步扩展到任意剖面，计及了多种运动形式，这在复杂工程问题中更能准确地预报砰击载荷。本发明建立在边界元理论上，与其他数值方法相比，有效地降低了计算的维度，显著提高了计算效率，为复杂几何结构多自由度入水载荷问题提供保障。
附图说明
[0069]
图1是一种任意二维剖面三自由度入水载荷确定的边界元方法流程图
[0070]
图2是任意剖面三自由度入水模型示意图
[0071]
图3是边界元系数求解的局部坐标系
[0072]
图4是剖面入水过程的自由面
[0073]
图5是剖面入水的压力时历
[0074]
图6是剖面入水的压力分布
具体实施方式
[0075]
为了使本发明的技术目的、技术方案和有益效果更加清楚，下面结合附图和具体实施例对本发明的技术方案做出进一步的说明。
[0076]
本发明提供了一种任意二维剖面三自由度入水载荷确定的边界元方法，计算流程如图1所示。以二维船体外飘剖面多自由度入水为具体实施例，该方法包括以下步骤：
[0077]
步骤1：建立三自由度入水砰击的数学模型；
[0078]
如图2所示，二维船体外飘剖面以三自由度入水时包含横向速度u，垂向速度v，转动速度w。整个流体域d的基本边界有物面s
s
，自由面s
f
，远方边界 s
c
和底部边界s
b
。
[0079]
步骤2：建立计算模型的控制方程和边界条件；
[0080]
假定流体无旋、无粘且不可压缩，则整个流体域d满足laplace方程
[0081][0082]
自由面s
f
上满足运动学方程和动力学方程
[0083][0084][0085]
在物面表面s
s
有
[0086][0087]
式中，(n
x
,n
y
)为物面单位法向量，指向流域外部，(x,y)为所示坐标系下物面的坐标。
[0088]
在远方控制面s
c
和底部s
b
上，流域未发生扰动，有
[0089][0090]
步骤3：建立入水初始时刻的理论解；
[0091]
初始入水阶段，自由面变化剧烈，直接利用边界元方法进行求解时，会导致计算量的急剧增大。通常的解决办法是引入一个近似的方法来模拟入水初始阶段的自由表面。对于二维入水砰击问题，引入外域解来给出任意时刻自由面满足的方程，即
[0092]
[0093]
求解上述常微分方程，得到自由表面表达式如下
[0094][0095]
式中，c(t)为半浸湿宽度，当从t0时刻开始采用边界元方法时，初始自由面形状由上式给出。
[0096]
步骤4：边界离散建立未知量求解方程组；
[0097]
根据势流理论的推导结果，流域任意点p处的速度势可用边界上的速度势及其法向导数表示，即
[0098][0099]
式中，n
q
为q点处边界法向指向流域外部；r
pq
为p，q两点的距离。将整个流域边界离散为若干单元，边界条件转化为单元本身的边界条件，采用常值单元将上式离散为
[0100][0101]
式中，n为划分的单元总个数；s
j
为第j个单元边界；c
ij
和d
ij
为影响系数。
[0102]
对于每个单元的影响系数在该单元的局部坐标系下计算，局部坐标系原点在单元中心，x轴在单元上，y轴垂直于单元且指向流域内部，如图3所示，则
[0103][0104][0105][0106]
将式(9)写成矩阵形式，有
[0107][0108]
式中，c，d分别为元素c
ij
和d
ij
组成的矩阵；e为单位阵，维数均为n
×
n。根据边界条件分离已知量和未知量，上式改写为
[0109][0110]
采用gauss消去法求解上述方程组可得到物面上的速度势和自由面上的法向速度，之后对时间步进求解。
[0111]
步骤5：自由面更新、光顺与重构；
[0112]
对于自由面的更新，假设t时刻自由面速度势和速度均已知(初始时刻均为 0)，给定步长dt，则下个时刻的坐标及速度势可按下式进行更新
[0113][0114]
式中，自由面x向速度和y向速度可根据速度势的切向导数φ
s
和法向导数φ
n
求出，即
[0115][0116]
式中，φ
s
为切向导数，分为三种情况进行讨论。
[0117]
当单元为第一个单元时，φ
s
记为
[0118][0119]
当单元为第2个至第n
‑
1个单元时，φ
s
记为
[0120][0121]
当单元为第n个单元时，φ
s
记为
[0122][0123]
将式(15)
‑
(18)代入式(14)，即可得到t dt时刻的自由面，此时得到的自由面并非光滑，需要进行光顺，避免后续时间步进中发散，本发明应用五阶三次平滑公式对自由面进行光顺。假设f
i
为更新后的y向坐标，则第一个节点处，f
i
更新为
[0124][0125]
第二个节点处，f
i
更新为
[0126]
[0127]
第三个至第n
‑
2个节点处，f
i
更新为
[0128][0129]
第n
‑
1个节点处，f
i
更新为
[0130][0131]
第n个节点处，f
i
更新为
[0132][0133]
式中，φ
s
为切向导数；h为单元长度；y
i
为节点的y向坐标。图4给出了两个时刻下自由面的形式。
[0134]
步骤6：入水载荷的求解；
[0135]
由bernoulli方程，物体表面的压力可写为
[0136][0137]
对于速度势对时间的导数由时间差分得到，当第一个时间步时，记为
[0138][0139]
当第n个时间步时，记为
[0140][0141]
图5给出了两侧对称位置的压力时历，图6给出了不同时刻下剖面所受压力分布。将所求压力沿物体表面积分，可得到剖面所受力和力矩，即
[0142][0143]
式中，f
x
为垂向力，f
y
为水平力，m为所受力矩。
[0144]
最后所应说明的是：上述实施例仅用于说明而非限制本发明的技术方案，任何对本发明进行的等同替换及不脱离本发明精神和范围的修改或局部替换，其均应涵盖在本发明权利要求保护的范围之内。