基于FEM－SBFEM的水下结构声固耦合分析方法

李上明

(中国工程物理研究院总体工程研究所 四川绵阳 621999)

浸没在无限水域中的水下结构声固耦合问题广泛存在于工程领域，如水下潜艇声振分析、海洋平台抗震分析、海底管道抗震分析、坝体抗震分析等。该类问题的难点之一为准确高效模拟无限水域。模拟无限水域的方法通常包括 FEM［1］、传递边界法(TBC)［2－6］、边界元 (BEM)［7］等。其中，FEM 与TBC方法主要思想为将无限水域截断成有限水域和吸收边界条件(即阻抗边界)。这些方法主要面临离散区域大、自由度多及只能近似模拟无限水域无反射特性等不足。其中，FEM模拟无限域时，其离散域大、单元多，而传递边界法虽然在离散域及单元数量上小于FEM，但其公式复杂，公式普适性相对较弱，特殊问题特殊处理。而BEM则在无限水域的边界上进行离散，无限水域的无反射特性利用BEM的基本解进行精确表述，具有离散域小，精度高等特点，是一种能有效模拟无限水域的方法，但其奇异积分严重影响了其计算效率。

作为模拟无限域的一种替代方法，比例边界有限元法［8－9］(SBFEM)是一种半解析力学数值模拟方法，在模拟无限域时展现出离散网格少、精度高、不含边界元法中出现的奇异积分、无需基本解等特点，已逐渐用来模拟无限声学水域［10］。

针对浸没在无限水域中的水下结构，文献［10］利用FEM模拟水下结构，SBFEM模拟整个无限水域，水下结构与无限水域的耦合作用通过交错迭代算法进行求解。在整个求解过程中，水下结构和无限水域分别利用不同的公式进行交错迭代求解，当当前时刻计算结果满足某一收敛条件时，则进行下一时刻求解。该求解方法通常成为弱耦合迭代算法。利用弱耦合迭代算法及FEM－SBFEM耦合算法，文献［11］、文献［12］实现了坝库耦合瞬态分析及冲击波作用下的水下结构瞬态分析。其中，FEM模拟水下结构或坝体，FEM－SBFEM耦合模拟无限水域。数值结果表明弱耦合迭代算法有效获取了水下结构的动力响应，但因交错迭代导致了更多的时间消耗，并且收敛条件的定义对计算效率有较大影响。此外，交错迭代算法难以应用于频域求解，如模态及谐响应等问题的求解。

本文将基于FEM和SBFEM理论，建立水下结构FEM公式与无限水域SBFEM公式的耦合公式，形成水下结构与无限水域耦合的FEM－SBFEM强耦合公式，并对整个声固耦合系统用同一控制方程(FEM－SBFEM强耦合方程)进行求解。

1 含无限水域的声固耦合系统

图1为浸没在无限水域中的水下结构声固耦合系统。在模拟该系统的数值方法中，水下结构通常用有限元(FEM)模拟，无限水域通过无限水域截断面，截断成有限水域和剩余的无限水域。有限水域用FEM模拟，而剩余的无限水域用比例边界有限元(SBFEM)模拟。水下结构为满足连续介质力学理论的连接介质，而无限水域为满足压力波动方程的声学水域。详细结构及声学水域控制方程及其有限元理论可参考文献［13］。

图1 含无限水域的声固耦合系统Fig.1 Acoustic－solid coupling system involving infinite water

2 水下结构与有限水域耦合方程

在商业有限元软件中，有限水域与水下结构的声固耦合算法相对成熟，常采用直接耦合法(即强耦合法)。直接耦合法的声固耦合方程为

式中M，C，K为结构与声学流体区域的总体质量、阻尼、刚度矩阵，上标P表示与声学流体相关的物理量，上标fs表示与声固耦合相对应的物理量。它们可用标准有限元程序获取，详情参见文献［13］。{u}表示结构节点位移，{P}为结构振动产生的声学流体节点压力。{Fu}为结构节点等效力列向量，{FP}为有限水域节点等效加速度列向量。将式(1)重新写成分块形式为:

式中下标1和2分别代表除了无限水域截断边界之外的所有自由度对应的变量和无限水域截断边界上自由度对应的变量。式(2)中，除了{}之外，其余物理量均可由标准有限元法获取，而{}由下节SBFEM理论获取。

3 无限水域的SBFEM公式

文献［10］给出了基于SBFEM的、在无限水域截断边界上的压力{P2}的表达式。其利用SBFEM将无限水域从无限水域截断面到无限远的区域离散成多个半无限扇形水域形式，网格形式见图2，即当用SBFEM离散无限水域时，维数减一，只需离散无限水域截断面，无限水域离散截断面的每个网格均代表相对应的半无限扇形水域(三维问题则为半无限楔形水域)。离散无限水域时，网格应具有相同的相似中心(即所有半无限扇形的两边应相交于一点)。

图2 无限水域的SBFEM离散模型Fig.2 SBFEM discretizationmodel of infinite water

根据SBFEM理论，离散后无限水域截断面上的控制方程［10］满足

4 水下结构声固耦合公式

将式(3)代入式(2)得:

式(4)为水下结构声固全耦合公式。该公式可采用Newmark直接积分法进行水下结构的瞬态响应分析。

5 算例

5.1 内压作用下水下椭圆壳瞬态响应

分析浸没在无限海水域中、由两共焦椭圆构成的椭圆环在均布内压作用下的瞬态响应，考察检验式(4)的正确性。文中将无限海水分割成有限水域和“无限”水域，见图3，“无限”水域分别采用无限元法和比例边界有限元法离散，有限水域与椭圆环分别用四节点声学有限元和四节点平面应力有限元单元离散。采用无限元时，“无限”水域截断半径分别为长半长轴的1.2倍、7倍、20倍、40倍，而采用SBFEM时，“无限”水域截断半径只取为长半长轴的1.2倍。图3给出了1.2倍时的水下结构与有限水域的网格形式。均布内压在5 ms内从零逐渐变为单位压力，之后保持不变。共焦椭圆的半长轴分别为53.85 cm 和 52.83 cm，半短轴分别为 36.32 cm和34.79 cm，弹性模量为 200 GPa，密度为 7 849.7 kg/m3。海水的密度为1 023.1 kg/m3，海水中声音速度为 1 498.6 m/s。

图3 水下椭圆筒和无限水域的有限元网格Fig.3 Finite elementmesh ofsubmerged cylindrical ellipse and infinite water

图4给出了用对应静态位移规则化后的椭圆环半短轴处径向位移图。采用无限元离散“无限”水域时，随着截取半径的增大，结果逐渐收敛。20倍和40倍的结果完全重合，而1.2倍的结果则远偏离20倍的结果，7倍的结果则较接近20倍的结果，这是因为无限单元只能近似模拟无限水域无反射特性，用其模拟无限水域时要求保留较大的有限水域，即较大的截取半径(如6倍以上)。当采用SBFEM模拟无限水域时，因其在径向方向是解析解，无限水域截取半径的大小对结果影响不大。图4中给出了用SBFEM模拟无限水域截取半径为长半长轴1.2倍时的结果，其与用无限单元模拟无限水域时20倍、40倍的结果重合。这表明SBFEM离散无限水域比无限元更具优势。图4中1.2倍、7倍、20倍、40倍的结果采用Ansys软件，利用有限元与无限元耦合分析而得，而SBFEM 1.2倍结果由式(4)计算而得。

图4 椭圆环半短轴处径向位移Fig.4 Radial displacement at theminor axis of cylindrical ellipse

5.2 地震激励下重力坝瞬态响应

利用式(4)计算文献［11］中已分析的重力坝库系统(见图5)在图6水平向地震激励下的坝库耦合瞬态响应。分析中，不考虑材料和几何的非线性，忽略水库自由表面波，库底、库岸的吸收性及坝体的材料阻尼，并假定库底刚性。

