雷電沖擊作用下土體的變形特性分析

饒平平，項遠兵，崔紀飛，吳志林，歐陽昢晧

(上海理工大學 環境與建筑學院，上海 200093)

雷電是地球上普遍存在的自然現象，全球平均每年約發生14億次，頻率約為45 次/s[1]，其中，云地閃電是一種強大雷電流向大地釋放能量的過程，每次云地閃電都涉及109～1010J的能量，盡管其中大部分能量在空氣傳導中被消耗以產生雷聲、熱空氣、光和無線電波[2-3]，仍有約108J的能量釋放進大地中，產生的峰值電流可達數百千安，持續數十微秒，巨大能量和強烈電磁場活動會引發顯著的沖擊效應、電磁效應、光效應以及熱效應，促使土體發生各種物理和化學過程[4-6]。近年來，雷電的相關研究已越來越細致深入，如部分電氣工程學者對接地裝置的雷電沖擊特性開展研究，通過室內試驗與數值模擬的方法對含接地極的土體雷電沖擊特性進行了大量研究，重點分析了土體電性參數、雷電電流波形及其幅值、接地極結構等因素的影響[7-8]，同時，通過X射線透射成像技術拍攝了土中的沖擊放電影像[9-10]。航天材料學者們重點研究了用于飛機表層防護的碳纖維增強復合材料(CFRP)的雷擊熱效應[11-12]和力學損傷[13-16]。

然而，至今國內外以巖土體為對象，針對雷擊對其造成的影響研究仍較少。實際上，雷擊足以對土體造成破壞， Yokoyama[17]提到雷擊會導致一些地方機場停機坪出現塌陷或地坑；Wakasa等[18]通過室內試驗模擬雷擊對巖石造成的破壞，指出巨大的雷電能量足以使巖石或土體失穩及破壞；Fajingbesi等[19]采用實驗與數值模擬探究了雷擊大地后土體電導率的變化以及土中放電電流分布；Castro等[20]通過實驗分析了閃電與火山巖相互作用過程中的物理和化學變化，并使用電弧焊設備再現了天然閃電熔巖的地球化學和紋理特征；Chen等[21-22]采用理論分析和野外觀測指出雷擊對巖石中閃電熔巖的形成起著重要作用；Rao等[23-25]利用沖擊電流發生器和自行設計的試驗裝置對不同條件(不同含水量、密實度、含鹽量和溫度)下的土體進行了模擬雷擊試驗，同時，采用電熱耦合模型分析了雷電作用下巖土體的熱效應[26]。另外，在防雷裝置設計工作中，雷電沖擊作用下，土中接地極等防雷裝置的雷電沖擊散流效應和防雷效益受其周圍土體的影響極大，與土體相關參數息息相關，然而鑒于學科間的差異性和局限性，目前雷電作用下巖土工程和電氣工程的交叉融合研究仍十分欠缺。因此，探究雷擊對土體的熱力等效應、雷電致災機制等研究可為降低地下雷電災害損失、指導防雷系統設計等工作提供一定指導，對巖土工程防雷減災具有重要意義。

本文從雷電沖擊波理論與電磁場理論出發，通過理論分析得到雷電過程中對土體引起的雷電沖擊壓力與脈沖電磁力，并采用有限元方法對土體在雷擊沖擊作用下的動態變形特性以及應力分布特征進行求解分析，確定土中變形等動態發展特征，在此基礎上深入探究土體和雷電流等因素對雷擊土體的動態變形響應影響規律。

1 雷擊的力學沖擊描述

1.1 基本假定

1)假設雷電通道垂直于土體表面，且土體表面為電絕緣平面，無電流流出。

2)土體為連續均質的各向同性材料，忽略蠕變和松弛效應。

3)雷擊沖擊作用時間很短，假定雷擊期間土體不發生固結。

1.2 雷電沖擊壓力

