郭 君，張阿漫，楊文山，李 佳

(哈爾濱工程大學(xué) 船舶工程學(xué)院，哈爾濱 150001)

水下爆炸載荷作用下的艦船結(jié)構(gòu)毀傷一直是國(guó)內(nèi)外學(xué)者研究的熱點(diǎn)［1－3]，與之相對(duì)應(yīng)的水下沖擊流固耦合問題主要是求解瞬態(tài)沖擊波激勵(lì)下結(jié)構(gòu)的響應(yīng)，其研究?jī)?nèi)容涉及到結(jié)構(gòu)、材料、流體力學(xué)等諸多領(lǐng)域。當(dāng)結(jié)構(gòu)遭受水下沖擊波載荷作用時(shí)，可壓縮瞬態(tài)波在流體中傳播并撞擊到結(jié)構(gòu)上，在結(jié)構(gòu)濕表面處會(huì)發(fā)生繞射，如果結(jié)構(gòu)是彈性的，且周圍水的靜壓力足夠低，散開的稀疏波可能會(huì)導(dǎo)致子區(qū)域內(nèi)產(chǎn)生氣穴，我們稱之為氣穴現(xiàn)象［4]。產(chǎn)生氣穴區(qū)域的水的體積模量比未產(chǎn)生氣穴區(qū)域的水的體積模量要小好幾個(gè)數(shù)量級(jí)，呈現(xiàn)非線性特性。由于結(jié)構(gòu)尺寸、入射波衰減長(zhǎng)度和氣穴區(qū)域在空間尺寸上差別很大，氣穴會(huì)對(duì)近自由面結(jié)構(gòu)如水面艦船產(chǎn)生很大影響，如果一個(gè)45 kg的HBX－1炸藥在水下爆炸，產(chǎn)生氣穴區(qū)域，產(chǎn)生的氣穴區(qū)域在水平上最長(zhǎng)為489 m，深度上最大為9 m，其中船體尺寸為:寬9.1 m、高4.6 m、吃水 3.1 m。

為處理氣穴現(xiàn)象，Newton［5]提出了一種非常有效的計(jì)算方案來處理氣穴聲學(xué)流體，該方案采用位移勢(shì)作為流體域有限元方程中的主要變量，位移勢(shì)是一個(gè)標(biāo)量。由于傳統(tǒng)流體有限元方程中采用的變量為位移矢量，位移矢量表達(dá)式中的變量是位移勢(shì)表達(dá)式中的未知量數(shù)目的三倍，而且不能自動(dòng)對(duì)流體運(yùn)動(dòng)實(shí)行無旋處理，因此，與位移矢量相比，位移勢(shì)的選取具有相當(dāng)大的優(yōu)勢(shì)。Newton方案的第二個(gè)特點(diǎn)是應(yīng)用了顯式時(shí)間積分方法。后來，F(xiàn)elippa和 DeRuntz［6]對(duì) Newton的方法進(jìn)行了擴(kuò)展，對(duì)流固耦合方程進(jìn)行三維表述，并對(duì)結(jié)構(gòu)和流體進(jìn)行有限元處理。然而，有限元方法采用的單元通常為低階單元，在解決波的傳播問題時(shí)會(huì)出現(xiàn)數(shù)值發(fā)散，要想得到準(zhǔn)確的數(shù)值結(jié)果需要十分精細(xì)的網(wǎng)格，尤其當(dāng)有氣穴現(xiàn)象出現(xiàn)時(shí)，網(wǎng)格細(xì)化更是尤為重要，從而導(dǎo)致計(jì)算量很大。為此，基于Felippa和DeRuntz的三維表述和結(jié)構(gòu)有限元處理，本文引進(jìn)譜單元對(duì)流體進(jìn)行處理，以探索一種更為有效的水下爆炸流固耦合求解方法。文中采用基于勒讓德多項(xiàng)式的譜單元［7－8]，這是一種高階單元，采用高階單元可以不需要細(xì)化網(wǎng)格便獲得收斂性。此外，譜單元還可以很好地顯示波的傳播特性［9－10]。譜單元方法綜合了有限元方法的幾何靈活性和譜方法的計(jì)算高精度性，在保證精度的同時(shí)還能提高計(jì)算速度。

