水下沖擊波作用下結構損傷的數值預報誤差分析

張藝凡 宗 智 張文鵬

大連理工大學 工業裝備結構分析國家重點實驗室 運載工程與力學學部，遼寧 大連 116024

水下沖擊波作用下結構損傷的數值預報誤差分析

張藝凡 宗 智 張文鵬

大連理工大學 工業裝備結構分析國家重點實驗室 運載工程與力學學部，遼寧 大連 116024

準確模擬結構在水下爆炸沖擊波作用下的損傷是水下爆炸數值仿真中的難點，有必要研究數值預報誤差及其產生的原因。采用ABAQUS提供的聲固耦合方法，研究方板模型水下爆炸數值計算的誤差。將幾何模型劃分為3個不同的有限元模型，分析了網格大小對計算結果的影響，以及聲固耦合方法的計算誤差。通過方板模型仿真、實驗室水池試驗和方板實船水下爆炸實驗結果的比對分析表明：數值計算結果與實驗值間的誤差約在30%以內。總波和散波兩種流固耦合計算方法的誤差及其產生原因分析表明總波計算公式存在夸大空化對結構影響的可能。

水下沖擊波；塑性變形；誤差分析；流固耦合

1 引言

水下爆炸會產生很大的沖擊載荷使附近的水中結構物產生損傷。因此，從二戰以來，水下爆炸產生的結構損傷問題成為各國海軍關注的焦點問題。

在過去60年間，對于結構在水下沖擊波作用下的彈性響應［1］，出現了許多成熟的理論和計算方法。比如，Huang［2］提出的聲學耦合理論方法，Geers［3－4］和 Deruntz［5］提出的 DAA 方法都可 以精確地預測水下爆炸所引起的結構彈性響應。通過商業軟件，比如ABAQUS和LS－DYNA也可以實現結構在水下爆炸載荷作用下動態彈性響應的計算。

然而，一旦結構變形超過了彈性進入塑性，無論是理論值還是數值模擬往往與實驗結果吻合得不好。早在二戰時期，Kirkwood［6］就研究過圓板受水下沖擊波載荷的塑形變形，提出了結構的損傷機理，但是計算得到的結果遠遠小于實驗值。Zong［7－8］基于 Kirkwood 的理論，對損傷機理進行了改進，大大提高了理論預測的精度。吳成［9］研究了方形板的損傷。朱錫等［10］研究了板架在水下爆炸作用下結構損傷。但是，系統地用數值方法預報結構在水下沖擊波作用下產生的損傷的研究還不夠。特別是對于任意給定的一個結構，使用數值方法對其損傷（或者塑性變形）進行預報時，往往很難得到令人滿意的結果。

其原因在于其中涉及的流固耦合作用尚不明確。本文使用聲固耦合方法，通過對兩個實驗結果和一個實船試驗結果的對比分析，考察數值模擬預報水下爆炸造成的損傷方面存在的誤差，希望對數值方法的使用提供有益的參考。

2 計算方法簡介

本文在進行數值仿真計算時，計算方法采用ABAQUS提供的聲固耦合方法。關于這方面的介紹在許多論文中均有涉及，對此不再贅述。需要強調的是ABAQUS中提供的兩種載荷計算方法：不考慮空化的散波公式法和可以考慮空化的總波公式法。在線性假設下，水中沖擊波可以分為三個部分：不考慮結構存在時，由爆源傳播到空間某一點的壓力，稱為沖擊波；假設結構是剛體，入射沖擊波打到結構上產生的反射波；輻射波是指在沒有入射波和反射波時，結構在水中振動產生的波。在這種假設下來建立求解方法的就是散波公式法，如圖1所示。

圖1 水中沖擊波示意圖Fig.1 Schematic of underwater shock

當結構的變形速度較大時，結構前的壓力就會降到蒸氣壓力以下，從而發生空化現象。

圖2 水中空化示意圖Fig.2 Schematic ofwater cavitation

考慮到這種空化現象的計算方法就是總波公式法?？栈F象改變了附近流場的壓力分布，從而使得水域的加載方式有所不同。

通過這兩種方法計算得到載荷，就可以通過耦合的方法來計算結構的損傷。流固耦合計算方法在流體區域求壓力時，滿足速度邊界條件；而在求解結構位移時，沖擊波壓力作為力的邊界條件加入到結構的控制方程中。

水中沖擊波入射波的計算在遠場往往采用如下經驗公式：

對于TNT炸藥峰值Pm和衰減常數τ采用如下公式獲得：