雷擊放電過程中，由地表發出的回擊電流可在幾微秒內通過焦耳熱的形式使空氣急劇升溫[27]，形成可見的雷電等離子體通道。此過程中，在這極短時間內不足以使空氣中雷電離子體通道內的粒子(電子、離子等)密度顯著降低，這種溫度的快速上升導致等離子體通道內壓力大幅增加并產生超壓，進而形成以超音速徑向膨脹的圓柱形沖擊波。假設沖擊波內部等離子體的流動是等熵流動的且沖擊波內空氣為比熱恒定的理想氣體，則雷電沖擊波內的控制方程為等離子體流動的運動方程、連續性方程和熱力學流動方程：

(1)

(2)

(3)

式中：u、ρ和p分別為t時刻坐標為r處的速度、密度和壓力；κ=cp/cV為空氣的質量熱容比，對于理想氣體，κ=1.4。

在沖擊波波陣面上，其以一定的速度移動，是一個強斷面，前后的介質狀態是不同的，根據Rankine-Hugoniot條件，由前后介質間的質量、動量和能量守恒關系，可以得到速度、密度和壓力邊界條件：

(4)

(5)

(6)

式中：u1、ρ1和p1分別為波陣面后的速度、密度和壓力；ρ0、p0和c0為波陣面前無擾動空氣的密度、壓力和聲速，其中，ρ0=1.293 kg/m3；U=dR/dt為波陣面運動速度。

根據Lin[28]、Karch等[29]的研究，可采用自相似方法對式(1)～(6)進行求解，沖擊波前沿半徑R(t)可由下式給出：

(7)

雷電沖擊波的壓力p=(η,t)可由下式描述：

p(η,t)=0.180(ρ0E0)1/2f(η)t-1

(8)

式中：η=r/R(t)為無量綱半徑；f(η)為無量綱壓力函數，描述了從起始點η=0到沖擊波前沿η=1間的柱狀壓力分布，由于壓力的自相似性，對任意時刻t都是有效的。根據Lin[28]的研究，其沖擊波內各處的無量綱壓力分布如圖1所示，將各f(η)點擬合可得

圖1 f(η)-η擬合曲線

f(η)=5.751×10-6e11.74η+0.434 8

(9)

至此，雷電沖擊壓力可表示為

p(η,t)=0.180(ρ0E0)1/2(5.751×

10-6e11.74η+0.434 8)t-1

(10)

1.3 雷電脈沖電磁力

雷擊發生時伴有強大的雷電流，其會在周圍產生強烈的電磁場，土中將產生電磁力。Lee等[16]采用簡化的二維麥克斯韋方程推導了雷電流在CFRP平面板內產生的脈沖電磁力表達式，適用于二維平面結構，但其并不適用于具有三維空間結構的土體，因此有必要進行雷電流在三維土體中的脈沖電磁力分析。如圖2所示，取土中任一縱截面，以雷擊點為坐標原點，建立柱坐標系，將雷電流看作位于原點處的點電流源，假設雷電流僅在土中傳導，即土體上表面為電絕緣平面。根據電磁場理論[30]，土中將產生如圖2虛線的半球形等勢面，其各點電勢V為

圖2 雷擊電磁力計算示意

(11)

式中：I為雷電電流，σ為土體電導率。

由J=-σ×?V可得

(12)

由于土體表面為電絕緣平面，無電流穿過，則有

J(r,z)=0,z<0

(13)

根據麥克斯韋方程，其各點磁場強度H的散度等于電流密度J，即?×H=J，其在柱坐標系內可寫為

(14)

由于雷電流引起的磁場強度H呈對稱分布，其與θ無關，將式(12)、(13)代入式(14)，?×H=J可化簡為

(15)

(16)

求解式(15)和(16)，可得

(17)

(18)

式中：C1、C2為待定常數，i、e分別表示土體內部與外部空氣區域。

代入邊界條件

(19)

(20)

可得

C1=1，C2=I/(2π)

(21)

則式(17)、(18)可寫成

(22)

(23)

