龐崇安，劉玉鋒，王 震

(1.浙江同濟(jì)科技職業(yè)學(xué)院，杭州311231；2.浙江大學(xué) 建筑工程學(xué)院，杭州310058；3.浙大城市學(xué)院，杭州310015)

我國(guó)水利工程發(fā)達(dá)，大壩數(shù)目眾多，作為關(guān)系到民生、國(guó)防、環(huán)境、能源等重要領(lǐng)域的構(gòu)筑物，其毀壞的后果的是不堪設(shè)想的，中國(guó)古代歷史上就有多次通過(guò)決堰毀壩來(lái)達(dá)到戰(zhàn)略目的的戰(zhàn)役，可見(jiàn)大壩的安危影響之重大?，F(xiàn)代戰(zhàn)爭(zhēng)強(qiáng)調(diào)精確打擊，恐怖襲擊或意外爆炸作用導(dǎo)致的失事?lián)p失巨大，工程界對(duì)大型公共結(jié)構(gòu)的抗爆性能也日益重視，大壩抗爆性能的重要性也不言而喻[1,2]。針對(duì)混凝土壩在空中和水下爆炸荷載作用下的作用效應(yīng)與毀傷模式，國(guó)內(nèi)外學(xué)者開(kāi)展了多項(xiàng)研究。Kalateh基于歐拉-拉格朗日有限元法研究了空中爆炸沖擊荷載作用下混凝土壩的動(dòng)態(tài)響應(yīng)[3]，并對(duì)大壩的不同高度和起爆距離進(jìn)行了參數(shù)化分析。文獻(xiàn)[4]比較了水下和空中爆炸時(shí)大壩的動(dòng)態(tài)響應(yīng)認(rèn)為水下爆炸沖擊荷載破壞效應(yīng)更為嚴(yán)重，容易造成大壩結(jié)構(gòu)的破壞。文獻(xiàn)[5]用數(shù)值模擬研究蓄水位和炸藥位置等因素對(duì)沖擊波傳遞、大壩動(dòng)態(tài)響應(yīng)、混凝土損傷因子等方面的影響。針對(duì)水下接觸爆炸對(duì)大壩的動(dòng)態(tài)響應(yīng)和破壞特征，Xu[6]、Wang[7]、Li等研究認(rèn)為初始靜水壓力[8]、初始應(yīng)力場(chǎng)對(duì)混凝土重力壩的沖擊波傳播、爆炸振動(dòng)和破壞模式具有不可忽略的影響。Ren等基于離心機(jī)實(shí)驗(yàn)[9]，建立三維數(shù)值模型并進(jìn)行參數(shù)化分析，研究了混凝土重力壩在非接觸式水下爆炸攻擊下的破壞和失效情況。Moradloo等研究了拱形混凝土壩在水下爆炸作用下的動(dòng)態(tài)響應(yīng)和損傷特性[10]，結(jié)果表明離爆源越近，爆炸的損害越早，而對(duì)于一定數(shù)量的炸藥，大壩上游面位移隨爆炸深度的增加而增加。文獻(xiàn)[11]總結(jié)了爆炸荷載分類和目前常用的技術(shù)路線，并指出各種爆炸條件需要關(guān)注不同問(wèn)題。然而，由于混凝土大壩在水下爆炸作用下的破壞模式較為復(fù)雜，有必要針對(duì)混凝土大壩在水下爆炸荷載作用下的動(dòng)力響應(yīng)作進(jìn)一步的深入研究。

利用動(dòng)力顯式有限元軟件ANSYS/LS-DYNA對(duì)水下爆破荷載作用下的大壩進(jìn)行數(shù)值模擬，考慮其流固耦合效應(yīng)，研究分析了考慮重力和不考慮重力兩種模型下大壩在爆破荷載下的動(dòng)力響應(yīng)，以期為實(shí)際大壩工程動(dòng)力分析及抗爆設(shè)計(jì)提供參考。