式中，G為炸藥質量，kg；R為爆炸距離，m。該公式直接由實驗值獲得，具有很高的精度。同時，結構在空氣中的損傷計算也很成熟，具有較高的精度。因此，水中沖擊波造成結構損傷的數值模擬計算誤差最有可能來自于載荷計算方法和流固耦合方法中的誤差。

本文將通過3個具體的例子來考慮數值模擬的誤差問題。

3 數值模擬的收斂性分析

數值模擬計算必須進行收斂性驗證。收斂性是指數值模型計算網格越細，數值結果越接近真實值。一個數值方法只有具有收斂性才可以使用。

為考察收斂性，我們取一塊方板模型進行水下爆炸數值計算，結構參數和計算工況如表1所示。

表1 計算模型尺寸參數Tab.1 Parameters of calculation model

對于數值計算問題，建立的有限元模型會對最終的結果產生較大影響。本文計算模型中，主要影響參數為流場大小和有限元網格尺寸。流場大小的影響分以下3個方面：重力影響、阻尼影響和慣性影響。通常比較關心的是慣性影響，此時部分流場將參與船體的總振動，該部分流場質量稱之為附連水質量，與結構本身質量為同一量級，在計算時附連水質量不可忽視。

在進行水下爆炸數值模擬分析時，為獲得較準確的模擬數據，必須保證流場區域足夠大。但是，由于現有計算軟件對數值模型單元數的限制，流場區域不可能取為無限大。表2［11］給出了隨著流場尺寸與結構尺寸比值的變化，相應附加質量率的變化值。本文中，流場尺寸表示流場邊緣到板中心的長度，結構尺寸表示板邊緣至板中心的長度，附加質量率表示無限流場的附加水質量與有限元模型流場的附加水質量的比值。

表2 有限元模型中流場大小與附加質量率的關系Tab.2 The relationship of add ed mass ratio and finite element size of flow field

為提高精度并合理安排計算周期，本文中流場尺寸與結構尺寸之比取6。幾何模型如圖3所示，其中，z方向為垂直于板面的方向，x方向為沿板長度方向。

圖3 幾何模型示意圖Fig.3 Geometrymodel

為分析網格大小對計算結果的影響，對上述幾何模型劃分為3個不同的有限元模型進行計算。模型中板結構采用殼單元，流場單元采用六面體聲學單元，流場大小為4.2m×2.4 m×2m，表3中給出了流場單元的尺寸，流場網格大小如圖4所示。板單元的大小與相對應的流場單元大小相同。

表3 有限元模型流場單元尺寸Tab.3 Size of finite elements of flow field

圖4 不同尺寸流場網格示意圖Fig.4 Different sizes of finite element flow field

數值計算時采用總波公式（在ABAQUS中通過關鍵字Totalwave實現）進行顯式動態分析計算。板材料為鋼，在計算中材料考慮率相關的影響，采用目前普遍接受的Cowper－Symonds提出的本構方程形式：

水域與結構接觸面設置聲固耦合條件，以傳遞壓力和位移。關鍵字設置為：

其中，Water－Surface為與結構接觸的水域，Structure－Surface為與水域接觸的結構面。

為消除水域邊界的影響在模型邊界上定義無反射條件，表示沖擊波可以穿過該表面而不發生反射。本文計算模型的水域邊界形狀均為平面，分別對模型邊界面采用如下關鍵字進行設置：

*SIMPEDANCE，NONREFLECTING＝PLANAR

Water－1

其中Water－1表示水域模型邊界面。通過如上所述的設置后，提交任務進行計算，施加的沖擊波載荷如圖5所示。

圖5 沖擊波載荷示意圖Fig.5 The curve of shock pressure

計算得到3個模型中心點的撓度，如圖6所示。

圖6板中心點撓度示意圖Fig.6 The displacement of the center of plate

通過對計算結果的分析可以看出，模型model－1 與 model－2 之間的誤差為 1.65%，model－2 與model－3之間的誤差在 1.24%。因此，在本算例中，當流場網格密度達到11 161個／m3時，數值計算即可認為收斂。計算時，為排除計算時間步長對結果的影響，采用ABAQUS提供的自動匹配時間步長設置以保證計算精度。

4 實驗室試驗的對比研究

首先考慮兩塊圓板的損傷，通過對具體實驗工況的數值模擬，將計算得到的結果與Kirkwood的實驗結果進行比較。這兩項實驗在水池中進行，炸藥當量很小。板的尺寸和爆炸工況如表4所示。