土中電磁力可表示為電流密度J與磁感應強度B的叉乘，則在z≥0土中區域，電磁體積力可計算為

Fzez=Jr(r,z)er×Beθ=Jr(r,z)er×uru0Heθ=

(24)

Frer=Jz(r,z)ez×Beθ=Jz(r,z)ez×uru0Heθ=

(25)

式中：u0=4π×10-7H/m為空氣磁導率；ur為相對磁導率，對于土體，ur=1。

為驗證電磁力推導的正確性，將本文電磁力結果與文獻[16]所得雷電作用在CFRP平面板內產生的脈沖電磁力結果進行對比。由于文獻[16]僅得出電磁力的豎向壓力表達式，只適用于二維平面分析，而本文土體中電磁力基于三維空間坐標，可進行退化對比，將式(25)沿z方向積分，可得本文脈沖電磁力的豎向分量表達式

(26)

式(26)結果與文獻[16]豎向壓力電磁力結果相同，因此，可以驗證本文雷電作用下土中脈沖電磁力算法的正確與合理性。

1.4 標準雷電流波形

標準雷電流波形的主要參數為雷電流峰值Im、波前時間T1和半峰值時間T2，如圖3所示，雷電流峰值Im表示雷電全過程中的最大電流值，波前時間T1為電流峰值的10%～90%時間間隔的1.25倍，描述雷電流從開始到峰值的時間，半峰值時間T2為電流隨時間衰減到峰值50%的時間。雙指數波形是國際雷電防護標準IEC62305[31]和建筑物防雷設計規范GB 50057—2010規定的雷電流標準波形，其描述的波形特征和實際測量波形相近，可以很好地反映雷電流的特征值[32]，且形式簡單，便于計算，本文采用雙指數波形模型作為雷電流表達式，即

圖3 雷電電流波形

I(t)=λIm(e-αt-e-βt)

(27)

式中：I(t)為雷電流的瞬時值，Im為雷電流峰值，λ為波峰修正系數，α為波前系數，β為波尾系數。

2 有限元模型的建立

基于有限元方法的COMSOL Multiphysics軟件具有強大的多物理場仿真功能，以靈活開放、易用高效著稱，能夠很方便地定義出本文由多變量控制的荷載函數。采用COMSOL Multiphysics中的固體力學模塊建立雷擊作用下土體的動力特性瞬態模型。

2.1 土體本構模型

雷擊沖擊作用下，土體將產生不可恢復的塑性變形，應采用彈塑性本構模型，COMSOL Multiphysics中提供了匹配Mohr-Coulomb準則的Drucker-Prager外接圓(DP1)準則[33]，本文采用該模型，其形式為

(28)

式中：I1、J2分別為應力張量第一不變量和偏應力第二不變量；α、k為與土體黏聚力c和內摩擦角φ相關的常數，α和k的不同取值在π平面內表示不同的圓，對于DP1準則，α和k表示為

(29)

2.2 模型建立

在假定土體為各向同性材料的條件下，雷電沖擊作用呈軸對稱分布，為簡化模型與計算量，可采用二維軸對稱模型。如圖4所示，土體幾何尺寸為1.0 m×0.5 m，在左側邊界上設置為軸對稱邊界條件，由于雷擊作用是一個瞬態過程，模型計算時還應防止土中壓力波和剪切波在邊界處發生反射，造成計算結果反復震蕩，影響計算精度，將模型下側以及右側邊界設置為低反射邊界以截斷計算域，可以很好地減少邊界反射造成的影響。采用國際雷電防護標準IEC62305[31]推薦的T1/T2=10/350 μs的標準雷電流波形進行雷電作用模擬，雷電流峰值Im=100 kA，其雷電流與土體參數如表1所示。

表1 雷電流和土體參數

圖4 有限元計算模型