1 譜單元數(shù)值模型

1.1 控制方程

在推導(dǎo)產(chǎn)生氣穴的聲學(xué)流體的控制方程時(shí)，首先采用雙線性模型［11]來考慮流體中的氣穴現(xiàn)象。雙線性模型是一種簡(jiǎn)單合理的流體構(gòu)成模型，在這種模型中，當(dāng)壓縮率(負(fù)體積應(yīng)變)為正時(shí)，體積模量為水作為聲學(xué)介質(zhì)時(shí)的值;當(dāng)壓縮率為負(fù)時(shí)，體積模量為零。雙線性狀態(tài)方程如下:

式中，Ba是不存在氣穴時(shí)的體積模量是流體的總壓力，pv是溫度T0時(shí)水的蒸汽壓力，ρ0是溫度T0飽和流體的密度是瞬時(shí)密度是流體質(zhì)點(diǎn)的位置矢量，t表示時(shí)間。假定在整個(gè)討論域內(nèi)T0、ρ0和ρv是常量，定義，其中c是無氣穴聲學(xué)流體中的聲速，S=(ρ/ρ0－1)是總的壓縮率，通過這些定義可將狀態(tài)方程改寫為:

假設(shè)流體是無旋的，且忽略粘性的影響，從而保證了流場(chǎng)用勢(shì)來進(jìn)行表述。用位移勢(shì)來表示流場(chǎng)的位移為，其中和分別為流體質(zhì)點(diǎn)的總位移和平衡位移為位移勢(shì)是流體質(zhì)點(diǎn)的整體位置矢量。

根據(jù)流體的連續(xù)性方程和歐拉運(yùn)動(dòng)方程，并引用位移勢(shì)，可得到流體的平衡方程如下:

其中peq=－ρ0(Vg+h)是平衡壓力。

式(2)與式(3)便構(gòu)成了流體的控制方程組，利用關(guān)系式peq=pv+ρ0(Vg+h)，并定義壓縮率s=ρ0(SSeq)，位移勢(shì) ψ =ρ0Ψ，則根據(jù)式(2)和式(3)，控制方程組可寫為:

式中為動(dòng)態(tài)壓縮率為動(dòng)態(tài)位移勢(shì)為流體節(jié)點(diǎn)坐標(biāo)，peq為水的平衡壓力，c為水中的聲速，pv為水的汽化壓力，由于pv比平衡壓力小很多，故可忽略不計(jì)，上標(biāo)點(diǎn)表示變量對(duì)時(shí)間的導(dǎo)數(shù)。

1.2 譜單元空間離散

與有限元離散相同，采用等參元對(duì)整個(gè)流體域進(jìn)行離散，將其離散成ne個(gè)六面體單元，則每個(gè)單元的幾何形狀可表示為:

φ仍然采用三線性函數(shù)。

然而，與有限元離散不同的是，在自然坐標(biāo)系內(nèi)，將壓縮率s和位移勢(shì)ψ離散為:

式中:φ為由N階多項(xiàng)式構(gòu)成的基函數(shù)，s為每個(gè)單元內(nèi)(N+1)3個(gè)節(jié)點(diǎn)分別對(duì)應(yīng)的壓縮率，ψ為每個(gè)單元內(nèi)(N+1)3個(gè)節(jié)點(diǎn)分別對(duì)應(yīng)的位移勢(shì)，ξ，η，ζ為單元的自然坐標(biāo)(－1≤ξ，η，ζ≤1)。

