基于自由落體的落石沖擊土層的理論及LS-DYNA模擬研究

王 星，周天躍，師江濤,2，夏永旭

(1. 長安大學(xué) 公路學(xué)院，西安 710064； 2. 陜西交通建設(shè)集團，西安 710086)

滾石災(zāi)害已經(jīng)成為目前中國的3大地質(zhì)災(zāi)害問題之一，尤其是雨季時隧道洞口處的滾石災(zāi)害情況更為嚴(yán)峻.各種滾石防護結(jié)構(gòu)便應(yīng)需而設(shè)，而防護結(jié)構(gòu)設(shè)計的核心即是要準(zhǔn)確解析并計算出滾石的沖擊力.

國內(nèi)外學(xué)者針對滾石沖擊力的計算做出了大量的研究：其中傳統(tǒng)的Hertz接觸理論及Thornton彈塑性理論被廣泛引用[1-2].文獻(xiàn)[3]深入分析了物體間在靜態(tài)接觸、滑動、滾動、沖擊等作用下的力學(xué)響應(yīng)規(guī)律，并給出了相應(yīng)的計算公式.文獻(xiàn)[4]通過采用室內(nèi)試驗的方法，得出了滾石沖擊力的半經(jīng)驗半理論算法.文獻(xiàn)[5]以室內(nèi)試驗及動量定理為基礎(chǔ)，給出了一種滾石沖擊力的半理論半經(jīng)驗算法.文獻(xiàn)[6]從現(xiàn)場測試的角度出發(fā)，獲得了滾石沖擊力及沖擊深度的量綱為1的計算表達(dá)式.文獻(xiàn)[7]研究了滾石直接沖擊混凝土板時的力學(xué)響應(yīng)情況.文獻(xiàn)[8]對滾石沖擊力進行了研究，并得出了滾石沖擊力呈現(xiàn)出脈沖式變化規(guī)律的結(jié)論.文獻(xiàn)[9]同樣證明了滾石沖擊過程屬于一種脈沖式的碰撞過程.文獻(xiàn)[10]以理論計算為手段，分別推導(dǎo)出了滾石沖擊力的空腔膨脹算法以及能量守恒算法.文獻(xiàn)[11]采用實驗及動量定理的方法，歸納出了滾石沖擊力的計算方法.文獻(xiàn)[12]同樣在動量定理的基礎(chǔ)上總結(jié)出了滾石沖擊力計算方法.文獻(xiàn)[13]研究了滾石沖擊荷載下墊層材料的動力響應(yīng)情況.文獻(xiàn)[14]分析了滾石沖擊柔性防護網(wǎng)的力學(xué)響應(yīng)規(guī)律.文獻(xiàn)[15]結(jié)合理論算法推導(dǎo)出了滾石沖擊防護結(jié)構(gòu)墊層的沖擊力.

針對滾石沖擊力，國內(nèi)外學(xué)者研究成果較多，意義深遠(yuǎn).然而目前國內(nèi)的相關(guān)設(shè)計規(guī)范及其他滾石沖擊力的算法，其所得結(jié)果多為滾石的平均沖擊力，同時日本、瑞士學(xué)者所給出的經(jīng)典算法均是基于試驗所得出的半理論半經(jīng)驗算法，且算式中各字母變量的意義尚不十分明確.

為此，本文作者首先基于滾石沖擊力脈沖式的變化特點，推導(dǎo)出滾石最大沖擊力的脈沖式算法；其次再考慮滾石沖擊墊層土體的彈塑性特點，得出滾石最大沖擊力的彈塑性算法；再次通過采用LS-DYNA軟件模擬分析了滾石的沖擊力及侵徹深度隨時間的變化規(guī)律，并歸納出滾動沖擊力的LS-DYNA算法.最終將所獲得的3種算法與已有算法進行對比分析，以期為滾石防護設(shè)計提供參考.

1 滾石沖擊力的脈沖算法

文獻(xiàn)[11]提出了一種滾石沖擊力的計算方法，定義為關(guān)寶樹算法，其沖擊力P及沖擊時間tc的表達(dá)式分別為

(1)

(2)

式中：m為滾石的質(zhì)量；a為滾石沖擊加速度；Q為滾石的重量；g為重力加速度；h為墊層材料的厚度；He為滾石的下落高度.

