基于完全耦合算法的绕水翼流固耦合特性研究

吴 钦，黄 彪，王国玉，曹树良

（1.清华大学 热能工程系，北京 100084；2.北京理工大学 机械与车辆学院，北京 100081）

基于完全耦合算法的绕水翼流固耦合特性研究

吴 钦1，黄 彪2，王国玉2，曹树良1

（1.清华大学 热能工程系，北京 100084；2.北京理工大学 机械与车辆学院，北京 100081）

基于完全耦合算法对绕二维NACA0009水翼流固耦合特性进行了数值模拟研究。采用Theodorsen模型和Munch模型对刚性和弹性水翼的水弹性响应进行了数值计算，分析了流体与结构的相互作用关系，研究了影响结构水弹性响应和流固耦合特性的因素。研究结果表明：考虑了流体黏性的Munch模型与基于势流理论的Theodorsen模型对气动弹性响应的数值计算结果基本一致，而Theodorsen模型由于没有考虑流体黏性在一定程度上低估了结构的水弹性响应。结构的惯性、阻尼和刚度力矩与流体的相应附加载荷均处于同一数量级，故流体与结构的相互作用不可忽略，尤其对于弹性水翼，流体的惯性、附加阻尼作用增大，流固耦合算法的数值稳定性对流固耦合特性的计算结果影响将更大。外部激励频率为非共振频率时，结构的刚度作用是影响水弹性响应的主要因素，外部激励频率为共振频率时，流体的附加阻尼和附加刚度作用减弱，除结构的刚度作用外，流体与结构的惯性作用对水弹性响应和流固耦合特性的影响也较大。

完全耦合算法；水弹性；流固耦合

0 引 言

流固耦合振动现象在工程和自然界中广泛存在，对这种现象的研究在机械、航空、航天、海洋、建筑和生物等领域都具有十分重要的意义。人们对于流固耦合振动现象的早期认识源于机翼及叶片的气动弹性问题[1-4]。在与水动力学相关的水力机械与船舶工程领域中，复杂水力环境下的水弹性响应会导致系统动力特性明显改变，并与流速、压力、湍动能、涡湍粘性等流动参数的变化紧密关联，随流态变化而异，这显著加剧了水动弹性问题流固耦合研究的理论难度[5-6]。随着对水力机械的安全稳定性要求日益提高，国内外学者对绕水翼流固耦合特性及其影响因素进行了大量的研究工作。Gowing[7]针对两种不同水弹性自适应复合材料水翼展开实验研究，结果表明叶片变形可以有效减小叶片攻角，从而减小叶尖载荷，推延叶尖初生空化。Young[8]建立了一种低阶基于势流理论的三维边界元方法用来求解空泡边界和压力波动，结合有限元软件确定叶片的动态响应，建立了边界元与有限元的流固耦合计算方法。Ausoni[9]应用数值计算与实验相结合的方法研究了绕水翼空化流动和流固耦合特性对流场涡结构产生机理的影响，结果表明，在翼型非共振条件下，初生空化数与雷诺数的平方根成正比，而在翼型共振条件下，涡脱落频率与结构的特征频率一致，且涡结构的空间相干性具有准二维性质。

流固耦合问题的重要特征是流体和固体结构之间的相互作用，包括固体结构在流体载荷作用下产生的变形或运动、固体的变形和运动对流场的影响，因此，流固耦合问题的求解需要同时考虑流场和结构场的求解及其耦合。耦合场的求解算法一般包括完全耦合法（或同步求解法）和分步求解法[10-14]。完全耦合法对流体和结构建立统一的耦合方程，在一个时间步内对流体域和固体域中所有的未知量进行求解；分步求解法分别对流场和结构场选择合适的数值算法进行独立求解，流固耦合界面数据通过反复迭代求解并在界面间反复传递直至获得收敛解。虽然完全耦合法需要自编程序实现，但由于其建立的流场与结构场强烈的相关性而作为解决流固耦合问题的有效方法被广泛应用。Ryzhakov等[15]利用完全耦合法基于拉格朗日体系建立了流场和结构场控制方程，对水流冲击弹性平板和注水气球弹性变形等算例进行了分析研究，结果表明这种方法能精确求解出流体和结构的响应，且数值收敛性较好。Michler[16]针对活塞与流体相互作用的一维模型问题，从算法稳定性、计算精度和计算效率等方面对比了不同流固耦合算法在流固耦合问题求解中的应用情况，结果表明完全耦合法无条件稳定且计算精确度相对较高。

