基于面光滑有限元法的三维声传播问题研究

尤小健，李 威，柴应彬

(1.武汉第二船舶设计研究所,湖北 武汉 430064； 2.华中科技大学 船舶与海洋工程学院，湖北 武汉 430074； 3.船舶和海洋水动力湖北省重点实验室，湖北 武汉 430074；4.高新船舶与深海开发装备协同创新中心，上海 200240)

基于面光滑有限元法的三维声传播问题研究

尤小健1，李 威2,3,4，柴应彬2

(1.武汉第二船舶设计研究所,湖北 武汉 430064； 2.华中科技大学 船舶与海洋工程学院，湖北 武汉 430074； 3.船舶和海洋水动力湖北省重点实验室，湖北 武汉 430074；4.高新船舶与深海开发装备协同创新中心，上海 200240)

为解决标准有限元方法“过刚”的特性，采用一种基于面的光滑有限元法(FS-FEM)来求解三维声学问题，在FS-FEM模型当中，应用四节点的四面体单元来离散三维声场域，通过与四面体单元面相关的光滑域并运用伽辽金弱形式来建立线性离散系统的控制方程，在面光滑域上进行的梯度光滑操作能够“软化”原来“过刚”的FEM模型，因此能够减少数值离散误差。数值仿真结果表明：运用相同的网格密度，FS-FEM能够得到比FEM更加精确的计算结果。

面光滑有限元法；梯度光滑技术；声学；数值方法

0 引 言

充液管道在机械、船舶和石油化工等许多工程领域有着极为广泛的应用背景，因此对充液管道内声传播及声辐射特性问题的研究具有重要的学术价值和工程意义，国内外许多学者在这个研究领域做了大量的工作[1-2]。

对于一些简单的结构，其声学问题的解析解是可以得到的。但是对于现实当中大多数的复杂工程结构，往往没有解析解。此时通常用各种各样的数值方法来得到其近似解，其中有限元方法和边界元方法是求解声学问题当中发展的最成熟，应用最广泛的数值方法[3-5]。在运用有限元方法进行声场数值模拟时，由于有限元模型总体刚度与系统的真实刚度相比总是表现的“偏刚”，这就不可避免地就会出现数值计算误差[6]，而且这种数值计算误差会随着计算频率的提高而不断变大。为了尽量减少这种数值误差，提高声学数值仿真的计算精度，通常的做法就是不断加大有限元网格的密度，从而使有限元模型的刚度更加接近模型的真实刚度。但是随着网格密度的增加，计算成本也大大增加，尤其对于复杂的三维声学问题，这个问题就显得更加突出。

为了解决上述有限元模型表现的“过刚”的问题，国内外许多学者在这方面做了大量的工作。近年来，Liu等[7-10]提出的光滑有限元方法(SFEM)是将由Chen等[11-12]提出的无网格分区应力光滑技术与传统标准的有限元结合在一起的一种新型数值方法。该方法通过对分区光滑域的光滑操作，将原来“过刚”的有限元计算模型进行适当“软化”，使其能够更加接近计算模型的真实刚度，因而能够得到更加精确的数值计算结果。本文应用一种基于面的光滑有限元方法(FS-FEM)计算三维充液管道内部的声传播问题，计算结果表明，与标准有限元法相比，FS-FEM具有更高的计算精度和更好的收敛性，具有良好的工程应用前景。

1 声学光滑有限元计算模型

在理想流体介质中，忽略流体粘性的影响，简谐小振幅振动声波满足以下的波动方程：

(1)

式中：▽2为Laplace算子；p为声场中任意一点的声压；c为声速；t为时间。

在声场当中任意一点的质点振动速度和该点的声压有以下关系：

▽p+jρωv=0，

(2)

式中：ρ为流体介质的密度；ω为圆频率；v为质点的振动速度。

在声场计算中，考虑3种声学边界条件：Dirichlet边界条件，Neumann边界条件和Robin边界条件，这3类控制方程的表述如下[9,13]：

p=pD，

(3)

▽p·n=-jρωvn，

(4)

▽p·n=-jρωAnp。

(5)

式中vn和An分别为质点振动的法向速度和边界上的导纳系数。

采用加权值法并应用上述边界条件，可得式(1)的Galerkin弱积分形式[9]:

(6)

式中Γ1和Γ2分别为式(4)给定的速度边界条件和式(5)给定的声阻抗边界条件。声压变量由式(7)得到：

(7)

式中：Ni为形函数；pi为节点声压。