從而有滾石沖擊過程中的平均沖擊力Pa為

(3)

由于Pa產(chǎn)生的沖量與實際沖擊力P(t)產(chǎn)生的總沖量It相等，則有

(4)

式中：2t0為滾石沖擊時間.

滾石沖擊過程屬于一種脈沖式行為，且沖擊力隨時間基本呈現(xiàn)出近似正弦變化的趨勢[8-10,16]，故考慮P(t)近似滿足圖1中的正弦變化規(guī)律.

圖1 滾石沖擊力的變化曲線Fig.1 Variation curve of rockfall impact force

則有滾石實際沖擊力P(t)與滾石最大沖擊力Pmax之間的關(guān)系為

P(t)=Pmaxsin(ω0t)

(5)

根據(jù)正弦曲線的性質(zhì)可得

(6)

式中：T為周期；ω0為角速度.

聯(lián)立式(4)～(6)可得

(7)

聯(lián)立式(3)及式(7)，便可求出滾石的最大沖擊力Pmax為

(8)

考慮到滾石在沖擊墊層土體后會產(chǎn)生一定的回彈，假定滾石的初始速度為v初，反彈的速度為v彈，v初與v彈的方向相反，則有恢復(fù)系數(shù)κ為

κ=v彈/v初

(9)

根據(jù)動量定理，則有考慮恢復(fù)系數(shù)影響后，對滾石最大沖擊力Pmax進行修正的滾石沖擊力P修的表達(dá)式為

P修t=mv初-(-mv彈)

(10)

聯(lián)立式(9)～(10)，可以得出

P修t=(1+κ)mv初=(1+κ)Pmaxt

(11)

從而有

P修=(1+κ)Pmax

(12)

考慮恢復(fù)系數(shù)κ的影響，聯(lián)立式(8)及式(12)，則可得出滾石沖擊力的脈沖算法的表達(dá)式，定義為脈沖算法

(13)

參數(shù)κ可參考文獻(xiàn)[17]的實驗結(jié)果進行取值.

進而求取滾石侵徹深度的理論計算值.假定v(t)為隨沖擊時間變化的動態(tài)滾石沖擊速度，則結(jié)合牛頓第二定律有

(14)

聯(lián)立式(5)～(6)、(14)，則有

(15)

對式(15)進行積分后，則有

(16)

(17)

式中：C為常數(shù).

此處存在邊界條件t=2t0,v=0，則有

(18)

從而可得出滾石的沖擊深度為

(19)

對式(19)進行積分后，則有

(20)

聯(lián)立式(2)、(13)、(20)，可得

(21)

2 滾石沖擊力的彈塑性算法

滾石在沖擊下部墊層土體后，最先與滾石接觸的正下方土體區(qū)域?qū)⒔?jīng)歷：接觸-彈性變形-彈性變形的臨界狀態(tài)-塑性變形，而距離滾石較遠(yuǎn)的區(qū)域則一直保持著彈性變形的狀態(tài).在沖擊行為發(fā)生后，滾石與墊層土體之間將會形成球冠接觸面，為了對滾石沖擊力進行計算，需通過轉(zhuǎn)換為滾石與墊層土體之間的接觸面積大小進行考慮并計算，最終形成圖2所示的接觸面積的區(qū)域劃分.圖2中，滾石與墊層土體之間形成接觸半徑為fp的塑性區(qū)，而彈性接觸面積外緣的半徑為fr.

圖2 墊層的彈塑性區(qū)域分布Fig. 2Regional distribution of elastoplasticity of cushion

滾石球體沖擊墊層土體時的彈性接觸壓力Pe與變形量Se之間的關(guān)系式為

(22)

式中：Ra為等效半徑；Ea為等效彈性模量.

Ra與Ea計算方式為

(23)

式中：υ1、υ2分別為滾石、墊層土體的泊松比；E1、E2分別為滾石、墊層土體的彈性模量；R1、R2分別為滾石、墊層土體的半徑.

滾石沖擊墊層土體后會形成一個球冠接觸面，其接觸面積Ac為

Ac=πf2=πR1S

(24)

式中：f為滾石的接觸半徑；S為對應(yīng)于接觸半徑f下的侵徹深度.

從而有

f2=R1S

(25)

土體在彈性變形階段的極限屈服應(yīng)力py及對應(yīng)的接觸半徑fy之間滿足關(guān)系式[18]

