劉云龍,汪 玉,張阿漫
(1.武漢第二船舶設計研究所,湖北 武漢 430205;2.哈爾濱工程大學船舶工程學院,黑龍江 哈爾濱150001;3.中國人民解放軍92857部隊,北京100161)
圓柱殼作為潛艇典型的結構形式,常被用作潛艇結構強度機理性研究的對象。在很多情況下對圓柱殼結構抗沖擊性能深入研究的成果可以推廣到實際潛艇上。因此,本文中以雙層圓柱殼結構為研究對象,對水下爆炸載荷下圓柱殼結構的鞭狀運動特性進行分析,旨在為潛艇抗沖擊設計和研究提供參考。水下爆炸問題具有瞬態強非線性特點,其涉及的流固耦合問題一直是艦船抗沖擊研究的熱點和難點。早在20世紀60年代,研究人員就嘗試通過解析的方法求解一些簡單的沖擊問題。比較有代表性的是H.Huang等早期的工作,其中得到了一些球殼以及圓柱殼等簡單結構在平面階躍波[1-2]及球面波[3]作用下的彈性響應。盡管此類方法對于實際復雜結構并不適用,但在后來的研究中通常用來驗證理論和數值模型的精度和有效性。隨著計算機性能的提高和邊界元的發展,流固耦合問題的數值求解技術得到了很大的提高。對現代水下沖擊流固耦合研究影響較大的是T.L.Geers在1971年提出的雙漸近(double asymptotic approximation,DAA)法[4]。該方法分別采用平面波假設和勢流假設對水下沖擊問題的前期和后期響應進行近似,而中間頻段則通對兩者進行匹配得到。隨后的研究中,T.L.Geers[5]對DAA法[4]進行了改進,考慮了結構曲率的影響以及前后期的非線性匹配問題,形成了二階雙漸近(second-order doubly-asymptotic approximation,DAA2)法。DAA2法已被廣泛應用于水下沖擊防護研究。針對潛艇結構聲學覆蓋層問題,劉建湖[6]對DAA2法進行改進,考慮波在聲學覆蓋層中的反射和透射特性,將該過程分為聲學反射階段和DAA2階段分別求解,發現用傳統DAA2法處理聲學覆蓋層問題時低估了沖擊波載荷的作用。隨后,孫士麗[7]、王詩平等[8]通過計入非線性的二階速度項,對DAA2法進行進一步改進,得到了非線性雙漸近(non-linear second-order doubly-asymptotic approximation,NDAA2)法。宗智等[9]結合DAA2法和有限元法,計算了氣泡脈動載荷下簡單船體結構的動態響應和應力集中現象。董海等[10]針對潛艇的鞭狀運動問題,將DAA2法和船體梁模型結合在一起,采用DAA2法計算流體載荷,用變截面梁模型計算結構響應,對結構的時域和頻域響應均進行了深入分析。然而,由于DAA2方程是基于無限域三維波動方程推導的,因此對內流場問題不適用。盡管M.A.Sprague等、T.L.Geers等對經典DAA2法進行修改并推廣到內域問題11-13,但在實際計算中經典DAA2法仍有較多限制。因此,本文中對雙層圓柱殼結構舷間水對應的內流場采用聲固耦合求解,而對外流場采用DAA2法求解,兩者在非耐壓殼處耦合,形成適用于雙層圓柱殼結構的流固耦合分析方法,并據此分析水下爆炸載荷下雙層圓柱殼結構的鞭狀運動特性。

圖1 DAA2法后期時域近似示意圖Fig.1 Scheme of later approximation in time domain of DAA2
關于DAA2法,對前期近似采用在平面波假設下的流固耦合模型。在該假設下,任意時刻時邊界上的流體動壓力不受附近其他點的影響,推導過程相對簡單,本文不再贅述。對后期近似,有不同的推導方法。王詩平等[14]通過將延遲積分方程進行拉普拉斯變換,并在拉氏空間內對方程進行簡化,最后通過拉氏逆變換變回時域空間,為DAA2法提供了新的切入點。采用常規格林函數法,重新推導DAA2法后期近似公式。該推導方法更簡便,物理意義更明確,旨在為DAA2法后續研究提供新的視角。假設流場速度勢φ滿足波動方程:

式中:cf為流體中的聲速。若取格林函數,則有:

