水下爆炸數值模擬聯合算法求解結構穩態響應

李亞偉， 杜志鵬， 李 俊， 李曉彬， 張 磊

(1.武漢理工大學 交通學院， 湖北 武漢 430063； 2.海軍研究院， 北京 100161)

0 引 言

加筋板架是船體結構最基本的組成部分[1]，加筋板結構的設計與水面艦船的好壞有直接關系，關于加筋板架的研究[2]至關重要。JIANG等[3]考慮流固耦合和應變率效應，提出一種簡化的理論分析方法，用于計算板架水下爆炸動塑性響應。吳有生等[4]研究在爆炸載荷作用下艦船板架的變形與破損，并在試驗中得到了較好的驗證。GUPTA等[5]對船體板架結構用ABAQUS/Explicit進行仿真，計算板架彈塑性變形與破壞。牟金磊等[6]用MSC.Dytran對不同藥量、不同爆距下的板架毀傷模式進行仿真研究，并進行試驗驗證。王龍侃等[7]用ABAQUS/Explicit對加筋板架進行水下爆炸響應仿真計算，得到加筋板架的變形與破損，并對比仿真結果與試驗結果。綜合上述分析發現，以往基于ABAQUS的水下爆炸數值模擬方法都是基于瞬態顯式求解器設置的，而顯式求解器在求解加筋板架的穩態響應方面存在一定的缺陷，即加筋板架變形在最終塑性變形位置處無限振蕩，從而可能導致仿真結果被低估或誤判(取值誤差)。鑒于此，提出采用“瞬態顯式-穩態隱式”聯合數值仿真方法求解加筋板架動力響應。該方法基于ABAQUS/Explicit-Standard算法，針對瞬態顯式和穩態隱式靜力平衡迭代建立聯合仿真分析方法。該方法在求解結構穩態響應時能夠有效避免瞬態顯式分析結果中存在的缺陷問題。

1 穩態響應

在水下爆炸沖擊作用下，加筋板架遭受沖擊載荷后產生彈塑性變形，并在彈塑性變形最大時發生回彈，然后在最終塑性變形位置振蕩并無限趨近于最終塑性變形。從變形回彈到無限趨近于最終塑性變形的過程稱為穩態響應，具體如圖1所示。

圖1 水下爆炸加筋板架變形時歷曲線

1.1 瞬態顯式求解穩態響應

采用瞬態顯式中心差分法對穩態響應求解，將板架結構的穩態響應處理成自由振動問題，引入適當的阻尼項來消除引起的動態振蕩，使板架在一段時間后達到穩定狀態[8]。施加黏性阻尼系統的動態平衡方程為

(1)

式中：M、C、K分別為結構質量、阻尼、剛度矩陣；U為節點位移向量。對于一個簡單地施加臨界阻尼的振蕩系統，系統位移為

u=αe-ω t+βte-ω t

(2)

式中：ω為自然角頻率；α、β為與邊界條件、幾何、材料性質有關的系數；t為時刻。

1.2 隱式算法求解穩態響應

采用隱式求解算法計算板架的穩態響應過程。此時板架所受載荷已遠小于爆炸沖擊載荷，幾乎可忽略不計。因此，隱式算法將板架穩態響應看作是一個擬靜態的響應過程，采用增量形式的虛功方程求解整個板架的非線性方程組，并采用靜力平衡迭代的方法求解板架的穩態響應，即隱式算法求解穩態響應的過程。

隱式算法計算得到的穩態響應結果可靠性強、計算效率高，而采用瞬態顯式算法求解穩態響應必須引入合適的阻尼系數，且引入的系數又缺乏明確的物理意義。因此，本文提出一種聯合“瞬態顯式-穩態隱式”算法的有限元仿真分析方法，該方法基于ABAQUS中的瞬態顯式動力求解與隱式靜力分析。

2 加筋板架瞬態-穩態響應分析

2.1 加筋板架模型

根據文獻[1]和文獻[7]中的加筋板架試驗模型，分別對加筋板架進行仿真研究，并將仿真結果與試驗結果進行對比驗證。加筋板架模型如圖2所示。所有模型材料均采用Q235鋼。采用有限元軟件ANSA對加筋板架及流場水域進行前處理。其中，模型2加筋板架及流場水域的有限元模型如圖3所示。

在仿真分析時，炸藥選用TNT裝藥，炸藥均放置在加筋板架正下方，加筋板架邊界條件為四周剛性固定，試驗工況如表1所示。采用Geers-Hunter的水下爆炸雙重漸近模擬計算水下爆炸載荷。由于只關注沖擊波載荷對結構的影響，因此只需將爆炸載荷加載時間設置為小于氣泡脈動周期即可。

表1 2種模型水下爆炸計算工況

圖2 加筋板架模型

圖3 模型2加筋板架及流場水域的有限元模型

2.2 瞬態顯式仿真結果與試驗結果分析

采用瞬態顯式法對船體板架結構有限元模型進行計算，計算結果如圖4所示，為模型1與模型2在特定時刻的變形云圖。由圖4可知：在水下爆炸載荷作用下，加筋板架發生大幅度的塑性變形，且不存在缺口，即加筋板架的破壞模式為塑性大變形。