运用广义梯度光滑技术对声场当中的光滑域进行分片光滑操作，可得如下光滑的Galerkin弱积分形式：

(8)

由式(8)可得到如下的声学有限元计算模型[9]：

(9)

式中：

(10)

(11)

为光滑刚度矩阵；

(12)

(13)

上述光滑刚度矩阵不同于标准有限元方法中的刚度矩阵，它通过在每个光滑域上而不是在每个单元内部进行数值积分得到。这种操作是对原来“过刚”的有限元模型的一种修正。

2 基于面的光滑域构造

对任意一个三维问题，每个光滑域的构造是基于每个单元的面，构造过程是：首先将三维问题域离散为4节点的四面体单元，然后将相邻四面体的中点及每个面上的3个顶点顺序相连，这就形成了基于面的光滑域。因此，整个问题域当中有多少个面，就会形成多少个基于面的光滑域[9,14]，具体如图1所示。

图1 基于四面体网格的光滑域Fig.1 The smoothing domain based on tetrahedral mesh

在运用广义梯度光滑技术对原有的梯度场进行光滑操作时实际上是对每一个光滑域上的速度梯度场进行光滑，光滑后的速度由式(14)得到：

(14)

式中Vk为基于面k的光滑域的体积。将式(7)代入式(14)可得：

(15)

式中n为与每个光滑域相关的节点数。

对于三维问题

(16)

(17)

式中：Vk为每个光滑域的体积；Ns为每个光滑域边界被分成的份数；Ng为每份积分域当中高斯点的个数；w为与之相应的权系数。

进一步可得每个光滑域Ωk的单元刚度矩阵：

(18)

和标准有限元类似，将每个基于面的光滑域上的单元刚度矩阵总体组装，就可得到系统的总体刚度矩阵。

3 数值计算结果

三维充液管道计算模型如图3所示，管道的长度L=1m，半径R=0.05m，管道内部区域Ω充满密度ρ=1 000kg/m3, 声速c=1 500m/s的流体介质。管道的边界条件为：在左端面施加一个vn=10sinωtm/s的速度边界条件，管道的右端和管壁为刚性边界条件。

图2 三维充液管道模型Fig.2 The three dimensional tube model filled with liquid

3.1 模态分析

为了验证本文所提出计算方法的有效性与可行性，首先对图3所示的计算模型进行模态分析，计算过程中不考虑阻尼的影响。计算网格如图3所示，由于对称性，图3只给出了1/4的网格模型，该网格模型包含894个节点和3 141个单元，单元的平均尺寸为h=0.025m。根据声学网格单元尺寸所要满足的“拇指法则”，即每个计算频率所对应的波长至少应包含6个单元[9,17]，该网格模型下可分析的上限频率为10 000Hz。表1给出了FEM和FS-FEM在相同网格模型下前10阶模态的计算结果，以及解析解结果。

图3 三维充液管道1/4网格模型Fig.3 A quarter mesh of the three dimensional tube model filled with liquid

从表1可看出：在低阶模态时，FEM和FS-FEM都可得到较为理想的计算结果，但与FEM相比，FS-FEM的计算精度更高。随着模态阶数的增高，FEM和FS-FEM的计算误差都逐渐增大，但是，FS-FEM的计算结果仍然要优于FEM。这说明通过前面所述的“梯度光滑技术”得到的FS-FEM计算模型有效地“软化”了原来“过刚”的FEM计算模型，因而能够得到更加精确的计算结果。

表1 三维充液管道模态计算结果Tab.1 The results of the three dimensional tube eigenfrequencies

3.2 声传播问题

分别运用FEM和FS-FEM在前述的边界条件下对图2所示的管道模型的声传播问题进行计算和分析，计算网格仍为图3所示的网格模型。为了考察声传播频率对计算结果的影响，这里计算了f=2 000Hz，f=4 000Hz和f=7 000Hz等3种不同频率下管道内的声压幅值分布情况，结果如图4～图6所示，另外该声场问题的解析解为[9]：

(19)

(20)

为了便于比较，图中也给出了解析解的计算结果。

从图4～图6可看出，在低频区域(1 000Hz)，FEM和FS-FEM的计算结果都与解析解吻合得很好；在高频区域，FEM和FS-FEM结果的误差都增加，但FS-FEM的结果要优于FEM，这就进一步验证了FS-FEM的精确性和有效性。

图4 声传播频率f=2 000 Hz声压幅值分布图Fig.4 Acoustic pressure distribution along the tube with frequency f=2 000 Hz