式中:R=p-q為流場中從源點q指向場點p的矢量。于是,令式(1)×G-式(2)×φ,得:

對式(3)在整個流域V內積分并應用格林第2公式可得:

式中:S為流場V的邊界,n為其單位法向量。式(4)即為求解波動方程的DBEM(domain boundary element method)積分方程。對式(4)左端積分,可采用常規邊界元法離散實現,而右側體積分項,由于積分域是自由場中的無界空間,只能進行近似。根據DAA2法的思想,考慮t cf?L時對方程進行近似,其中L為結構的特征尺度,此時被結構擾動的流場可認為是一個半徑為t cf的球體,結構的影響可看作球體中心的點源和偶極子的貢獻。由于點源和偶極子的誘導速度勢分別按照進行衰減,因此略去偶極子影響所產生的二階項,結構對流場的貢獻相當于球體中心一個點源的影響。于是根據波動方程三維空間的時域基本解,流場中任意點的速度勢可以近似表示為:

式中:m (t)為源強的任意函數。于是式(4)右端的體積分可化為:

式中:C為待定常數。為確定C的取值,將其化為積分形式代入式(4),當t=0時,顯然C=0。記u·n=將式(6)代入式(4),有:

可見式(7)即為DAA2后期近似積分方程。通過平面波假設,沖擊問題的前期近似可表示為[5,14]:

式中:χ為邊界的局部曲率。將前期近似和后期近似在頻域內匹配,整理可得DAA2邊界積分方程[5,14]:

水下爆炸氣泡脈動是一個強非線性的物理化學過程,針對不同的問題已提出了很多氣泡模型。其中T.L.Geers等基于DAA2法提出的單個水下爆炸氣泡脈動模型[15]考慮了流場的可壓縮性和氣泡內部氣體的聲學特性,因此在氣泡各次脈動過程中,能夠真實反映氣泡最大半徑和氣泡脈動壓力峰值的衰減。另外考慮了氣泡中心的上浮,并通過以往水下爆炸氣泡實驗數據對上浮中的阻力進行了匹配,得到結果與實驗結果吻合良好。在隨后的研究中,T.L.Geers等對Geers-Hunter模型進行了一些優化調整[16]。采用Geers-Hunter模型計算得到氣泡的脈動和遷移運動后,則可通過以下方法計算流體中任意一點的流體速度勢[17]:


圖2 水下爆炸氣泡脈動載荷計算示意圖Fig.2 Scheme of the calculation of underwater explosion bubble pulsation load
式中:右側第1項為氣泡脈動等效的點源所誘導的速度勢,第2項為氣泡遷移等效的偶極子所誘導的速度勢,后面2項為自由面鏡像的氣泡所誘導的速度勢。其中e1和e2分別為源強和偶極子強度;kf為自由面的反射因數,由于本文中分析深水中結構對水下爆炸的動態響應,因此不考慮自由面影響,kf取值為0;r1和r2分別為氣泡中心和鏡像氣泡中心的距離;θ1和θ2分別為2連線與鉛垂線的夾角,見圖2。
然后通過非定常伯努利方程得到流場任意位置壓力。氣泡等效的變強度點源和偶極子強度需要滿足氣泡表面的連續性條件,可分別表示為:

式中:a為氣泡等效半徑,vm為氣泡垂向遷移速度,d為水深。
對FEM(finite element method)求解器選擇ABAQUS/Explicit,可通過時間顯式積分的方法計算結構在外載荷作用下的響應。由于不需要迭代過程,計算一般不存在收斂性問題,而時間積分的誤差對于爆炸沖擊這種瞬態分析來講通常為小量,不會對結果產生明顯影響。
圖3中的數字表示時間增量步,在同一個時間增量步中,FEM計算得到的結構變形、速度以及加速度通過接口程序傳遞到BEM程序中,然后根據DAA2法計算得到結構濕表面的流體動壓力,并在下一個時間增量步中傳遞給ABAQUS的FEM求解器,更新FEM中的載荷。重復以上流程,實現BEM和FEM之間的數據傳遞和沿時間的顯式推進。

圖3 ABAQUS用戶子程序數據流示意圖Fig.3 Data flow in the user-subroutine of ABAQUS
FEM同BEM的數據傳遞過程中存在一個延遲問題,即FEM求解器在計算時采用的是上一個時間增量步的壓力數據,對此可采用前2個增量步的載荷數據外插的方法得到當前增量步的載荷:

式中:下標表示增量步數。對于BEM和FEM間的數據傳遞過程,本文中采用非共節點的網格映射方法,一方面減小BEM的計算規模,另一方面盡可能地保持較高的精度。
通過和P.Zhang等[18]針對充水球殼經典算例推導的解析解進行對比,驗證本文計算模型的有效性。該算例旨在考察厚度為1mm、內部充滿水的彈性鋼質球殼在1MPa平面階躍波作用下的響應。
圖4為不同時刻球殼內流場的聲壓云圖。從圖4可知:在相互作用初期,內流場透射沖擊波幅值與入射波基本一致,且波陣面以平面向前推進;隨著時間的增加,由于結構中應力波的波速大于流場波速,因此遠離殼體的球殼中部的沖擊波明顯慢于外側靠近球殼部分流體,從而使波陣面不再保持平面,見圖4(b);當透射波到達球殼下部時,又產生了一定的反射,透射壓力疊加,形成高壓區,見圖4(c)和(d)。

圖4 不同時刻內流場聲壓云圖Fig.4 Acoustic pressure of the inner field at different times
圖5分別是不同方法計算結果中球殼上迎爆點和背爆點在沖擊波方向的速度對比。可以看到,在圖5(a)中內外域均為ABAQUS/聲固耦合計算結果的高頻部分衰減過快,因此在后期的速度曲線較光滑。在圖5(b)中,背爆點產生了2次明顯的加速過程:第1次發生在t=0.6ms左右,為應力波通過殼體傳遞到背爆面所致;第2次發生在t=1.3ms左右,為內流場透射壓力到達背爆面所致。3種方法的計算結果基本一致,表明本文中所建立的計算模型可用于水下沖擊問題的工程和理論研究。

圖5 數值模型對比驗證Fig.5 Comparison among different methods
計算模型如圖6~7所示。計算模型長約75m,非耐壓殼直徑約8m,耐壓殼直徑約6m,肋距0.6m,一階濕模態頻率約為2.6Hz。對外流場采用DAA2法模擬,對內流場采用聲固耦合法模擬,兩者在非耐壓殼處滿足連續性條件。由沖擊動力學和結構動力學的相關理論可知,舷間流場的存在并未改變細長圓柱殼的縱向體質量和剛度的分布情況,因此對總體低階振動特性不會產生明顯影響。但舷間流場對高頻沖擊波的透射作用會顯著改變結構的局部高頻響應特征,尤其是前期沖擊波作用階段更明顯。

圖6 雙層圓柱殼總體響應計算示意圖Fig.6 Configuration of the total response simulation of double cylindrical structure

圖7 雙層圓柱殼計算模型Fig.7 Calculation models for double cylindrical shell structures
對艦艇在水下爆炸作用下的總體響應已進行了大量的理論和實驗研究,通常認為當氣泡脈動頻率接近或略小于艦艇一階模態頻率,即氣泡脈動周期略大于結構固有周期時,鞭狀總體響應最明顯。氣泡脈動周期可通過庫爾估算公式獲得[1]:

式中:K 為材料常數,對于 TNT,K≈2.11s·kg-1/3·m5/6;W 為裝藥質量,kg;d 為水深,m;10.3m 為大氣壓對應的壓力水頭。
為了深入研究不同工況參數對鞭狀運動的影響,定義周期參數γT=Tb/Ts和爆距參數γR=R/L,作為衡量水下爆炸氣泡誘導鞭狀運動的量綱一參量。Tb為氣泡脈動周期,可通過式(14)計算;Ts為結構一階濕模態對應的固有頻率;R為爆距;L為結構特征長度,對加筋圓柱殼,本文中L取殼體長度。
2.2.1 周期比的影響
首先考慮周期參數對鞭狀運動的影響,工況設置見表1,其中TNT炸藥質量均為1t,爆距均為20m,來分析不同水深d對鞭狀運動的影響。事實上,方位角對結構鞭狀運動的影響:一方面,體現在氣泡向上遷移時不同方位角下的爆距改變趨勢不同;另一方面,則是由于氣泡遷移引起的流場非對稱性導致的。方位角為90°時氣泡向上遷移對結構鞭狀運動響應影響最小。因此,為減少變量數量,本文的分析均是針對方位角為90°的情況。

