基于SPH的充液刚体刚-液耦合动力学建模

张诗琪，刘锦阳; Zhang Shiqi; Liu Jinyang

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

摘要

本文考虑刚-液耦合效应，基于SPH方法建立充液刚体动力学模型，该动力学模型不仅适用于微幅晃动问题，还适用于大幅晃动问题.首先建立了非惯性系下液体晃动的方程，其中,边界条件采用镜像虚粒子法.通过计算各流体粒子的绝对加速度，得到粒子惯性力对刚体质心的主矢和主矩，用达朗贝尔原理计算得到向质心简化的晃动力和晃动力偶矩，并引入刚体的约束方程，建立了考虑刚-液耦合的充液刚体的动力学方程.首先对给定周期性运动的充液刚体的液面振荡进行数值仿真，验证了建模方法的有效性;然后通过对受周期性激励力作用的充液刚体的运动和液面振荡进行数值仿真，分析了不同参数下系统的刚-液耦合特性.

关键词

充液刚体; 刚-液耦合; 动力学建模; 大幅度晃动

引言

对于液体微幅晃动研究，一种常用的方法是基于线性速度势理论和微幅晃动假设，计算液体晃动的解析解^[

1]；另一种方法是采用等效摆模型得到了微幅晃动下的液体等效方程，并计算出液体的晃动稳定性边界^{[参考文献 2
百度学术}2].随着数值计算方法和技术的进步，不少学者利用数值模拟方法研究充液系统的大幅晃动问题.常用的数值仿真方法是基于网格思想的有限元法（FEM），有限体积法（FVM）和边界元法（BEM）等.Goudarzi^{[参考文献 3
百度学术}3]等用非线性有限元方法对周期性载荷与地震激励同时作用下液舱的晃动问题进行研究，通过实验对比，验证了计算结果的准确性.石怀龙等^{[参考文献 4
百度学术}4]基于绝对节点坐标方法描述的实体有限单元进行液体网格划分，研究了液罐车中的液体晃动问题，模拟了液罐车横向晃动和纵向晃动时液体自由表面的形状变化.

对于大幅晃动问题，用网格方法进行计算，有时需要重新划分边界网格，给数值计算带来困难.通过采用无网格方法，可避免上述基于网格的传统数值方法的缺点^[

5].SPH方法是一种无网格方法，即用粒子代替网格，粒子带有液体的物理属性，从而解决了网格求解存在的问题.SPH粒子的参数计算只受与其邻近的粒子影响，且只取决于SPH粒子当下的分布情况，因此可避免网格畸变的问题.SPH方法的这些特点使其可用来研究液体大晃动问题，甚至飞溅现象.

最初，SPH方法主要用于线性问题的研究，例如，Monaghan^[

5-7]采用SPH方法模拟平板剪切流、Poiseuille流等问题.随着流动问题研究的逐渐成熟，SPH方法开始被应用于液体晃动领域.郑兴^{[参考文献 8
百度学术}8]对外部液舱运动已知的液体晃动问题进行了研究，当液体做微幅晃动时，液体自由表面波高存在理论解，将理论解与SPH结果进行对比，证明了SPH研究液体晃动问题的可行性.刘富等^{[参考文献 9
百度学术}9]研究了基于SPH方法的三维液体晃动数值模拟，成功模拟出波浪翻卷和破碎现象.Rafiee等^{[参考文献 10
百度学术}10]用实验方法和SPH方法研究了在外部激励下矩形液舱晃动问题.李大鸣等^{[参考文献 11
百度学术}11]用SPH方法研究了二维水槽液体晃动问题，分析了充液率、防荡隔板对液体晃动的影响，并且与理论结果进行了对比.Negrut等^{[参考文献 12
百度学术}12]和Hu等^{[参考文献 13
百度学术}13，14]研究了复杂形状的液舱中的液体晃动问题.上述模拟均为已知液舱的运动规律时液体晃动问题的研究.

对于刚-液耦合问题，樊伟^[

15]等通过速度势函数，用伽辽金法离散，对含两个车轮的液罐车的刹车问题进行了研究，模拟了充液多刚体系统的耦合动力学特性，上述刚-液耦合模型仅适用于液体微幅晃动时的刚-液耦合问题研究.岳宝增等^{[参考文献 16-19}16-19]利用任意拉格朗日-欧拉( ALE) 有限元方法，对航天器刚体的运动与燃料晃动的耦合问题进行了研究，借助于网格移动算法实现数值仿真，不仅适用于小晃动问题，也适用于大晃动问题.