(26)

py與屈服應(yīng)力系數(shù)λ、材料屈服強度χM之間的關(guān)系，以及λ與υ2之間的關(guān)系為[19]

(27)

聯(lián)立式(25)～(27)，便可以求出滾石的彈性極限沖擊深度Sel為

(28)

聯(lián)立式(22)及式(28)，即可算出滾石的極限彈性沖擊力Pel為

(29)

結(jié)合圖2可得出相應(yīng)的力學(xué)平衡關(guān)系式為

(30)

式中：fp為土體墊層的塑性區(qū)半徑.

從而彈塑性沖擊條件下的力學(xué)平衡關(guān)系式為

(31)

式中：fmax為滾石沖擊過程中的最大接觸半徑；fr為彈性接觸半徑；Pep為彈塑性沖擊力.

聯(lián)立式(21)、(25)及式(31)便可得出滾石沖擊力的彈塑性算法的表達(dá)式，定義為彈塑性算法

Pep=Pel+πpyR1

(32)

3 LS-DYNA動態(tài)數(shù)值模擬

采用LS-DYNA軟件模擬滾石沖擊墊層土體的動態(tài)過程，并選取典型計算工況進行呈現(xiàn)：墊層的材料為砂土，尺寸為5 m(長)×5 m(寬)×1 m(高)，滾石半徑為0.5 m，沖擊速度為20 m/s.滾石與墊層土體間設(shè)置為“面-面接觸”，墊層土體中心1.5 m×1.5 m的區(qū)域進行了網(wǎng)格加密.墊層土體的下底面設(shè)置了3個方向的全約束，滾石的沖擊時間設(shè)置為0.06 s，計算步設(shè)置為200步，此次計算時間約為46 min.材料參數(shù)具體見表1.

表1 墊層土體材料參數(shù)

LS-DYNA的數(shù)值計算模型及滾石沖擊墊層土體后0.015 s、0.045 s及0.06 s時的Von mises應(yīng)力云圖如圖3所示，根據(jù)圖3中的應(yīng)力等級色條便可得出土體墊層各部分的應(yīng)力分布情況，同時也可判別相應(yīng)的彈塑性區(qū)域.

圖3 LS-DYNA數(shù)值模擬結(jié)果Fig.3 LS-DYNA numerical simulation results

滾石的沖擊力隨時間的變化曲線如圖4所示.由圖4可見，在滾石沖擊0.004 s后，滾石的沖擊力即達(dá)到峰值，而在此后的0.004～0.048 s期間，滾石的沖擊力在逐步的震蕩過程中減小至0，滾石沖擊過程完成.

圖4 沖擊力的變化曲線Fig. 4Variation curve of impact force

滾石沖擊過程中侵徹深度隨時間的變化曲線如圖5所示.可見在0～0.034 s期間，滾石的侵徹深度將會隨著沖擊時間的增加而呈現(xiàn)出近拋物線的增加趨勢，且在0.034 s時形成最大沖擊坑深0.301 m，而在0.034～0.048 s期間，產(chǎn)生了約2.6 cm的回彈，并最終形成了0.275 m的沖擊坑深.

圖5 侵徹深度的變化曲線Fig. 5Variation curve of penetration depth

滾石沖擊結(jié)束后，墊層土體結(jié)構(gòu)中軸線上的各節(jié)點的最終位移如圖6所示.

圖6 沖擊坑的位移曲線Fig.6 Displacement curve of impact crater

可見沖擊坑橫剖面的左右兩側(cè)形狀對稱，沖擊坑中心所形成的最大位移變形為0.275 m，并在坑的周圍形成了約8 cm的隆起，墊層土體發(fā)生變形的主要區(qū)域是坑中心處半徑為0.5 m圓的范圍.

4 計算結(jié)果對比分析

為對各算法下滾石沖擊力的計算結(jié)果進行對比分析，選取滾石半徑為0.5 m，滾石質(zhì)量為1 308.25 kg，墊層厚度取為1 m，其余的計算參數(shù)見表1及表2，計算工況中的各計算參數(shù)參考文獻(xiàn)[20-21]進行選取，各種算法的計算表達(dá)如下[22].

表2 滾石沖擊力計算工況表

1)Hertz算法：

(33)

2)瑞士算法：