另外，定義由位置變量r、z以及時間變量t控制的雷電沖擊壓力與脈沖電磁力函數表達式，并將雷電沖擊壓力作為邊界載荷施加在模型上側邊界，將雷電脈沖電磁力作為體載荷添加在模型左上方矩陣區域內。經對雷電沖擊壓力和電磁力表達式分析可知，其影響范圍主要集中于雷擊點附近，沖擊力大小隨距離增加而迅速減小，本文設置雷電沖擊壓力邊界范圍為0.35 m，雷電電磁力作用矩形的尺寸為0.35 m×0.35 m，以模擬雷電對土體的力學沖擊作用。采用映射方式劃分四邊形網格，并對荷載作用區域進行局部加密處理，以提高雷擊計算結果的精度，而外部影響較小區域網格相較疏松，以提高模型計算速度。瞬態計算時間步長采用Δt=0.000 1 ms，計算時間共3 ms。

3 雷擊動力特性分析

圖5和6分別為雷擊點以下(即r=0)不同深度處各點的豎向位移時程曲線和Mises應力時程曲線。由圖5可以看出，位移變化較大各點均表現為前期迅速變化，而后變形速度逐漸減緩，最終約在1.5 ms附近達到穩定狀態，這主要是因為雷電沖擊作用載荷與爆炸作用類似，呈現出急劇變大而后又迅速減小的特點。同時可以看出，豎向變形隨深度的增加衰減得很快，在雷擊點即z=0 mm的最終穩定位移最大，達到30.51 mm，而雷擊點以下各點的位移迅速減小，在z=-66.7 mm處，其穩定位移僅為0.97 mm，說明此處受影響很小，可認為此時雷擊有效影響深度為66.7 mm。

由圖6可知，雷擊作用下各深度處的Mises應力在雷擊初期陡增后隨時間推移又不斷衰減，整體表現出一個主峰值和多個較小峰值，約在1.5 ms附近衰減至0，而后保持不變，這與圖5中豎向位移隨時間變化類似，其應力與位移保持同步性。Mises應力曲線主峰值隨深度增大也變化明顯，在雷擊點處其值為12.43 MPa，而在z=-16.7 mm處減小為3.26 MPa。此外，從圖6局部放大圖可以看到，不同深度處的Mises應力最初上升時間約為0、0.005、0.02、0.08、0.15、0.20、0.25 ms，隨深度增加逐漸延后，反映了應力波在土中的分布特征。

圖5 不同深度處的位移時程曲線

圖6 不同深度處的Mises應力時程曲線

圖7為同時考慮雷電沖擊壓力和電磁力、僅考慮雷電沖擊壓力、僅考慮電磁力3種情況下土體表面最終穩定時的豎向位移分布。可以看出，3種情況下各處的豎向位移沿雷擊中心點呈對稱分布，并在中心點處取得最大值，隨著離雷擊點距離r增大豎向變形不斷減小，并在距離30 mm附近存在土體少量隆起，這是由于在強大的沖擊作用下土體內部會產生側向擠壓作用，使得距雷擊點周圍一定距離處的土體有向上運動趨勢，造成隆起變形。在隆起區域外側的土體表面處，豎向變形逐漸減小為0，此區域的土體位移影響很小。3種情況下表面土體位移隨r變化規律相同，但其各點處數值明顯不同，對于雷擊點處，即豎向位移最大處，其變形值分別為30.51、27.46和2.96 mm，比較這3種情況各自對土體變形的影響可知，電磁力造成的影響僅能達到雷電沖擊壓力的10.78%，但從圖7中同時考慮這兩種效應的曲線可以看出，電磁力對最終土體表面位移仍有較大影響。

圖7 土體表面豎向位移

4 影響因素分析

4.1 土體參數的正交分析

雷電作用下的土體變形是個多因素問題，為探究不同土體力學參數下地表的變形特性及主次影響因素，采用正交分析法[34]進行多因素分析，選取的主要影響因素為彈性模量E、泊松比μ、黏聚力c和內摩擦角φ，每個參數都有4個不同的水平，各因素的水平如表2所示。