表4 計算模型尺寸參數Tab.4 Parameters of calculation model（Steel 1）

對于第一塊板進行計算時共劃分了105個六面體流場單元和300個殼結構單元。水中聲速取為1 500m／s。分別采用散波公式和總波公式計算Steel1的損傷情況。數值計算模型如圖7所示。

圖7 有限元模型Fig.7 Finite elementmodel

載荷的計算采用ABAQUS提供的水下爆炸載荷計算公式進行加載。載荷參數具體設置如下：

將數值計算得到損傷結果與實驗值進行對比，如圖8所示。

圖8 Steel1中心點位移示意圖Fig.8 The displacement of the center of Steel 1

圖為Steel 1中心點的塑性位移情況。

從圖8中可以看出，散波計算的誤差比較小，為16.8%；總波計算的誤差為28.9%，比散波高一倍左右。在ABAQUS提供的兩種聲固耦合計算方法中，散波公式將水壓力簡化為線性問題；而對于空化這類非線性問題只能采用總波公式［12］。因此，通常分析認為總波計算能夠考慮空化現象，理應更加準確，但在本例中總波公式計算產生的誤差更大，分析其原因，很可能在于不恰當的考慮了空化效應。

為說明這一點，觀察Steel 1附近流場區域空化現象的演化過程。如圖9所示，從圖中可以清楚的看到沖擊波的傳播路徑。

為更清楚的觀察到空化區域的擴展，對流場壓力的顯示數值進行調整，僅取板周圍的流場壓力云圖進行分析，如圖10所示。當計算至0.08 ms時，與板直接接觸的部分流場出現空化現象，隨著計算時間的推移空化區域迅速擴大。

圖9 Steel 1流場壓力云圖Fig.9 Plot ofwater pressure of Steel 1

圖10 Steel 1空化區域云圖Fig.10 Plotof cavitations of Steel 1

圖中淺灰色的區域為空化區域。通過簡單的計算可以得到，流場內空化區域的傳播速度約為1 500 m／s，遠遠大于實際中空化的傳播速度（約10～20m／s），空化效應對結構的影響被放大。

另外，我們還計算了另外一塊圓形板Steel 2在水下沖擊波作用下的響應情況，幾何參數和計算工況見表5，圖11分別給出了總波和散波計算所得到的塑性變形時程曲線。

表5 計算模型尺寸參數Tab.5 Parameters of calculation model（Steel 2）

同樣的散波計算的誤差為6.33%，比總波計算的誤差13.8%小的多，其原因是由于總波計算法過分地夸大了空化的作用。

比較圖8和圖11可以看出，Steel 2的計算結果誤差比較小。產生的原因可能為Steel 2的相對變形（變形量 ／板的直徑）近似為 0.55，而 Steel 1的相對變形量為0.11。由于Steel 2的相對變形大，塑性變形成為整個變形過程中的主要因素，比較容易捕捉；Steel 1的相對變形較小，彈性變形的影響變得不可忽視。

圖11 Steel 2中心點位移示意圖Fig.11 The displacement of the cent er of Steel 2

5 實船試驗對比研究

另外，對一塊實船上的方板以及該實船的水下爆炸作用下的試驗結果與數值計算進行了比較。板的幾何尺寸和爆況如表6所示。

表6 計算模型尺寸參數Tab.6 Parameters of calculation model（Ship 1）

計算共劃分了260個板單元，流場大小為4.8m×12.78m×6m， 流場單元總數為 556 800。圖12所示為總波計算和散波計算的損傷模擬結果。

圖12 Ship 1中心點位移示意圖Fig.12 The displacement of the cent er of Ship 1

Ship 1總波計算和散波計算結果與實驗值的誤差分別為25.05%和22.73%。兩者數值相近。說明在實際應用中，總波和散波公式的精度區別不大。而在實驗室條件下，總波的公式結果都比散波大（Steel 1和Steel 2），說明總波計算公式有夸大空化的可能，但這一結論還需更多的研究來驗證。

6 參數分析

為比較各個計算參數變化對結果的影響，對屈服極限、藥量和爆炸距離分別進行了參數分析。對Ship1模型分別計算表6所給出的6種工況下的塑性位移情況。

計算得到板的塑性位移的變化量如表7所示。

表7 計算工況Tab.7 Calculation condition s

表8 位移變化率計算結果Tab.8 Calculation results of various displacement ratio

通過表格中的數據分析我們可以看出，當計算參數值各變化10%時，對計算結果均會產生一定的影響，其中爆炸距離變化對計算結果產生的影響最大。在實船實驗中，材料、爆炸距離等參數都會產生一定的測量誤差，使得數值計算的結果與實驗值的誤差增大。