譜單元方法和有限元方法的最大區(qū)別在于φi的選擇，在有限元法中，φi取為與幾何變換相同的三線基函數(shù)，而在譜單元方法中，φi取為由勒讓德多項(xiàng)式構(gòu)成的一維拉格朗日插值函數(shù)［12]，這也是目前國(guó)內(nèi)外應(yīng)用最廣泛最有效的譜單元基函數(shù):

式中 PN為 N階勒讓德多項(xiàng)式，ξi為第 i個(gè) Gauss-Lobatto-Legendre積分點(diǎn)，其值為下式的根:

可以證明，表達(dá)式(6)滿足下列關(guān)系式:

式中δij為克羅內(nèi)克δ。

用標(biāo)準(zhǔn)伽遼金法對(duì)控制方程式(4a)進(jìn)行離散，用φ左乘式(4a)，在流體域上對(duì)方程積分，并應(yīng)用格林第一公式可得:

式中Q，H，b分別為:

2 平板的水下爆炸流固耦合研究

2.1 計(jì)算模型

為驗(yàn)證所建立的數(shù)值模型的準(zhǔn)確性，需將數(shù)值解與準(zhǔn)確解進(jìn)行對(duì)比。這里采用的計(jì)算模型為經(jīng)典的Bleich-Sandler平板模型，如圖4所示，平板位于半無限長(zhǎng)流體柱上，遭受平面步指數(shù)沖擊波作用，平板速度響應(yīng)的解析解可由特征法得到［11]。

圖1 Bleich-Sandler平板模型Fig.1 Bleich-Sandler plate model

計(jì)算模型的物理屬性為:大氣壓力為Patm=0.101 MPa，重力加速度為g=9.81 m/s2，流體密度為 ρ=998 kg/m3，流體中聲速為 c=1 450 m/s，流體長(zhǎng)度為3.81 m，由100個(gè)等邊長(zhǎng)的六面體單元構(gòu)成，平板密度為ρ=5 698 kg/m3，邊長(zhǎng)為0.038 1 m，由一個(gè)單元構(gòu)成，平面步指數(shù)波的壓力峰值為p0=0.71 MPa，延遲時(shí)間為τ=0.1 ms，在沖擊波波陣面距離平板一個(gè)流體網(wǎng)格時(shí)開始積分計(jì)算［13]。

2.2 計(jì)算結(jié)果及討論

平板中心的速度時(shí)歷曲線如圖2所示，實(shí)線顯示的是由譜單元法得到的數(shù)值解，圓點(diǎn)和三角形點(diǎn)則分別表示的是由特征法得到的考慮氣穴和不考慮氣穴時(shí)的平板中心速度響應(yīng)的解析解。分別計(jì)算了人工阻尼β 取作0.25、0.50、0.75 和 1.0 的四種情況，臨界時(shí)間步根據(jù)2.3節(jié)中的表述進(jìn)行計(jì)算，計(jì)算時(shí)間步取作臨界時(shí)間步的0.75倍。

圖2給出了人工阻尼β取不同值時(shí)平板中心的速度隨時(shí)間變化的曲線，對(duì)計(jì)及氣穴和不計(jì)氣穴的兩種情況進(jìn)行了計(jì)算，并分別與Bleich-Sandler的解析解進(jìn)行對(duì)比。從圖2可以看出，由譜單元方法得到的數(shù)值結(jié)果與Bleich-Sandler解析解吻合較好，同時(shí)人工阻尼β對(duì)結(jié)構(gòu)的響應(yīng)幅值和響應(yīng)時(shí)間幾乎沒有影響，圖2(a)中曲線的微小振蕩是由流體網(wǎng)格的瞬時(shí)擾動(dòng)所引起，隨著β的增加，振蕩逐漸消失。從圖2還可以看出，氣穴的存在不會(huì)改變速度的峰值，但會(huì)延長(zhǎng)速度的衰減時(shí)間，從而延長(zhǎng)載荷對(duì)結(jié)構(gòu)的作用時(shí)間，對(duì)結(jié)構(gòu)產(chǎn)生更大的破壞。

