基于波有限元法的流固耦合结构波传导问题

倪广健， 林杰威

(天津大学 内燃机燃烧学国家重点实验室，天津　300072)



基于波有限元法的流固耦合结构波传导问题

倪广健， 林杰威

(天津大学 内燃机燃烧学国家重点实验室，天津300072)

摘要：采用波有限元方法研究流固耦合结构中的波传导问题。该方法以有限元法为基础，首先建立研究对象的有限元法模型，得到模型的动态刚度矩阵。通过对动态刚度矩阵的重新排列组合得到研究对象的传递矩阵，求解传递矩阵的特征值问题可以得到分别代表自由波传递的波数和波模。该研究首先分析独立流体结构和固体结构中的振动问题，并比较了采用波有限元法和理论方法求解得到的固体结构中波数分布情况，验证了模型的正确性。随后采用波有限元法分析流固耦合结构中的波传导问题。波有限元法的应用并不局限于所给出的均匀或周期性结构，还可将其应用于缓慢变化的非均匀结构。

关键词：波传导；流固耦合；有限元；波有限元

流固耦合结构在工程领域中很常见并且扮演着重要的角色。一般来讲，当流体结构与固体结构接触发生相互耦合作用时，耦合面的流场会受到固体结构位移的影响，反过来固体结构的振动情况也受到流体压强作用而产生变化。采用传统的理论分析来研究复杂系统特别是流固耦合系统是非常困难的，但是数值方法(例如有限元法)则为研究人员提供了良好的研究手段。有限元法将研究对象所在的连续求解域离散为一系列单元，采用适当的形函数来近似代表每个单元在求解域上待求的未知场函数。形函数通常由未知场函数及其导数在单元各节点的数值插值函数来表示[1]。这样一来，一个连续的无限自由度问题就可以采用离散的有限自由度问题来近似表示。如果待求解问题所涉及的最短波长内含有至少6个单元，则可以认为该有限元模型可以近似代表实际的连续系统[2]。

有限元法的优势在于模拟复杂结构，但是传统的有限元法只可以计算得到系统的整体响应，其中的详细物理意义则无从所知。波方法则可以从波传导和波衰减的角度解释系统的动态特性，为系统的动态响应赋予相应的物理含义。

对于均匀波导结构来讲，其几何和材料属性在其长度上均保持不变，类似的结构包括产生轴向振动的均匀直杆，产生弯曲振动的均匀直梁及具有恒定曲率的弯曲梁等等。波方法(Wave approach)在分析弹性结构，尤其是均匀结构的动态特性中应用广泛。Mace[3]根据波的反射、透射和传播情况研究了均匀直梁的振动情况。Miller等[4]研究了一维波导(杆和梁)中的能量流问题，指出两个方向相反的衰减波界面处的能量可以进行传递。Easwaran等[5]研究了传递矩阵和阻抗矩阵之间的关系，分别分析了具有对称系统、交互系统和保守系统的矩阵属性。Zhong等[6]根据传递矩阵的辛特性建立了一套高效、精确的特征值求解方法。传递矩阵的特征向量描述了结构中波运动情况，特征值则描述了波在通过波导时，幅值和相位的变化。Mead[7]总结了连续周期性结构中的波传导情况和一些传递矩阵的应用实例。Ichchou等[8]采用波有限元法研究了流固耦合结构中的损伤对波传导的影响。Ni等[9]采用波有限元法研究了薄锥壳结构中的波传导问题，并利用求解得到的波特征向量对传统有限元法计算得到的结果进行解耦分析，给出了不同波对壳结构振动的贡献程度。Ni等[10-11]还采用波有限元法对人耳耳蜗中的波传导情况进行了分析，得到了波数及各个波贡献度在空间域上的分部情况。根据其模型对称性假设，他们使用的结构为只涉及单层的流固耦合。

本文将采用相同的有限元和波方法相结合的方式来研究频域范围内流固耦合系统内的波传导情况。

1流固耦合系统自由振动

1.1流体系统

为了研究流体中的声波振动，本文采用8节点6面体声学有限元单元，如图1所示。该单元的长、宽、高为2a1×2a2×2a3， (ξ1,ξ2,ξ3)代表三个方向上的无量纲坐标。

图1　8节点6面体声学单元Fig.1 8-node hexahedral acoustic element

(1)

(2)