7 結 論

1）對于實驗室試驗，由于可控性好，散波公式法精度明顯高于總波公式計算法，說明總波公式可能夸大了空化效應的影響。

2）根據實船試驗結果對比分析，總波公式和散波公式計算誤差相差不大。

3）水下爆炸損傷的數值預報是一個非常復雜的問題，在現有的計算水平下計算誤差達到30%還是令人滿意的。從這個意義上講，本文所討論的聲固耦合計算方法可以應用于實際工程預報中。同時，為了進一步提高精度，還需要做更多的計算方法研究工作。

［1］ JONESN.Structural impact［M］.Cambridge： Cambridge University Press，1989：333－383.

［2］ HUANG H，Wang Y F.Transient interaction of spherical acoustic waves with a cylindrical elastic shell ［J］.Journal of the Acoustical Society of America， 1970，48：228－235.

［3］ GEERS T L.Residual potential and approximatemethods for three－dimensional fluid－structural interaction problem［J］.Journal of Acoustic Society of America，1971，49（5B）：1505－1510.

［4］ GEERS T L.Doubly Asymptotically approximations for transient analysis of submerged structures ［J］.Journal of Acoustic Society of America，1978，64（5）：1500－1508.

［5］ DERUNTZ J A.The underwater shock analysis code and its applications［C］//Proceedings of 60th Shock and Vibration Symposium，Virgina，USA，1989.

［6］ KIRKWOOD JG，RICHARDSON JM.The plastic deformation of circular diaphragms under dynamic loading by an underwater explosion wave［R］.OSRD 4200，1944.

［7］ ZONG Z，LAM K Y， LIU G R.Probabilistic risk p rediction of submarine pipelines subjected to underwater shock［J］.Journal of Offshore Mechanics and Arctic Engineering，1999，121（4）：251－254.

［8］ ZONG Z，LAM K Y.Viscoplastic response of a circular plate to an underwater explosion shock［J］.Acta Mechanica，2001，148（1/4）：93－104.

［9］ 吳成，金儼，李華新.固支方板對水中爆炸作用的動態響應研究［J］.高壓物理學報，2003，17（4）：275－282.

WU C，JIN Y，LIH X.A study on square plate dynamic response under underwater explosion ［J］.Chinese Journal of High Pressure Physics，2003，17（4）：275－282.

［10］ 朱錫，白雪飛，黃若波，等.船體板架在水下接觸爆炸作用下的破口試驗［J］.中國造船，2003，44（1）：46－52.

ZHU X，BAI X F，HUANG R B，et al.Crevasse experiment research of platemembrance in vessels subjected to underwater contact explosion ［J］.Shipbuilding of China，2003，44（1）：46－52.

［11］ ABAQUSusermanual ver.6.9EF［S］.ABAQUS， Inc..

［12］ ABAQUS analysis user’smanual，23.3.1 ［S］.ABAQUS，Inc.

Error Analysis of Numerical Prediction of Structure Damages Subjected to Underwater Shock

Zhang Yi－fan Zong Zhi ZhangWen－peng
Faculty of Vehicle Engineering and Mechanics，State Key Laboratory of Structural Analysis for Industrial Equipment， Dalian University of Technology，Dalian 116024，China

The fas t development of computing technology saw wide applications of numerical simulations in underwater explosion in recent years，but it′s still difficult to numerically and exactly predict the damages of amarine structure subjected to underwater shock.It is true that the numerical results and experimental results of a certain underwater explosion problems are in good agreement，but in themost general cases， the problem that how good the numerical results agree with experimental results remain unsolved.In this paper， through three examples， we try to reveal the numerical errors compared with experimental results.All of the three examples show that the numerical errors can be controlled within 30 percent.Two numerical methods， that is， total wave and scatter wave formulations， are employed in this paper.It seems that the totalwavemethod overestimates the effect of cavitations.

underwater shock； plastic deformation； error analysis； fluid structure interaction

U661.41

A

1673－3185（2011）06－38－07

10.3969/j.issn.1673-3185.2011.06.008

2010-08-02

國家自然科學基金（50921001）；973項目（2010CB832700）

張藝凡（1986－），女，碩士研究生。研究方向：流固耦合，水下爆炸。E-mail：yifanzhang929＠gmail.com

宗 智（1964－ ），男，教授，博士生導師。 研究方向：水下爆炸，流固耦合。E-mail：zongzhi＠dlut.edu.cn