3 譜單元方法和有限元方法的對(duì)比研究

如前面部分提到，譜單元法相比有限元法具有很多優(yōu)點(diǎn)，這節(jié)將分別從精度和效率上對(duì)此進(jìn)行討論。

3.1 計(jì)算模型

如圖3所示，計(jì)算模型為位于流體柱上的一個(gè)兩自由度的質(zhì)量彈簧振子，遭受平面步指數(shù)沖擊波作用，圖中質(zhì)量塊m1表示船體外板，質(zhì)量塊m2表示船體內(nèi)部的結(jié)構(gòu)和設(shè)備，兩質(zhì)量塊由剛度系數(shù)為k的彈簧連接，m1和m2的位移分別用u1(t)和u2(t)表示，速度分別用V1(t)和V2(t)表示。

計(jì)算模型的物理屬性為:大氣壓力為Patm=0.101 MPa，重力加速度為 g=9.81 m/s2，質(zhì)量塊由邊長(zhǎng)為0.3 m的正方形板表示，m1=76.9 kg，m2=384.5 kg，彈簧剛度系數(shù)為k=94 870 kg/s2，流體密度為ρ=1 026 kg/m3，流體中聲速為c=1 500 m/s，流體長(zhǎng)度為3 m，壓力峰值為 p0=16.2 MPa，延遲時(shí)間為 τ=0.42 ms，沖擊波波陣面距m1的距離為dinc=0.3 m時(shí)開始進(jìn)行積分計(jì)算。

圖2 平板中心速度隨時(shí)間的變化曲線Fig.2 The flat center speed changing with time curve

3.2 計(jì)算結(jié)果及討論

有限元方法作為一種數(shù)值方法，已經(jīng)被廣泛應(yīng)用于水下沖擊問題的求解中，并證明具有一定的精度［6]，針對(duì)3.1 節(jié)的物理模型，我們采用由高度細(xì)化的有限元網(wǎng)格(24 000個(gè)等長(zhǎng)單元)產(chǎn)生的有限元解作為基準(zhǔn)解，討論譜單元方法(SEM)和有限元方法(FEM)在不同細(xì)化程度時(shí)的精度。

圖3 計(jì)算模型Fig.3 Calculation Model

分別采用譜單元方法和有限元方法計(jì)算m1的速度響應(yīng)，并分別與基準(zhǔn)解對(duì)比。不同細(xì)化程度時(shí)質(zhì)量塊m1的速度隨時(shí)間的變化曲線如圖4所示。

圖4給出了分別采用譜單元方法和有限元方法在不同細(xì)化程度下算得的質(zhì)量塊m1的速度響應(yīng)曲線，分別將有限元解和譜單元解與基準(zhǔn)解進(jìn)行對(duì)比分析。誤差系數(shù)見表1。

表1 有限元解及譜單元解與基準(zhǔn)解之間的誤差系數(shù)cTab.1 Finite element solution and the spectral element solution and the error c between the coefficient of benchmark solutions

圖4 結(jié)構(gòu)的速度響應(yīng)對(duì)比曲線Fig.4 The comparison curve of the velocity response of the structure

從表1給出的誤差系數(shù)可以看出，譜單元方法和有限元方法所產(chǎn)生的誤差均隨著細(xì)化階數(shù)的升高逐漸減小，收斂于基準(zhǔn)解。然而，在細(xì)化階數(shù)相同(系統(tǒng)所具有的自由度相同)的條件下，譜單元方法得到的數(shù)值解更接近基準(zhǔn)解，因此，與有限元方法相比，譜單元方法具有更高的精度。

從表1還可以看出，在誤差相同的條件下，譜單元方法需要的自由度較少。例如:若將c=0.15設(shè)為允許誤差，用譜單元方法需要不足496個(gè)自由度便能達(dá)到允許誤差，而用有限元方法需要的自由度數(shù)則接近1 025。因此，譜單元方法在精度提高的同時(shí)可以大大減少計(jì)算所需要的自由度數(shù)，具有更高的計(jì)算效率。