式中：ρf为流体的平均体积密度，κ为流体的可压缩性。

对汉密尔顿原理中的积分式关于时间求导，可以得到以声压p来表示的流体运动方程[13]。如果流体系统的某一个面(r面)为流固耦合系统的交界面，则第r面上的由声压p和与该面正交的虚位移δur产生的耗散功可以表示为

δWd,r=∫SrpδurdSr

(3)

式中：ξr1,2,3代表ξr1、ξr2和ξr3，ξr1和ξr2为第r耗散面尺寸，面积Sr=ar1×ar2。ξr3代表耗散面的位置。作用在流体上的虚功可以表示为

δWq=∫VpδqdV

(4)

式中：q(ξ1,2,3)为流体单位体积下的体位移分布律。

1.2固体系统

本文所涉及的固体结构为图2所示的矩形薄固体结构。该单元由4个节点组成，每个节点包含3个自由度。节点的各自由度可以表示为w、θx=∂w/∂y和θy=-∂w/∂x，分表代表垂直平板平面的位移w、绕x轴的转角θx和绕y轴的转角θy。其中w以沿z轴正向为正，转角按右手螺旋法则用矢量表示，矢量以沿x,y轴正方向为正。单元尺寸为2a×2b，(ξ,η,ζ)代表局部坐标系，ξ=x/a，η=y/b。

图2　四节点矩形单元Fig.2 4-node rectangular element

该薄板单元的动能和势能可以表示为[2]

(5)

(6)

式中：Iz=h3/12为单位宽度上的截面二阶矩，h为板厚度，ρs为材料密度，E为材料杨氏模量，χ为曲率向量χ(ξ,η,t)，D为薄板弯曲弹性系数矩阵。

如果存在耗散因素或有外力作用，则相对应的虚功可以表示为

(7)

δWf=∫ApzδwdA

(8)

式中：μ为单位面积上的黏性阻尼系数，pz(ξ,η,t)为作用在单元表面上的单位面积的横向分布力，δw(ξ,η,t)为竖直方向虚位移(虚挠度)。

1.3流固耦合系统

流固耦合系统的振动特性可以根据汉密尔顿原理推导而得。在汉密尔顿原理中，对于下层的流体系统，如图3所示，该积分必须包含固体结构作用在流体上的虚功δWar。第r个耦合面上由声压-pz,r和虚位移δwr所产生的虚功可以写为

δWar=-∫Srpz,rδwrdSr

(9)

式中：ξr1,2,3代表ξr1、ξr2和ξr3，ξr1和ξr2为第r耦合面尺寸，且面积Sr=ar1×ar2。ξr3代表耦合面的位置。对于固体结构，该哈密尔顿积分也必须包含流体作用在该板上的虚功(与流体系统的虚功互补)δWsr。此虚功由声压pr虚位移δur产生，表示为

δWsr=∫AprδwrdSr

(10)

式中：ξr1=ξ，ξr2=η，ar1=a及ar2=b。

此处所涉及的运动方程均只考虑自由振动，即不存在耗散力和声源体位移引起的虚功。根据前面的分析，下层流体系统的动能、势能、薄板对流体的虚功，可作为流体系统的总能量和总的虚功。这样针对流体系统的汉密尔顿积分应满足

(11)

式(11)除了应该满足流体系统边界条件外，当时间t=t1及t=t2时，仍需满足δφ(t)=0。

类似流体系统，薄板系统的总能量和总虚功可以用薄板系统的动能、势能、流体作用在薄板上的虚功来表示，如下

(12)

图3　流固耦合结构示意图Fig.3 Schematic diagram of a fluid-structure coupled system

根据式(11)和(12)可推导得到图3所示的流固耦合耦合系统的运动方程：

(13)

式中：下标“t”代表上层流体，“b”代表下层流体。本文为简化表示，假设上下层流体尺寸和材料均相同。模型所涉及的材料参数和几何尺寸详见表1。

表1　模型几何及材料参数

图4　流固耦合系统固有频率及振型Fig.4 Natural frequencies and vibration modes of the fluid-structure coupled system

本文先计算了耦合系统的低阶模态频率和振型，用于解释流固耦合作用的一些基本特性。如图4所示，耦合系统低阶的固有频率和模态振型是由固体结构主导的，耦合系统中的固体结构模态振型与独立的固体结构相同。受固体结构主导的振型，其固有频率要低于单独的固体结构在该阶振型下对应的固有频率，这是由于流体的介入，相当于给固体结构增加了质量，因此其固有频率会下降。