目前，国内外在水动弹性力学的流固耦合数值计算研究中多基于理想有势流动，考虑振荡效应、流体流动等因素对物体振动特性的影响较少，确定动态参振质量、动态刚度和动态阻尼的理论方法尚不完善。本文基于完全耦合算法，利用自编程软件对二维NACA0009水翼模型在流固耦合作用下的水弹性响应进行了数值计算，通过对比不同流体介质和不同结构参数下的水弹性响应，分析了流体与结构的相互作用关系，并基于计算所得弹性系统产生的附加质量、附加阻尼和附加刚度，对绕水翼流固耦合特性进行了研究。

1 控制方程及其求解

1.1 基本控制方程

利用有限元对流场中的连续质量结构进行瞬态分析，结构动力学方程定义为：

其中：［Ms］、［Cs］和 ［Ks］分别为结构的质量矩阵、阻尼矩阵和刚度矩阵，｝分别为结构的位移、速度和加速度，FEX是结构的外部激振力，FHE是流体对结构的附加作用力。

计算采用单自由度的二维NACA0009水翼模型，考虑翼型仅在外部激励作用下绕弹性轴发生振荡运动，忽略翼型弦长方向和展向的变形。图1给出了计算模型示意图，图中U为无穷远来流速度，b为半弦长，弦长c=2b=0.1m，弹性轴位于翼型1/2弦长处，与结构的重心、压心重合，θ代表翼型的扭转变形角度。因此，（1）式可简化为：

其中：Iθ为结构关于弹性轴定义的转动惯量，Cθ为结构扭转变形的阻尼系数，Kθ为结构扭转刚度系数，θ、θ˙和θ¨分别为结构的扭转变形角度、角速度和角加速度，M0sinωe（ ）t是幅值为M0、频率为ωe的外部激励力矩，Mfluid是流体对结构的附加作用力矩。本文分别采用Theodorsen模型[17]和Munch模型[18]对流体的附加作用力矩Mfluid进行数值计算，分别记为

图1 单自由度二维NACA0009水翼模型Fig.1 NACA0009 hydrofoil with one degree offreedom

Theodorsen模型[17]基于不可压缩、无黏流动，流体对结构的附加作用力矩MTfluid表述为：

其中：MfT、CfT和KfT分别为根据Theodorsen模型确定的与结构变形加速度、速度和位移有关的附加水动力载荷，即结构的附加质量、附加阻尼和附加刚度，ρf为流体密度，k=ωeb/U为折合频率，Theodorsen函数C（）k是折合频率的复函数，表述为：

C（k）是由于考虑了自由涡的作用而引起的修正项，H21（k），H20（k）分别为第二类一阶和零阶Hankel函数[19-20]。

为了考虑流体黏性对流体作用力的影响，Munch[18]基于绕NACA0009水翼不可压缩、湍流黏性流动的大量数值计算与实验结果进行曲线拟合，提出Munch非定常动力模型，将流体对结构的附加作用力矩定义为：

1.2 流固耦合算法

本文采用完全耦合算法（Fully Coupled，以下简称FC）对流固耦合场进行数值求解。FC算法的优点是计算无条件收敛，计算稳定性和精确度较高，由于同时求解结构场和流场方程组，计算效率较高，适用于流体和结构之间具有强烈非线性特性的耦合问题求解。

图2给出了FC算法流程图。FC算法对流体和结构建立统一的耦合方程，同时对流场和结构场进行数值求解计算，根据解得的结构变形进行流场和结构场的网格更新，对下一离散时间步重复上述求解过程。

将（2）式在时域内进行离散可得：

基于Theodorsen模型和Munch模型，将（3）式和（5）式分别代入（6）式，其中分别定义为因此（6）式可分别表述为：

图2 完全耦合算法流程图Fig.2 Flowchart of the FC algorithm

其中：下标n表示时间增量，采用二阶Crank-Nicholson方法[21]对方程进行离散求解。