4 結(jié)論

本文基于流體雙線性構(gòu)成關(guān)系考慮流體中可能發(fā)生的氣穴現(xiàn)象，采用基于勒讓德多項(xiàng)式的譜單元方法建立了水下爆炸流固耦合動(dòng)力學(xué)模型，應(yīng)用經(jīng)典的平板沖擊問題對(duì)模型進(jìn)行了驗(yàn)證。分別采用譜單元方法和有限元方法對(duì)彈簧－平板模型進(jìn)行水下爆炸流固耦合的求解，基于廣義誤差系數(shù)對(duì)譜單元方法和有限元方法進(jìn)行對(duì)比研究，得到以下結(jié)論:

(1)采用譜單元方法對(duì)Bleich-Sandler平板水下爆炸流固耦合問題進(jìn)行計(jì)算，得到的數(shù)值解與Bleich-Sandler的基本解吻合良好，譜單元方法能較好地應(yīng)用于水下爆炸流固耦合問題的求解中。

(2)流體中產(chǎn)生氣穴時(shí)，未產(chǎn)生氣穴的流體區(qū)域?qū)⒔Y(jié)構(gòu)與產(chǎn)生氣穴的流體區(qū)域分開，形成流體堆積。氣穴區(qū)域閉合時(shí)流體中形成正壓力沖擊作用，產(chǎn)生二次加載現(xiàn)象，從而對(duì)結(jié)構(gòu)產(chǎn)生更嚴(yán)重的破壞，因此在計(jì)算中應(yīng)予以考慮。

(3)人工阻尼對(duì)結(jié)構(gòu)響應(yīng)具有光順作用，但基本不會(huì)影響結(jié)構(gòu)的響應(yīng)幅值和響應(yīng)時(shí)間。隨著人工阻尼的增加，由流體網(wǎng)格瞬時(shí)擾動(dòng)引起的結(jié)構(gòu)響應(yīng)的振蕩逐漸消失。

(4)與有限元方法相比，在同樣網(wǎng)格細(xì)化的條件下，譜單元方法具有更高的精度;在誤差相同的條件下，譜單元方法所需的自由度數(shù)大大小于有限元方法所需的自由度數(shù)。因此，與有限元方法相比，譜單元方法在提高精度的同時(shí)，能大量節(jié)省計(jì)算時(shí)間。

［1] 牟金磊，朱 錫，黃曉明，等.水下爆炸氣泡載荷在加筋板塑性變形中的作用［J] .振動(dòng)與沖擊，2010，29(5):74－77.

［2] 李海濤，朱 錫，趙小龍，等.箱形梁在水下近距非接觸爆炸作用下的整體毀傷研究［J] .振動(dòng)與沖擊，2010，29(3):158－161.

［3] Ramajeyathilagam K，Vendhan C P，Bhujanga Rao V.Nonlinear transient dynamic response of rectangular plates under shock loading［J] . International Journal of Impact Engineering，2000，24:999 －1015.

［4] Kennard E H.Cavitation in an elastic liquid［J] .Physical Review，1943，63(5/6):172－181.

［5] Newton R E.Finite element analysis of shock-induced cavitation，1980. Preprint 80 － 110， ASCE Spring Convention.

［6] Felippa C A，DeRuntz J A.Finite element analysis of shockinduced hull cavitation［J] .Computer Methods in Applied Mechanics and Engineering，1984，44:297－337.

［7] 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:2273－2299.

［8] Maday Y，Patera A T.Spectral element methods for the incompressible Navier-Stokes equations［J] .The American Society of Mechanical Engineers，New York，1989.

［9] 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.

［10] 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.

［11] Bleich H H，Sandler I S.Interaction between structures and bilinear fluids［J] .International Journal of Solids and Structures，1970，6:617－639.

［12] 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:2273－2299.