2流固耦合系统内的自由波

2.1波有限元法

波有限元法首先基于传统有限元法建立起所研究的波导的一个分段，获得其动态刚度矩阵并进行矩阵变换，得到该分段的传递矩阵，通过求解传递矩阵的特征值问题获得分别代表波经过该段波导所产生的幅值和相位变化的波数值(特征值)，以及相对应的波模(特征向量)。当考虑所研究的流固耦合系统上长度为Δ的一个分段时，该分段的运动方程可以表示为：

Dq=f

(14)

式中：D为该分段动态刚度矩阵，q为自由度向量，f为力向量。式(14)的矩阵形式为

(15)

式中：下标“L”和“R”表示分段的左侧和右侧。波有限元方法就是利用式(15)中的元素来建立特征值问题方程。对于均匀波导结构，下列关系成立：

(16)

式中：“·T”代表矩阵转置。针对每个分段，传递矩阵可以定义为[14]

(17)

文献[14]中所述的周期性条件表示了该分段上的位移以及力的关联关系

(18)

式中：λ表示通过该分段波的幅值和相位变化。因此，通过长度为Δ的流固耦合系统分段的自由波运动情况可以用如下的特征值问题表示出来。

(19)

式中传递矩阵T可以采用式(15)中动态刚度矩阵的各个元素表示为

(20)

假设波导分段一侧所具有的自由度数量为n，那么整个传递矩阵T的大小则为2n×2n。这样根据式(19)就可以得到n对特征值和特征向量。式(19)中的特征值λ与长为Δ的波导分段中的波相关，可以表示为

λj=e-ikjΔ(j=1,2…,n)

(21)

图5显示了板结构中的波数分布情况。可以看出，采用波有限元法计算得到的结果和理论结果[17]非常接近。波有限元方法的计算误差可以概括为有限元离散误差和惯性项的舍入误差。当波在长为Δ的波导中传播时产生的相位变化变大时，有限元的离散误差也随之变大。这是因为有限元模型是系统的一种近似模拟，不可避免的会产生数值误差。通常情况下，每个波长上至少要包含6个或更多单元才能精确表示系统的运动[1]。

图5　板结构波数分布的理论结果(实线)和波有限元结果(虚线)Fig.5 Wavenumber distribution in the thin plate structure calculated analytically (solid line) and numerically using the wave finite element method (dashed line)

2.2流固耦合系统内的波传导

根据前述的流固耦合系统运动方程，可以利用波有限元方法对如图3所示的流固耦合系统进行分析，研究其内波的传递情况。流体系统除了与固体结构耦合的面其余表面均为声学刚性，固体结构在y=0和y=W处的边界条件分别为简支和固支。该分段上共含有70个自由度，因此会产生70个不同的波(正向35个，反向35个)。在此仅列出系统中四个显著的正向传播的波的波数分布情况，如图6所示。

图6　流固耦合结构中的波数分布情况。其中实线代表波数的实部，虚线代表波数的虚部。Fig.6 Wavenumberdistribution in the fluid-structure coupled system, in which real parts are denoted by solid line and imaginary parts are denoted by dashed line

为了解释各个波的含义，接下来将分析每个波所携带的动能和势能，及波传递的功率流。长为Δ的固体结构的平均动能和势能可以表示为[15, 18]

(22)

(23)

Es=Ek,s+Ep,s

(24)

式中：λj和wj分别为第j个波相关的特征值和位移向量，K和M为结构的刚度矩阵和质量矩阵，Es固体结构单位长度上的总能量密度。类似的，流体中的平均动能、势能和功率流可以表示为[19]

(25)

(26)

Ef=Ek,f+Ep,f

(27)

式中：pj代表第j个波相关的压强，ρf流体密度，H和Q为刚度和质量矩阵，Ef流体单位长度上的总能量密度。固体结构[15]和流体[19]的功率流可以表示为

(28)

式中：fs和ff分别代表波有限元法计算中固体结构和流体的内力。