1 材料模型與參數(shù)

1.1 炸藥

炸藥采用LS-DYNA中的* MAT_HIGH _EXPLOSIVE _BURN，炸藥狀態(tài)方程采用JWL方程

(1)

式中：A、B、R1、R2、ω為常數(shù)；V為體積；E為能量。文中具體參數(shù)如表1所示。其中:ρ為密度;D為起爆速度;PCJ為爆壓。

表 1 TNT炸藥模型參數(shù)

1.2 水介質(zhì)

對(duì)高壓下的水體采用Gruneisen方程

(γ0+αμ)E

(2)

式中：P為爆轟壓力；E為單位體積內(nèi)能；μ=ρ/ρ0-1；C為聲音在水中傳播速度；S1～S3、γ0為狀態(tài)方程系數(shù)，a為γ0的一階矯正系數(shù)。本文中水介質(zhì)的具體參數(shù)如表2所示。其中ρ為某個(gè)時(shí)刻的密度，E0為初始內(nèi)能，V0為相對(duì)體積。

1.3 空氣介質(zhì)

對(duì)氣體一般采用多方狀態(tài)方程

P=(C0+C1μ+C2μ2+C3μ3)+

(3)

式中：C0～C6為常數(shù)；V為相對(duì)體積；E為單位體積內(nèi)能，具體取值如表3所示。其中:ρ為密度;E0為初始內(nèi)能;V0為初始體積。

表 2 水介質(zhì)參數(shù)

表 3 空氣介質(zhì)參數(shù)

1.4 混凝土

混凝土采用HJC(Holmquist-Johnson-Cook)本構(gòu)模型[12]，在LS-DYNA中為*MAT_JOHNSON_HOLMQUIST_CONCRETE，該模型反映混凝土在大變形、高應(yīng)變情況下的特征，并能模擬其損傷過(guò)程?；炷恋牡刃?qiáng)度是關(guān)于靜水壓力、混凝土應(yīng)變率和損傷因數(shù)的函數(shù)，其中混凝土材料的損傷累積主要來(lái)自于塑性體積應(yīng)變、等效塑性應(yīng)變和靜水壓力的影響，具體參數(shù)見(jiàn)表4。其中:ρ為密度;G為剪切模量;A為標(biāo)準(zhǔn)化的凝聚力強(qiáng)度;B為標(biāo)準(zhǔn)化的壓力硬化強(qiáng)度;C為應(yīng)變率系數(shù);N為壓力硬化指數(shù);Fc為準(zhǔn)靜態(tài)單軸抗壓強(qiáng)度;T為最大靜水拉伸壓力;ETS0為準(zhǔn)靜態(tài)閾值應(yīng)變率;EFMIN為斷裂塑性應(yīng)變;SFMAX為標(biāo)準(zhǔn)最大強(qiáng)度;PC為壓碎應(yīng)力;UC為壓碎體積應(yīng)變;PL為鎖定壓力;UL為鎖定體積應(yīng)變;D1～D2為損傷常數(shù);K1～K3為壓力常數(shù);FS為失效塑性應(yīng)變。

表 4 混凝土HJC模型參數(shù)

1.5 巖石

巖石采用雙線性彈塑性模型，LS-DYNA中材料類別為*MAT_PLASTIC_KINEMATIC。具體參數(shù)如表5。其中:ρ為密度;E為彈性模量;PR為泊松比;SIGY為屈服強(qiáng)度;ETAN為強(qiáng)化模量，強(qiáng)化系數(shù)β=1代表采用等向強(qiáng)化模型;FS為破壞有效塑性應(yīng)變。

表 5 巖石參數(shù)

2 建立有限元模型