表2 不同水平下的土體因素

根據正交分析法的設計原則，由于每個因素有4個水平，選用L16(45)的正交表進行研究計算，該表可最多安排5個參數，將所選的4個因素安排在正交表前4列，剩余的第5列不安排因素，作為空白列檢驗誤差水平。各方案設計與地表最大豎向變形計算結果如表3所示，共有16種方案，相對于全面計算方案的256種，工作量顯著減少。

根據表3中的結果進行極差分析，結果如表4和圖8所示。通過極差比較，各因素對最大變形的影響從主到次的排序為內摩擦角φ、彈性模量E、泊松比μ、黏聚力c。由圖8可以看出，隨著各因素水平的增加，雷擊下的地表最大變形均呈不同程度的減小趨勢，其中，內摩擦角的影響較為明顯，而黏聚力的影響相對較小。

表3 正交分析方案和計算結果

表4 最大變形的極差分析

圖8 各影響因素水平與最大變形均值趨勢

對表3結果進行方差分析，探究各因素顯著性水平，結果見5所示。選取顯著性水平α分別為0.01、0.05和0.1，由于各因素的自由度為3，誤差項的自由度為15-12=3，根據F分布可知F0.01(3,3)=29.46，F0.05(3,3)=9.28，F0.1(3,3)=5.39。若計算的F值大于F0.01時，表示該因素對地表最大變形的影響高度顯著，記為“***”；若計算的F值大于F0.05但不大于F0.01時，表示該因素的影響顯著，記為“**”；若計算的F值大于F0.1但不大于F0.05時，表示該因素對地表最大變形有一定影響，記為“*”；若計算的F值小于F0.05時，表示該因素的影響不顯著。表5的結果表明：對于雷電沖擊作用下的地表最大變形，內摩擦角φ的影響高度顯著，這是因為土體變形與本構模型有關，本文采用的DP1準則在主應力空間和π平面上的屈服面半徑與內摩擦角密切相關，內摩擦角越小時，主應力空間和π平面上的屈服面半徑就越小，彈性范圍越窄，此時相同應力狀態下土體更易進入塑性狀態，土體變形會明顯增加，因此，內摩擦角的影響最為顯著。

表5 最大變形的方差分析

4.2 內摩擦角φ的影響

由4.1節正交分析法可知，雷電沖擊作用下，土體內摩擦角φ對土體變形的影響最為顯著，本節重點考察內摩擦角的影響規律，采用內摩擦角φ=10°、20°、30°、40° 4種情況分別進行模擬計算，其他雷電流和土體參數同表1一致。

圖9為不同內摩擦角φ條件下的土體表面豎向位移結果。從宏觀上看，不同內摩擦角條件下的表面整體變形規律與前述相同，即沿雷擊中心點呈對稱分布，且在雷擊點處變形最大，而后兩側位移逐漸減小，并伴有少量隆起，最終最外側土體幾乎不受影響；從數值上看，內摩擦角的改變對土體變形影響較大，其值越小時，雷擊點的豎向位移越大，且土體周圍向上隆起程度也明顯變大，當φ=10°時，其最大豎向位移為48.78 mm，最大隆起高度為3.99 mm，而當φ=40°時，最大豎向位移僅為12.98 mm，而此時周圍土體幾乎沒有隆起現象；從整體變形曲線來看，內摩擦角較大時，其整體變形曲線仍為深凹狀，而隨著內摩擦角的減小逐漸呈扁平狀，此時雷擊點周圍一定寬度內的豎向變形大致相同。

圖9 不同內摩擦角下的土體表面豎向位移

