跨介质航行器弹性舵翼空化流固耦合仿真分析

杜特专，黄仁芳，王 畅

(1.中国科学院力学研究所流固耦合系统力学重点实验室，北京 100190;2.中国科学院大学工程科学学院，北京 100049)

0 引言

跨介质航行器高速入水过程中，在与水接触的短时间内航行器头部产生短时高幅值冲击波，冲击波在航行器内传播，诱导航行器头部塑性变形。入水之后航行器与水介质相互作用，诱导与航行器侧面接触的水介质空化演化与脱落，航行器在头部流体阻力及水介质空化演化的共同作用下航行。与此同时，航行器在时空分布的水动力载荷作用下，结构发生变形，并产生长时的振动效应。航行器入水过程中速度较高，翼面、舵面等关键结构表面能够产生空泡。空泡的非定常演化及溃灭会在结构表面形成复杂的压力脉动，促使结构发生显著的动态响应。同时，航行器在穿越自由表面过程、浸水特性的变化会显著影响结构自身的动态特性，与通常的机翼气动弹性问题具有很大的区别，给耦合分析带来很大的难度。

随着跨介质高速航行器轻量化和高速的需求，空化流固耦合问题逐渐成为重要问题。针对包含空化的流固耦合问题，已有研究主要集中于柔性水翼和复合材料螺旋桨等对象。Young[1-2]利用边界元方法模拟空化流动，结合有限元分析结构的变形，模拟了水轮机叶片在空化和无空化条件下的耦合效应，并利用实验进行了验证，分析了复合材料螺旋桨的性能。Amromin等[3]利用梁模型求解水翼振动，分析了附着空化条件下水翼弹性振动的特征。Wu等[4-6]利用实验和数值模拟研究了柔性NACA66水翼的空化流激振动幅值与不同空化类型的关系。Ducoin等[7]运用商业流体力学计算软件CFX 建立了2自由度的动网格模型来模拟柔性水翼的变形，分析了柔性水翼动力学的规律。任怀勋[8]基于商业软件ANSYS及脚本控制对水翼空化问题进行了流固耦合计算，分析了三维水翼在空化流动作用下的流体动力和结构载荷特性。Campbell等[9]采用开源软件OpenFOAM对大变形的柔性叶片进行数值模拟，与实验结果吻合良好。已有研究工作为空化流固耦合提供了有效的研究方法，加深了研究人员对于空化流动作用下结构振动特征的认识，为进一步深入研究空化流固耦合效应奠定了良好的基础。

当前对于超空泡流动条件下结构与水汽界面的相互作用的研究还较少，翼/舵与水气界面的撞击产生的冲击力可能激起较大的结构振动，给结构安全性带来了很大的挑战。本文基于径向基函数方法建立了流固耦合计算方法，对航行器舵翼在超空化流动下的动态响应问题进行数值仿真，分析姿态角(攻角和侧滑角)对空泡形态以及舵翼的振动响应。

1 计算模型

流固耦合仿真根据求解方法可以分为3种：一是以固体求解器为主，流动过程采用简化模型处理；二是以流体求解器为主，结构振动采用简化的模型进行处理；三是流体和结构分别采用不同的求解器，利用数据接口达成耦合计算。当前对于空化流固耦合问题的分析主要以后两种方法为主。本文采用第2类方法，以流体求解器为主建立松耦合的计算模型。带空化的流固耦合计算模型主要包括流体计算模型、结构模型和流固耦合接口，下面将分别介绍这3个模型(接口)及其控制方程。

1.1 流体模型

流动包含了多相界面、水汽之间的相变以及湍流流动等问题，具体控制方程如下

(1)

(2)

(3)

(4)

ρ=αvρv+(1-αv)ρl

(5)

μ=αvμv+(1-αv)μl

(6)

式中，下标l和v分别代表液态水和水蒸气。

本文采用基于有限体积法的SIMPLE方法求解流体方程组，时间采用二阶有界中心差分格式，压力、动量和能量的插值格式分别采用PRESTO、有界中心差分格式、二阶迎风格式。

1.2 结构模型

本文关注的舵翼截面为三角型，底边长3.2 mm，高18 mm，展长40 mm。舵翼的材料为304不锈钢，弹性模量E=194 GPa，密度为7 930 kg/m3。在结构计算中，将舵翼简化成悬臂梁，并采用模态叠加法进行求解。悬臂梁固有频率ωi和对应的质量归一化的模态φi(x)如下

(7)

φi(x)=cosβix-coshβix-ϑi(sinβix-sinhβix)

(8)

(9)

其中，β1l=1.875,β2l=4.694,β3l=7.855。结构振动方程采用4阶龙格-库塔方法进行求解。

1.3 耦合接口及动网格方法

流场和结构的耦合求解是通过流固界面上的数据交换来实现的，但是在流固界面上流场网格和结构网格并不是一一对应的，这就需要通过插值来传递位移和受力的信息。本文的插值方法基于径向基函数(RBF)[10]，其插值的一般形式为

首先，培养学生勤于思考的习惯。爱因斯坦说“学习知识要善于思考、思考、再思考，我就是靠这个学习方法成为科学家的”，孔子也有名言“学而不思则罔”，所以在学习中思考发挥着至关重要的作用。思考的过程是对新学知识的吸收、消化、升华的过程。因此，提高阅读的效果，就需要培养学生勤于思考的阅读习惯。

(10)

(11)

(12)

式(10)具有N+4个未知系数，将N个结构节点的位置和位移代入式(10)可以建立N个方程组，为了使该方程组有解，还需补充4个方程