考虑到SPH方法在大幅晃动计算方面具有优势，本文基于SPH方法建立大幅晃动情况下刚-液耦合充液刚体的动力学模型.首先建立了在非惯性系下液体晃动的连续性方程和动量方程，固壁边界条件采用镜像虚粒子法.计算出各流体粒子的绝对加速度，得到流体粒子的惯性力对刚体质心作用的主矢、主矩，通过达朗贝尔原理，给出了液体对刚体向质心简化的晃动力、晃动力偶矩.以此为基础，建立了充液的单个刚体动力学方程.在此基础上，通过液体微幅晃动的算例，将本文方法的计算结果与理论解^[

1]和文献[3]的实验结果进行对比，验证了SPH方法研究液体小幅晃动和大幅晃荡的准确性，然后通过液舱受周期性外力时的刚-液耦合问题，探讨了不同参数情况下液体晃动对刚体运动的影响.

1 SPH方法基本理论

通过核函数近似和粒子近似这两个步骤，可以推导出物理量及其导数的SPH形式.物理量 $f (r)$ 经过核函数近似^[

5]，表达式为

f (r) \approx \int_{Ω} f (x) W (r - x, h) d x

(1)

其中，W为光滑核函数，h为光滑长度.本文采用高斯核函数作为光滑函数，表达式为：

W (R, h) = α_{d} e^{- R^{2}}

(2)

其中， $R = \frac{|r - x|}{h}$ ， $α_{d}$ 为常系数，在二维问题中取 $\frac{1}{π h^{2}}$ .

对核函数近似的表达式进一步做粒子近似.问题域经过SPH方法被离散为一系列粒子，粒子之间没有相互连接.如图1所示，与粒子a距离小于kh的粒子称为粒子a的相邻粒子，粒子a在周围相邻粒子的作用下自由运动，其中，k为与光滑函数有关的参数.以粒子a所在位置为圆心，以kh为半径的圆形区域称为粒子的支持域.

图1 粒子a支持域

Fig.1 The support domain of particle a

对于粒子a，用粒子的函数值叠加求和代替场函数积分，经过粒子近似后的表达式为：

<f (r_{a})> = \sum_{b = 1}^{N} \frac{m_{b}}{ρ_{b}} f (r_{b}) W_{a b}

(3)

其中，粒子b为粒子a支持域内的任意粒子，N为支持域内粒子总数， $m_{b}$ , $ρ_{b}$ 和 $r_{b}$ 分别为粒子 $b$ 的质量、密度和位置矢量， $r_{a}$ 为粒子 $a$ 位置矢量， $W_{a b} = W (r_{a b}, h)$ 为光滑核函数， $r_{a b}$ 为粒子 $a$ 与粒子 $b$ 的距离.同样，场函数梯度经过核近似和粒子近似后，表达式为

<\nabla \cdot f (r_{a})> = \sum_{b = 1}^{N} \frac{m_{b}}{ρ_{b}} f (r_{b}) \cdot \nabla_{a} W_{a b}

(4)

其中， $\nabla_{a} W_{a b} = \frac{r_{a b}}{r_{a b}} \cdot \frac{\partial W_{a b}}{\partial r_{a b}} = \frac{(r_{b} - r_{a})}{r_{a b}} \cdot \frac{\partial W_{a b}}{\partial r_{a b}}$ .

2 基于SPH方法的刚-液耦合多体系统动力学方程

2.1　非惯性系下流体控制方程

在SPH方法中，密度对各物理量变化影响很大，因此对粒子密度的计算十分重要.粒子a的密度通过SPH方法近似后，求解表达式为^[

8]：

\frac{d ρ_{a}}{d t} = \sum_{b = 1}^{N} m_{b} v_{a b}^{β} \frac{\partial W {}_{a b}}{\partial r_{a}^{β}}

(5)

其中， $α$ 和 $β$ 表示坐标轴方向， $v$ , $ρ$ 和 $m$ 分别为粒子的绝对速度、密度和质量.

对流体的动量方程，使用SPH方法进行离散，得到求解表达式为：