2 结果与讨论

本文针对NACA0009水翼模型，编程实现流固耦合控制方程的求解，对绕翼型流固耦合特性进行分析。表1给出了NACA0009刚/弹性水翼模型的物质属性，其中：ρs为结构密度，ξs为结构阻尼系数，ωθ为结构在真空中的特征频率，ωn为结构在水中的特征频率。采用的计算工况为来流速度5 m/s，雷诺数5×105。

表1 刚/弹性水翼模型物质属性Tab.1 Model parameters for the case of rigid and flexible hydrofoils subject to forced pitching motion

图3分别给出了在不同激励条件和流体介质中，数值计算的刚/弹性水翼俯仰角度随时间的变化情况，当ωe/ωn=0.01时，表征外部激励频率远小于共振频率，当ωe/ωn=1时外部激励等于共振频率。

对比图3（a）、（b）与图3（c）、（d）可知，当流体介质为空气时，采用Theodorsen模型和Munch模型数值计算得到的水翼俯仰角度随时间的变化情况基本一致；当流体介质为水时，采用Theodorsen模型和Munch模型对翼型水弹性响应的数值计算结果出现差异，且共振激励频率下不同模型的数值计算结果差异比非共振激励频率下更为明显。造成这种差异的原因在于，空气黏性较小，因此，考虑了流体黏性的Munch模型与基于势流理论的Theodorsen模型的数值计算结果基本一致，而水的黏性不可忽略，当流体介质为水时，Theodorsen模型由于不考虑流体黏性而在一定程度上低估了结构的弹性变形量，从而导致Theodorsen模型和Munch模型的数值计算结果差异较大。

图3 在不同激励条件和流体介质中，数值计算得到的刚性与弹性水翼俯仰角度随时间的变化Fig.3 Comparison of the predicted pitching motion obtained using Theodorsen model and Munch model for the steel/flexible hydrofoil in air/water atωe/ωn=0.01 andωe/ωn=1

对比图3（a）、（c）与图3（b）、（d）可以看出，当结构所受外部激励的频率为非共振频率时，刚性和弹性水翼的俯仰角度随时间的变化情况基本一致，均随外部激励发生周期性变化；当结构所受外部激励的频率为共振频率时，振幅明显大于非共振激励频率下的振幅，且刚性和弹性水翼发生共振时的振幅和频率相差较大，水翼在水中的弹性变形量呈明显发散趋势。这反映了当结构所受外部激励的频率接近系统固有频率时，系统振幅显著增大的共振现象，同时结构的材料属性对水弹性响应的振幅和频率影响较大。

图4给出了分别采用Theodorsen模型和Munch模型数值计算得到的刚/弹性水翼附加惯性力矩、阻尼力矩和刚度力矩随时间的演变情况。在本小节的研究中，设定流体介质为水，外部激励频率为非共振频率，即ωe/ωn=0.01。表2给出了采用不同模型数值计算得到的流固惯性力矩、阻尼力矩和刚度矩系数之比。

由图4（a）-（b）可知，对于相同结构，基于Theodorsen模型和Munch模型计算的惯性力矩随时间的演变规律相同，这是因为，由表2可知，基于Theodorsen模型和Munch模型计算所得附加惯性力矩系数（附加质量）相等，即结合（3）式和（5）式可知，Theodorsen模型与Munch模型对附加惯性力矩系数项的定义相同，附加惯性力矩系数仅与流体密度和结构几何参数有关，与折合频率、雷诺数等无关。同时，附加惯性力矩与结构惯性力矩随时间的演变呈同相分布，说明在水弹性响应的影响下，由流体与结构组成的弹性系统内部惯性作用加强。此外，由图4（a）-（b）不难发现，初始时刻，刚性水翼的附加惯性力矩小于结构的惯性力矩，而弹性水翼的附加惯性力矩大于结构的惯性力矩。结合表2可知，刚性和弹性水翼的流固惯性力矩比值分别为0.245和1.718，刚性水翼的密度较大，流固密度比较小，因此流体对刚性结构的惯性作用比流体对弹性结构的惯性作用小。