图5 声传播频率f=4 000 Hz声压幅值分布图Fig.5 Acoustic pressure distribution along the tube with frequency f=4 000 Hz

图6 声传播频率f=7 000 Hz声压幅值分布图Fig.6 Acoustic pressure distribution along the tube with frequency f=7 000 Hz

3.3 收敛性分析

众所周知，运用有限元等数值方法求解声学问题时不可避免地会出现数值误差，数值误差越小，数值方法的收敛性就越好。为了验证FS-FEM的收敛性，本文给出了声学数值计算结果的全局部误差因子[15-16]：

(16)

3.3.1 单元尺寸对全局误差的影响

网格划分的质量对数值方法的计算结果有着重要的影响，一般情况下，网格单元的平均尺寸越小，网格越精细，计算结果的误差就越小。这里给出了4种不同的网格，单元平均尺寸分别为h=0.025m，h=0.02m，h=0.015m和h=0.012m。计算频率f=4 000Hz，该计算频率满足上述4种声学网格所决定的上限频率，图7给出了FEM和FS-FEM在不同单元尺寸网格下的全局误差。从图7可看出：随着单元平均尺寸的减小，网格密度的增大，FEM和FS-FEM的全局误差因子都迅速减小，可见这2种方法都是收敛的，但在同一单元尺寸下，FS-FEM的全局误差要小于FEM，这说明FS-FEM的收敛性要优于FEM。

图7 全局误差与单元平均尺寸之间的关系Fig.7 The global error against the mesh size

3.3.2 计算频率对全局误差的影响

在用数值方法求解声学问题时，计算频率的影响也是至关重要的，通常在低频区域，很多数值方法都可得到较为理想的数值结果，随着频率的提高，数值方法的误差就会变大，有时甚至不能得到可靠的结果。本小节主要讨论计算频率对全局误差及收敛性的影响。网格计算模型如图4所示，单元平均尺寸为0.025m，为满足声学网格中每个波长至少要包含6个单元的要求，这里的分析频率为0～9 000Hz，计算结果如图8所示。

图8 全局误差与声传播频率之间的关系Fig.8 The global error against the frequency

从图8可看出，在低频区域，FEM和FS-FEM的全局误差都较小，随着频率的提高，FEM和FS-FEM的全局误差都逐渐增大，但在整个计算频率范围内FS-FEM的全局误差要小于FEM，这进一步说明了FS-FEM的计算精度更高，收敛性更好。另外从图中可看出，全局误差在系统的各阶固有频率附近会急剧增大，说明在系统的固有频率附近，FEM和FS-FEM的计算误差都较大，这与先前其他学者所得的结果相吻合[15-16]。出现这种现象的原因是：在求解式(9)时，当声波的传播频率和系统的固有频率相近时，导致式(9)的系数矩阵特性接近于一个奇异阵，这样在计算过程当中任何微小的数值误差都会导致求解声压时出现较大的误差。

4 结 语

本文采用基于面的光滑有限元法(FS-FEM)对三维充液管道内的声传播问题进行研究，通过对计算结果的分析和比较可以得到如下的结论：

1)在三维充液管道的模态预报中，无论是低阶模态还是高阶模态，FS-FEM都能得到比FEM更加精确的结果。

2)在声场计算中，当声传播频率较低时，FEM和FS-FEM都可以得到较为理想的结果，但随着声传播频率的提高，FS-FEM的计算精度要明显高于FEM。

3)随着网格单元平均尺寸的减小，FEM和FS-FEM的全局误差都迅速减小，但FS-FEM的全局误差始终要小于FEM,表明FS-FEM在收敛性方面要优于FEM。

[1] 刘敬喜，李天匀，刘土光，等.弹性介质中充液管道的波衰减特性[J].华中科技大学学报(自然科学版),2003,31(10):90-92.

[2]PER-ANDERSH,GÖRANS.Dynamicfiniteelementanalysisoffluid-filledpipes[J].ComputerMethodsinAppliedMechanicsandEngineering,2001,190:3111-3120.

[3]CHENGL,ROBERTD,WHITEKG.Three-dimensionalviscousfiniteelementformulationforacousticfluid-structureinteraction[J].ComputerMethodsinAppliedMechanicsandEngineering,2008,197:60-72.

[4]CHENGH,CHENJ,ZHANGYB,etal.Amulti-domainboundaryelementformulationforacousticfrequencysensitivityanalysis[J].EngineeringAnalysiswithBoundaryElements,2009,33:815-821.