圖4 加筋板架0.15 s時刻爆炸沖擊變形云圖

圖5所示為模型2在試驗中加筋板架的整體塑性變形圖，其最終變形撓度為15.79 mm。圖4(b)為模型2在特定時刻的數值模擬塑性變形云圖，其最終變形撓度為17.50 mm。

圖5 模型2加筋板架0.15 s時刻整體塑性變形圖

根據對數值仿真結果的分析，可以得出這2種加筋板架模型的變形值。表2給出了2種模型其加筋板架的仿真值與試驗實測值。從表2可以發現，這兩種模型的仿真誤差都在10%以上。

表2 2種模型的加筋板架仿真值與試驗值對比

2.3 瞬態顯示-穩態隱式聯合仿真結果與試驗結果分析

(1) ABAQUS隱式算法對加筋板架穩態響應求解具體分析步驟如下。

① 首先，需導入沖擊變形后的加筋板架。在環境欄Module后面選擇part，進入part分析模塊，并導入瞬態計算完成的.odb文件，選擇分析步中的最后一個增量步，即可導入變形的加筋板架。

② 裝配加筋板架并創建分析步。環境欄切換至Assembly模塊，將加筋板架加載到裝配體，然后切換到Step模塊，打開非線性，因為內應力突然釋放，為穩定收斂，設置自動穩定因子，分析步初始增量尺寸值可以設為0.1。

③ 施加邊界條件并引入瞬態變形結果。切換至Load模塊，對沿x方向的邊界施加關于與x軸(坐標軸1)垂直的平面反對稱邊界條件，對沿y方向平行的邊界施加關于與y軸(坐標軸2)垂直的平面反對稱邊界條件。在菜單欄Predefined Field中創建一個初始狀態，輸入之前瞬態變形生成的job文件，則瞬態變形后的加筋板架(變形網格、硬化材料屬性、應力應變等)會替換本模型中的加筋板架。

④ 創建并提交分析。在求解結束后，即可進入后處理，查看穩態變形結果[9-12]。

(2) 采用瞬態顯式-穩態隱式聯合算法計算加筋板架穩態響應。

以模型2的仿真分析為例，采用聯合算法計算模型2加筋板架的穩態響應問題。具體方法如上所述：首先將瞬態顯式算法計算完成的結果文件導入ABAQUS中；然后創建并設置分析步，分析步采用穩態隱式靜力分析；最后提交作業分析，并查看最終仿真結果。

圖6(a)所示為模型2在0.02 s時刻的塑性變形云圖，其最終變形撓度為17.50 mm；圖6(b)為加筋板架在瞬態顯式計算完成后，采用穩態隱式算法計算完成后的穩態響應變形云圖。穩態前的位移參考值被重置為0，故穩態分析的位移即為穩態變形值。由圖6(b)可知，穩態最大變形量為0.68 mm。

圖6 水下爆炸作用下加筋板架瞬態-穩態響應變形云圖(模型2)

為更清晰地觀察加筋板架的撓度變化，圖7(a)與圖7(b)分別給出加筋板架采用“瞬態顯式-穩態隱式”聯合算法和單獨采用瞬態顯式方法的加筋板架中心點處垂向位移時歷曲線對比圖。

圖7 采用兩種不同算法的加筋板架垂向位移時歷曲線(模型2)

由圖7(a)和圖7(b)可知，單獨采用瞬態顯式算法得到的仿真結果中可能存在取值誤差。“瞬態顯式-穩態隱式”聯合算法仿真結果不存在取值誤差缺陷。表3給出了聯合算法的仿真結果與試驗結果的對比分析值。

表3 “瞬態顯式-穩態隱式”聯合算法加筋板架仿真值與試驗值對比

由表3可知：采用聯合“瞬態顯式-穩態隱式”算法得到的仿真結果與試驗值的誤差約6.5%，與單獨采用瞬態顯式算法的仿真結果相比，誤差降低約39.8%，說明在采用聯合“瞬態顯式-穩態隱式”算法仿真計算時，在排除取值誤差的缺陷后，得到的加筋板架穩態變形的仿真結果更接近試驗值。

3 結 論

采用ABAQUS“瞬態顯式-穩態隱式”聯合算法，計算船體加筋板架在水下爆炸作用下的穩態響應，求解加筋板架在水下爆炸沖擊載荷下的穩態變形特性，并將數值仿真結果與試驗結果進行對比分析，結論如下：

(1) 采用“瞬態顯式-穩態隱式”聯合算法得到的結果與試驗結果相差在10%以內，說明該算法的仿真方法可用于水下爆炸作用下結構的動力響應數值模擬分析。

(2) 采用聯合“瞬態顯式-穩態隱式”算法仿真得到的穩態變形值比采用瞬態顯示算法更接近試驗值。由于聯合算法在穩態響應分析時釋放了瞬態響應階段的殘余應力和應變，因此此時的穩態變形更接近試驗值。這說明，在瞬態響應分析的基礎上，穩態隱式靜力平衡分析能夠提高加筋板架數值仿真的變形精度。