表1 柱殼總體響應工況設置Table 1 Case parameters for total responses of cylindrical shell structures
由于鞭狀運動主要是一階模態形式下的運動,因此采用船體梁中點撓度來評估圓柱殼鞭狀運動的幅值。考慮到圓柱殼在沖擊載荷下產生的剛體位移,將圓柱殼撓度定義為:

式中:Uc為圓柱殼中部的垂向位移;Ul和Ur分別為圓柱殼兩端的垂向位移。
先以工況1為例,圓柱殼到達鞭狀運動幅值各波峰波谷時的變形和應力云圖如圖8所示。

圖8 工況1鞭狀運動應力云圖Fig.8 Mises stress contours of the whipping response in case 1
在工況1下,質量為1t的TNT在水下60m深處爆炸的氣泡一次脈動周期為610ms,約為圓柱殼結構一階固有周期的1.60倍。從圖8可以看出,工況1的圓柱殼在水下爆炸的前期的沖擊波作用下已被激起了一階鞭狀運動,但幅值不大。而鞭狀運動過程中在240和450ms出現了2個相鄰的波峰,間隔時間較短,3個鞭狀運動幅值沒有明顯增大。

圖9 工況1鞭狀運動幅值和沖擊載荷曲線Fig.9 Curves of whipping response to shock pressure in case 1
圖9為工況1中鞭狀運動幅值和沖擊載荷曲線。其中,綠虛線為沖擊波和氣泡脈動載荷,實線為計算得到的三維圓柱殼結構的鞭狀運動幅值曲線。該工況下結構最大幅值僅有0.1m左右;另外,在第2次鞭狀運動周期出現2個緊鄰的波峰。當沖擊波載荷過后,圓柱殼開始第1次鞭狀運動,并很快達到第1個波峰。在回落過程中,受氣泡脈動載荷負壓區的作用,對圓柱殼結構輸入了較大的能量,使第1個鞭狀運動波谷的幅值大于第1個波峰的幅值,約為0.17m。由于工況1的氣泡脈動周期較長,當圓柱殼結構第2個鞭狀運動峰值出現時,氣泡一次脈動壓力尚未到來,因此在這段時間內,圓柱殼的鞭狀運動主要是由本身的有阻尼震蕩和脈動壓力的負壓區共同主導的,幅值僅有0.10m,較第1個波谷有明顯的減小。氣泡一次脈動載荷峰值在鞭狀運動第2個波峰的回落過程中出現,從而使結構調轉運動方向,形成了另一個波峰。在這種水深較小的工況下,氣泡脈動周期較長,使氣泡一次脈動壓力峰值到來的相對較晚,若恰巧出現在鞭狀運動的回落階段,則會產生類似的鞭狀運動的雙峰現象,即存在2個間隔很小的波峰。在這種情況下,通常不會引起劇烈的鞭狀運動。
對于工況3,1t的TNT在水下80m深處爆炸的氣泡一次脈動周期為490ms,約為圓柱殼結構一階濕模態對應的周期的1.3倍。工況3鞭狀運動幅值和沖擊載荷時歷曲線見圖10。

圖10 工況2鞭狀運動幅值和沖擊載荷曲線Fig.10 Curves of whipping response to shock pressure in case 2