\begin{array}{l} \frac{d v_{a}^{α}}{d t} = - \sum_{b = 1}^{N} m_{b} (\frac{p_{a}}{ρ_{a}^{2}} + \frac{p_{b}}{ρ_{b}^{2}}) \frac{\partial W {}_{a b}}{\partial r_{a}^{α}} \\ + \sum_{b = 1}^{N} m_{b} (\frac{μ_{a} ε_{a}^{α β}}{ρ_{a}^{2}} + \frac{μ_{b} ε_{b}^{α β}}{ρ_{b}^{2}}) \frac{\partial W {}_{a b}}{\partial r_{a}^{β}} + \frac{F_{o, a}^{α}}{m_{a}} \end{array}

(6)

其中， $p$ 和 $μ$ 分别为粒子的压强和动力粘度系数， $F_{o}^{α}$ 表示作用在粒子上的重力，切向应变率 $ε_{a}^{α β}$ 可表示为

\begin{array}{l} ε_{a}^{α β} = - \sum_{b = 1}^{N} (\frac{m_{b}}{ρ_{b}} v_{b a}^{β} \frac{\partial W {}_{a b}}{\partial r_{a}^{α}} + \frac{m_{b}}{ρ_{b}} v_{b a}^{α} \frac{\partial W {}_{a b}}{\partial r_{a}^{β}}) \\ - (\frac{2}{3} \sum_{b = 1_{}}^{N} \frac{m_{b}}{ρ_{b}} v_{b a}^{} \cdot \nabla_{a} W_{a b}) δ^{α β} \end{array}

(7)

