改进的湿模态法在流固耦合中的应用

刘瑞骏， 郝志勇， 郑 旭， 谈江林

(浙江大学 能源工程学院，杭州 310027)

改进的湿模态法在流固耦合中的应用

刘瑞骏， 郝志勇， 郑 旭， 谈江林

(浙江大学 能源工程学院，杭州 310027)

针对流固耦合运动方程，指出了流体刚度矩阵的奇异性将导致现有的湿模态法无法直接使用，并指出了在湿模态法中使用分离零频项的方法存在的矛盾。湿模态法的问题根源在于为了简化计算而提出的不合理的前提假设，即同时忽略流体的可压缩性及自由液面波动的影响；据此，提出了直接求解法和不动点法两种解决方法；前者适合流体规模较小的计算，而后者更适合大规模流体的耦合计算；应用不动点法，编制了Matlab程序，对某油底壳在不同盛油状态时的结构特性进行了仿真计算，对比表明随着盛油量从0逐步增加，结构模态频率呈先减小后增大的趋势并趋于平缓；加入流体后，结构模态的振幅大幅下降。除了第1阶模态，流体较少时的湿模态振型与干模态振型区别较大，但随着流体的增加，振型逐渐趋于稳定，与干模态振型接近。

流固耦合；湿模态；直接求解法；不动点法

与流体相接触的结构在发生振动时，其周围流场也会同时产生变化。这种流场的变化反过来会使得结构所受到的流体动力发生变化，形成反馈的流体-结构相互作用的流固耦合问题。在海洋、船舶、航空、水利等许多工程领域都存在此类流固耦合的振动问题。

当下，对于与流体的相对运动不大的结构，计算其流固耦合的结构特性时，通常可以把流体视为无黏、无旋的理想流体，采用的方法主要有虚质量法、附加质量法、边界元法以及湿模态法[1-4]。Lewis等假定半沉入水中的圆柱体周围的水流状态与全沉入水中的相同，利用已知速度势得到附加水质量的计算公式[5]。柳瑞锋等[6]对比了虚质量法和附加质量法在计算船体低阶湿模态时的误差，结果表明两种方法的精度均满足工程要求，且虚质量法的计算过程更简单，结果也较精确。吴绍亮等[7]使用流体边界元法对船体局部结构的流体附加质量进行计算，与经验公式法的对比表明两者对于简单结构均有足够的工程计算精度，但对于复杂结构，经验公式法的精度不如流体边界元法。

为了提高计算结构动力响应时的效率，模态叠加法[8]是不错的选择，此时就需要利用固有模态矩阵进行坐标变换。当不考虑流体的影响时，计算所得的模态称为干模态，而考虑了周围流体影响的则称为湿模态。湿模态法是最常使用的一种计算流固耦合振动特性的方法，现有的文献中已大量报导了湿模态法的应用实例[9-15]。但是，实际应用过程中，此方法的使用却存在一个不可忽视的问题，即流体刚度矩阵的奇异性。王勖成[16]提出使用分离流体零频的方法来消除奇异性，方法本身是合理的，但本文将指出其使用过程中的矛盾，并提出湿模态法难以求解的根本原因及其解决方法。此外，本文将利用内燃机上常用的某个油底壳结构，计算不同盛油量对结构模态频率及振型的影响。

1 湿模态法

1.1 一般表达形式

对于流固耦合的结构，其流体方程及边界条件为

(1) 流体域内

(1)

式中：p为流体动压力；c为流体内声速。

对于不可压缩流体，c→∞，则有

2p=0

(2)

(2) 流固耦合面上

(3)

(3) 固定界面上

(4)

(4) 流体自由表面上

(5)

式中：z为重力方向；g为重力加速度。

为简化后续计算，当不考虑流体自由表面的影响时，式(5)将改写为

(6)

(5) 流体无限远处

(7)

式中：r为足够远处边界的法向；t为时间。

对于不可压缩流体，有

(8)

使用Galerkin法离散化后，直角坐标系内流场内任一点的压力分布p(x,y,z,t)可近似表示为

(9)

式中：l为流体单元节点数；N为流体单元形函数。

由有限元分析可知

∭ΩNi(2p*)dΩ=0

(10)

式中，Ω为流体单元体积。若令

(11)

(12)

式中，k为流体域总节点数。则式(10)可扩展为

∭ΩN(2p*)dΩ=0

(13)

