張博利

（中國民航大學(xué)基礎(chǔ)實驗中心，天津 300300）

近十多年來，全球范圍內(nèi)的恐怖爆炸襲擊事件不斷上升，已構(gòu)成對城市乃至國家的重要威脅。政府辦公和重要的基礎(chǔ)設(shè)施已成為恐怖襲擊的重要目標(biāo)。機場作為重要的基礎(chǔ)設(shè)施，現(xiàn)已成為恐怖襲擊的重點對象。爆炸恐怖襲擊是最常見而且容易實施的重要手段，像汽車炸彈、郵件（包裹）炸彈、人體炸彈、固定箱包炸彈等方式。由于汽車炸彈相對攜帶的爆炸物多一些，即使在與機場一定距離處爆炸，也可對機場進行毀滅性的破壞。然而其他三種炸彈，可通過嚴(yán)格的安全檢查措施加以控制。由于近距離爆炸作用于建筑物的沖擊波荷載非常大，可導(dǎo)致建筑物嚴(yán)重的破壞。因此對于機場來說，對汽車炸彈的防范成為重中之重。機場防恐防爆的主要目的是航站樓不會倒塌，減小結(jié)構(gòu)產(chǎn)生的整體破壞，控制產(chǎn)生次生殺傷的飛片。由于航站樓結(jié)構(gòu)設(shè)計的特殊性，不適合在四周設(shè)置防護墻。安裝防撞桿可以起到與設(shè)置防撞墻同樣的作用，并且防撞桿易于成型，堅實實用。本文對防撞桿進行安全性分析，并利用LS-DYNA軟件對它進行碰撞模擬實驗，以期為防撞桿設(shè)計提供參考。

Ls-dyna軟件是世界上最著名的通用顯示非線性動力分析程序，能模擬真實世界的各種復(fù)雜幾何非線性（大位移、大轉(zhuǎn)動和大應(yīng)變）、材料非線性和接觸非線性問題，特別適合求解各種二維、三維非線性結(jié)構(gòu)的高速碰撞、爆炸和金屬成形等非線性動力沖擊問題。

1 非線性有限元理論

圓薄壁鋼管在沖擊作用下的屈曲過程是一個動態(tài)的大位移和大變形的過程，接觸狀態(tài)和沖擊載荷都影響著屈曲的全過程，系統(tǒng)具有集合非線性和材料非線性等多重非線性。一般的線性有限元方法都基于線性的小位移系統(tǒng)，本例題所研究的系統(tǒng)模擬計算則需要采用動態(tài)大變形非線性有限元方法。

1.1 物體變形狀態(tài)的描述

由于物體在變形前后發(fā)生很大的變形，因此必須考慮物體受力前后在空間位置上的改變。

為了描述物體變形前后兩種不同狀態(tài)，須引入?yún)⒖紶顟B(tài)和變形狀態(tài)。在某一瞬時，物體在空間所占據(jù)的區(qū)域稱為物體的構(gòu)型（Configuration）。令在時間t=0時，物體的初始構(gòu)型為V0，參考于一固定的直角坐標(biāo)系{xi}，物體的任一質(zhì)點P的位置可由一向徑p（x1，x2，x3）或其質(zhì)點坐標(biāo) x1，x2，x3確定。構(gòu)型 V0稱為物體的參考狀態(tài)。在后來某一瞬間t，物體被移動到空間的另一位置，其構(gòu)型為V，這時的狀態(tài)稱為變形狀態(tài)。描述這一變形狀態(tài)，用另一直角坐標(biāo)系（yi）。初始構(gòu)型中的P點，變形后被移動到空間的P點，可由一向徑p（x1，x2，x3）確定。如果令坐標(biāo)系（xi）和（yi）重合，則在二維情況下同一個質(zhì)點變形前后的關(guān)系

其中，yi是x的單值連續(xù)函數(shù)，且雅可比行列式（Jacobi）的值不等于 0。

如果ui為質(zhì)點xi軸方向的位移，那么

或?qū)懗上蛄啃问剑?/p>

在描述物體變形前后的不同狀態(tài)時，有如下兩種方法。

1.2 拉格朗日（Lagrange）描述