(34)

3)日本算法：

(35)

4)隧道手冊算法：

(36)

(37)

(38)

式中：ME為墊層土體彈性模量；ε為拉梅常數(shù)；v0為滾石初始沖擊速度；vc為滾石沖擊墊層的壓縮波速.

碎石土、砂土、黏土墊層情況下墊層各算法所得出的滾石沖擊力如圖7～圖9所示.通過對圖7～圖9進行分析，可得出如下結(jié)論：

圖7 碎石土墊層的“沖擊速度-沖擊力”曲線Fig.7 “Impact speed - impact force” curves of gravel soil cushion

圖8 砂土墊層的“沖擊速度-沖擊力”曲線Fig. 8 “Impact speed - impact force” curves of sandy soil cushion

圖9 黏土墊層的“沖擊速度-沖擊力”曲線Fig. 9 “Impact speed - impact force” curves of clay cushion

1)Hertz算法所得出的滾石沖擊力值無疑是偏大的，這主要是由于該算法是基于完全彈性理論所得出.瑞士算法結(jié)果略微偏大，因為該算法考慮了墊層材料的彈性模量的影響，并取彈性模量值的2/3次方作為因子，同時瑞士算法明顯滿足P碎石土>P砂土>P黏土，而當(dāng)墊層土體為黏土材料時，瑞士算法沖擊力已小于日本算法，由此可見，瑞士算法結(jié)果整體較大，同時墊層土體的彈性模量對瑞士算法結(jié)果影響較大.

2)日本算法所得結(jié)果整體較為偏大，且沖擊力值整體較為穩(wěn)定，是因該算法考慮將滾石重力值的2/3次方及下落高度的3/5次方作為滾石沖擊力的計算因子，而拉梅常數(shù)值一般取為106.

3)關(guān)寶樹算法所得出的計算結(jié)果整體偏小，這是由于其算法所得出的沖擊力為整個沖擊過程中的平均沖擊力值，且實驗過程中選用的是5 N及8 N的試驗錘，或有可能存在一定的局限性，當(dāng)試驗中的滾石質(zhì)量達(dá)到13 082.5 N時，關(guān)寶樹算法可能會存在一定的偏差，從而導(dǎo)致計算結(jié)果偏低.

4)隧道手冊算法所計算出的沖擊力值普遍偏小，是因為在該算法中滾石沖擊力時間的計算是通過2倍的墊層厚度與壓縮波速之比所得出的，墊層厚度直接影響到?jīng)_擊時間，故而沖擊時間勢必存在一定的偏差，最終致使沖擊力與實際存在偏差.

5)LS-DYNA數(shù)值模擬結(jié)果與日本算法、瑞士算法、彈塑性算法及脈沖算法之間的吻合度較好，實現(xiàn)了各算法之間相互驗證.脈沖算法整體與日本算法、瑞士算法及數(shù)值模擬算法保持一致，但仍略有偏小.彈塑性算法與日本算法、瑞士算法、數(shù)值模擬算法及脈沖算法整體一致性較好，從而證明了該算法的可靠性，這或可能是由于該算法充分考慮了滾石沖擊過程中墊層土體的彈塑性變化情況，從而得出了較為準(zhǔn)確的計算公式.

5 滾石沖擊力的LS-DYNA算法

由于LS-DYNA數(shù)值模擬計算結(jié)果為滾石的最大沖擊力，而隧道手冊算法所得出的計算結(jié)果為滾石的平均沖擊力，且隧道手冊算法充分考慮了墊層土體的厚度、彈性模量、密度、泊松比及滾石的重量及下落高度等因素.故而嘗試定義將滾石的最大沖擊力與平均沖擊力之比作為沖擊力的放大系數(shù)，且給出放大系數(shù)的建議值，并通過計算獲取碎石土、砂土及黏土墊層的滾石沖擊力放大系數(shù)ζ與滾石沖擊速度v0之間的關(guān)系，建立v0-ζ散點圖，并進而擬合出相應(yīng)的滾石沖擊力放大系數(shù)曲線，具體見圖10.

圖10v0-ζ曲線Fig.10v0-ζ curves

對于碎石土、砂土及黏土墊層，其滾石沖擊力的放大系數(shù)曲線表達(dá)式為

(39)

