基于場分離的水下爆炸流固耦合計算方法*

郭 君，楊文山，姚熊亮，張阿漫，任少飛

(哈爾濱工程大學船舶工程學院，黑龍江 哈爾濱150001)

結構在水下爆炸中的動態響應問題涉及結構、材料、流體力學等諸多領域，是研究的熱點和難點。C.A.Felippa 等［1］用有限元方法對水下爆炸進行了數值模擬，此后，譜單元法作為一種計算精度更高的算法發展起來。譜單元法沿用了有限元法的基本框架，用由勒讓德多項式構成的高階拉格朗日插值函數代替有限元法中的三線性函數作為基函數，經過研究比較發現，譜單元法相比有限元法在計算精度和計算效率上有著更大的優勢［2］。

當入射波在流體網格中傳播時，數值計算的誤差會使實際入射波大小和數值計算中的不一致，誤差的累積會導致計算結果的巨大誤差。研究表明，隨著入射波在流體網格中傳播距離的增加，總體結構動力響應誤差呈增大的趨勢［3］。為此，本文中采用場分離技術，將總場分為數值已知的入射場、平衡場和數值未知的散射場，流體網格中僅進行散射場的計算，入射場和平衡場作為已知量施加到結構上進行動響應分析，從而避免了入射波在網格中的傳播，有利于計算精度的提高。

1 散射場數值模型

1.1 場分離

根據壓力產生的原因及特征可以把總壓分成3 部分:平衡壓、入射壓和散射壓。平衡壓由大氣壓力和靜水壓力相加得到，入射壓由在無界均勻流體中傳播的聲學沖擊波壓力得到，散射壓由結構和自由面引起，由控制方程迭代計算得到。相應于3 種壓力定義3 種場來計算水下爆炸流固耦合問題，這3 種場為:平衡場、入射場和散射場。在進行水下爆炸流固耦合計算時，僅用散射場中的物理量進行流體物理量迭代計算，入射場和平衡場作為已知量參與計算，從而避免了因數值誤差導致的入射波在流體中傳播的失真，使計算精度大大提高，相應于以往采用的總場模型，把這種數值模型稱為散射場模型［1］。

1.2 散射場控制方程

總場計算時流體運動控制方程組如下［4-5］

式中:s(X，t)為總壓縮率，ψ(X，t)為總位移勢，X 為流體節點坐標，pe為水的平衡壓力，c 為水中的聲速，pv為水的汽化壓力，可忽略不計。

根據場分離的介紹，式(1)中總壓縮率可分為s=si+ss，其中si為入射壓縮率，ss為散射壓縮率。總位移勢可分為ψ=ψi+ψs，其中ψi為入射位移勢，ψs為散射位移勢。s 和ψ 均為減掉平衡場后的值，因此等式中沒有se和ψe項，其中se為平衡壓縮率，ψe為平衡位移勢。

將s=si+ss，ψ=ψi+ψs代入式(1)，得到散射場控制方程［4-5］

本文中初始入射壓力pi的表達式為

式中:p0和τ 分別是壓力峰值和衰減時間，H()是海維西特階躍函數。入射位移勢ψi可通過對pi二次積分得到。

1.3 散射場空間離散

采用等參元對流體域進行離散，將其離散成n 個六面體單元，則每個單元的幾何形狀可表示為

式中:X、Y、Z 為整體坐標系下的節點坐標，φ 為形函數，這里采用三線性函數，上標T 代表轉置，ξ、η、ζ為自然坐標系下的節點坐標。

在自然坐標系內，散射壓縮率ss和散射位移勢ψs可離散為［4-5］

式中:ss，ijk( t) 為每個單元內(N+1)3個節點分別對應的壓縮率，ψs，ijk( t) 為每個單元內(N+1)3個節點分別對應的位移勢，φi為由N 階多項式構成的基函數，在譜單元方法中，φi取為由勒讓德多項式構成的一維拉格朗日插值函數［6］

式中:PN為N 階勒讓德多項式，ξi為第i 個Gauss-Lobatto-Legendre 積分點。

用標準伽遼金法對控制方程進行離散，用φ 左乘式(1)，在流體域上對方程積分，并應用格林第一公式可得

式中:Γ 是流體域的表面，n 是Γ 的法線矢量，背離流體為正。將式(5)代入式(7)可得

1.4 散射場時間步控制

由于數值解不穩定，在流體積分過程中可能會出現數值振蕩，因此采用人工阻尼β 來抑制數值振蕩的出現。此外，為控制計算過程的穩定性，必須嚴格控制積分時間步長，對于流體方程，臨界時間步長［7］

式中:c 為水中聲速，λmax為單元的廣義本征問題(H-λQ)z=0 的最大本征值，可由下式得到

式中:Hij和Qii分別為H 和Q 中的元素。

1.5 誤差系數

采用T.L.Geers 提出的廣義誤差系數分析瞬態響應數值解與基本解之間的誤差［8］，廣義誤差系數

式中:c(t)是瞬態響應數值解，b(t)是對應時刻的基本解，t1≤t ≤t2是所考核的時間區域。

2 散射場數值計算及結果分析

2.1 計算模型

如圖1 所示，計算模型為位于流體柱上的一個兩自由度的質量彈簧振子，遭受平面沖擊波作用，圖中質量塊m1表示船體外板，質量塊m2表示船體內部的結構和設備，兩質量塊由剛度系數為k 的彈簧連接，m1和m2的位移分別用u1(t)和u2(t)表示，速度分別用v1(t)和v2(t)表示［9-10］。