图4 刚/弹性水翼的结构惯性、阻尼和刚度力矩（黑色曲线）及流体附加惯性、附加阻尼和附加刚度力矩（红色曲线为基于Munch模型的数值计算结果，蓝色曲线为基于Theodorsen模型的数值计算结果）随时间的演变情况（ωe/ωn=0.01）Fig.4 Comparison of the predicted fluid to solid inertial moments,damping moments and stiffness moments for the steel/flexible hydrofoil in water atωe/ωn=0.01

由图4（c）-（d）可知，采用Theodorsen模型数值计算得到的附加阻尼力矩小于采用Munch模型计算得到的附加阻尼力矩。结合表2可知，采用不同模型数值计算得到的刚性水翼流固阻尼力矩比值分别为0.186（Theodorsen模型）和8.914（Munch模型），而对于弹性水翼，数值计算得到的流固阻尼力矩比值则分别为0.905（Theodorsen模型）和18.429（Munch模型）。这是因为Theodorsen模型基于无黏势流假设，低估了流体的阻尼作用力。此外，采用Munch模型计算的附加阻尼力矩明显大于结构的阻尼力矩，两者同相分布，说明水弹性响应进一步增强了弹性系统的阻尼作用，其中流体附加阻尼作用强于结构阻尼作用。

由图4（e）-（f）可知，采用Theodorsen模型和Munch模型计算的刚度矩随时间的演变情况基本一致，结合表2可知，采用不同模型计算得到的流固刚度矩比值也很接近。这是由于结构的刚度系数较大，流固刚度矩比值较小，因此，流体的附加刚度作用对水弹性响应影响较小。与惯性力矩和阻尼力矩不同的是，流体附加刚度矩与结构刚度矩呈反相分布，这说明水弹性响应在一定程度上削弱了弹性系统的刚度作用。

综合图4（a）-（f）可以看出，当外部激励频率为非共振频率时，结构的刚度力矩远大于结构的惯性力矩、阻尼力矩以及流体的相应附加载荷，因此结构的刚度作用是影响结构水弹性响应的主要因素。同时，结构的惯性、阻尼和刚度力矩与流体的相应附加载荷均处于同一数量级，因此流体与结构的相互作用不可忽略，尤其对于弹性水翼，由于流体的惯性、附加阻尼作用增大，流固耦合算法的数值稳定性将对弹性水翼流固耦合特性的数值计算结果产生更大的影响。

表2 采用不同模型数值计算得到的流固惯性力矩、阻尼力矩和刚度力矩系数之比Tab.2 Comparison of the fluid-solid ratio of the inertial,damping and stiffness moments coefficients

为了进一步研究振荡频率对结构水弹性响应和流固耦合特性的影响，图5对比了外部激励频率为共振频率（ωe/ωn=1）时，采用Theodorsen模型和Munch模型计算得到的刚/弹性水翼附加惯性力矩、阻尼力矩和刚度力矩随时间的演变情况。

图5 刚/弹性水翼的惯性力矩、阻尼力矩和刚度矩（黑色曲线）及流体附加惯性、附加阻尼和附加刚度力矩（红色曲线为基于Munch模型的数值计算结果，蓝色曲线为基于Theodorsen模型的数值计算结果）随时间的演变情况（ωe/ωn=1）Fig.5 Comparison of the predicted fluid to solid inertial moments,damping moments and stiffness moments for the steel/flexible hydrofoil in water atωe/ωn=1