由于篇幅限制，本文仅列出不同波中的功率流计算结果。图7可以看到，波1与梁结构中的弯曲波相似。该波能够在耦合结构中很好的传播。该波波数的虚部在各研究频率下均为0，而且在该波中，当频率高于约550 Hz时，固体结构的振动处于主导地位，且固体结构的总能量密度与流体总能量密度的比值总是大于1，这就是为什么该波与梁结构中的弯曲波相似。在势能方面，流体的势能大于固体结构所具有的势能。由于图7中流体和固体结构的功率流皆为正，所以该波为传播方向为正向。

图7　波1中固体和流体的归一化功率流分布，固体功率流(实线)，流体功率流(虚线)。Fig.7 Normalized power flow in solid (solid line) and fluid (dashed line) associated with wave 1

如图8所示，波2的截断频率出现在1 500 Hz 左右，低于该频率，波2无法传播。当频率高于该阶段频率是，波2开始传播，流体在波中占主导地位，流体的动能快速衰减(此时该波为衰减波)，高于该截断频率，波2开始传播。当频率高于截断频率时，固体结构开始在波2中占据主导地位。由于功率流为正，所以波2的传播方向同样是正向。

图8　波2中固体和流体的归一化功率流分布，固体功率流(实线)，流体功率流(虚线)。Fig.8 Normalized power flow in solid (solid line) and fluid (dashed line) associated with wave 2

如图9所示，波3与波2相类似，其处截断频率为4 000 Hz左右。低于该截断频率，该波无法传播。高于该截断频率，固体结构开始在波中占主导地位，且波3开始正向传播。

图9　波3中固体和流体的归一化功率流分布，固体功率流(实线)，流体功率流(虚线)。Fig.9 Normalized power flow in solid (solid line) and fluid (dashed line) associated with wave 3

波4为一个纯粹的声波(纵波)，其波数为一个实数，该波中固体结构不含任何能量，因此该波中固体的功率流始终为0。

组速度是描述波的一个重要参数，实际计算中可以采用能量法[15, 18]，有限差分法[21]或特征值微分法[22]来求解波的组速度。本文采用能量法来计算上述4个的组速度。在能量法中，组速度可以表示为功率和能量的比值[18]

(29)

式中：Pj和Ej为第j波所含的平均功率和能量。图10显示了图6中4个波各自的组速度分布情况。

图10　不同波的组速度分布Fig.10 Group velocity of each wave

波4的组速度为一个在整个分析频率范围内为一个恒定值，大小约为1 500 m/s，等于水中声波的传播速度。波1的组速度随频率的增加而增加，波3和波4的组速度在截断频率之下为0，即没有传播，当频率高于截断频率时，这两个波开始传播，这和两个波的功率流分布趋势相符。

3结论

本文利用有限元法对流固耦合系统中的自由振动和波进行了研究，阐述了耦合作用对二者固有频率和振型的影响。随后采用波有限元方法对流固耦合系统中的波传导情况进行了分析，包括不同波所的功率流，以及在该波传播过程中固体系统和流体系统对其的贡献度等等。

本文采用了相同的波有限元方法来研究频域范围内流固耦合系统内的波传导情况，但耦合结构相比Ni等[10-11]所建立的模型更加复杂。本文所提出的模型中，固体结构沉浸在两层流体之内，更接近于耳蜗的实际结构[12]。在此基础上，本文所建立的模型可以进一步推广应用于耳蜗建模或其它流固耦合领域。

利用波有限元方法研究耦合结构中的波传导情况是一个新的尝试。在无法采用解析方法分析复杂耦合结构中的波传导时，波有限元可以作为一个很好的选择来进行数值计算、模拟，了解结构中的波传导及各个波所携带的能量，这对耦合系统的振动控制非常重要。

参 考 文 献

[ 1 ] Petyt M. Introduction to finite element vibration analysis[M].Cambridge: Cambridge University Press, 1990.

[ 2 ] Fahy F, Gardonio P. Sound and structural vibration: Radiation, transmission and response[M]. 2nd ed.Oxford, UK: Elsevier Academic Press, 2007.

[ 3 ] Mace B R. Wave reflection and transmission in beams[J]. Journal of Sound and Vibration,1984, 97:237-246.

[ 4 ] Miller D W, Von Flotow A. A travelling wave approach to power flow in structural networks[J]. Journal of Sound and Vibration, 1989,128: 145-162.