圖11 工況5鞭狀運動幅值和沖擊載荷曲線Fig.11 Curves of whipping response to shock pressure in case 5
從圖10可看出:與工況1相同,在從第1個波峰回落的過程中,由于同時受到氣泡脈動載荷負壓區的能量輸入,使第1個鞭狀運動波谷的幅值大于第1個波峰的幅值;而氣泡一次脈動載荷又恰好位于鞭狀運動的第2次上升階段,同樣以做正功的形式向圓柱殼結構輸入能量,使鞭狀運動幅值繼續增大,達到0.28m;在接下來的第3和第4個周期中,氣泡脈動載荷同樣為系統輸送能量,從而使峰值逐個增大,分別達到0.35和0.41m;而第5個周期時,脈動壓力已經衰減到很小,所做的功不足以抵消整個系統的能量耗散,使鞭狀運動幅值開始回落。
在工況5中,周期參數約為1.1,所對應圓柱殼鞭狀運動最大幅值僅有0.37m,遠小于工況3。圖11中共有2個氣泡脈動壓力峰值出現在鞭狀運動上升階段,鞭狀運動最大運動幅值出現在第3個運動周期。第1個周期的波峰同樣來自沖擊波載荷,氣泡前2次脈動壓力峰值均在鞭狀運動上升階段,使運動幅值逐次增加。而氣泡第3次脈動壓力峰值剛好處于圓柱殼鞭狀運動第3個波谷,此時鞭狀運動速度為零,即此次氣泡脈動載荷基本對系統不做功。同時由于系統能量的耗散,使第4次的運動幅值開始減小。隨后的氣泡脈動壓力峰值開始出現在結構鞭狀運動的下降階段,從而加快了鞭狀運動衰減的速度。從圖12可知,隨著水深的增大,圓柱殼結構鞭狀運動周期逐漸減小。綜合所有工況,分析周期比與鞭狀運動最大幅值Um的關系,將其繪于圖13。

圖12 不同工況下,鞭狀運動時歷曲線Fig.12 Whipping response curves in different cases

圖13 鞭狀運動最大幅值隨周期比的變化Fig.13 The maximum amplitude of the whipping response varied with period ratio
從圖13可知,隨著周期比的增大,鞭狀運動幅值逐漸增大,水深為80~90m,即1.19<γT<1.30時,鞭狀運動幅值達到最大,隨后迅速衰減。綜合以上分析可知,當1.19<γT<1.30時,可是使連續多個氣泡脈動峰值出現在圓柱殼鞭狀運動幅值的上升階段,從而連續向系統輸送能量,達到共振。當γT從1.19減小時,氣泡脈動峰值處于鞭狀運動上升階段的個數逐漸減小,從而使鞭狀運動最大幅值隨之減小。當γT從1.30增大時,氣泡第一次脈動峰值直接后移至鞭狀運動幅值下降階段,對系統的能量輸入變為負,在鞭狀運動的第2個峰值附近緊接著出現第3個峰值,從而使鞭狀運動最大幅值迅速減小。
2.2.2 爆距比的影響
除水深外,爆距也會明顯影響鞭狀運動的形式和幅值。為分析爆距對圓柱殼結構鞭狀運動的影響,進行一系列的計算,工況設置見表2,TNT炸藥的質量均為1t,水深均為80m,氣泡脈動周期為490ms。不同工況下圓柱殼鞭狀運動的幅值曲線見圖14~15。從圖14~15可知:隨著爆距的增大,鞭狀運動幅值迅速減小。其次,爆距越小的工況,鞭狀運動第1個周期越長,如R=10m的工況第1個鞭狀運動周期約為550ms,超過圓柱殼結構固有周期380ms約45%,因此使后續周期相對于其他工況后延了約1/4個周期。而對于R≥20m的工況,鞭狀運動每個周期的長度基本都是均勻的。圖14給出了80m水深下各工況最大鞭狀運動幅值隨爆距的變化趨勢。從圖14可以看出,隨著爆距的增大,圓柱殼結構鞭狀運動幅值迅速減小。這種減小一方面來自于水下爆炸載荷隨爆距的倒數關系,另一方面則是隨著爆距的增大,爆炸載荷在結構縱向的分布逐漸均勻,使結構更多傾向于剛體平動,而非鞭狀運動形式的彎曲響應。

表2 柱殼總體響應工況設置Table 2 Case parameters for total response

圖14 爆距對鞭狀運動幅值的影響Fig.14 Effect of explosion distance on the amplitude of the whipping response