圖10為不同內摩擦角φ下雷擊點處即變形最大處的豎向位移時程曲線。可以看出，在不同內摩擦角下，雷電作用初期土體變形基本一致，而后隨著時間的推移，整體位移呈現先增大后減小而后逐步趨于穩定的趨勢，且內摩擦角越小時其位移達到穩定所需時間越長，這是由于當摩擦角較小時，土體屈服面半徑較小，彈性范圍較窄，更容易進入塑性狀態，在雷電沖擊作用數值達到峰值后土體依然處于塑性屈服狀態，其變形將繼續變大，因此，達到穩定耗時最長。另外，在雷電沖擊下土體先出現彈塑性變形，而后隨著沖擊作用減小直至消失，其彈性變形逐漸恢復，最終穩定時僅剩余塑性變形，這在內摩擦角較大時更加明顯。這是因為摩擦角較大時土體屈服面半徑較大，彈性范圍較廣，產生的彈性變形更多，摩擦角較大時回彈現象更明顯。如從φ=40°條件下的豎向位移時程曲線可以明顯看出，在0.12 ms后出現顯著的彈性位移恢復現象，整個雷擊作用中，該點最大豎向位移為16.78 mm，最終穩定位移為12.95 mm，恢復了22.82%，而對于其他內摩擦角條件下，其回彈量不足5%，這也是致使圖9中呈扁平狀的重要原因。

圖10 不同內摩擦角φ下的位移時程曲線

4.3 雷電流峰值Im的影響

雷電流峰值Im作為衡量雷電能量的重要指標，與雷電沖擊壓力和電磁力的大小直接相關。為探究雷電流峰值對土體雷擊變形的影響，保持表1其他參數不變，分別選用Im=50、100和150 kA進行模擬計算，并分別考量在此過程中脈沖電磁力的影響，各情況下的表面豎向位移結果如圖11所示。可以看出，其表面變形規律與前所述相同，并呈深凹狀。此外，雷電流峰值越大時其變形量越大，并且從僅考慮雷電沖擊壓力與同時考慮雷電沖擊壓力和電磁力對比結果中發現，同時考慮這兩種沖擊作用時土體變形量更大，其雷擊點豎向位移分別比僅考慮雷電沖擊壓力時大7.45%、11.11%和14.53%，說明雷電流峰值越大時電磁力對土體變形影響越大，在大電流作用下應著重考慮其影響。

圖11 不同雷電流峰值Im下的土體表面豎向位移

5 結 論

應用雷電沖擊波理論，并基于麥克斯韋方程推導了雷電在土中產生的電磁力表達式，以此建立雷電沖擊壓力和電磁力作用下的土體變形特性有限元瞬態模型，主要結論如下：

1)雷擊沖擊作用下土體先產生彈塑性變形，隨著雷擊的結束其彈性變形逐漸恢復，最終僅剩余塑性變形，其整體變形規律表現為雷擊初期表面土體急劇變形，而后變形速率逐漸減緩，部分土體回彈后趨于穩定，此外，土中應力分布特征規律與應力波傳播規律類似。

2)通過正交分析可以得到土體各力學參數對雷擊效應的影響大小，內摩擦角影響較為顯著，彈性模量影響次之，泊松比和黏聚力的影響相對較小。隨著內摩擦角和彈性模量的減小，其變形均表現為增大的趨勢，并伴有一定隆起，且隨著內摩擦角的增大，表面豎向位移曲線由深凹狀逐漸過渡為扁平狀。

3)雷電流峰值是衡量雷電荷載大小的重要因素，隨著雷電流峰值的不斷增大，更大的雷電沖擊力將引起更大的土體變形，并且雷電脈沖電磁力對整體變形的影響貢獻也逐漸提升，在進行較大雷電流作用對土體的沖擊效應研究中應予以重視。

4)在工程防雷設計中，就特定因素而言，可增大布置在防雷裝置周圍的土體內摩擦角等力學參數以減小土體變形，減少防雷裝置擾動影響，保證防雷裝置的雷電沖擊散流效應和防雷效益。

5)本文研究成果對建筑及電力工程地基基礎防雷接地系統設計和改造以及保障工程建設區域安全具有一定意義。另外，評估雷擊后土體力學性狀也為從巖土工程角度研究防雷減災提供了新思路。