式中：ζ碎石、ζ砂土、ζ黏土分別為碎石土墊層、砂土墊層、黏土墊層時的沖擊力放大系數(shù).

鑒于式(39)中3條擬合曲線之間存在較小差異，故在綜合考慮滾石沖擊過程中的多方因素后，得出最終統(tǒng)一的沖擊力放大系數(shù)ζe的表達(dá)式

ζe=0.655+0.241ln(v0+2)

(40)

聯(lián)立式(36)～(38)、(40)，并考慮滾石恢復(fù)系數(shù)的影響，從而得出滾石沖擊力為

(41)

在考慮滾石沖擊角度及滾石自重因素的影響后，得出滾石最大沖擊力的LS-DYNA算法的最終表達(dá)式，定義為LS-DYNA算法

(42)

式中：θ為滾石沖擊速度與墊層土體表面的夾角；β為滾石重力的方向與墊層表面的夾角.

為驗證式(42)的準(zhǔn)確性，將其計算結(jié)果與Hertz算法、日本算法等8種算法的計算結(jié)果進行對比，并分別考慮碎石土、砂土及黏土墊層的情況，對比結(jié)果如圖11～13所示.

圖11 碎石土墊層的驗證結(jié)果Fig. 11Verification results of gravel soil cushion

圖12 砂土墊層的驗證結(jié)果Fig.12 Verification results of sandy soil cushion

圖13 黏土墊層的驗證結(jié)果Fig. 13Verification results of clay cushion

可見LS-DYNA算法結(jié)果與日本算法、瑞士算法、數(shù)值模擬結(jié)果、彈塑性算法及脈沖算法結(jié)果具備較好的一致性，脈沖算法及隧道手冊算法結(jié)果仍略微偏小，而Hertz算法仍明顯偏大，而關(guān)寶樹算法無疑偏小，與之前的對比驗證結(jié)果保持一致，其較好地驗證了式(42)的準(zhǔn)確性.

在得出滾石最大沖擊力值后，此時便可結(jié)合文獻(xiàn)[12]計算出滾石沖擊力作用在防護結(jié)構(gòu)頂板上的壓力q，見圖14.

圖14 防護結(jié)構(gòu)頂板上的沖擊壓力Fig. 14Impact pressure on roof of protective structure

由圖14可知q為

(43)

假定滾石沖擊力作用于防護結(jié)構(gòu)頂板上部的墊層的面積為滾石球體的等效截面，即落石等效直徑為2R1的區(qū)域.其中在隧道手冊算法中中假定θ為30°，而在路基規(guī)范中假定θ為35°，此處由文獻(xiàn)[22]可得

θ=45°-φ/2

(44)

式中：φ為墊層土體的內(nèi)摩擦角.

6 結(jié)論

1)基于關(guān)寶樹算法以及滾石沖擊力的脈沖式變化規(guī)律，同時考慮了滾石沖擊恢復(fù)系數(shù)的影響，得出了滾石沖擊力以及侵徹深度的脈沖算法.結(jié)合滾石沖擊過程中墊層土體的彈塑性變化過程，建立了滾石沖擊力的力學(xué)平衡關(guān)系式，并基于侵徹深度與接觸半徑之間的關(guān)系，推導(dǎo)出了滾石沖擊力的彈塑性算法.

2)通過采用LS-DYNA軟件分析研究了滾石沖擊墊層土體的微觀力學(xué)響應(yīng)機理，得出了滾石沖擊力及侵徹深度隨沖擊時間的變化規(guī)律，并給出了最終塑性沖擊坑的形狀.

3)結(jié)合具體的滾石沖擊計算工況，對各算法進行了對比分析.研究結(jié)果表明：完全彈性的Hertz算法結(jié)果明顯偏大，日本算法、瑞士算法、數(shù)值模擬算法、脈沖算法、彈塑性算法結(jié)果吻合性較好，但相互之間存在少量差異，脈沖算法及隧道手冊算法略微偏小，關(guān)寶樹算法則明顯偏小.

4)以數(shù)值模擬結(jié)果及隧道手冊算法為基礎(chǔ)，建立了v0-ζ之間的擬合曲線，并在考慮κ、θ、滾石自重影響的基礎(chǔ)上，得出了滾石最大沖擊力的LS-DYNA算法，最終結(jié)合8種算法結(jié)果對LS-DYNA算法的適用性進行了驗證.