引 言
流体晃荡现象普遍存在于航空、航天、船舶与海洋工程等领域内，是一种常见的流体运动现象，通常发生在部分装载的液舱中．当外界激励频率与液舱内流体自由液面固有频率相接近时，由晃荡引起的冲击荷载使舱壁结构产生较大的结构响应，进而可能导致结构破坏．在船舶与海洋工程领域中，海上载液货轮的摇晃会引起液舱内流体晃荡，如果液舱内的流体晃荡剧烈，则可能导致倾覆．在航空航天领域内，燃料储箱中的液体晃荡会影响飞行姿态，可能会引发事故[1]．因此，晃荡问题是多领域学者们关注的热点问题．
由于流体晃荡是一种非常复杂的流体运动现象，具有高度的非线性及随机性，其理论研究还不完善，且随着计算机技术的飞速发展，数值模拟方法逐渐成为研究晃荡问题的重要手段．朱仁庆等[2]基于CFD技术，提出了数值预报 LNG 船薄膜型液舱晃荡压强的方法．宁德志等[3]采用高阶边界元方法，对纵摇容器中的液体晃荡问题进行了数值模拟研究．祁江涛等[4]采用VOF方法，并结合动网格技术对部分装载液舱的晃荡进行了仿真计算．卫志军等[5]基于SPH方法，对二维矩形舱的晃荡问题进行了研究．于强等[6]基于光滑粒子流体动力学理论，对航天器贮箱内的液体大幅度晃动进行了动力学分析．在液舱内部增设隔板是常见的阻晃方法，有许多国内外学者对阻晃隔板进行了研究．Faltinsen 等[7]提出了一种在液舱中心加入旋转隔板的抑制液体晃荡方法，多孔隔板垂直于流体运动方向，流体穿过多孔板会导致能量耗散．宁德志等[8]对三维液舱内不同位置、吃水高度和半径的浮子抑制晃荡的效果进行了研究．朱小松等[9]采用光滑粒子水动力学方法，研究了二维液舱内T型和I型隔板对抑制矩形液舱内流体晃荡的影响，分析了隔板高度和宽度对晃荡的影响．黄志涛等[10]基于光滑粒子法，对罐车的行驶稳定性、晃荡抑制措施进行了研究．Wang等[11]首次提出等几何边界元法(IGABEM)来研究水平矩形、圆形或椭圆环形二维多孔挡板储罐的流体晃荡问题．张友林等[12]基于MPS-FEM耦合方法开发了流固耦合求解器，利用该求解器对带刚性挡板和弹性挡板的液舱晃荡问题进行了研究．Xue 等[13]基于数值研究的方法，对环形挡板的阻晃效应进行了研究．上述学者的研究均是针对液舱内固定的隔板形状来进行的，而Guan等[14]采用 VOF 方法对液舱内单个隔板的晃荡问题进行了数值计算，并结合拓扑优化技术得到了使壁面受力最小的隔板形态．因此，采用拓扑优化技术研究阻晃隔板形状对晃荡的抑制效果是可行的，本文针对液舱内固定形状的双隔板以及拓扑优化后的双隔板的晃荡问题进行了数值计算，并分析了增设双隔板后流场的运动学和动力学特性．
本文的主要内容如下：首先给出了本文所研究的晃荡问题的数学模型；其次分析了液舱内增设双隔板后的运动学和动力学特性；最后分析了拓扑优化后双隔板的减晃效果和结论．
1 数 学 模 型
1.1 控制方程
本文应用的控制方程是均质不可压缩流体的三维非定常Navier-Stokes方程：
div v=0,
(1)
2v，
(2)
式中，v表示流体的速度矢量，Fb表示单位质量流体上的质量力矢量，p表示压力，μ表示黏性系数．
1.2 自由液面处理
本文采用流体体积法，即VOF方法来处理自由液面．VOF方法是由Hirt和Nichols首先提出的一种模拟界面的方法[15]．VOF方法的基本原理是通过定义一个流体体积函数F来确定自由液面，针对某一计算单元，体积分数会出现以下三种情况：
1) F=0，代表该单元内全部为气体；
2) F=1，代表该单元内全部为液体；
3) 0<F<1，代表该单元内既有气体又有液体，存在气液交界面．
晃荡过程中的自由液面就是当F介于0～1之间的情况，在找出自由液面的位置后，可以通过F的梯度来确定自由液面的发展方向，从而近似确定自由液面的形状．
1.3 界面追踪
界面追踪是上述VOF方法中的关键步骤，界面追踪包括多种方法：最小方差法、ELVIRA法、PLIC方法(界面重构法)等[16]．本文采用PLIC方法来追踪界面．对于PLIC界面重构技术，在每个时间步长，对每个计算网格，若指定该网格的体积分数和界面法向量，就可以在该网格中构造出一个法向量一致的平面，该平面会随着流动而传播，从而必须求出每种流体流入相邻网格中的体积、质量和动量通量．更新了整个区域中流体流入的体积分数、质量场和动量场，就可对下一步的流动进行数值模拟．因此，界面追踪过程可简化为以下三个步骤：首先，在交错网格下，利用经典中心差分格式，估计界面的法向量；其次，利用界面法向量和体积函数F对界面位置进行确定和重构；最后，利用Lagrange方法跟踪计算界面的传播．
1.4 最优化方法及其求解过程
本文所采用的优化算法是一种基于求解常微分方程的方法．优化问题的目标是选定目标参数x1，x2，x3，…，xn,使函数ω值尽可能小，优化问题如下：
(3)
式中m1，m2，m3是满足关系0≤m1≤m2≤m3的整数，fi是一个实数，gi是一个连续方程．其迭代求解过程如下：
1) 提供初始近似解
初始目标值
和步长h>0;
2) 判断约束性质；
3) 建立并求解一阶常微分方程组的初值问题，找到搜索方向；
4) 进行线性搜索以减小目标值
当优化目标参数的值超出可行域时，对优化参数进行修正，使其返回可行域；
5) 调整搜索步长h，并判断是否进行下一次迭代．
本文利用拓扑优化技术对隔板形状进行优化设计，以得到壁面受力最小时的隔板形态．在本文使用的算法中，以隔板透水率为自变量，壁面所受的动水压力合力为目标函数．
自变量隔板透水率Ω为每个网格透水部分面积与单个网格总面积之比．Ω的取值范围在0～1之间，当Ω=0时，表示隔板不透水，即隔板上面不存在“洞”；当Ω=1时，表示隔板完全透水，即不存在隔板．在优化程序中，我们将每个隔板划分为6×6的网格，每个网格对应一个透水率参数，故在本次模拟中需要对两个隔板上的72个参数进行优化计算．由上述可知，透水率Ω的优化约束条件如下：
0≤Ω(m,n)≤1, m=1，2，…，6, n=1，2，…，6,
(4)
式中Ω(m,n)表示隔板上不同位置处的透水率参数值．在本次模拟计算中，为了简化拓扑优化程序的计算过程，人为将Ω的值界定在0.3～0.7之间，即当Ω=0.3时，隔板不透水，当Ω=0.7时，隔板透水．
目标函数壁面动水压力合力A的计算公式如下：
A=|∑fforce(ijk)|,
fforce(ijk)=p(ijk)-ρghs； i=1； j=1，2，3，…， jmax； k=1，2，3，…，kmax,
(5)
式中A为壁面动水压力的合力，i， j，k分别表示x，y，z方向上的网格数，p为压力,ρ表示液体密度，g表示重力加速度，h表示距离液面高度，s表示单位网格面积．
在晃荡过程中实现优化的具体过程为：首先设定初始透水率Ω和优化开始时间tstart，以Ω=0.5为初始值，此时隔板处于半透水状态，同时记录此时目标函数A0的数值；其次当液体晃荡一段时间之后(t>tstart)，每增加一个时间步长Δt时都会进行优化计算，优化后的透水率会替代上一时刻的透水率，再调用流体晃荡模块计算此时目标函数值A1；最后判断此时目标函数A1与A0的大小，若A1是最优，则继续计算下一时刻直到晃荡结束，反之则重复优化程序, 其优化流程图如图1所示．
图1 优化算法流程图
Fig. 1 The flow chart for the optimization algorithm
2 液舱中增设双隔板后流场的运动学和动力学特性
2.1 研究对象
本文的物理量均是无量纲的．模拟液舱的长宽高分别为12×4.8×1.5，液舱内液面高度为1．激励函数为
u=0.8cos(50t),
(6)
其中液舱沿其长度方向(x轴方向)进行周期性运动，运动周期为T=0.126．沿液舱的长度方向均匀设置两个隔板，即在x=4位置处设置第一块隔板，在x=8位置处设置第二块隔板，隔板的高度与液面高度相同，即隔板高度为1．液舱的边界条件设置为固定边界无滑移．本文对液舱内固定形状的双隔板和拓扑优化的双隔板进行了数值计算．液舱增设双隔板后的示意图如图2所示．
2.2 增设双隔板后流场的运动学特性
液体晃荡是一种具有大幅自由液面运动的强非线性现象，尤其是在共振区晃荡、大幅值运动时，非线性效应更加明显．为了研究液体在晃荡过程中的运动学特性，我们所选取的激励振幅为小振幅，激励频率远离共振区．我们对晃荡的第一个周期进行分析，液舱内液体的流场图如图3所示．
由液舱内液体的流场图可以看出，液舱内隔板处的流速方向与外界激励保持一致，这表示隔板固定在液舱上，随液舱一起运动．当t=0时，液舱内液体开始向x轴正方向运动，从图3(a)、(c)、(e)中可以看出：隔板两侧的液体受隔板的影响而导致液体的流动方向发生变化；当t=0.03时，液体的流速降低，其流动方向开始发生改变；当t=0.06时，液体沿x轴负向运动；当t=0.09时，液体流速再次降低，其流动方向发生改变；当t=0.12时，液舱内液体重新变为沿x轴正方向运动，至此液体晃荡完成一个周期的变化．从图3可以看出，液体流速的方向随着时间进行周期性变化．在刚开始晃荡的时候，液体的运动总趋势是随着晃荡驱动的方向刚运动．在左右边界及隔板处，液体的运动产生了较大的速度变化，当向x轴正向运动时，隔板左侧边界大部分流体向下运动，右侧边界流体大部分向上运动;当向x轴负向运动时，则出现了相反的情况，隔板两侧液体速度的变化可以解释在不同区域内液体表面一高一低的现象．隔板两侧液体表面一高一低的现象在文献[17]的物理模型实验中得到了证实．
图2 双隔板位置示意图
Fig. 2 Schematic diagram of double baffle positions
(a) t=0.01
(b) t=0.03 (c) t=0.06
(d) t=0.09 (e) t=0.12
图3 液舱内流场矢量图
Fig. 3 The velocity vectors of liquid flow fields in the tank
2.3 增设双隔板后流场的动力学特性
液体晃荡易对结构造成砰击，若晃荡剧烈则可导致结构破坏，因此有效抑制液体晃荡荷载是工程和科学界研究的重点，需要对液体晃荡过程中的动力学特性进行研究．我们对前两个周期内液舱两侧壁面以及隔板受到的液体冲击压力进行分析，晃荡过程中液舱内的动压图如图4所示．
(a) t=0.03 (b) t=0.06
(c) t=0.09 (d) t=0.12
(e) t=0.15 (f) t=0.18
(g) t=0.21 (h) t=0.24
图4 晃荡过程中液舱内的动压图
Fig. 4 Dynamic pressure diagrams in the tank during sloshing
注 为了解释图中的颜色，读者可以参考本文的电子网页版本，后同．
图4中红色区域表示高压，蓝色区域表示低压．如图4(a)～(d)所示，当t=0.03时，晃荡幅度达到正向最大，此时液舱内隔板所受的液体冲击压力明显，隔板左侧直接受到液体的冲击，压力区域的变化最大；当t=0.06时，隔板两侧冲击压力趋于零；当t=0.09时，晃荡幅度达到负向最大，隔板右侧直接受到液体的冲击，压力区域的变化最大；当t=0.12时，隔板两侧冲击压力重新趋于零．图4(e)～(f)是第二周期内液舱内的动压图，由此可知隔板的动压力也是周期性变化的．从图4(f)和图4(h)可以看出，由于后期晃荡自由液面的非线性，在液体并未对隔板直接冲击的情况下，隔板两侧以及液舱壁面处仍有压力存在．从以上计算结果显示，当晃荡位移达到极值时会出现压力极值；随着时间的增加，由于晃荡自由液面的非线性，即使不对隔板造成直接冲击，隔板两侧仍有压力存在．
3 隔板拓扑优化后对晃荡的抑制效果
3.1 优化后隔板形态
在优化程序中，我们将每个隔板划分为6×6的网格，每个网格对应一个透水率参数，故在本次模拟中需要对两个隔板上的72个参数进行优化计算．在每一个时间步内，均可以得到72个透水率参数的数值，进而得到当前时刻下使壁面受力最优的透水率参数．由于隔板的透水率在每一时刻都是变化的，因此隔板形态也在时刻变化．因此我们将每一时刻的透水率进行平均处理，求出晃荡过程中不同时刻透水率的平均值，再将透水率平均值赋值到晃荡程序中进行计算，查看此时的壁面受力情况．取平均值之后的平均透水率数值如图5所示，隔板透水率云图如图6所示．在图6中，X是沿着液舱宽度方向，即垂直于液舱激励方向(y轴)，Y是沿着液舱高度方向(z轴)．平均隔板形态如图7所示．
(a) 第一块隔板平均透水率 (b) 第二块隔板平均透水率
(a) The average permeability of the first baffle(b) The average permeability of the second baffle
图5 平均透水率数值
Fig. 5 Average water permeability values
(a) 第一块隔板平均透水率云图 (b) 第二块隔板平均透水率云图
(a) The average permeability contour of the first baffle(b) The average permeability contour of the second baffle
图6 平均透水率云图
Fig. 6 The contour of the average permeability of the baffles
(a) 第一块隔板形状示意图 (b) 第二块隔板形状示意图
(a) The shape diagram of the first baffle (b) The shape diagram of the second baffle
图7 平均隔板示意图
Fig. 7 Average partition diagrams
图5～7分别表示了隔板透水率平均值、隔板透水率云图以及隔板形态示意图．由图5可知，隔板上最下层网格的透水率参数均为0.302，即经过优化计算后的隔板最底层基本不透水，其他网格上的透水率参数各不相同．从透水率云图和隔板形态示意图可以看出，第一块隔板呈现条状分布，隔板由三个条状板组成，第二块隔板的形态较为复杂，不具有规律性，这也符合拓扑优化随机“掏洞”的特点．
3.2 受力结果分析
为了检验拓扑优化技术的应用是否可以减小壁面所受到的冲击作用力，我们将上一小节得到的透水率参数平均值重新赋值到程序中进行计算．模拟计算得到的壁面冲击合力与未优化前的冲击合力对比如图8所示．图中点划线表示无隔板时的壁面合力，虚线表示隔板未优化前的壁面合力，实线表示隔板优化后的壁面合力．由图8可知，增设双隔板之后的壁面合力有时大于无隔板情况，有时小于无隔板情况；双隔板拓扑优化之后的壁面合力明显小于未优化前和无隔板情况下的壁面合力．结果表明，对隔板进行拓扑优化可以降低流体晃荡对液舱壁面的冲击作用．因此，对隔板进行拓扑优化是降低晃荡冲击荷载的一个重要思路．
此外，我们还与文献[14]的单隔板拓扑优化计算结果进行了对比，其中单隔板位于液舱底部中间位置处，即在x=6处设置单隔板．单隔板与双隔板的壁面合力对比如图9所示．图中实线表示双隔板优化后的壁面合力，虚线表示单隔板优化后的壁面合力．从图中可以明显看出，单隔板和双隔板经过拓扑优化之后，双隔板情况下的壁面受力明显低于单隔板情况下的壁面受力，则经过优化之后的双隔板比单隔板对液体晃荡冲击荷载的抑制效果更佳．
图8 双隔板拓扑优化前后合力的对比图 图9 双隔板与单隔板拓扑优化后的合力对比图
Fig. 8 Comparison of forces before and Fig. 9 Comparison of forces between double and after topology optimization single baffles after topology optimization
4 结 论
本文对液舱晃荡进行了数值模拟计算．通过分析液舱增设双隔板后的运动学和动力学特性发现：液体的流速和隔板所受的动压力随着液舱的往复运动周期性变化；在晃荡过程中，隔板两侧产生了较大的速度变化，与物理模型实验中隔板两侧自由液面一高一低的现象相符合；在后期晃荡过程中，由于自由液面的非线性及复杂性，即使液体没有直接冲击隔板，隔板处也有动水压力的作用．通过引入拓扑优化技术对双隔板进行优化计算，发现拓扑优化后的双隔板可以降低流体晃荡对液舱壁面的冲击作用．因此，对隔板进行拓扑优化计算，为船舶与海洋工程领域和航空航天领域中的减晃问题研究提供了一个重要思路．
[1] 朱仁庆. 液体晃荡及其与结构的相互作用[D]. 博士学位论文. 无锡: 中国船舶科学研究中心, 2002.(ZHU Renqing. Time domain simulation of liquid sloshing and its interaction with flexible structure[D]. PhD Thesis. Wuxi: China Ship Scientific Research Center, 2002.(in Chinese))
[2] 朱仁庆, 马海潇, 缪泉明, 等. LNG船液舱晃荡压强预报[J]. 船舶力学, 2013, 17(1/2): 42-48.(ZHU Renqing, MA Haixiao, MIU Quanming, et al. Prediction of pressure induced by liquid sloshing for LNG carrier[J]. Journal of Ship Mechanics, 2013, 17(1/2): 42-48.(in Chinese))
[3] 宁德志, 宋伟华, 滕斌. 纵摇容器中液体晃荡的非线性数值模拟[J]. 船舶力学, 2017, 21(1): 15-22.(NING Dezhi, SONG Weihua, TENG Bin. Nonlinear numerical simulation of liquid sloshing in a container subjected to pitch excitation[J]. Journal of Ship Mechanics, 2017, 21(1): 15-22.(in Chinese))
[4] 祁江涛, 顾民, 吴乘胜. 液舱晃荡的数值模拟[J]. 船舶力学, 2008, 12(4): 574-581.(QI Jiangtao, GU Min, WU Chengsheng. Numerical simulation of sloshing in liquid tank[J]. Journal of Ship Mechanics, 2008, 12(4): 574-581.(in Chinese))
[5] 卫志军, 张文首, 王安良, 等. 基于SPH方法的二维矩形舱液体晃荡数值研究[J]. 大连理工大学学报, 2014, 54(6): 597-603.(WEI Zhijun, ZHANG Wenshou, WANG Anliang, et al. Numerical investigation of liquid sloshing in a 2D rectangular tank based on SPH method[J]. Journal of Dalian University of Technology, 2014, 54(6): 597-603.(in Chinese))
[6] 于强, 王天舒. 航天器贮箱内液体大幅晃动动力学分析[J]. 中国科学: 物理学 力学 天文学, 2019, 49(2): 127-134.(YU Qiang, WANG Tianshu. Dynamic analysis of large-scale amplitude liquid sloshing in the spacecraft[J]. Scientia Sinica: Physica, Mechanica & Astronomica, 2019, 49(2): 127-134.(in Chinese))
[7] FALTINSEN O M, TIMOKHA A N. Sloshing[M]. Cambridge: Cambridge University Press, 2009.
[8] 宁德志, 苏朋, 张崇伟, 等. 三维液舱内浮子式减晃荡结构的水动力特性[J]. 哈尔滨工程大学学报, 2019, 40(1): 154-161.(NING Dezhi, SU Peng, ZHANG Chongwei, et al. Hydrodynamic characteristics of an anti-sloshing floating body structure in a 3D tank[J]. Journal of Harbin Engineering University, 2019, 40(1): 154-161.(in Chinese))
[9] 朱小松, 滕斌, 吕林, 等. 液舱内不同结构形式对晃荡的影响分析[J]. 水道港口, 2011, 32(4): 297-304.(ZHU Xiaosong, TENG Bin, LÜ Lin, et al. Analysis on the effect of sloshing in different liquid tank[J]. Journal of Waterway and Harbor, 2011, 32(4): 297-304.(in Chinese))
[10] 黄志涛, 杨瑜, 邵家儒, 等. 罐车防晃结构SPH模拟研究[J]. 应用数学和力学, 2020, 41(7): 760-770.(HUANG Zhitao, YANG Yu, SHAO Jiaru, et al. Numerical simulation of sloshing-mitigating structures in tank trucks with the SPH method[J]. Applied Mathematics and Mechanics, 2020, 41(7): 760-770.(in Chinese))
[11] WANG W Y, ZANG Q S, WEI Z J, et al. An isogeometric boundary element method for liquid sloshing in the horizontal eccentric annular tanks with multiple porous baffles[J]. Ocean Engineering, 2019, 189: 106367.
[12] 张友林, 陈翔, 万德成. 基于MPS-FEM耦合方法对比研究刚性与弹性挡板对液舱晃荡的抑制作用[J]. 应用数学和力学, 2016, 37(12): 1359-1377.(ZHANG Youlin, CHEN Xiang, WAN Decheng. An MPS-FEM coupled method for the comparative study of liquid sloshing flows interacting with rigid and elastic baffles[J]. Applied Mathematics and Mechanics, 2016, 37(12): 1359-1377.(in Chinese))
[13] XUE M A, LIN P. Numerical study of ring baffle effects on reducing violent liquid sloshing[J]. Computers & Fluids, 2011, 52: 116-129.
[14] GUAN H, XUE Y F, WEI Z J, et al. Numerical simulations of sloshing and suppressing sloshing using the optimization technology method[J]. Applied Mathematics and Mechanics(English Edition), 2018, 39(6): 845-854.
[15] HIRT C W, NICHOLS B D. Volume of fluid (VOF) method for the dynamics of free boundaries[J]. Journal of Computational Physics, 1981, 39(1): 201-225.
[16] KAUFMAN E H, LEEMING D J, TAYLOR G D. An ODE-based approach to nonlinearly constrained minimax problems[J]. Numerical Algorithms, 1995, 9(1): 25-37.
[17] WEI Z J, FALTINSEN O M, LUGNI C, et al. Sloshing-induced slamming in screen-equipped rectangular tanks in shallow-water conditions[J]. Physics of Fluids, 2015, 27(3): 03210.