作為本文研究對(duì)象的大壩，壩體高85 m，壩頂寬15 m，壩底寬60 m，壩坡高65 m。大壩左側(cè)為水，水深75 m，水面以上至壩頂高度為空氣，此部分流體域?qū)挾热?30 m。大壩以下為基巖，基巖深度取170 m，寬300 m，左邊界與流體域齊平。炸藥位于水體表面，與大壩接觸處。大壩和基巖直接共節(jié)點(diǎn)，即認(rèn)為屬于剛接，流體域和固體域之間是兩個(gè)分離的體塊，不需設(shè)置接觸，但為使固體發(fā)生大變形后流體能填補(bǔ)空檔，在流體域深入固體部分4 m，重疊部分用水介質(zhì)填充。將基巖的左面、右面和底面節(jié)點(diǎn)的三個(gè)平動(dòng)自由度約束。計(jì)算模型的幾何如圖1(a)所示。本文采用掃略網(wǎng)格方法劃分模型，在炸藥附近將網(wǎng)格劃分得較細(xì)，往外每隔一段距離將網(wǎng)格大小放大，共放大兩次，最后得到的有限元模型如圖1(b)所示。在厚度方向，嘗試了數(shù)種劃分方式，典型方式包括厚度1 m，劃分為一個(gè)單元和厚度2 m，劃分為四個(gè)單元，兩者的結(jié)果見(jiàn)圖2?？梢园l(fā)現(xiàn)2 m厚度模型由于爆炸產(chǎn)生的缺口明顯大于1 m模型，其應(yīng)力極值(超過(guò)100 MPa)也遠(yuǎn)大于1 m模型，甚至遠(yuǎn)高于混凝土的正常極限強(qiáng)度。且由于極值過(guò)大，云圖主要集中在缺口附近，對(duì)其他區(qū)域的應(yīng)力分布不能很好地顯示。綜合以上考慮，并結(jié)合相關(guān)文獻(xiàn)[13-15]，將建模方法確定為厚度方向1 m，劃分為一個(gè)單元。

圖 1 大壩結(jié)構(gòu)模型

圖 2 兩種厚度下爆炸后16ms時(shí)應(yīng)力云圖

對(duì)于固體計(jì)算常用的Lagrange算法而言，一個(gè)單元內(nèi)只能是一種物質(zhì)，即單元是附在物質(zhì)上的。若炸藥、空氣單元采用Lagrange法，很容易在爆炸過(guò)程發(fā)生嚴(yán)重的畸變，往往使計(jì)算中斷。對(duì)于ALE算法，一個(gè)單元中可包含不同的材料，并能在空間網(wǎng)格中完成物質(zhì)的輸送。建模時(shí)，炸藥、水和空氣三種材料可采用歐拉網(wǎng)格建模，單元使用多物質(zhì)ALE算法。

在k文件中，添加關(guān)鍵字*CONTRAINED_LAGRANGE_IN_SOLID定義流固耦合，耦合方式采用罰函數(shù)的方法，需要提醒的是，此關(guān)鍵詞卡中前兩個(gè)參數(shù)分別指流體域和固體域，需要事先定義PART組，可在前處理中定義或修改*SET_PART_LIST完成。

計(jì)算的能量控制參數(shù)在前處理中用EDENERGY命令設(shè)置，本報(bào)告的計(jì)算模型只打開(kāi)石墻能(Stonewall Energy)和滑移能(Sliding Interface)，此外，體積黏性系數(shù)取默認(rèn)值，無(wú)需特別設(shè)置。用EDCTS控制時(shí)間步迭代大小，此值過(guò)大會(huì)影響計(jì)算收斂，有高能炸藥時(shí)缺省為0.67，本文中取0.6。在k文件中添加關(guān)鍵字*CONTROL_ALE來(lái)控制流固耦合時(shí)的相關(guān)參數(shù)設(shè)置，采用默認(rèn)值。最后在k文件中添加關(guān)鍵字*INITIAL_DETONATION來(lái)定義起爆點(diǎn)和起爆時(shí)間，將起爆點(diǎn)位置定義在炸藥域的中心，無(wú)重力時(shí)起爆時(shí)間定為0，有重力時(shí)起爆時(shí)間為1.2 s。

