水下爆炸數(shù)值模擬聯(lián)合算法求解結(jié)構(gòu)穩(wěn)態(tài)響應(yīng)

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

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

0 引 言

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

1 穩(wěn)態(tài)響應(yīng)

在水下爆炸沖擊作用下，加筋板架遭受沖擊載荷后產(chǎn)生彈塑性變形，并在彈塑性變形最大時(shí)發(fā)生回彈，然后在最終塑性變形位置振蕩并無(wú)限趨近于最終塑性變形。從變形回彈到無(wú)限趨近于最終塑性變形的過(guò)程稱為穩(wěn)態(tài)響應(yīng)，具體如圖1所示。

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

1.1 瞬態(tài)顯式求解穩(wěn)態(tài)響應(yīng)

采用瞬態(tài)顯式中心差分法對(duì)穩(wěn)態(tài)響應(yīng)求解，將板架結(jié)構(gòu)的穩(wěn)態(tài)響應(yīng)處理成自由振動(dòng)問(wèn)題，引入適當(dāng)?shù)淖枘犴?xiàng)來(lái)消除引起的動(dòng)態(tài)振蕩，使板架在一段時(shí)間后達(dá)到穩(wěn)定狀態(tài)[8]。施加黏性阻尼系統(tǒng)的動(dòng)態(tài)平衡方程為

(1)

式中：M、C、K分別為結(jié)構(gòu)質(zhì)量、阻尼、剛度矩陣；U為節(jié)點(diǎn)位移向量。對(duì)于一個(gè)簡(jiǎn)單地施加臨界阻尼的振蕩系統(tǒng)，系統(tǒng)位移為

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

(2)

式中：ω為自然角頻率；α、β為與邊界條件、幾何、材料性質(zhì)有關(guān)的系數(shù)；t為時(shí)刻。

1.2 隱式算法求解穩(wěn)態(tài)響應(yīng)

采用隱式求解算法計(jì)算板架的穩(wěn)態(tài)響應(yīng)過(guò)程。此時(shí)板架所受載荷已遠(yuǎn)小于爆炸沖擊載荷，幾乎可忽略不計(jì)。因此，隱式算法將板架穩(wěn)態(tài)響應(yīng)看作是一個(gè)擬靜態(tài)的響應(yīng)過(guò)程，采用增量形式的虛功方程求解整個(gè)板架的非線性方程組，并采用靜力平衡迭代的方法求解板架的穩(wěn)態(tài)響應(yīng)，即隱式算法求解穩(wěn)態(tài)響應(yīng)的過(guò)程。

隱式算法計(jì)算得到的穩(wěn)態(tài)響應(yīng)結(jié)果可靠性強(qiáng)、計(jì)算效率高，而采用瞬態(tài)顯式算法求解穩(wěn)態(tài)響應(yīng)必須引入合適的阻尼系數(shù)，且引入的系數(shù)又缺乏明確的物理意義。因此，本文提出一種聯(lián)合“瞬態(tài)顯式-穩(wěn)態(tài)隱式”算法的有限元仿真分析方法，該方法基于ABAQUS中的瞬態(tài)顯式動(dòng)力求解與隱式靜力分析。

2 加筋板架瞬態(tài)-穩(wěn)態(tài)響應(yīng)分析

2.1 加筋板架模型

根據(jù)文獻(xiàn)[1]和文獻(xiàn)[7]中的加筋板架試驗(yàn)?zāi)Ｐ?，分別對(duì)加筋板架進(jìn)行仿真研究，并將仿真結(jié)果與試驗(yàn)結(jié)果進(jìn)行對(duì)比驗(yàn)證。加筋板架模型如圖2所示。所有模型材料均采用Q235鋼。采用有限元軟件ANSA對(duì)加筋板架及流場(chǎng)水域進(jìn)行前處理。其中，模型2加筋板架及流場(chǎng)水域的有限元模型如圖3所示。

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

表1 2種模型水下爆炸計(jì)算工況

圖2 加筋板架模型

圖3 模型2加筋板架及流場(chǎng)水域的有限元模型

2.2 瞬態(tài)顯式仿真結(jié)果與試驗(yàn)結(jié)果分析

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

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

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

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

根據(jù)對(duì)數(shù)值仿真結(jié)果的分析，可以得出這2種加筋板架模型的變形值。表2給出了2種模型其加筋板架的仿真值與試驗(yàn)實(shí)測(cè)值。從表2可以發(fā)現(xiàn)，這兩種模型的仿真誤差都在10%以上。

表2 2種模型的加筋板架仿真值與試驗(yàn)值對(duì)比