[ 5 ] Easwaran V, Gupta V H, Munjal M L. Relationship between the impedance matrix and the transfer matrix with specific referenceto symmetrical reciprocal and conservative systems[J]. Journal of Sound and Vibration, 1993,161:515-525.

[ 6 ] Zhong W X, Williams F W. On the direct solution of wave propagation for repetitive structures[J]. Journal of Sound and Vibration, 1995,181: 485-501.

[ 7 ] Mead D M.Wave propagation in continous periodic structures: research contributions from southampton[J].Journal of Sound and Vibration, 1996,190: 495-524.

[ 8 ] Ichchou M N, Mencik J M, Zhou W. Wave finite elements for low and mid-frequency description of coupled structures with damage[J]. Comput Meth Appl Mech Eng, 2009,198: 1311-1326.

[ 9 ] Ni G, Elliott S J. Wave interpretation of numerical results for the vibration in thin conical shells[J]. Journal of Sound and Vibration, 2014,333: 2750-2758.

[10] Elliott S J, Ni G, Mace B R, et al.A wave finite element analysis of the passive cochlea[J].The Journal of the Acoustical Society of America, 2013,133: 1535-1545.

[11] Ni G, Elliott S J.Wave finite element analysis of an active cochlear model[J]. The Journal of the Acoustical Society of America, 2013,133: 3428.

[12] De Boer E. Mechanics of the cochlea: modelling efforts[M]//Dallos P, Popper A N, Fay R R(Eds.) The cochlea, Springer, New York, 1996:258-317.

[13] Craggs A. The transient response of a coupled plate-acoustic system using plate and acoustic finite elements[J]. Journal of Sound and Vibration, 1971,15: 509-528.

[14] Brillouin L. Wave propagation in periodic structures[M]. Mineola, NY:Dover Publications, 2003.

[15] Mace B R, Duhamel D, Brennan M J, et al. Finite element prediction of wave motion in structural waveguides[J]. The Journal of the Acoustical Society of America, 2005,117: 2835-2843.

[16] Duhamel D, Mace B R, Brennan M J. Finite element analysis of the vibrations of waveguides and periodic structures[J].Journal of Sound and Vibration, 2006,294: 205-220.

[17] Graff K G. Wave motion in elastic solids[M]. London:Oxford University Press, 1991.

[18] Cremer L, Heckl M, Petersson B A T. Structure-borne sound: structural vibrations and sound radiation at audio frequencies[M]. Berlin:Springer, 2005.

[19] Maess M, Herrmann J, Gaul L.Finite element analysis of guided waves in fluid-filled corrugated pipes[J]. The Journal of the Acoustical Society of America, 2007,121: 1313-1323.

[20] Biot M A. General theorems on the equivalence of group velocity and energy transport[J]. Physical Review, 1957,105: 1129.

[21] Stroud K A, Booth D J. Advanced engineering mathematics[M]. New York:Palgrave MacMillan, 2003.

[22] Finnveden S. Evaluation of modal density and group velocity by a finite element method[J]. Journal of Sound and Vibration, 2004,273: 51-75.

Wave propagation in a fluid-structural coupled system based on wave finite element method

NIGuang-jian,LINJie-wei

(State Key Laboratory of Engines, Tianjin University, Tianjin 300072, China)

Abstract:Wave propagation in a fluid-structural coupled system was presented here using a combined finite element method and wave approach (the wave finite element method). The method was based on the finite element descriptions of the model. The eigenvalue problem of the model transfer matrix derived from the dynamic stiffness matrix of the model was solved to give eigenvalues and eigenvectors, they determined free wave propagation. Free vibration problems of the separated structure and independent fluid structure were analyzed firstly. Wavenumber distributions in a plate strip were calculated using both the wave finite element method and the analytic method. The results showed that the model and the method are accurate. Numerical examples of wave propagation in a fluid-structural coupled system was then presented. Moreover, it was shown that the application of the wave finite element method is not limited to uniform or periodic structure, but can be extended to non-uniform structures with slowly varying properties.

Key words:wave propagation; fluid-structural coupling; finite element; wave finite element

中图分类号:O327

文献标志码：A

DOI:10.13465/j.cnki.jvs.2016.04.033

通信作者林杰威 男，博士后，1984年生

收稿日期：2014-09-12修改稿收到日期：2014-12-12

基金项目：高等学校博士学科点专项科研基金(20130032130005)

第一作者 倪广健 男，博士后，1981年生