其中，狄拉克函数 $δ^{α β} = \{\begin{matrix} 1, α = β \\ 0, α \neq β \end{matrix}$ .

对于做平面运动的刚体，在刚体质心上建立连体基.通过加速度合成公式，可求出粒子绝对加速度表达式为：

\frac{d v_{a}^{α}}{d t} = a_{e, a}^{α} + a_{c, a}^{α} + \frac{{}^{B} d v_{r, a}^{α}}{d t}

(8)

其中， $a_{e, a}^{α}$ 和 $a_{c, a}^{α}$ 分别表示粒子a的牵连加速度与科氏加速度.粒子牵连加速度与刚体质心加速度、刚体角速度和角加速度相关，粒子科氏加速度与刚体角速度和粒子相对速度相关. $\frac{{}^{B} d v_{r, a}^{α}}{d t}$ 为粒子相对刚体质心上的连体基的加速度.将上式代入流体的动量方程中，可得到

\begin{array}{l} \frac{{}^{B} d v_{r, a}^{α}}{d t} = - \sum_{b = 1}^{N} m_{b} (\frac{p_{a}}{ρ_{a}^{2}} + \frac{p_{b}}{ρ_{b}^{2}}) \frac{\partial W {}_{a b}}{\partial r_{a}^{α}} \\ + \sum_{b = 1}^{N} m_{b} (\frac{μ_{a} ε_{a}^{α β}}{ρ_{a}^{2}} + \frac{μ_{b} ε_{b}^{α β}}{ρ_{b}^{2}}) \frac{\partial W {}_{a b}}{\partial r_{a}^{β}} \\ + \frac{F_{o, a}^{α}}{m_{a}} + \frac{F_{i, a}^{α}}{m_{a}} \end{array}

(9)

其中， $F_{i, a}^{α}$ 为作用在粒子a上的惯性力，表达式为

F_{i, a}^{α} = - m_{a} (a_{e, a}^{α} + a_{c, a}^{α})

(10)

SPH方法在计算过程中，可能会出现粒子发散、非物理性穿透等现象，为了解决这一问题，可采用XSPH方法^[

8]，计算式如下：

\frac{d r_{a}^{α}}{d t} = v_{a}^{α} + ς \sum_{b = 1}^{N} \frac{m_{b}}{ρ_{b}} r_{a b}^{α} W_{a b}

(11)

其中， $ς$ 为 $0 \leq ς \leq 1$ 的常系数.XSPH方法充分考虑了相邻粒子对所求粒子的影响，使相邻粒子之间速度较为接近，从而使得流场运动比较整齐，避免了粒子之间的穿透.

SPH方法中流体粒子的运动主要取决于粒子之间的压力梯度，因此如何求解粒子压强是十分重要的.如果通过泊松方程来求解粒子压强，则需要很小的时间步长，因此，可通过引入人工压缩率，将不可压缩流体视为弱可压缩，得到密度与体积相关的状态方程^[

7]如下：

P = P_{0} [{(\frac{ρ}{ρ_{0}})}^{γ} - 1]

(12)

其中，一般取 $γ = 7$ ,P₀根据不同问题选取不同的值， $ρ_{0}$ 为流体的初始密度.通过状态方程可以看出，即使粒子压强发生很大变化，粒子密度的变化都很小，一般来说，保证粒子密度变化在1%即满足弱可压缩性质.

对于液体晃动问题的固壁边界条件，采用镜像虚粒子法^[

8]进行模拟.对于邻近边界的粒子，镜像对称分布虚粒子，如果固壁边界是不可穿透且不可滑移的，则虚粒子与流体粒子速度大小相等、方向相反；如果固壁边界仅需满足不可穿透条件，则虚粒子与流体粒子在平行于边界的方向上速度相等，垂直边界的速度相反.镜像虚粒子法可大大减少边界对流体的扰动，从而保证计算精度，特别是靠近边界的流场.

与基于网格的数值方法相同，SPH方法中用得比较多的积分格式有龙格库塔法、蛙跳法和预估校正法.由于蛙跳法的计算效率高，因此本文选取蛙跳法.

在SPH方法中，为了计算任意粒子a的物理量及其导数，需要搜索粒子a的相邻粒子.常用的相邻粒子搜索法有基于直观思想的全配对搜索法、链表搜索法和树形搜索法.由于全配对搜索法十分简单直接，故本文采用全配对搜索法.

2.2　液体与刚体之间的相互作用力和力偶

研究刚-液耦合问题时，本文给出了刚体的晃动力和晃动力偶矩的计算方法.通过计算各SPH粒子的绝对加速度，得到各粒子的惯性力对刚体质心的主矢、主矩，根据达朗贝尔原理，计算出向质心简化的晃动力和晃动力偶矩.

$F_{C}^{'} = {[\begin{matrix} F_{C}^{1} & F_{C}^{2} \end{matrix}]}^{T}$ 和 $M_{C}$ 分别为刚体边界作用于粒子的接触力向刚体质心简化后，得到的合力在连体基下的列阵和合力偶矩.通过达朗贝尔原理，可得表达式：

F_{C}^{α} - \sum_{a = 1}^{N} m_{a} \frac{d v_{a}^{α}}{d t} + \sum_{a = 1}^{N} F_{o, a}^{α} = 0, α = 1,2

(13)

M_{C} - \sum_{a = 1}^{N} m_{a} (- y_{a}^{'} \frac{d v_{a}^{1}}{d t} + x_{a}^{'} \frac{d v_{a}^{2}}{d t}) + \sum_{a = 1}^{N} M_{o, a}^{} = 0

(14)

其中，式(13)第二项表示粒子惯性力矢量和在刚体的连体基下的分量，第三项表示粒子重力的矢量和在刚体连体基下的分量.式(14)第二项表示粒子惯性力对刚体质心的力矩，( $x_{a}^{'}$ , $y_{a}^{'}$ )表示粒子在刚体连体基下的坐标，第三项表示粒子的重力对刚体质心的力矩之和.

根据式(13)和(14)，可计算出 $F_{C}^{'}$ 和 $M_{C}$ .设 $F_{L} = {[\begin{matrix} F_{L}^{1} & F_{L}^{2} \end{matrix}]}^{T}$ 为粒子作用于刚体的接触力的主矢在惯性坐标系下的列阵， $M_{L}$ 为接触力对刚体质心的主矩，定义 $F_{L}$ 为液体对刚体的晃动力， $M_{L}$ 为晃动力偶矩，根据牛顿第三定律，可得表达式：

F_{L} = - A F_{C}^{'}

，

M_{L} = - M_{C}

(15)

其中， $A$ 为连体基相对惯性基的方向余弦阵.

2.3　充液刚体的动力学方程

通过2.2中的方法计算出刚体受到的晃动力和晃动力偶矩，可建立充液刚体的平面运动动力学方程如下：

\begin{array}{l} m {\ddot{x}}_{C} = F_{L}^{1} + F_{E}^{1} + F_{N}^{1} \\ m {\ddot{y}}_{C} = F_{L}^{2} + F_{E}^{2} + F_{N}^{2} \\ J_{C} \ddot{ϕ} = M_{L} + M_{E} + M_{N} \end{array}

(16)

其中，m为充液刚体质量，J_C为刚体关于质心的转动惯量， $F_{E} = {[\begin{matrix} F_{E}^{1} & F_{E}^{2} \end{matrix}]}^{T}$ 和 $M_{E}$ 分别为刚体受到主动力的主矢在惯性基下的列阵和关于刚体质心的主矩， $F_{N} = {[\begin{matrix} F_{N}^{1} & F_{N}^{2} \end{matrix}]}^{T}$ 和 $M_{N}$ 分别为刚体受到的理想约束力的主矢在惯性基下的列阵和关于质心的主矩.刚体的动力学方程可用矩阵形式写为

M \ddot{q} = {\hat{F}}_{S}^{} + {\hat{F}}_{E}^{} + {\hat{F}}_{N}^{}

(17)

其中， $q = {[\begin{matrix} x_{C} & y_{C} & φ \end{matrix}]}^{T}$ 为刚体的位形坐标阵， $M$ 为增广质量阵， ${\hat{F}}_{L}^{}$ 、 ${\hat{F}}_{E}^{}$ 和 ${\hat{F}}_{N}^{}$ 分别为增广晃动力阵、增广主动力阵和增广理想约束力阵，分别可表示为

M = [\begin{matrix} m & 0 & 0 \\ 0 & m & 0 \\ 0 & 0 & J_{C} \end{matrix}]

{\hat{F}}_{L}^{} = [\begin{matrix} F_{L}^{1} \\ F_{L}^{2} \\ M_{L} \end{matrix}], {\hat{F}}_{E}^{} = [\begin{matrix} F_{E}^{1} \\ F_{E}^{2} \\ M_{E} \end{matrix}], {\hat{F}}_{N}^{} = [\begin{matrix} F_{N}^{1} \\ F_{N}^{2} \\ M_{N} \end{matrix}]

(18)

设充液刚体的运动学约束方程为 $Φ = 0$ ，则加速度约束方程为：

\ddot{Φ} = Φ_{q} \ddot{q} - γ = 0

(19)

其中， $Φ_{q} = \partial Φ / \partial q$ 为雅可比矩阵， $γ$ 为加速度约束方程的右项，可表示为

γ = - {(Φ_{q} \dot{q})}_{q} \dot{q} - 2 Φ_{q t} \dot{q} - Φ_{t t}

(20)

其中， $Φ_{q t} = \partial Φ_{q} / \partial t$ ， $Φ_{t t} = \partial^{2} Φ / \partial t^{2}$ .

在研究刚-液耦合问题时，使用SPH方法离散的流体连续方程、动量方程，以及刚体的动力学方程、约束方程需要耦合求解.

在初始时刻，通过对上述方程联合求解，消去 $F_{L}$ ， $M_{L}$ ， $F_{N}$ ， $M_{N}$ ，得到初始时刻粒子的密度和相对速度对时间的导数，以及刚体质心加速度和角加速度.引入 $\dot{q}$ ，将刚体的二阶动力学方程扩充为一阶常微分方程组，用蛙跳法分别对粒子的密度和相对速度，刚体的质心加速度和角加速度，以及刚体的质心速度和角速度进行时间积分，得到下一个时刻粒子的密度和相对速度，刚体的质心速度和角速度，以及刚体的质心位置坐标和姿态角，通过时间循环，可以得到以上这些变量随时间的变化曲线.

3 仿真算例

3.1　液体晃动仿真

充液矩形液舱宽度2l=2m，高度h=2m，水深d=1m,水的密度ρ=1000kg/m³.在静水面中央建立非惯性系，如图2所示.已知液舱沿x轴方向运动规律为 $s = s_{0} (1 - c o s Ω t)$ ，其中s₀=0.002m为位移幅值， $Ω$ 为液舱横向运动频率， $Ω_{0} = \sqrt[]{\frac{g π}{2 l} t a n h \frac{π d}{2 l}}$ 为液舱一阶固有频率，且 $Ω / Ω_{0} = 0.8$ .

图2 充液矩形液舱

Fig.2 Rectangular tank filled with water

液舱内部流体建模如图3所示.液舱内流体粒子初始呈规则的矩形分布，粒子初始间距dx=dy=0.025m，问题域内共有 $80 \times 40 = 3200$ 个粒子.施加镜像粒子边界模型条件，光滑长度恒定.液体的初始压强按静压分布，根据状态方程得到流体粒子初始密度，时间步长为dt=0.0002s，采用XSPH方法，取 $ς = 0.3$ .计算了2s内无粘流体的晃动情况.

图3 初始时刻模型

Fig.3 Initial model

t=1s和t=2s时的自由表面的波形曲线，如图4(a)和图4(b)所示.将本文方法的计算结果与文献[

1]给出的小幅晃动下液面高度解析解的计算结果进行对比，可以看出，当液舱运动幅值较小时，内部流体做微幅晃动，本文方法的计算结果与理论解基本一致，说明了本文方法研究充液液舱内液体微幅晃动的可行性.

(a) t=1s

(b) t=2s

图4 本文结果与解析解对比

Fig.4 Comparison of the present result and exact solution

为了验证本文方法在大幅晃动时的准确性，给出以下算例.

矩形水箱长和高分别为0.96m和1m，水深为0.192m，在正弦激励 $x (t) = D s i n (ω t)$ 的作用下运动，其中D=0.005m， $ω$ =4.387rad/s.边界条件采用镜像粒子法，计算得到15s内液舱左壁晃动波高，并与文献[

1]的解析解和文献[3]中的实验结果进行比较，如图5所示.可以看出，在大幅晃动情况下，本文方法结果与文献[3]实验结果吻合很好，说明本文方法研究充液液舱内液体大幅晃动的准确性.由于文献[1]的解析解是在小幅晃动假设下给出的，与实验结果具有一定误差.

图5 液舱左壁晃动波高时变曲线

Fig.5 Time history of the wave height of the left wall

3.2　刚-液耦合问题

下面通过本文SPH方法研究水箱受周期性激励力平动时的刚-液耦合问题.如图6所示，矩形水箱几何尺寸和静水面深度与3.1中相同，充液矩形水箱受到周期性激励力F的作用，在光滑的地面上平动，水箱的运动与内部液体的晃动相互影响.

图6 充液水箱示意图

Fig.6 Diagram of tank filled with water

对于充液水箱的平动问题，位形坐标阵为 $q = {[\begin{matrix} x_{C} & y_{C} & ϕ \end{matrix}]}^{T}$ ，其中 $(x_{C}, y_{C})$ 为水箱在惯性基下质心的坐标， $ϕ$ 为水箱的姿态角.

水箱动力学模型如图7所示：

图7 水箱简化的动力学模型

Fig.7 Simplified dynamic model of tank

水箱平面运动的动力学方程为：

m {\ddot{x}}_{C} = F - F_{L x}^{'}

(21)

m {\ddot{y}}_{C} = F_{N} - m g - F_{L y}^{'}

(22)

J_{C} \ddot{ϕ} = M_{N} - M_{L}^{'}

(23)

其中，m为水箱质量，F为周期性激励力， $F = F_{0} s i n (π t)$ ，根据牛顿第三定律， $F_{L x}^{'}$ 、 $F_{L y}^{'}$ 和 $M_{L}^{'}$ 为流体作用于水箱的晃动力和晃动力偶矩， $F_{N}$ 和 $M_{N}$ 为地面作用于水箱的支持力和力矩， $J_{C}$ 为水箱绕质心的转动惯量，由于地面光滑，故不存在地面摩擦力.

水箱的加速度约束方程为：

\ddot{Φ} = [\begin{matrix} {\ddot{y}}_{C} \\ \ddot{ϕ} \end{matrix}] = 0

(24)

设F₀=30N,m=1000kg，通过仿真计算得出水箱位移和速度的时间历程如图8和图9所示，其中，实线表示考虑液体晃动的情况，而虚线表示液体全充满时无液体晃动的情况.

图8 水箱位移时间历程 (F₀=30N, m =1000kg)

Fig.8 Displacement of tank (F₀=30N, m =1000kg)

图9 水箱速度时间历程 (F₀=30N, m =1000kg)

Fig.9 Velocity of tank (F₀=30N, m =1000kg)

从图8和图9可以看出，如果不考虑液体晃动的影响，相当于外力作用于一个质量为水箱质量和箱内液体质量之和的刚体上，水箱速度呈周期性变化；如果考虑液体晃动的影响，水箱速度不再呈周期性变化，振动的幅值低于不考虑液体晃动的情况，这是因为一部分平动能量转变为液体晃动的能量.此外可以看出，由于刚体质量较大，周期性激励力较小，只引起微幅晃动，液体晃动对水箱运动速度和位移的影响很小.

设F₀=30N,m=100kg，水箱位移和速度的时间历程如图10和图11所示.通过数值对比可以看出，不考虑液体晃动时，水箱的速度呈周期性变化；考虑液体晃动时，速度振荡的幅值低于不考虑液体晃动的情况.当周期性激励力幅值保持不变时，随着水箱质量减小，两种模型速度的差异更加显著.t=1s和t=2s时的自由表面波形如图12-13所示，可以看出由于激励力F较小，液体晃动的幅值较小，但是由于水箱质量较小，液体晃动对水箱运动的影响不容忽视.与速度相比，两种模型的水箱位移的差异较小.

图10 水箱位移时间历程 (F₀=30N, m =100kg)

Fig.10 Displacement of tank (F₀=30N, m =100kg)

图11 水箱速度时间历程 (F₀=30N, m =100kg)

Fig.11 Velocity of tank (F₀=30N, m =100kg)

图12 t= 1s时自由表面波形(F₀=30N, m=100kg)

Fig.12 Free surface when t=1s (F₀=30N, m=100kg)

图13 t = 2 s时自由表面波形 (F₀=30N, m =100kg)

Fig.13 Free surface when t = 2s (F₀=30N, m=100kg)

设F₀=500N,m=100kg，水箱位移和速度的时间历程如图14和图15所示，t=1s和t=2s时的自由表面波形如图16-17所示.从图中可以看出，当水箱质量不变时，随着激励力的幅值增大，液体振荡幅值显著增大，考虑和不考虑液体晃动两种情况下水箱速度的差异非常显著，呈现非线性特征，导致水箱位移也出现了明显的差异.仔细观察自由表面波形可以发现，在大幅晃动情况下，中点处波形的斜率不为零，而微幅晃动情况下，中点处波形的斜率为零（如图12-13所示）.

图14 水箱位移时间历程 (F₀=500N, m =100kg)

Fig.14 Displacement of tank (F₀=500N, m =100kg)

图15 水箱速度时间历程 (F₀=500N, m=100 kg)

Fig.15 Velocity of tank (F₀=500N, m=100 kg)

图 16 t = 1 s时自由表面波形(F₀=500N, m=100kg)

Fig.16 Free surface when t=1s(F₀=500N, m=100kg)

图 17 t=2s时自由表面波形(F₀=500N, m=100kg)

Fig.17 Free surface when t=2s (F₀=500N, m=100kg)

4 结论

本文基于SPH方法建立了可以适用于大幅晃动的刚-液耦合充液刚体的动力学模型，不仅适用于平动刚体，也适用于定轴转动和平面一般运动刚体.首先建立了非惯性系下液体晃动的连续性方程和动量方程，得到液体对外部刚体的晃动力和晃动力偶矩的表达式，在此基础上建立充液刚体的动力学方程.将液体晃动的连续性方程和动量方程，以及充液刚体的动力学方程和加速度约束方程联立计算，得到不同时刻液体晃动变量和刚体运动变量.

首先研究了外部液舱运动已知时内部液体的晃动问题.在微幅晃动情况下，将本文SPH模拟结果与理论解对比，验证了本文动力学模型对于小幅晃动问题的适用性，然后在大幅晃动情况下，将本文SPH模拟结果与实验结果进行对比，验证了本文动力学模型计算大幅晃动问题的准确性.在此基础上对比分析了考虑和不考虑液体晃动的两种模型的计算结果，研究了考虑刚-液耦合时，液体晃动对受周期性激励力作用的平动充液水箱的动力学特性的影响.结果表明当激励力的幅值不变时，随着水箱质量减小，水箱内液体的晃动会对水箱的运动产生明显影响.当水箱的质量不变时，随着激励力的幅值增大，液体晃动的幅值增大，液体晃动对充液水箱的动力学特性的影响愈加显著.

参考文献

Wu G X.Second-order resonance of sloshing in a tank. Ocean Engineering, 2007, 34: 2345~2349 [百度学术]

李清，余本嵩，金栋平.圆柱形贮箱液晃系统稳定性边界分析. 动力学与控制学报, 2017, 15(5): 467~471 [百度学术]

Li Q, Yu B S, Jin D P. Stability boundaries of liquid-sloshing system for a cylindrical tank. Journal of Dynamics and Control, 2017, 15(5): 467~471（in Chinese） [百度学术]

Goudarzi M A, Sabbagh-Yazdi S R. Investigation of nonlinear sloshing effects in seismically excited tanks. Soil Dynamics and Earthquake Engineering, 2012, 43: 355~365 [百度学术]

石怀龙，王勇，邬平波.基于拉格朗日描述的罐车内液体晃动模拟.动力学与控制学报,2018,16(2):157~164 [百度学术]

Shi H L, Wang Y, Wu P B. Liquid sloshing simulating in railway tank car cased on total Lagrangian approach. Journal of Dynamics and Control, 2018,16(2):157~164（in Chinese） [百度学术]

Gingold R A, Monaghan J J. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Royal Astronomical Society, 1977, 181: 375~389 [百度学术]

Monaghan J J,Kos A.Solitary waves on a cretan beach. Journal of Waterway, Port, Coastal, and Ocean Engineering, 1999, 125:145~155 [百度学术]

Monaghan J J. Smoothed particle hydrodynamics. Reports on Progress in Physics, 2005, 68: 1703~1759 [百度学术]

郑兴.SPH方法改进研究及其在自由面流动问题中的应用[博士学位论文] .哈尔滨：哈尔滨工程大学, 2010（Zheng X. An investigation of improved SPH and its application for free surface flow[Ph.D Thesis].Harbin: Harbin Engineering University, 2010（in Chinese）） [百度学术]

刘富，童明波，陈建平.基于SPH方法的三维液体晃动数值模拟.南京航空航天大学学报,2010, 42(1):122~126 [百度学术]

Liu F, Tong M B, Chen J P. Numerical simulation of three-dimensional liquid sloshing based on SPH method. Journal of Nanjing University of Aeronautics & Astronautics , 2010, 42(1): 122~126（in Chinese） [百度学术]

Rafiee A, Pistani F, Thiagarajan K. Study of liquid sloshing: numerical and experimental approach. Computational Mechanics, 2011, 47(1): 65~75 [百度学术]

李大鸣，李玲玲，高伟，等.液舱内液体晃荡的SPH模型.哈尔滨工程大学学报, 2012, 33(1): 37~41 [百度学术]

Li D M, Li L L, Gao W, et al. A model of liquid sloshing in tanks based on the smoothed particle hydrodynamics(SPH) method. Journal of Harbin Engineering University, 2012, 33(1): 37~41（in Chinese） [百度学术]

Negrut D, Tasora A, Mazhar H, et al. Leveraging parallel computing in multibody dynamics. Multibody System Dynamics, 2012, 27(1): 95~117 [百度学术]

Hu W, Pan W X, Rakhsha M, et al. A consistent multi-resolution smoothed particle hydrodynamics method.Computer Method in Applied Mechanics and Engineering, 2017, 324:278~299 [百度学术]

Hu W, Tian Q, Hu H Y. Simulating coupled dynamics of a rigid-flexible multibody system and compressible fluid. Science China,Physics,Mechanics & Astronomy,2018, 61(4): 044711-1~15 [百度学术]

樊伟.考虑多场耦合的多体系统动力学[硕士学位论文].上海：上海交通大学, 2013 [百度学术]

Fan W. Multi-field coupling dynamics for multi-body system[Master Thesis].Shanghai: Shanghai Jiao Tong University,2013（in Chinese） [百度学术]

岳宝增. 液体大幅晃动动力学，北京：科学出版社，2012（Yue B Z. Dynamics of large-amplitude liquid sloshing, Beijing：Science press, 2012（in Chinese）） [百度学术]

岳宝增, 唐勇.球形贮箱中三维液体大幅晃动数值模拟, 宇航学报, 2016, 37(12): 1405~1410 [百度学术]

Yue B Z, Tang Y. Numerical simulation of three dimensional large-amplitude liquid sloshing in spherical containers. Journal of Astronautics, 2016, 37(12): 1405~1410（in Chinese) [百度学术]

Tang Y, Yue B Z, Yan Y L. Improved method for implementing contact angle condition in simulation of liquid sloshing under microgravity. International Journal for Numerical Methods in Fluids, 2018, 89: 123~142 [百度学术]

岳宝增，于嘉瑞，吴文军.多储液腔航天器刚液耦合动力学与复合控制.力学学报, 2017, 49(2): 390~396 [百度学术]

Yue B Z, Yu J R, Wu W J. Rigid and liquid coupling dynamics and hybrid control of spacecraft with multiple propellant tanks. Chinese Journal of Theoretical and Applied Mechanics, 2017,49(2):390~396（in Chinese） [百度学术]

基于SPH的充液刚体刚-液耦合动力学建模

摘要

关键词

引言

1 SPH方法基本理论