圖15 鞭狀運動最大幅值隨爆距比的變化Fig.15 The maximum amplitude of the whipping response varied with explosion distance ratio
從波動方程出發,推導了DAA2法的后期近似,并結合聲固耦合法初步解決了雙層圓柱殼內流場問題,采用ABAQUS接口耦合實現了圓柱殼水下爆炸沖擊響應的計算分析。通過簡單算例的驗證表明,本文中建立的數值模型可用于圓柱殼水下沖擊的工程和理論研究。據此分析了雙層圓柱殼結構在水下爆炸載荷下的總體響應特性,發現對于本研究的參數范圍,周期比為1.2到1.3之間時,能夠誘發最嚴重的鞭狀運動;當周期比從1.2減小時,鞭狀運動幅值逐漸減小;而當周期比從1.3增大時,鞭狀運動幅值迅速減小;圓柱殼結構鞭狀運動幅值隨爆距的增加呈平方關系減小。
[1]Huang H.Transient interaction of acoustic plane waves with a spherical elastic shell[J].The Journal of the Acoustical Society of America,1969,45(3):661-670.
[2]Huang H.An exact Analysis of the transient interaction of acoustic plane waves with a cylindrical elastic shell[J].Journal of Applied Mechanics,1970,37:1091-1099.
[3]Huang H,Lu Y P,Wang Y F.Transient interaction of spherical acoustic waves and a spherical elastic shell[J].Journal of Applied Mechanics,1971,38:71-74.
[4]Geers T L.Residual potential and approximate methods for three-dimensional fluid-structure interaction problems[J].Journal of the Acoustical Society of America,1971,49(5B):1505-1510.
[5]Geers T L.Doubly asymptotic approximations for transient motions of submerged structures[J].Journal of the A-coustical Society of America,1978,64(5):1500-1508.
[6]劉建湖.艦船非接觸水下爆炸動力學的理論與應用[D].無錫:中國船舶科學技術研究所,2002:89-112.
[7]孫士麗.瞬態載荷作用下大幅運動航行體流固耦合方法及應用研究[D].哈爾濱:哈爾濱工程大學,2009:15-34.
[8]王詩平,孫士麗,張阿漫,等.沖擊波和氣泡作用下艦船結構動態響應的數值模擬[J].爆炸與沖擊,2011,31(4):367-372.Wang Shi-ping,Sun Shi-li,Zhang A-man,et al.Numerical simulation of dynamic response of warship structures subjected to underwater explosion shockwaves and bubbles[J].Explosion and Shock Waves,2011,31(4):367-372.
[9]張弩,宗智,張文鵬,等.基于雙漸進方法的水下爆炸氣泡載荷作用下艦船的動態響應分析[J].振動與沖擊,2012,32(23):50-56.Zhang Nu,Zong Zhi,Zhang Wen-peng,et al.Dynamic response of a ship hull structure subjected to an underwater explosion bubble based on doubly asymptotic approximation method[J].Journal of Vibration and Shock,2012,32(23):50-56.
[10]董海,劉建湖,吳有生.水下爆炸氣泡脈動作用下細長加筋圓柱殼的鞭狀響應分析[J].船舶力學,2007,11(2):250-258.Dong Hai,Liu Jian-hu,Wu You-sheng.Whipping response analysis of slender stiffened cylindrical shell subjected to underwater explosion with bubble pulse[J].Journal of Ship Mechanics,2007,11(2):250-258.
[11]Geers T L,Hunter K S.An integrated wave-effects model for an underwater explosion bubble[J].The Journal of the Acoustical Society of America,2002,111(4):1584-1601.
[12]Geers T L,Park C K.Optimization of the G & H bubble model[J].Shock and Vibration,2005,12(1):3-8.
[13]Sprague M A,Geers T L.Response of empty and fluid-filled submerged spherical shells to plane and spherical,step-exponential acoustic waves[J].Shock and Vibration,1999,6:147-157.
[14]王詩平,孫士麗,張阿漫,等.可壓縮流場中氣泡脈動數值模擬[J].力學學報,2012,44(3):513-519.Wang Shi-ping,Sun Shi-li,Zhang A-man,et al.Numerical simulation of bubble dynamics in compressible fluid[J].Journal of Theoretical and Applied Mechanics,2012,44(3):513-519.
[15]Geers T L,Zhang P Z.Doubly asymptotic approximations for submerged structures with internal fluid volumes:Formulation[J].Journal of Applied Mechanics,1994,61(4):893-899.
[16]Geers T L,Zhang P Z.Doubly asymptotic approximations for submerged structures with internal fluid volumes:Evaluation[J].Journal of Applied Mechanics,1994,61(4):900-906.
[17]Wilkerson S A.Elastic whipping response of ships to an underwater explosion loading[D].Washington:George Washington University,1985:57-72.
[18]Zhang P,Geers T L.Excitation of a fluid-filled,submerged elastic shell by a transient acoustic wave[J].Journal of the Acoustical Society of America,1993,93(2):696-705.