3 大壩水下爆炸荷載動(dòng)力響應(yīng)

修改前處理器生成的k文件，放入LS-DYNA求解器中計(jì)算，可得到爆炸后大壩隨時(shí)間變化的響應(yīng)，除了觀察其應(yīng)力分布之外，另外提取幾個(gè)典型位置點(diǎn)，讀取其時(shí)程歷史數(shù)據(jù)，典型位置的分布如圖3所示。

圖 3 典型位置點(diǎn)分布

3.1 未考慮重力模型

取炸藥當(dāng)量1 t(炸藥當(dāng)量根據(jù)炸藥密度乘體積算得)位于近水面緊靠大壩處，大壩上的Von-Mises應(yīng)力云圖如圖4所示。

模型的Von-Mises應(yīng)力云圖變化為與炸藥相鄰的壩體先被炸出一個(gè)裂口，即一小塊混凝土達(dá)到破壞準(zhǔn)則而退出工作，應(yīng)力波從此處向外擴(kuò)散，7 ms時(shí)四種模型的應(yīng)力波均擴(kuò)散至整個(gè)壩頂區(qū)域，之后壩頂區(qū)的應(yīng)力水平維持在一個(gè)較高的水平，而向下傳播的應(yīng)力波逐漸擴(kuò)散并衰減。整個(gè)過(guò)程中，與缺口相鄰的上方和下方，以及壩坡轉(zhuǎn)角處是應(yīng)力最大的位置，這也是容易產(chǎn)生應(yīng)力集中的地方。爆炸初始大壩上的應(yīng)力極值可超過(guò)50 MPa，后續(xù)極值穩(wěn)定在10～30 MPa之間震蕩。

圖 4 無(wú)重力爆炸后不同時(shí)刻Von-Mises應(yīng)力云圖

3.2 考慮重力后的模型

LS-DYNA中，通過(guò)定義加載曲線施加重力，LS-DYNA是動(dòng)力分析軟件，用其來(lái)進(jìn)行重力計(jì)算需要等待重力荷載擴(kuò)散并穩(wěn)定，即所謂的“擬靜力”，在此過(guò)程中，必須用關(guān)鍵字*DAMPING_PART_MASS_SET引入阻尼，否則結(jié)構(gòu)會(huì)不斷震蕩而不能達(dá)到穩(wěn)定。重力作用下和爆炸后的應(yīng)力云圖如圖5所示。

圖 5 重力作用下爆炸后不同時(shí)刻Von-Mises應(yīng)力云圖

重力作用下，壩踵處有著最高的應(yīng)力，可達(dá)數(shù)兆帕，壩趾處也有一定應(yīng)力集中的現(xiàn)象，但程度小于壩踵。爆炸發(fā)生后沖擊波擴(kuò)散的趨勢(shì)仍然比較明顯，但爆炸使混凝土破碎而產(chǎn)生的缺口消失了，過(guò)程中大壩上的最高應(yīng)力水平和未施加重力的模型相差不大，爆炸初始時(shí)高達(dá)50 MPa以上，隨后衰減到10～20 MPa，說(shuō)明沖擊波本身的應(yīng)力遠(yuǎn)大于重力造成的應(yīng)力大小，但是壩踵區(qū)在應(yīng)力波到達(dá)之前就有與之相近的應(yīng)力水平，說(shuō)明重力在該處造成的應(yīng)力集中水平與爆炸的應(yīng)力波約處于一個(gè)量級(jí)。然而，應(yīng)力波傳到壩踵與重力疊加后，壩踵的應(yīng)力水平并未上升，約3～8 MPa，與僅在重力作用下該處的應(yīng)力水平(4～8 MPa)相比并無(wú)顯著變化，說(shuō)明重力產(chǎn)生的應(yīng)力與爆炸產(chǎn)生的應(yīng)力在該處并不是線性疊加的關(guān)系，且偏于安全，但由于壩踵是影響大壩整體安全穩(wěn)定的關(guān)鍵部位，仍有必要對(duì)其進(jìn)行較為保守的估計(jì)。