2.3 瞬態(tài)顯示-穩(wěn)態(tài)隱式聯(lián)合仿真結(jié)果與試驗(yàn)結(jié)果分析

(1) ABAQUS隱式算法對(duì)加筋板架穩(wěn)態(tài)響應(yīng)求解具體分析步驟如下。

① 首先，需導(dǎo)入沖擊變形后的加筋板架。在環(huán)境欄Module后面選擇part，進(jìn)入part分析模塊，并導(dǎo)入瞬態(tài)計(jì)算完成的.odb文件，選擇分析步中的最后一個(gè)增量步，即可導(dǎo)入變形的加筋板架。

② 裝配加筋板架并創(chuàng)建分析步。環(huán)境欄切換至Assembly模塊，將加筋板架加載到裝配體，然后切換到Step模塊，打開(kāi)非線性，因?yàn)閮?nèi)應(yīng)力突然釋放，為穩(wěn)定收斂，設(shè)置自動(dòng)穩(wěn)定因子，分析步初始增量尺寸值可以設(shè)為0.1。

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

④ 創(chuàng)建并提交分析。在求解結(jié)束后，即可進(jìn)入后處理，查看穩(wěn)態(tài)變形結(jié)果[9-12]。

(2) 采用瞬態(tài)顯式-穩(wěn)態(tài)隱式聯(lián)合算法計(jì)算加筋板架穩(wěn)態(tài)響應(yīng)。

以模型2的仿真分析為例，采用聯(lián)合算法計(jì)算模型2加筋板架的穩(wěn)態(tài)響應(yīng)問(wèn)題。具體方法如上所述：首先將瞬態(tài)顯式算法計(jì)算完成的結(jié)果文件導(dǎo)入ABAQUS中；然后創(chuàng)建并設(shè)置分析步，分析步采用穩(wěn)態(tài)隱式靜力分析；最后提交作業(yè)分析，并查看最終仿真結(jié)果。

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

圖6 水下爆炸作用下加筋板架瞬態(tài)-穩(wěn)態(tài)響應(yīng)變形云圖(模型2)

為更清晰地觀察加筋板架的撓度變化，圖7(a)與圖7(b)分別給出加筋板架采用“瞬態(tài)顯式-穩(wěn)態(tài)隱式”聯(lián)合算法和單獨(dú)采用瞬態(tài)顯式方法的加筋板架中心點(diǎn)處垂向位移時(shí)歷曲線對(duì)比圖。

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

由圖7(a)和圖7(b)可知，單獨(dú)采用瞬態(tài)顯式算法得到的仿真結(jié)果中可能存在取值誤差?！八矐B(tài)顯式-穩(wěn)態(tài)隱式”聯(lián)合算法仿真結(jié)果不存在取值誤差缺陷。表3給出了聯(lián)合算法的仿真結(jié)果與試驗(yàn)結(jié)果的對(duì)比分析值。

表3 “瞬態(tài)顯式-穩(wěn)態(tài)隱式”聯(lián)合算法加筋板架仿真值與試驗(yàn)值對(duì)比

由表3可知：采用聯(lián)合“瞬態(tài)顯式-穩(wěn)態(tài)隱式”算法得到的仿真結(jié)果與試驗(yàn)值的誤差約6.5%，與單獨(dú)采用瞬態(tài)顯式算法的仿真結(jié)果相比，誤差降低約39.8%，說(shuō)明在采用聯(lián)合“瞬態(tài)顯式-穩(wěn)態(tài)隱式”算法仿真計(jì)算時(shí)，在排除取值誤差的缺陷后，得到的加筋板架穩(wěn)態(tài)變形的仿真結(jié)果更接近試驗(yàn)值。

3 結(jié) 論

采用ABAQUS“瞬態(tài)顯式-穩(wěn)態(tài)隱式”聯(lián)合算法，計(jì)算船體加筋板架在水下爆炸作用下的穩(wěn)態(tài)響應(yīng)，求解加筋板架在水下爆炸沖擊載荷下的穩(wěn)態(tài)變形特性，并將數(shù)值仿真結(jié)果與試驗(yàn)結(jié)果進(jìn)行對(duì)比分析，結(jié)論如下：

(1) 采用“瞬態(tài)顯式-穩(wěn)態(tài)隱式”聯(lián)合算法得到的結(jié)果與試驗(yàn)結(jié)果相差在10%以內(nèi)，說(shuō)明該算法的仿真方法可用于水下爆炸作用下結(jié)構(gòu)的動(dòng)力響應(yīng)數(shù)值模擬分析。

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