x1，x2，x3和t作獨立變量來描述物體的運動（或變形），稱為物體描述（material description）或拉格朗日（Lagrange）描述，而 x1，x2，x3和 t稱為拉格朗日變量，上述式（1）～式（3）就是這種描述形式。此種描述是以變形前的狀態(tài)，即以己知的初始構(gòu)型為V0參考狀態(tài)，隨著運動的質(zhì)點研究質(zhì)點的運動狀態(tài)。在描述狀態(tài)中其中p代表坐標(biāo)為x1，x2，x3的質(zhì)點，是不隨著時間t變化的。vt表示速度表示加速度，M表示質(zhì)量，p表示質(zhì)點軌跡，ρ0（p）是構(gòu)型V0在p點處的密度。

1.3 歐拉（Euler）描述

取y1，y2，y3和t作為獨立變量來描述物體的運動（或變形），稱為空間描述（spatial description）或歐拉（Euler）描述，y1，y2，y3和 t為歐拉變量?？臻g描述是以當(dāng)前的構(gòu)型V為參考狀態(tài)，研究空間各個不同的質(zhì)點經(jīng)過空間某一定點y1，y2，y3時的狀態(tài)。在描述狀態(tài)中，，其中，p 代表坐標(biāo)為 y1，y2，y3的質(zhì)點，vi表示速度表示加速度，M 表示質(zhì)量，V表示質(zhì)點軌跡，ρ0（p）是構(gòu)型V在點p處的密度。

在碰撞模擬計算中接觸問題的處理是一個重點，接觸問題的處理實質(zhì)上就是接觸力的計算，對于接觸力的計算主要采用罰函數(shù)法與拉格朗日乘子法。罰函數(shù)法雖然只是近似的滿足接觸碰撞邊界條件，但是容易實現(xiàn)，并與顯示算法完全相容，因此在實踐中得到廣泛的應(yīng)用。本文求解接觸碰撞問題所用的LS－DYNA有限元軟件就是采用罰函數(shù)法確立接觸界面的。

2 有限元模型的建立

某種防撞桿的結(jié)構(gòu)具體參數(shù)如下：圓鋼管直徑D=300 mm，壁厚 t＝2.5 mm，長度 L=1 000 mm，一端完全固定，另一端施加載荷。模擬分析采用10 ms內(nèi)施加0.4 N的力，分析鋼管在整個過程中的變形及應(yīng)力分布情況。鋼管底端為固定端，約束所有的線位移自由度，鋼管頂端除UZ方向外，其余全部約束。由于鋼管表面在沖擊壓縮過程中，可能發(fā)生表面折疊和自相接觸現(xiàn)象，需定義為單面接觸算法ASSC，接觸面的摩擦系數(shù)為0.1，求解時間為0.015 s。

表1 材料參數(shù)Tab.1 Material parameters

表2 沖擊載荷數(shù)組Tab.2 Array of impacting loading

采用固定網(wǎng)格劃分后，在沖擊后鋼管0.008 24s時的變形和應(yīng)力情況，如圖1所示：

由于薄壁管在沖擊作用下，將發(fā)生極大的塑形變形，管件內(nèi)部的等效應(yīng)變和等效應(yīng)變率分布是不均勻的，且應(yīng)變集中在某些區(qū)域，這導(dǎo)致在網(wǎng)格中只有小部分單元格發(fā)生干涉或嚴(yán)重變形，而大部分單元格尤其是內(nèi)部單元格的質(zhì)量依然良好。因此當(dāng)檢測到網(wǎng)格中發(fā)生畸變或干涉后，可以先只對少量畸變或干涉單元進行局部調(diào)整。在有限元計算歐拉方法中，計算網(wǎng)格固定于幾何空間，材料通過流進流出，計算單元質(zhì)量、動量和能量的變化，從而確定物體的變形。歐拉方法能夠處理材料嚴(yán)重的扭曲變形，但很難跟蹤物質(zhì)的變形，不能識別材料界面，計算材料的彈塑性不夠精確。拉格朗日方法把計算固定在材料上，可以跟蹤物質(zhì)的運動，能夠清楚地識別材料的邊界和交界面，并且具有較高的計算速度和精度。然而對于有較大的網(wǎng)格畸變，材料破碎使得簡單連接區(qū)域變成多連接，以及網(wǎng)格重疊造成負(fù)質(zhì)量時，將導(dǎo)致時間步長縮小，求解精度降低，甚至停止計算，解決這個問題的好方法就是采用自適應(yīng)網(wǎng)格方法。