由于Theodorsen模型和Munch模型对附加惯性力矩系数定义相同，由图5（a）—（b）可以看出，不同模型数值计算得到的附加惯性力矩随时间的演变情况相同，且由于刚性水翼的密度远大于流体的密度，而弹性水翼的密度与流体的密度相当，因此刚性水翼的附加惯性力矩小于结构的惯性力矩，而弹性水翼的附加惯性力矩大于结构的惯性力矩。由于附加惯性力矩系数与折合频率无关，由表2可知，外部激励频率为共振频率时，数值计算得到的刚/弹性水翼流固惯性力矩系数之比分别为0.245和1.718，与外部激励频率为非共振频率时的计算结果一致。值得注意的是，基于无黏势流假设的Theodorsen模型是针对任意变形的翼型简谐运动推导结果，当激励频率较高，结构变形量较大时，采用Theodorsen模型的数值计算结果可靠性大大降低。由图5（c）—（f）不难发现，基于Munch模型的数值计算结果表明，结构阻尼力矩小于流体附加阻尼力矩，且二者同相分步；结构刚度矩大于流体附加刚度矩，且二者呈反相分布，这与外部激励频率为非共振频率时数值计算得到的水弹性响应趋势是一致的。由表2可知，与外部激励频率为非共振频率时的数值计算结果相比，共振激励频率下，采用Munch模型计算所得的流固阻尼力矩系数比和刚度距系数比明显减小，流体的附加阻尼和附加刚度作用减弱。此外，结合图3（d）可知，当结构所受外部激励频率为共振频率时，结构的变形量呈现发散趋势，流体的附加载荷也随时间的发展逐渐增大。综合图5（a）—（f）可知，当外部激励频率为共振频率时，弹性水翼的惯性、刚度力矩与流体附加惯性力矩大小相当，说明除结构的刚度作用外，流体与结构的惯性作用对结构水弹性响应和流固耦合特性的影响也较大。

3 结 论

本文基于完全耦合算法对二维NACA0009水翼模型进行了数值模拟研究，分析了刚性和弹性水翼在流固耦合作用下的水弹性响应，并对影响绕水翼流固耦合特性的因素进行了研究。主要结论如下：

（1）由于空气的黏度较小，采用Theodorsen模型和Munch模型对气动弹性响应的数值计算结果基本一致；而水的黏性不可忽略，不考虑流体黏性的Theodorsen模型在一定程度上低估了结构在水中的变形量。

（2）结构的惯性、阻尼和刚度力矩与流体的相应附加载荷均处于同一数量级，因此流体与结构的相互作用不可忽略，尤其对于弹性水翼，流体的附加惯性、附加阻尼作用增大，流固耦合算法的数值稳定性对弹性水翼流固耦合特性的数值计算结果影响也将增大。

（3）当外部激励频率为非共振频率时，结构的刚度作用是影响结构水弹性响应的主要因素，而当外部激励频率为共振频率时，流体的附加阻尼和附加刚度作用减弱，除结构的刚度作用外，流体与结构的惯性作用对结构水弹性响应和流固耦合特性的影响也较大。

[1]Todd O N,Thomas W S.Aeroelastic response ofa rigid wing supported by nonlinear springs[J].Journal of Aircraft,1998, 35(4)：616-622.

[2]So R M C,Jadic I,Mignolet M P.Fluid-structure resonance produced by oncoming alternating vortices[J].Journal of Fluids and Structures,1999,13(4)：519-548.

[3]Kamakoti R,Shyy W.Fluid-structure interaction for aeroelastic applications[J].Progress in Aerospace Sciences,2005,40：535-558.

[4]崔 鹏,韩景龙.新型运输机机翼的颤振特性分析[J].振动工程学报,2011,24(2)：192-198.Cui Peng,Han Jinglong.Flutter analysis of new transport-type wings[J].Journal of Vibration Engineering,2011,24(2)：192-198.

[5]叶永林,吴有生,邹明松,倪启军.基于水弹性力学的SWATH船结构振动与噪声分析[J].船舶力学,2013,17(4)：430-438.Ye Yonglin,Wu Yousheng,et al.Analysis of the structural vibration and noise radiation of a SWATH ship based on hydroelastic method[J].Journal of Ship Mechanics,2013,17(4)：430-438.

[6]Chae E J.Dynamic response and stability of flexible hydrofoils in incompressible and viscous flow[D].University of Michigan,2015.

[7]Gowing S,Coffin P,Dai C.Hydrofoil cavitation improvements with elastically coupled composite materials[C]//Proceeding of 25th American Towing Tank Conference,Iowa City,USA,1998.

[8]Young Y L,Kinnas S A.Numerical modeling of supercavitating propeller flows[J].Journal of Ship Research,2003,47(1)：48-62.

[9]Ausoni P,Escaler X,Avellan F.Cavitation influenced on von karman vortex shedding and induced hydrofoil vibrations[J]. ASME Journal of Fluids Engineering,2007,129：966-973.