对式(13)采用Green公式，可化为

(14)

将边界条件及式(9)代入后，式(14)变为

∭ΩN

(15)

式中，SI为流固耦合面积。

(16)

(17)

其中，

H=∭ΩNNTdΩ

(18)

因其表达式与结构刚度矩阵具有相同的形式，故而称为流体刚度矩阵。

(19)

式中，B为流固耦合矩阵。

另一方面，结构的运动方程表示为

(20)

式中，Ms、Cs、Ks分别为结构的质量矩阵、阻尼矩阵以及刚度矩阵；fp为流固耦合面上流体动压力产生的激励矢量；f0为其他外界激励矢量。

利用虚位移法可得

fp=-BTp

(21)

因此，对于无阻尼结构，流固耦合中结构的运动方程为

(22)

联立式(17)与式(22)，当考虑自由振动时，且q0=0，则有

(23)

求解式(23)，易得

(24)

(25)

其中，

Ma=ρBTH-1B

(26)

式中，Ma为流体附加质量。

利用一般多自由度线性系统的自由振动求解方法[17]，流固耦合中结构的模态与振型很容易通过式(25)获得。

1.2 湿模态中的问题

计算流体附加质量需要用到式(26)，其中需要计算H的逆矩阵。然而，形函数的性质决定了被离散物体本身的刚度矩阵必须是奇异的。因此，流体刚度矩阵H是不可逆的，必须将流体刚度矩阵的奇异性消除才能求解式(23)～式(26)。

1.3 奇异性的消除方法及其问题

王勖成采用分离流体零频的方法来消除H的奇异性，其中具体的公式推导本文不再复述。

(27)

式中，SF为流体自由表面面积。可见，Mf由两项构成，前一项由流体的可压缩性引起，后一项由自由液面的波动引起。而回到问题最初的假设，流体是不可压缩的且不考虑液面波动的影响，因此将假设后的边界条件代入上式可得Mf=0，故而m1=0，m1出现在分母上是不合理的。

所以，此方法的应用前提是流体的可压缩性和自由液面的影响不能同时忽略，即Mf≠0，此时式(23)必须改写为

(28)

其中式(28)中第二式为

(29)

然而此时，流体刚度矩阵H奇异与否已经不影响进一步的计算了。

2 求解湿模态的方法

通过“1”节已经了解到，当前提假设为流体不可压缩且不考虑自由液面波动的影响时，H为奇异矩阵，直接求解式(23)或是通过“1.3”节的分离零频的方法都是不可取的。因此湿模态法的根本问题在于为了简化计算而提出了错误假设。事实上，当同时假设了流体的不可压缩性，以及忽略自由液面的影响，流体域处于完全自由边界的条件下，在结构激励下会表现出类似结构刚体运动的响应形式。对此，本文将给出两种解决方法。

2.1 直接求解法

考虑流体的可压缩性以及自由液面波动的影响，则对于无阻尼系统，流固耦合方程应采用式(28)表达。采用复指数法[18]，令r=Rejωt，p=Pejωt并代入式(28)，则

(30)

求解式(30)的第二式可得

P=ρω2(H-ω2Mf)-1BR

(31)

代入式(30)的第一式得

[ω2Ms-Ks+ρω2BT(H-ω2Mf)-1B]R=0

(32)

令

A(ω)=ω2Ms-Ks-ρω2BT(H-ω2Mf)-1B

(33)

则式(32)有非零解的唯一条件是A(ω)的行列式为零，即

det[ω2Ms-Ks+ρω2BT(H-ω2Mf)-1B]=0

(34)

求解式(34)可求得所有的特征值及特征向量，进而获得模态频率及振型。

由于矩阵的行列式为零等价于至少存在一个特征值为零的特征向量，可令其为x。考虑以下问题

(35)

式(34)的解一定是这个问题的解。这个问题中没有行列式，只有矩阵向量的运算，是一个约束优化问题[19]，采用拉格朗日乘子法即可求解[20]。

考虑到A(ω)中包含了(H-ω2Mf)-1，此逆矩阵的阶数等于流体域的节点数，且ω作为未知参数被代入求逆矩阵，因此仅当流体域较小时，此方法较为适用。

2.2 不动点法