3.3 未施加重力和施加重力模型對(duì)比

圖6是無(wú)重力作用下典型位置的位移時(shí)程曲線和壓強(qiáng)時(shí)程曲線。由于含重力模型在爆炸前各單元就有初始位移，為了更直觀比較爆炸后的響應(yīng)的變化，將含重力模型的位移減去重力作用下的初始位移。

由圖6(a)可知，無(wú)重力模型中各單元位移總體隨時(shí)間增大，而有重力模型的單元位移則減小，說(shuō)明爆炸產(chǎn)生的位移與重力產(chǎn)生位移有一部分互相抵消，且這種抵消與沖擊波的傳遞方向無(wú)關(guān)。兩個(gè)模型下爆炸缺口上方的位移都大于缺口下方，而位于壩體左側(cè)中部和中下部?jī)蓚€(gè)距離炸藥較遠(yuǎn)的點(diǎn)的位移較小，可以忽略。含重力模型爆炸后離炸藥較近的典型位置單元位移變化量絕對(duì)值明顯小于無(wú)重力模型，說(shuō)明重力對(duì)爆炸產(chǎn)生的動(dòng)力響應(yīng)產(chǎn)生了一定的抵抗作用，且不是線性疊加關(guān)系。

圖 6 典型位置單元時(shí)程曲線

無(wú)重力模型中，四個(gè)典型位置的壓強(qiáng)極值分別約4 MPa、3 MPa、2 MPa和1.5 MPa，缺口兩側(cè)的典型位置有多個(gè)壓強(qiáng)峰值，說(shuō)明沖擊波產(chǎn)生過(guò)多次反射，但極值均發(fā)生在第一個(gè)峰值處，且發(fā)生時(shí)間與距離炸藥距離呈負(fù)相關(guān)。重力造成的初始?jí)簭?qiáng)相對(duì)爆炸來(lái)說(shuō)較小，但有重力模型的四個(gè)典型位置的壓強(qiáng)峰值都小于無(wú)重力模型，這從另一個(gè)方面說(shuō)明重力的存在減輕了爆炸的作用。

4 結(jié)論

本文對(duì)考慮重力和不考慮重力兩種模型的大壩水下爆炸沖擊荷載及其動(dòng)力響應(yīng)進(jìn)行數(shù)值模擬，得到以下主要結(jié)論：

(1)爆炸初始產(chǎn)生的缺口是沖擊波擴(kuò)散的起點(diǎn)，沖擊波會(huì)多次反射，缺口和壩坡轉(zhuǎn)角處是容易產(chǎn)生應(yīng)力集中處。爆炸初始時(shí)最大應(yīng)力可超過(guò)50 MPa，之后在10～30 MPa之間震蕩。

(2)未考慮重力時(shí)，各點(diǎn)位移總體上呈隨時(shí)間增長(zhǎng)的趨勢(shì)，離炸藥較近的點(diǎn)位移較大。不同點(diǎn)的壓力都在爆炸后沖擊波第一次反射時(shí)達(dá)到極值，且極值隨離炸藥距離的增大而減小，可將第一次反射壓力作為爆炸沖擊荷載。

(3)重力作用下壩踵應(yīng)力最大，考慮重力時(shí)爆炸后大壩上沒(méi)有明顯缺口，但應(yīng)力極值并無(wú)顯著變化。

(4)重力產(chǎn)生的位移和壓力對(duì)爆炸效應(yīng)起抵抗作用，但(尤其在離炸藥接近的區(qū)域)兩者并不是線性疊加關(guān)系。