[10]Hubner B,Walhorn E,Dinkler D.A monolithic approach to fluid-structure interaction using space-time finite elements [J].Computer Methods in Applied Mechanics and Engineering,2004,193：2087-2104.

[11]Farhat C,Zee K vander,Geuzaine Ph.Provably second-order time-accurate loosely-coupled solution algorithms for transient nonlinear aeroelasticity[J].Computer Methods in Appllied Mechanical Engineering,2006,195(17)：1973-2001.

[12]Wu Q,Huang B,Wang G.Experimental and numerical investigation of hydroelastic response of a flexible hydrofoil in cavitating flow[J].International Journalof Multiphase Flow,2015,74：19-33.

[13]Matthies H G,Steindorf J.Partitioned strong coupling algorithms for fluid-structure interaction[J].Computers and Structures,2003,81：805-812.

[14]Dowell E H,Hall K C.Modeling of fluid-structure interaction[J].Annual Review of Fluid Mechanics,2001,33：445-490.

[15]Ryzhakov P B,Rossi R,Idelsohn S R,Oñate E.A monolithic Lagrangian approach for fluid-structure interaction problems[J].Computational Mechanics,2010,46(6)：883-899.

[16]Michler C,Hulshoff S J,Van Brummelen E H,De Borst R.A monolithic approach to fluid-structure interaction[J].Computers&Fluids,2004,33(5)：839-848.

[17]Theodorsen T H.Generaltheory ofaerodynamic instability and the mechanism offlutter[R].NACA Rep.496,1935.

[18]Munch C,Ausoni P,Braun O,Farhat M,Avellan F.Fluid-structure coupling for an oscillating hydrofoil[J].Journal of Fluids and Structures,2010,26：1018-1033.

[19]Abramowitz M,Stegun I A.Handbook of mathematical functions[D].National Bureau of Standards,Applied Math,Series #55,Dover Publications,1965.

[20]Sears W R.Some aspects of non-stationary airfoil theory and its practical applications[J].Journal of the Aeronautical Sciences,1941,8：104-108.

[21]Crank J,Nicolson P.A practical method for numerical evaluations of partial differential equations of the heat-conduction type[J].Mathematical Proceedings ofthe Cambridge Philosophical Society,Cambridge Press,1947,43(01)：50-67.

Fluid structure interaction analysis of a hydrofoil based on fully coupled algorithm

WU Qin1,HUANG Biao2,WANG Guo-yu2,CAO Shu-liang1
(1.Department of Thermal Engineering,Tsinghua University,Beijing 100084,China;2.School of Mechanical Engineering,Beijing Institute of Technology,Beijing 100081,China)

A fully coupled algorithm for modeling of fluid structure interaction in viscous flow around the NACA0009hydrofoil is presented.The numerical simulations are performed by using the Theodorsen model and the Munch model to investigate the fluid structure interaction and characterize the factors that significantly affect the hydroelastic responses.It is revealed that Theodorsen’s approximation of the hydroelastic moment is based on inviscid,potential flow theory and hence it underestimated the hydroelastic responses to some extent.The fluid inertial,damping,and stiffness forces are not negligible compared to their structural counterparts.For the flexible hydrofoil,the fluid inertial and damping forces increase,suggesting that the numerical instability of fluid structure coupling algorithms may lead to a more significant impact on the numerical prediction of the fluid structure interactions.Meanwhile,the structural stiffness dominates the response of the hydrofoil when the external excitation frequency is set to be the non-resonance frequency,while as the external excitation frequency equals to the resonance frequency,the effect of the fluid damp-ing and stiffness is less significant and the structural response of this case is dominated by the fluid and solid inertial effects in addition to the structural stiffness.

fully coupled algorithm;hydroelastic;fluid structure interaction

TV131.2

A

10.3969/j.issn.1007-7294.2017.07.002

1007-7294（2017）07-0804-10

2017-01-03

国家自然科学基金资助项目（51679005）；北京市自然科学基金资助项目（3172029）；高等学校博士学科点专项科研基金资助项目（20131101120014）

吴 钦（1989-），女，博士，助理研究员，E-mail：wuqin919@163.com；

黄 彪（1985-），男，副教授，硕士生导师。