对于节点数较多的流体域，“2.1”节所采用的直接求解法就难以适用了。为了仍能通过式(23)～式(26)求解流固耦合方程，则可依然在流体不可压缩和忽略自由液面影响的前提假设下，于流体域内假定一个不动点，此点的流体动压力恒为零，等效于对流体进行了约束。实际处理过程中，可采用大数法，然后通过式(23)～式(26)即可求解耦合运动方程。

这个不动点最好选取在流固耦合面上且结构节点位移为零处，此处流体动压也为零，又由于静压对结构振动不产生影响，因此这样的节点最适合作为不动点。如果不存在这样的节点，则可在无限远处或自由表面选取一点，对于大型的、节点数较多的复杂流体域，由于约束的区域相对很小，因此同样可以满足精度要求。

3 仿真计算

采用Matlab编程的方法，对内燃机中常用的某油底壳进行仿真计算，如图1(a)所示。令最多盛油量时油液高度为h(h=69 mm)，则图1(b)～图1(e)分别代表油液高度为h/4、h/2、3h/4以及h时的状态。

图1 不同盛油量时的油底壳Fig.1 The oil pans with different oil mass

油底壳采用三角形壳单元建立，油液采用一阶四面体单元建立，具体参数如表1所示。

表1 计算模型参数

由于流体节点数较多，因此采用直接求解法较为困难，此处采用不动点法进行计算，约束节点为流体自由表面某一点。表2与图2为5种油底壳状态前5阶模态频率(去除刚体模态)对比。

表2 不同油液高度下前5阶模态频率

由图2可知，当盛油量逐渐增加时，油底壳的前5阶模态频率基本呈现先减小后增大的趋势，并趋于平缓，这与张亮等[21]中所得结论并不一致，可能是由于油底壳结构不同而导致的。由式(26)可知，附加质量矩阵取决于H及B，而由B的表达式可以看到对它产生影响的不止有流体、固体的形函数，还有坐标变换矩阵，它是由耦合面的大小、形状、方向等几何特征共同决定的。这也同样说明，流体越多，附加的质量越大，耦合结构模态频率越低的说法并不是一定成立的。

图2 不同油液高度下前5阶模态频率Fig.2 The first 5 order modal frequencies with different oil height

另一方面，从宏观上来说，模态频率取决于结构的质量与刚度，因此由图2也可以推断出，对于本文所使用的油底壳结构，当耦合系统中仅有少量流体时，流体的影响仅体现在附加的质量上，而当流体增多时，附加质量增加的同时，流体附加的刚度逐渐提高，因此导致了模态频率回升的现象。

为了进一步验证仿真计算的结论是可靠的，在图1的基础上将翻边上的螺栓孔约束，通过试验测取了油液高度为0、h/4、h/2、3h/4的油底壳约束模态，前10阶模态频率的对比，如图3所示。

由图3可知，随着油液高度的增加，油底壳的约束模态频率基本也呈现先减小后增大的趋势，与自由模态有较好的一致性。

图3 不同油液高度下前10阶试验模态频率Fig.3 The first 10 order test modal frequencies with different oil height

图4为图2中5种油底壳盛油状态前3阶结构振型。其中，各个节点的绝对振幅以结构形变表示，相对振幅的范围为0～1。

由图4可知，从绝对振幅上看，图4(d)～图4(o)所呈现的结构形变都远小于图4(a)～图4(c)；这表明，流体的加入在很大程度上抑制了结构模态的振幅。因此，结合表2与图2，当进行到后续的动力学及声学计算时，如果不将流体的影响计入，则将很有可能使油底壳的振-声峰值频率后移并放大振动及噪声，导致出现不合理的仿真计算结果。

纵向对比每一阶模态，除第1阶振型在各种状态下变化不大以外，第2阶、第3阶振型均在加入流体后产生了明显的变化，如图4(b)、图4(c)与图4(e)、图4(f)的对比所示。并且，随着盛油量的增加，振型趋于稳定，且较为接近不盛油状态时的干模态振型。

图4 5种油底壳状态前3阶结构振型Fig.4 The first 3 order structural modal shapes of the 5 oil pans

4 结 论

(1) 本文指出了湿模态法在应用时所遇到的问题，即流体刚度矩阵的奇异性。同时，指出了通过分离零频来消除流体刚度矩阵奇异性的方法中存在的矛盾。这些问题说明了，湿模态法不能被直接使用的根本原因是错误的前提假设，即同时忽略了流体的可压缩性以及流体自由表面的影响。