[5]WUHJ,LIUYJ,JIANGWK.Alow-frequencyfastmultipoleboundaryelementmethodbasedonanalyticalintegrationofthehypersingularintegralfor3Dacousticproblems[J].EngineeringAnalysiswithBoundaryElements,2013,37:309-318.

[6]BOUILLARDPH,SULEAUS.Element-freegarlekinsolutionsforhelmholtzproblems:formulationandnumericalassessmentofthepollutioneffect[J].ComputerMethodsinAppliedMechanicsandEngineering,1998,162:117-135.

[7]LIUGR,NGUYENTT,LAMKY.Anedge-basedsmoothedfiniteelementmethod(ES-FEM)forstaticfree,andforcedvibrationanalysis[J].JournalofSoundandVibration,2009,320:1100-1130.

[8]CUIXY,LIUGR,LIGY,etal.Analysisofplatesandshellsusingedge-basedsmoothedfiniteelementmethod[J].ComputationalMechanics,2010,45:141-156.

[9]HEZC,LIUGR,ZHONGZH,etal.Anedge-basedsmoothedfiniteelementmethod(ES-FEM)foranalyzingthree-dimensionalacousticproblems[J].ComputerMethodsinAppliedMechanicsandEngineering，2009,199:20-33.

[10]HEZC,LIUGR,ZHONGZH,etal.AcoupledES-FEM/BEMmethodforfluid-structureinteractionproblems[J].EngineeringAnalysiswithBoundaryElements,2011,35:140-147.

[11]CHENJS,WUCT,YOONS,etal.Astabilizedconformingnodalintegrationorgalerkinmeshfreemethods[J].InternationalJournalforNumericalMethodsinEngineering,2001,50:435-466.

[12]YOOJW,MORANB,CHENJS.Stabilizedconformingnodalintegrationinthenatural-elementmethod[J].InternationalJournalforNumericalMethodsinEngineering,2004,60:861-890.

[13] 何智成，李光耀，成艾国，等.光滑有限元的声学研究：时域和频域分析[J].振动与冲击,2012,31(16):122-127.

[14]HEZC,LIUGR,ZHONGZH,etal.Acouplededge-/face-basedsmoothedfiniteelementmethodforstructural-acousticproblems[J].AppliedAcoustics,2010,71:1955-1964.

[15] 夏百战，于德介，姚凌云.光滑有限元的声学研究：二维多流体域耦合声场的光滑有限元解法[J].声学学报,2012,37(6):601-609.

[16]HEZC,LIUGR,ZHONGZH,etal.Dispersionfreeanalysisofacousticproblemsusingthealphafiniteelementmethod[J].ComputationalMechanics，2010,46:867-881.

[17]HEZC,LIUGR,ZHONGZH,etal.Coupledanalysisof3Dstructural-acousticproblemsusingtheedge-basedsmoothedfiniteelementmethod/finiteelementmethod[J].FiniteElementsinAnalysisandDesign，2010，46：1114-1121.

Analysis of sound of propagation in three-dimensional fluid-filled pipes based on face-based on smoothed finite element method (FS-FEM)

YOU Xiao-jian1，LI Wei2,3,4，CHAI Ying-bin2

(1.Wuhan Second Ship Design and Research Institute,Wuhan 430064,China； 2.School of Naval Architecture and Ocean Engineering, Huazhong University of Science and Technology,Wuhan 430074,China； 3.Hubei Key Laboratory of Naval Architecture & Ocean Engineering Hydrodynamics (HUST) ,Wuhan 430074,China； 4.Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration (CISSE) ,Shanghai 200240,China)

To solve the “overly-stiff” property of the standard FEM, an face-based smoothed finite element method (FS-FEM) was adopted for analyzing three-dimensional acoustic problems.In the FS-FEM model, four-node tetrahedron elements are used to discretize the three-dimensional acoustic domain, the discretized linear system equations are established using the smoothed Galergin weak form with smoothing domains associated the surface of the tetrahedrons.The gradient smoothing operations over the faced-based smoothing domains provides a proper softening effect to the “overly-stiff” FEM model and hence can reduce the numerical dispersion error.The numerical simulation results show that the FS-FEM model can achieve more accurate results than the FEM model with the same mesh owing to the softened stiffness to the acoustic model.

face-based smoothed finite element method(FS-FEM);gradient smoothing technique;acoustics;numerical method

2015-01-04；

2015-02-02

尤小健(1975-)，男，高级工程师，研究方向为舰船总体设计。

O42

A

1672-7649(2015)12-0104-06

10.3404/j.issn.1672-7649.2015.12.021