自適應(yīng)網(wǎng)格方法是指在計算某些變化較為劇烈的區(qū)域，如大變形、激波面、接觸斷面和滑移面等，網(wǎng)格在迭代過程中不斷調(diào)節(jié)，將網(wǎng)格細(xì)化。做到網(wǎng)格點分布與物理解的耦合，從而提高解的精度和分辨率的一種技術(shù)。自適應(yīng)網(wǎng)格希望在物理解變動較大的區(qū)域網(wǎng)格自動密集，而在物理解變化平穩(wěn)區(qū)域網(wǎng)格相對稀疏，這樣在保持計算高效率的同時得到高精度的解。在LS-DYNA中，自適應(yīng)網(wǎng)格劃分方法分為二種：hadaptive方法和r-adaptive方法。h-adaptive方法是指單元變形較大時，將單元分為更小的單元以改善精度，僅適用于殼單元，主要用于金屬成型模擬、薄壁結(jié)構(gòu)受壓屈曲等問題。r-adaptive方法也稱重分區(qū)技術(shù)或者自適應(yīng)網(wǎng)格部分，是指重新移動和劃分單元，將網(wǎng)格節(jié)點重新排列，得到時步長內(nèi)最優(yōu)縱橫比。

LS-DYNA采用自適應(yīng)網(wǎng)格方法的目的在于使用有限的資源獲得最大的計算精度。用戶設(shè)置好了初始網(wǎng)格和自適應(yīng)劃分級別后，程序根據(jù)需要將某些單元進行分割。雖然這種方法不能完全解決過程中的誤差，但與固定網(wǎng)格相比，可以使用較少的單元和計算資源來獲得盡可能高的計算精度。

薄壁圓鋼管在沖壓過程中，采用的是h-adaptive方法。只是將固定網(wǎng)格通過裂變，細(xì)分成更小的單元。裂變是指某些單元由于精度需要細(xì)分為更小的單元，如圖2所示。裂變后，新的邊長尺寸是原來的1/2。本例采用的是3級網(wǎng)格細(xì)化。細(xì)化后，鋼管的在同一時間的壓縮情況和特定單元的應(yīng)力都將發(fā)生變化，如圖3所示。

在圖4和圖5中，A代表圓管壓縮過程中變形較小的單元，B代表圓管壓縮過程中變形較大的單元。通過二圖中A曲線的比較，在如圖4所示的固定網(wǎng)格劃分的情況下，Z方向的應(yīng)力變化不大。而在圖5自適應(yīng)網(wǎng)格劃分后，由于網(wǎng)格個數(shù)的增多，Z方向的應(yīng)力變化較大，真實反映了變形較小單元在受沖擊過程中的應(yīng)力隨時間的變化情況。在圖4和圖5中B曲線在0.005 s和0.009 s左右時的應(yīng)力值幾乎達到同一值，說明在采取自適應(yīng)網(wǎng)格后，在大變形的情況下，最大值和最低值的變化不大，只是對變形過程的反映更加真實化、準(zhǔn)確化。

3 結(jié)語

通過以上數(shù)據(jù)可以看到，自適應(yīng)網(wǎng)格方法的目的在于使有限的計算資源獲得最大的計算精度，當(dāng)設(shè)置好初始網(wǎng)格和自適應(yīng)劃分級別后，程序根據(jù)需要將某些單元進行分割，雖然這種方法并不能完全解決求解過程中的誤差，但與固定網(wǎng)格相比，可以使較少的單元和計算資源盡可能提高求解精度。自適應(yīng)網(wǎng)格技術(shù)在薄板沖壓變形模擬，薄壁結(jié)構(gòu)受壓屈曲以及機場防撞桿的受沖擊等大變形情況中，LS-DYNA程序通過自適應(yīng)網(wǎng)格剖分功能，使出現(xiàn)大彎曲變形的起皺部分自動加密網(wǎng)格，對計算精度有了顯著提高。因此，自適應(yīng)單元網(wǎng)格劃分技術(shù)是LS-DYNA處理單元大變形問題的重要手段。

[1]陸 勇，曹立波.對汽車撞柱的仿真研究[J].農(nóng)業(yè)裝備與車輛工程，2006（1）：28-31.

[2]田志敏，張想柏，杜修力.防恐怖爆炸重要建筑物的概念設(shè)計[J].土木工程學(xué)報.2007，40（1）：34-41.

[3]何 濤，楊 競.ANSYS10.0/LS-DYNA非線性有限元分析實例指導(dǎo)教程[M].北京：機械工業(yè)出版社，2007.

[4]尚曉江，蘇建宇.ANSYS/LS-DYNA動力分析方法與工程實例[M].北京：中國水利水電出版社，2006.

[5]LAN SHENGRUI，CRAWFORD JOHN E.Development of shallow footing anti-ram bollard system through modeling and testing[J].Transactions of Tianjin University，2006，19：46-50.