(2) 本文提出了两种可行的湿模态法的计算方法。其中，直接求解法使用的前提是流体的可压缩性和自由表面的影响不能被同时忽略，计算时可将行列式方程的问题转化为约束二次优化问题，并采用拉格朗日乘子法求解，但该方法仅适用于规模较小的流体域。另一种方法为不动点法，此方法仍在原先的前提假设下，选取流体中的一点，令其动压力恒为零，则可避免流体刚度矩阵无法求逆的问题并继续使用湿模态法求解。不动点法更适用于大型、复杂的流体域，以保证局部约束给整体刚度带来的影响较小。

(3) 本文采用不动点法对不同盛油量的某油底壳进行了模态计算。模态频率的对比表明，随着盛油量从0逐步增加，该油底壳的模态频率呈先减小后增大的趋势，并趋于平缓。通过与张亮等研究的结论对比，表明流体的增多并不一定使结构的模态频率单调递减，流体导致的模态频率上的变化还取决于流固耦合面的几何特征。对于该油底壳而言，当盛油量较少时，流体的影响以附加质量为主，随着盛油量逐渐增多，附加质量增加的同时，附加刚度的作用逐渐体现。

(4) 不同盛油状态时模态振型的对比表明，流体的加入很大程度上抑制了结构的模态振幅，结合频率的变化说明，在进行流固耦合系统的结构特性、动力学以及声学计算时有必要将流体的影响考虑在内，否则将导致峰值频率的后移以及振-声幅值的放大，以致不合理的计算结果。另外，纵向对比各阶模态可见，随着盛油量的增加，模态振型趋于稳定，并接近干模态振型。

[ 1 ] 陆鑫森. 高等结构动力学[M]. 上海: 上海交通大学出版社, 1992.

[ 2 ] 贾维新. 发动机结构噪声和进气噪声的数字化仿真及优化设计研究[D]. 杭州: 浙江大学, 2008.

[ 3 ] 郑运虎, 李颖. 立式圆柱薄壳容器的振动特性研究[J]. 西华大学学报(自然科学版)， 2016,35(1): 24-28.

ZHENG Yunhu, LI Ying. Research on the vibration characteristics of vertical cylindrical shell container [J]. Journal of Xihua University (Natural Science)， 2016,35(1): 24-28.

[ 4 ] 周勇, 车驰东. 水下航行器附连水质量的理论及实验研究[J]. 上海交通大学学报, 2016,50(2): 176-181.

ZHOU Yong, CHE Chidong. Theoretical and experimental study of added mass for underwater vehicles [J]. Journal of Shanghai Jiao Tong University, 2016,50(2): 176-181.

[ 5 ] 孙洋. 附加水质量对船体总振动固有频率影响的研究[D]. 大连: 大连理工大学, 2008.

[ 6 ] 柳瑞锋, 黄嵘, 周相荣, 等. 船体低阶湿模态计算方法对比研究[J]. 船舶工程， 2014(4): 25-28.

LIU Ruifeng, HUANG Rong, ZHOU Xiangrong, et al. Contrast study on calculation method for lower order wet mode of ship hull [J]. Ship Engineering, 2014(4): 25-28.

[ 7 ] 吴绍亮, 金咸定. 流固耦合计算方法在船舶局部结构中的应用[J]. 振动与冲击, 2003,22(4): 28-30.

WU Shaoliang, JIN Xianding. Computation method for fluid / structure interaction in local structures of ship [J]. Journal of Vibration and Shock, 2003,22(4): 28-30.

[ 8 ] 李东旭. 高等结构动力学[M]. 北京: 科学出版社, 2010.

[ 9 ] 黄晓明, 朱锡, 牟金磊, 等. 整体结构模型低阶湿模态仿真计算方法[J]. 舰船科学技术, 2011,33(5): 9-12.

HUANG Xiaoming, ZHU Xi, MU Jinlei, et al. Simulation and experimental investigation on transverse lower order wet mode of whole structure model [J]. Ship Science and Technology, 2011,33(5): 9-12.

[10] 王峥, 洪明, 刘城. 基于FEM/BEM的浸水结构振动及声辐射特性国内研究综述[J]. 船舶力学, 2014,18(11): 1397-1414.

WANG Zheng, HONG Ming, LIU Cheng. Domestic review of the submerged structure vibration and acoustic radiation characteristics based on FEM/BEM [J]. Journal of Ship Mechanics， 2014,18(11): 1397-1414.