图5 重力坝库系统Fig.5 Gravity dam －reservoir system

图6 斜面水平向地震加速度Fig.6 Ramped horizontal ground acceleration

坝体用平面应变四节点等参单元离散，近场用四节点声学流体单元离散，远场用两节点比例边界有限元单元离散。远场离坝底6 m。其网格图见图7。图中颜色相对较浅的为坝体单元，而颜色较深的为声学流体单元。

图7 重力坝及有限水域的网格图Fig.7 Mesh of gravity dam and finite water

坝体分析参数为:密度2 400 kg/m3，泊松比0.2，杨氏模量25 GPa。水的密度为1 000 kg/m3，水中声音速度为1 438.656 m/s。Newmark时间积分的时间增量Δt=0.002 s。坝底在图6斜面水平向地震加速度作用下的动水压力如图8所示。

图8中弱耦合的结果来自文献［11］，而强耦合结果由式(4)计算而得。两者结果几乎完全重合。弱耦合结果的正确性在文献［11］中已说明。

图8 斜面加速度下重力坝坝底水动压力Fig.8 Hydrodynamic pressure at the heel of gravity dam under ramped horizontal ground acceleration

6 结论

文中建立了水下结构声固强耦合公式，分析了水下结构的瞬态响应。数值算例表明了该公式的正确性。该公式具有以下特点:(1)利用该公式计算水下结构响应时，不需要进行水下结构与无限水域响应计算的交错迭代，能有效提高分析效率;(2)该公式继承了SBFEM方法模拟无限水域的优点，即离散区域小，自由度少，精确满足无限水域无穷远处的无反射特性;(3)该公式实现了SBFEM与FEM的无缝耦合;(4)该公式为水下结构声固耦合提供了一种新的求解方法。

［1］KALATEH F，ATTARNEJAD R.Finite element simulation of acoustic cavitation in the reservoir and effects on dynamic response of concrete dams［J］.Finite Elements in Analysis and Design，2011，47(5):543－558.

［2］GOGOI I，MAITY D.A non－reflecting boundary condition for the finite element modeling of infinite reservoir with layered sediment［J］.Advances in Water Resources，2006，29(10):1515－1527.

［3］BENJAMIN M，NAJIB B.Simplified evaluation of the vibration period and seismic response of gravity dam－water systems［J］.Engineering Structures，2010，32(8):2488－2502.

［4］SAMIIA，LOTFIV.Application of H－W boundary condition in dam－reservoir interaction problem［J］.Finite Elements in Analysis and Design 2012，50(1):86－97.

［5］KHAZAEE A，LOTFIV.Application of perfectly matched layers in the transient analysis of dam－reservoirsystems［J］.Soil Dynamics and Earthquake Engineering 2014，60(2):51－68.

［6］WANG G H，ZHANG S R.Damage prediction of concrete gravity dams subjected to underwater explosion shock loading［J］.Engineering Failure Analysis 2014，39(1):72－91.

［7］SEGHIR A，TAHAKOURT A，BONNET G.Coupling FEM and symmetric BEM for dynamic interaction of dam－reservoir systems［J］.Engineering Analysis with Boundary Elements，2009，33(10):1201 －1210.

［8］WOLF JP，SONGCM.Finite－elementmodeling of unbounded media［M］.New York:Wiley，1996:11－15.

［9］WOLF JP.The Scaled Boundary Finite Element Method［M］.Chichester:Wiley，2003.51－54.

［10］FAN SC，LISM，YU G Y.Dynamic fluid－structure interaction analysis using boundary finite elementmethod－finite element method［J］.Journal of Applied Mechanics，ASME，2005，72:591－598.

［11］李上明.基于比例边界有限元法的坝库瞬态耦合分析［J］.华中科技大学学报:自然科学版，2011，39(9):108－111.

［12］李上明.基于比例边界有限元法和有限元法的水下结构瞬态分析方法［J］.工程力学，2013，30(11):42－46.

［13］王勖成.有限单元法［M］.北京:清华大学出版社，2003.523－533.