計算模型的物理屬性為:大氣壓力patm=0.101 MPa，重力加速度g=9.81 m/s2，質量塊由邊長為0.3 m 的正方形板表示，m1=76.9 kg，m2=384.5 kg，彈簧剛度系數k=94.87 t/s2，流體密度ρ=1.026 t/m3，流體中聲速c=1.5 km/s，流體長度為3 m，壓力峰值p0=16.2 MPa，延遲時間τ=0.42 ms，沖擊波波陣面距m1的距離為di=2.7 m時開始進行積分計算。

2.2 計算結果及討論

2.2.1 兩種數值模型的結構響應對比

有限元方法作為一種數值方法，廣泛應用于水下沖擊問題的求解中，并證明具有一定的精度［1］，針對2.1 節的物理模型，采用由高度細化的有限元網格(24 000 個等長單元)產生的有限元解作為基準解，以1.4 節中闡述的廣義誤差系數為對比參數，討論總場模型和散射場模型在不同細化程度時的精度。

通過自編程序，分別采用總場模型和散射場模型計算m1的速度響應，并分別與基準解對比。不同細化程度時質量塊m1的速度隨時間的變化曲線如圖2 所示。

為方便討論，根據圖2 中的數值結果，可以得到總場模型、散射場模型與基準解的誤差Ct、Cs如表1所示。從表1 可以看出，在相同細化階數的條件下，散射場模型的數值結果與總場模型相比誤差較小，說明散射場模型在計算水下爆炸流固耦合問題時具有更高的精度，此外，散射場模型和總場模型所產生的誤差均隨細化階數的升高逐漸減小，當N=4 時，兩者結果非常接近，均收斂于基準解。從圖2 可以看出，散射場模型和總場模型的結構響應均有微小振蕩，這是由數值誤差引起的，這種振蕩隨著細化階數的升高而不斷光順，在細化階數為4 時已經能得到較為光順的數值解。

圖1 計算模型Fig.1 Calculation model

表1 總場模型及散射場模型與基準解之間的誤差Table 1 Error between the response of TFM，SFM and benchmark

圖2 結構的速度響應對比曲線Fig.2 Contrast curves of structure’s velocity response

2.2.2 兩種數值模型的入射波失真對比

下面以N=3 為例討論連續入射波在非連續流體網格中傳播時的失真現象。圖3 給出了不同時刻入射波大小隨深度變化的曲線。由圖3 可以看出，入射波在網格中傳播時，隨著傳播距離的增加失真程度呈增大趨勢，這主要是傳播過程中入射波數值誤差不斷積累的緣故。由于散射場模型有效避免了入射波在傳播過程中產生的數值誤差，所以散射場模型相比總場模型具有更高的精度，上文中的結構響應誤差對比分析已經充分證明了這一點。

圖3 入射壓傳播圖Fig.3 Incident-wave pressure profile

3 結 論

根據水下爆炸中個各物理量產生的原因及特征將整個流體分為平衡場、入射場和散射場，通過理論推導建立了散射場數值模型，以水下爆炸典型問題——遭受水下爆炸沖擊波作用的平板彈簧模型為算例，經自編程序，分別用散射場模型和總場模型計算了結構響應，通過對比分析得到以下結論:

(1)由于數值誤差的不斷積累，隨著入射波在網格中傳播距離的增加，其失真程度呈增大趨勢。散射場模型有效避免了入射波在傳播過程中產生的數值誤差，相比總場模型具有更高的精度，可以較好地應用于水下爆炸流固耦合問題的求解中。

(2)散射場模型的誤差隨著細化階數的升高逐漸減小，且隨著模型的細化，數值結果收斂于基準解，因此可以通過提高細化階數來提高散射場模型的計算精度。

［1］ Felippa C A，DeRuntz J A.Finite element analysis of shock-induced hull cavitation［J］.Computer Methods in Applied Mechanics and Engineering，1984，44(3):297-337.

［2］ Komatitsch D，Vilotte J P.The spectral element method:An efficient tool to simulate the seismic response of 2D and 3D geological structures［J］.Bulletin of the Seismological Society of America，1998，88(2):368-392.

［3］ Shin Y S，Santiago L D.Surface ship shock modeling and simulation:Two-dimensional analysis［J］.Shock and Vibration，1998，5(2):129-137.

［4］ Priolo E，Carcione J M，Seriani G.Numerical simulation of interface waves by high-order spectral modeling techniques［J］.Journal of the Acoustical Society ofAmerica，1994，95(2):681-693.

［5］ Bleich H H，Sandler I S.Interaction between structures and bilinear fluids［J］.International Journal of Solids and Structures，1970，6(5):617-639.

［6］ Ronquist E M，Patera A T.A Legendre spectral element method for the Stefan problem［J］.International Journal for Numerical Methods in Engineering，1987，24(2):2273-2299.

［7］ Mulder W A.Spurious modes in finite-element discretizations of the wave equation may not be all that bad［J］.Applied Numerical Mathematics，1999，30:425-445.

［8］ Geers T L.An objective error measure for the comparison of calculated and measured transient response histories［J］.The Shock and Vibration Bulletin，1984，54(2):99-108.

［9］ Sprague M A，Geers T L.A spectral-element method for modelling cavitation in transient fluid-structure interaction［J］.International Journal for Numerical Methods in Engineering，2004，60(15):2467-2499.

［10］ Sprague M A，Geers T L.Computational treatments of cavitation effects in near-free-surface underwater shock analysis［J］.Shock and Vibration，2001，8(2):105-122.