[11] 武大江，梅志远，王永历. 充水密加筋夹层结构固有特性仿真及试验研究[J]. 振动与冲击, 2016,35(15): 134-139.

WU Dajiang, MEI Zhiyuan, WANG Yongli. Nature characteristics simulation and tests for water-filled multi-stiffened sandwich structures[J].Journal of Vibration and Shock,2016,35 (15):134-139.

[12] 吴健，李泽成，熊晨熙. 基于Abaqus的水下结构声辐射仿真方法[J]. 计算机辅助工程, 2015,24(6): 37-41.

WU Jian, LI Zecheng, XIONG Chenxi. Simulation method of sound radiation of underwater structure based on Abaqus [J]. Computer Aided Engineering， 2015,24(6): 37-41.

[13] LIU Ruijun, HAO Zhiyong, XU Wang, et al. A study of the influence of cooling water on the structural modes and vibro-acoustic characteristics of a gasoline engine[J]. Applied Acoustics, 2016, 104(1): 42-49.

[14] 陈冬冬，王辉，杨景玲，等. 油底壳流固耦合动力学特性分析[J]. 噪声与振动控制， 2014,34(6): 17-19.

CHEN Dongdong, WANG Hui, YANG Jingling, et al. Fluid-structure coupled model for modal analysis of an oil pan [J]. Noise and Vibration Control, 2014,34(6): 17-19.

[15] 冯威，袁兆成，刘伟哲. 用液固耦合方法研究柴油机油底壳辐射声场[J]. 内燃机学报, 2005,23(6): 536-540.

FENG Wei, YUAN Zhaocheng, LIU Weizhe. Study on radiated acoustic field of oil pan by liquid-solid coupled [J]. Transactions of CSICE, 2005,23(6): 536-540.

[16] 王勖成. 有限单元法[M]. 北京: 清华大学出版社, 2003.

[17] CLOUGH R W, PENZIEN J, GRIFFIN D S. Dynamics of structures[M]. [S.l.] ： Wiley, 2010.

[18] 程耀东. 机械振动学[M]. 杭州: 浙江大学出版社, 1988.

[19] 张永恒. 工程优化设计与MATLAB实现[M]. 北京: 清华大学出版社, 2011.

[20] 华东师范大学数学系. 数学分析[M]. 北京: 高等教育出版社, 2010.

[21] 张亮，袁兆成，黄震. 流固耦合有限元法用于油底壳模态计算[J]. 振动与冲击, 2003,22(4): 102-103.

ZHANG Liang, YUAN Zhaocheng, HUANG Zhen. Fluid-structure coupled finite element method applied in oil pan’s mode calculation [J]. Journal of Vibration and Shock, 2003,22(4): 102-103.

Theapplicationofanimprovedwetmodemethodinfluid-structureinteraction

LIURuijun,HAOZhiyong,ZHENGXu,TANJianglin

(CollegeofEnergyEngineering,ZhejiangUniversity,Hangzhou310027,China)

The singularity of the fluid stiffness matrix in the equations of motion was pointed out, which made the existing wet mode method unavailable. Meanwhile, the contradiction of the method isolating the zero-frequency item was indicated. The paper shows that the root of the problem of the wet mode method is the unreasonable presumption. The presumption simultaneously neglects the compressibility of the fluid as well as the influence of the free surface fluctuation which aims at the easier computation. Two solutions called the direct method and the fixed point method were put forward based on the above conclusion. The direct method is more suitable for the computation with a small fluid domain scale while the fixed point method is fit for a large one. In this paper, the structural characteristics of an oil pan with different oil mass were calculated by a Matlab program based on the fixed point method. The comparison demonstrates that the structural modal frequencies decrease firstly, then increase and finally become stable with the oil mass increasing from 0. The deformation magnitude of the modal shapes decrease considerably when the fluid is added. For the same order of mode, the wet modal shapes differ a lot from the dry modal shapes with less fluid addition except the 1st mode. However, the wet modal shapes tend to be stable and the same as the dry modal shapes with the continuous increase of oil mass.

fluid-structure interaction; wet mode; direct method; fixed point method

2016-05-17 修改稿收到日期： 2016-09-18

刘瑞骏 男，博士生，1990年生

郝志勇 男，教授，博士生导师，1955年生

TB122

A

10.13465/j.cnki.jvs.2017.22.031