(13)

通过这N+4个方程可以求出系数γi和αi，则得到流体节点的位移为

Δxa=Aasλ

(14)

其中,

(15)

(16)

前面提到，流场网格节点数通常较多，对于复杂工程问题，把所有耦合界面上的流体节点全都取为控制点将会造成计算量过大。对此，可采用贪婪算法减少物面插值基点数目，从而提高计算效率。该方法在空化流固耦合计算中能够得到较好的结果[11]。

1.4 几何模型及网格划分

本文研究对象的主体是圆台和圆柱的组合体，在主体尾部对称分布两个舵翼(图1)。

图1 航行器外型及舵翼几何外型Fig.1 The geometry of vehicle and the rudder wing

这里不考虑航行器的耦合运动，只分析均匀来流条件下舵翼的流激振动。航行器主体和Z轴负向的舵翼设置为刚体，Z轴正向舵翼设置为弹性体。计算域网格采用分块结构化网格划分(图2、图3),能够更好地捕捉水气界面的形状和演化。近壁面边界层网格保证y+≈30。

图2 航行器网格划分Fig.2 The grids for the vehicle

图3 舵翼网格划分Fig.3 The grids around the rudder wing

2 计算结果

本文模拟了不同来流攻角和侧滑角，用以分析不同航行状态下空泡与三角型截面翼的相互作用。其中攻角在XOY平面内，侧滑角在ZOX平面内，具体计算工况如表1所示。

表1 计算工况Tab.1 Calculation conditions

从计算结果可以看出，在来流攻角为2°～6°时，航行器及舵翼大部分都被包在超空泡内部，仅有锥段末端与水直接接触(图4)。随着攻角的增大，航行器的沾湿面积也不断增大，但是舵翼仍然处于空泡内。结构振动的模态位移如图5所示，由于在给定的参数范围内舵翼始终被包裹在空泡内部振动，因此相当于在空气中振动，受力较小，振动位移也较小，其振动频率基本上与干模态基本相同，约为700 Hz。随着攻角的增大结构振动的平衡位置也不断增大(图5)。

图4 空泡形态(2°攻角) Fig.4 The supercavity shape (angle of attack is 2°)

(a)不同攻角下的振动时间曲线

(b)2°攻角下舵翼振动频谱图5 舵翼1阶振动特征 (工况1～3)Fig.5 The first-order vibration displacement of rudder wing (cases 1～3)

当来流攻角为8°时，由于超空泡向Y轴正向偏移，舵翼的自由端(翼梢)会穿出空泡壁，切割空泡界面，并在尾流区引起了空泡界面的褶皱(图6红色箭头处)。此时，舵翼翼梢的压力面为沾湿面，吸力面仍然处于空泡之中，舵翼在液态水和空泡内水蒸气的共同作用下振动，其振动特征发生了较大的变化。从图7(a)中可以看出，舵翼在切割自由面时受到较大的弯矩，产生较大的模态位移①。随后舵翼被空泡完全包围，受力减小，振动以为也随之降低，直到第2次切割空泡界面，产生第2个峰值②。经过两个周期后，舵翼一直处于切割空泡界面的状态③，在附加质量作用下其振动频率小于干模态的700 Hz(图7(b))。 图8给出了相同条件下弹性舵翼和刚性舵翼的受力对比，可以看到在t=0.17 s时刻刚性舵翼的受力有1个峰值，而弹性舵翼则峰值较小。这是由于弹性翼振动变形对空化形态造成影响，导致舵翼沾湿面积较小，从而降低其所受水动力。

图6 舵翼切割空泡(8°攻角)Fig. 6 The supercavity around rudder wing(angle of attack is 8°)

(a)8°攻角下的振动时间曲线

(b)8°攻角下舵翼振动频谱图7 舵翼1阶振动特征(工况4)Fig.7 The first-order vibration displacement of rudder wing(case 4)

图8 8°攻角条件下的刚体和弹性体的受力曲线Fig.8 The lift curves of the rigid wing and the elastic wing under 8° angle of attack

图9 舵翼一阶振动特征(工况5)Fig.9 The first-order vibration displacement of rudder wing(case 5)

对于攻角和侧滑角均为2°的航行工况，舵翼的沾湿面积相对于工况4来说更小，仅有翼梢顶端一小部分切割空泡壁。每个振动的峰值对应着舵翼切割空泡壁的过程，而低谷则表明舵翼完全处于泡内。该工况条件下，结构振动始终无法达到一个稳定的状态(图9)。这说明在该条件下结构的受力变化最为剧烈，比小攻角条件下的受力大了1个量级，可能对航行的稳定性带来较大隐患。

3 结论

本文在FLUENT的基础上进行二次开发，基于径向基函数实现了空化流激振动的耦合数据传递及网格变形，利用梁模型模拟舵翼的振动，建立了空化流固耦合计算方法。利用该方法模拟了包含航行器影响的舵翼流激振动过程，获得了不同攻角条件下空泡形态和结构振动特征。计算结果表明，在计算空化数条件下，当攻角较小时(不大于6°)，舵翼完全处于空泡内部，其振动与空气中的振动特征基本一致。当攻角为8°时，舵翼的自由端切割空泡界面会导致舵翼一开始的振动幅值比小攻角条件大1个量级，随后在一个较高的平衡位置上稳定振动。对于攻角和偏航角均为2°时，其振动无法达到稳定状态，这是因为舵翼的自由端始终在切割自由面和完全被气泡包围两种状态中切换。