三維斷裂彈性動力學的改進無單元Galerkin法

景永強 彭妙娟

摘要：為提高斷裂彈性動力學問題數值計算的精度，避免出現病態或奇異方程組，基于改進的移動最小二乘法建立三維彈性動力學問題的積分弱形式，采用罰函數法施加位移邊界條件，引入隱式時間積分并且結合三維斷裂力學的形函數考慮裂紋尖端的奇異性，探究將改進的無單元Galerkin（improved element-free Galerkin， IEFG）法用于斷裂彈性動力學問題的數值計算。通過懸臂梁、柱和矩形板等3個算例，討論節點分布、影響域比例參數、罰因子和時間步長等參數對計算精度的影響，證明IEFG法用于求解三維斷裂彈性動力學問題的正確性和有效性。

關鍵詞：

彈性動力學; 斷裂力學; 裂紋; 改進的移動最小二乘法; 形函數; 隱式時間積分

中圖分類號：O241.82; TB115.1

文獻標志碼：B

Improved element-free Galerkin method for

three-dimensional fracture elasto-dynamics

JING Yongqiang， PENG Miaojuan

（Department of Civil Engineering， School of Mechanics and Engineering Science， Shanghai University， Shanghai 200444， China）

Abstract：

To improve the accuracy of the numerical calculation of fracture elasto-dynamics， and to avoid the ill conditioned or singular equations， the improved element-free Galerkin（IEFG） method is applied to the numerical calculation of fracture elasto-dynamics. Based on the improved moving least square method， the integral weak form of the three-dimensional elasto-dynamics problem is established. The penalty function method is used to set the displacement boundary conditions. The singularity of crack tip is considered by introducing the implicit time integration and the shape function of three-dimensional fracture mechanics. Three examples（include the cantilever beam， the column and the square plate） are selected to analyze the influence of computing parameters（include the node distribution， the influence domain scale parameter， the penalty factor and the time step） on the calculation accuracy. It is proved that the IEFG method is correct and effective for solving three-dimensional fracture elasto-dynamics problems.

Key words：

elasto-dynamics; fracture mechanics; crack; improved moving least square method; shape function; implicit time integration

0 引 言

斷裂彈性動力學問題是斷裂力學與彈性動力學耦合的問題，傳統的解析方法難以求解，目前可采用的數值方法主要有有限元法和無單元法等。采用有限元法求解斷裂力學問題時，在裂紋的發生擴展過程中需要重新劃分網格，計算精度較低。無單元法無須劃分網格，只需要節點信息。對于含有裂紋的構件，為提高計算精度，無單元法可以在裂紋處布置加密節點，避免對網格和單元依賴的問題，在求解斷裂彈性動力學問題時具有一定的優勢。

BELYTSCHKO等[1]采用無單元Galerkin法研究三維幾何和材料非線性的動力學問題。CHEN等[2]采用無單元Galerkin法研究材料高速撞擊發生彈塑性大變形的問題。LEE等[3]采用無單元Galerkin法對裂紋的奇點進行不連續性建模，提出一種改進的裂紋分析技術。RAO等[4]提出一種改進的無單元法，用于均勻、等向、非線性彈性二維固體裂縫的斷裂分析。

無單元Galerkin法利用移動最小二乘法構造試函數，因此具有較高的計算精度，但計算過程中形成的方程組有時是病態的或奇異的。陳美娟等[5]提出改進的移動最小二乘法，基函數采用正交函數，可有效避免病態或奇異方程組的產生，不需要求解逆矩陣，可提高計算效率。ZHANG等[6]提出改進的無單元Galerkin（improved element-free Galerkin， IEFG）法，并用于解決三維勢問題、波動方程、斷裂力學和瞬態熱傳導問題。CHENG等[7]提出針對三維彈性力學、彈性動力學等問題的IEFG法。SUN等[8]研究共平面方形裂紋，對等向彈性體中二維和三維共平面方形裂紋的相互作用進行數值模擬。PENG等[9]在改進的移動最小二乘法的基礎上，提出針對三維黏彈性問題的IEFG法。ZHANG等[10]基于改進的移動最小二乘法，在裂縫尖端使用基函數，提出一種可以解決二維斷裂問題的IEFG法。魏慶節等[11]基于改進的移動最小二乘法建立三維彈性動力學問題的形函數，結合三維彈性動力學的Galerkin積分弱形式，采用罰函數法施加位移邊界條件，引入隱式時間積分，并提出三維彈性動力學的IEFG法。蔡小杰等[12]基于改進的移動最小二乘法，提出針對彈塑性大變形問題的IEFG法，討論權函數、影響域比例參數、罰因子、節點數和迭代步數等參數對計算精度的影響。鄒詩瑩等[13]在IEFG法的基礎上，將能反映裂紋尖端附近應力奇異性的特征項引入改進的移動最小二乘法的基函數中，將斷裂力學與IEFG法相結合，用于研究線彈性斷裂力學問題，并對含反射裂縫的機場復合路面層狀體系結構進行數值分析。

劉學聰等[14]基于擴展有限元法提出一種新的裂紋尖端加強函數，組合傳統的函數基，繼承傳統附加函數的特性，使節點的奇異附加自由度減少為2個，可減小總矩陣的規模，提高計算效率。石路楊等[15]引入裂紋交叉匯合加強函數，分析多裂紋的交叉匯合過程，在裂紋附近區域使用廣義形函數，并引入線性遞增函數消除混合單元，建立求解多裂紋擴展的擴展有限元法，可有效提高裂紋附近的分析精度。楊永濤等[16]基于包含裂紋尖端增強函數的數值流形方法，采用Newmark隱式動力學算法進行時間積分，重點研究動力載荷作用下裂紋動態應力強度因子的求解方法。張洪武等[17]基于精細網格的傳統有限元分析方法，針對非均質飽和多孔介質彈塑性動力問題，提出一種廣義耦合擴展多尺度有限元法。

分析以上研究發現，在求解斷裂力學問題時IEFG法優于有限元法。目前，IEFG法在三維彈性力學、彈性動力學、彈塑性大變形、黏彈性和線彈性斷裂力學等方面的研究較多，針對斷裂彈性動力學問題，IEFG法的研究很少。本文基于改進的移動最小二乘法，建立三維彈性動力學問題的積分弱形式，采用罰函數法施加位移邊界條件，引入隱式時間積分，并且結合三維斷裂力學的形函數考慮裂紋尖端的奇異性。采用改進的移動最小二乘法可避免病態或奇異方程組，保證IEFG法的計算精度。三維斷裂彈性動力學IEFG法不僅要考慮三維彈性動力學中加速度和阻尼的影響，還要考慮斷裂力學裂紋尖端的奇異性。本文結合改進的移動最小二乘法、彈性動力學和斷裂力學的無單元Galerkin法，建立三維斷裂彈性動力學的IEFG法的離散方程，并推導相應的計算公式，采用懸臂梁、柱和方形板等3個數值算例證明三維斷裂彈性動力學IEFG法的正確性和有效性。

1 三維斷裂彈性動力學基本方程

三維線彈性動力學的控制方程為

LTσ+b=ρu··+μu·， x∈Ω

（1）

式中：L為微分算子矩陣;σ為應力向量;b為外力向量;ρ為質量密度;u為位移向量;μ為阻尼系數;x為點的集合。L、σ、b和u分別定義為

L=x1000x2000x3x2x10x30x10x3x2

（2）

σ=[σ11 σ22 σ33 σ12 σ13 σ23]T

（3）

b=[b1（t） b2（t） b3（t）]T

（4）

u=[u1（t） u2（t） u3（t）]T

（5）

式中：xi（i=1，2，3）為點x緊支域內的點;σij（i，j=1，2，3）為不同方向的應力;t為時間。

三維線彈性動力學控制方程的邊界條件為

ui-i=0， i=1，2，3， x∈Γu

（6）

njσij-ti=0， i，j=1，2，3， x∈Γt

（7）

式中：i為邊界Γu上的已知位移;nj為邊界Γt的外法線方向的余弦;ti為邊界Γt上的已知面力。Γu∪Γt=Γ，Γ為求解域Ω的邊界。

三維線彈性動力學控制方程的初始條件為

ui（x，t0）=ui0（x）， x∈Ω（8）

u·i（x，t0）=vi0（x）， x∈Ω （9）

式中：vi0和ui0分別為點x的初始速度和初始位移;t0為初始時間。

三維線彈性動力學控制方程的本構方程為

σ=Dε

（10）

其中：

ε=Lu

（11）

D=E（1-v）（1+v）（1-2v）1v1-vv1-v000v1-v0v1-v000v1-vv1-v00000001-2v2（1-v）0000001-2v2（1-v）0000001-2v2（1-v）

（12）

式（11）即為三維線彈性動力學控制方程的幾何方程。

為得到彈性動力學的Galerkin弱形式，采用罰函數法施加位移邊界條件，即

∫ΩδuT·ρu··dΩ+ ∫ΩδuT·μu·dΩ+∫ΩδεT·σdΩ+

∫ΓtδuT·bdΩ-∫ΓtδuT·tdΩ+

α∫ΓtδuT·S（u-）dΓ=0

（13）

式中：α為罰因子;

=[1 2 3]T （14）

t=[t1 t2 t3]T（15）

將式（10）和（11）代入式（13），計算得到罰函數法的可行域為

S=S1000S2000S3

（16）

∫ΩδuT·ρu··dΩ+∫ΩδuT·μu·dΩ+∫Ωδ（Lu）T·（Lu）dΩ+∫ΩδuT·bdΩ-∫ΓtδuT·tdΩ+α∫ΓtδuT·S（u-）dΓ=0（17）

當x1、x2、x3有位移約束時，對應的S1、S2、S3等于1，否則等于0。

2 三維斷裂彈性動力學的IEFG法

利用改進的移動最小二乘法建立逼近函數。由文獻[13]中裂紋尖端的應力場表達式可知，裂紋尖端附近的位移與特征項r（r為點x到裂紋尖端的距離）成正比，應力具有1/r階奇異性。為有效計算裂紋尖端場的應力數值解，在基函數中加入擴展函數項，由改進的移動最小二乘法得到形函數。基于三維彈性動力學IEFG法，在基函數中加入文獻[13]中裂紋尖端位移場表達式中的重要項和重要的梯度項，使IEFG法的計算結果更精確。

結合文獻[13]得到三維斷裂力學IEFG法的試函數為

pT（x）=1 x1 x2 x3 rcos θ2 rsin θ2

rsin θ2sin θ rcos θ2sin θ

（18）

PI（x）=[w（x-x1）p（x1） w（x-x2）p（x2） …

w（x-xn）p（xn）]

（19）

式中：θ為裂紋的擴展角;n為影響域覆蓋點x的節點數;xn為點x的緊支域內的節點;w（x-xn）為具有緊支集特性的權函數。

類似文獻[11]中移動最小二乘法的形函數推導，可得到三維斷裂彈性動力學的逼近函數為

uα，h（x）=nI=1pT（x）A-1（x）PI（x）uα（xI）

（20）

式中：A-1（x）為由pi構成的對角矩陣。

三維斷裂彈性動力學的形函數為

ΦI（x）=pT（x）A-1（x）PI（x）

（21）

將式（21）代入式（20），可得

uα，h（x）=nI=1ΦI（x）uα（xI）

（22）

采用局部基函數擴展法，在基函數中加入特征項r，即

pT（x）=[1 x1 x2 x3 r]

（23）

將求解域離散為有限個節點，節點總數為M。利用改進的移動最小二乘法建立逼近函數，可以得到位移的逼近函數為

ui（x，t）=nI=1ΦI（x）uI（xI，t）， i=1，2，3

（24）

即 u（x，t）=nI=1ΦI（x）uI

（25）

將式（24）和（25）代入式（17）可得

∫ΩδnI=1ΦIuITρnJ=1ΦJu··JdΩ+∫ΩδnI=1ΦIuITμnJ=1ΦJu·JdΩ+∫ΩδnI=1LuITDnJ=1ΦJLuJdΩ-

∫ΩδnI=1ΦIuITbdΩ-∫ΓtδnI=1ΦIuITtdΓ+∫ΓtδnI=1ΦIuITSnJ=1ΦJu··JdΓ+∫ΓtδnI=1ΦIuITSu··dΓ=0

（26）

對式（26）進行積分運算，可以得到三維斷裂彈性動力學的離散方程為

MU··（t）+CU·（t）+（K+Kα）U（t）=F（t）+Fα

（27）

其中：

M=（Mij），Mij=∫ΩΦiρΦjdΩ， i，j=1，2，…，M

（28）

U=[uT1（t） uT2（t） … uTM（t）]T

（29）

C=（Cij），Cij=∫ΩΦiμΦjdΩ， i，j=1，2，…，M

（30）

K=（Kij），Kij=∫ΩBTiDBjdΩ， i，j=1，2，…，M

（31）

Bi=Φi，2000Φi，2000Φi，3Φi，2Φi，10Φi，30Φi，10Φi，3Φi，2， i=1，2，…，M

（32）

Kα=（Kij，α），Kij，α=α∫ΓtΦiSΦjdΓ，

i，j=1，2，…，M

（33）

F=（Fij），Fij=∫ΩΦibdΩ+∫ΓtΦjtdΓ，

i，j=1，2，…，M

（34）

Fα=（Fi，α），Fi，α=α∫ΓtΦiSdΓ，i=1，2，…，M

（35）

3 隱式時間積分

隱式時間積分方法可實現求解目標在大時間步長內的無條件穩定，被廣泛應用于各種動力學問題中。本文采用Newmark-β法對運動方程進行時間離散。假設時刻t的時間域[0，T]被離散為n個時間步（即Δt=T/n），位移Ut及其對應導數U·t、U··t的值均已知，對U·t+Δt和U··t+Δt進行二階Taylor展開可得

Ut+Δt=Ut+ΔtU·t+12（1-β2）Δt2U··t+12β2Δt2U··t+Δt

（36）

U·t+Δt=U·t+（1-β1）ΔtU··t+β1ΔtU··t+Δt

（37）

U··t+Δt=2β2Δt2（Ut+Δt-Ut）-2β2ΔtU·t-1β2-1U··t

（38）

式中：β1和β2均為Newmark參數，分別取β1=3/2，β2=8/5。

令α1=2β2Δt2、α2=2β2Δt、α3=1β2-1，可得

U··t+Δt=α1（Ut+Δt-Ut）-α2U·t-α3U··t

（39）

將式（39）代入式（37）可得

U·t+Δt=β1α2（Ut+Δt-Ut）+1-2β1β2U·t+

1-β1β2ΔtU··t

（40）

根據隱式時間積分方法，式（27）可寫成

MU··t+Δt+CU·t+Δt+（K+Kα）Ut+Δt=Ft+Δt+Fα

（41）

將式（39）和（40）代入式（41）可得

（α1M+β1α2C+K+Kα）Ut+Δt=

Ft+Δt+Fα+M（α1Ut+α2U·t+α3U··t）+

Cβ1α1Ut+2β1β2-1U·t+β1β2-1ΔtU··t

（42）

4 數值算例

選擇懸臂梁、柱和矩形板等3個算例，驗證IEFG法的正確性和有效性。3個算例的權函數均為三次樣條權函數，將計算結果與Abaqus軟件有限元法計算結果進行比較。

為分析IEFG法的精確度，定義均方差為

eI=1NNi=1（uIEFG，i-uAbaqus，i）

（43）

式中：uIEFG，i為IEFG法的計算結果;uAbaqus，i為Abaqus的計算結果;N為節點數。

為討論Abaqus解的準確性和收斂性，定義方差

eA=1NNk=1（uJ，k-IJ，k）2+（uI，k-IJ，k）2

（44）

式中：uI，k和uJ，k分別為節點分布I和J下第k個節點的數值解;IJ，k為uI，k和uJ，k的平均值。

4.1 受突加線性分布載荷作用的懸臂梁

選取某典型懸臂梁，上、下表面自由，一端固定，另一端受突加線性分布載荷作用，見圖1，其中a=120 mm、b=120 mm、c=360 mm、r=30 mm。材料參數分別取ρ=2 405 kg/m3、E=20.69 GPa、ν=0.15，g取9.8 m/s2。在x3=-60 mm和x3=60 mm處施加突加載荷P=1×103 N。

利用有限元軟件Abaqus對該算例進行計算，有限元網格采用20節點二次六面體單元C3D20R，分別取9×9×15、15×15×31、15×15×51、17×17×65和33×33×81等5種節點分布，懸臂梁自由端Q點（x1=360 mm、x2=-60 mm）位移數值解的方差見圖2。由此可知，隨著網格的細化，方差逐漸趨于0，說明Abaqus的數值解收斂。與IEFG方法的數值解進行對比時選擇節點分布為33×33×81，Abaqus中懸臂梁的網格劃分結果見圖3。

根據三維斷裂彈性動力學IEFG法的離散方程，采用與Abaqus計算時同樣的節點分布，利用MATLAB對懸臂梁Q點位移進行數值計算，不同節點分布時懸臂梁位移的均方差結果見圖4。由此可知，隨著節點數增加，IEFG法的均方差逐漸減小，計算精度逐漸提高。

影響域比例參數dmax在1.5～3.5范圍內取不同值時，懸臂梁位移均方差的變化見圖5。由此可知，當dmax在2.5～3.0范圍內時，懸臂梁位移的均方差較小，說明此時數值解的計算精度較高。

給定影響域比例參數dmax，取不同時間步長Δt，計算懸臂梁位移的均方差，結果見圖6。由此可知，隨著時間步長的增大，均方差逐漸增大。當時間步長為1.0×10-6 s時均方差較小，說明此時數值解的計算精度較高。

給定影響域比例參數dmax和時間步長，取不同罰因子α，計算懸臂梁位移的均方差，結果見圖7。由此可知，隨著罰因子的增大，均方差呈類似拋物線分布。當罰因子α=2.5×1015時均方差較小，說明此時數值解的計算精度較高。

取dmax=2.7、Δt=1.0×10-6 s、α=2.5×1015，Abaqus中節點分布取33×33×81，計算圖1中自由端Q點的位移u1、u3和應力σ隨時間的變化，見圖8～10。由此可知，IEFG法與有限元法的計算結果吻合，證明IEFG法解決三維斷裂彈性動力問題的正確性和有效性。

4.2 受突加載荷作用的柱

選取單邊帶有半圓型裂紋的柱形試件，裂紋位于試件中部，裂紋半徑r=5 cm，試件下端固定，上端施加單向突加載荷作用，見圖11。突加載荷P=20 kPa，材料和幾何參數取E=210 GPa、ν=0.35、a=30 cm、b=30 cm、c=120 cm。

利用Abaqus對柱位移進行計算，有限元網格采用20節點二次六面體單元C3D20R，取9×9×17、9×9×31、17×17×64、21×21×71和41×41×81等5種節點分布，計算得到Q點（x1=-15 cm、x3=120 cm）位移的方差，見圖12。隨著節點數增加，Abaqus數值解的方差趨于0，說明Abaqus的數值解收斂。與IEFG法的數值解進行對比時選擇節點分布為41×41×81，Abaqus中柱的網格劃分結果見圖13。

根據三維斷裂彈性動力學IEFG法的離散方程，取dmax=2.8、Δt=1.0×10-6 s、α=2.5×1015，利用MATLAB對柱位移進行數值計算，圖11中Q點的位移u1、u3和應力σ隨時間的變化見圖14～16。由此可知，IEFG法與有限元法計算結果吻合，證明IEFG法的正確性和有效性。

4.3 受突加載荷作用的矩形板

選取帶有半圓型裂紋的矩形板，下表面固定，上表面受突加均布載荷作用，見圖17，其中r=0.03 m、a=3 m、b=3 m、c=1 m。矩形板的彈性模量E=210 GPa、泊松比ν=0.15，突加載荷P=100 MPa。

利用Abaqus計算矩形板位移，有限元網格采用20節點二次六面體單元C3D20R，取21×21×11，31×31×21，41×41×31，61×61×41和81×81×61等5種節點分布，圖17中矩形板Q點（x1=1.5 m、x3=1.0 m）位移的方差見圖18。隨著節點數增加，方差趨于0，說明Abaqus的數值解收斂。與IEFG法數值解對比時選擇節點分布為81×81×61，Abaqus中矩形板的網格劃分結果見圖19。

根據三維斷裂彈性動力學IEFG法的離散方程，利用MATLAB對矩形板位移進行數值計算，取dmax=2.8、Δt=1.0×10-6 s、α=2.5×1015，圖17中Q點的位移u1、u3和應力σ隨時間的變化見圖20～22。IEFG法與有限元法計算結果吻合，證明IEFG法的正確性和有效性。

5 結束語

基于改進的移動最小二乘法，提出三維斷裂彈性動力學的IEFG法。通過懸臂梁、柱和矩形板3個算例，對比分析節點分布、影響域比例參數、罰因子和時間步長等參數對計算精度的影響。本文提出的三維斷裂力學的IEFG法計算結果與有限元法計算結果吻合，證明該方法的正確性和有效性。

參考文獻：

[1] BELYTSCHKO T， KRYSL P， KRONGAUZ Y. A three-dimensional explicit element-free Galerkin method[J]. International Journal for Numerical Methods in Fluids， 1997， 24（12）： 1253-1270. DOI： 10.1002/（SICI）1097-0363（199706）24：12〈1253：：AID-FLD558〉3.0.CO;2-Z.

[2] CHEN Y P， ESKANDARIAN A， OSKARD M， et al. Meshless analysis of high-speed impact[J]. Theoretical and Applied Fracture Mechanics， 2005， 44（3）： 201-207. DOI： 10.1016/j.tafmec.2005.09.007.

[3] LEE S H， YOON Y C. An improved crack analysis technique by element-free Galerkin method with auxiliary supports[J]. International Journal for Numerical Methods in Engineering， 2003， 56（9）： 1291-1314. DOI： 10.1002/nme.611.

[4] RAO B N， RAHMAN S. An enriched meshless method for non-linear fracture mechanics[J]. International Journal for Numerical Methods in Engineering， 2004， 59（2）： 197-223. DOI： 10.1002/nme.868.

[5] 陳美娟， 程玉民. 改進的移動最小二乘法[J]. 力學季刊， 2003， 24（2）： 266-272. DOI： 10.3969/j.issn.0254-0053.2003.02.019.

[6] ZHANG Z， WANG J F， CHENG Y M， et al. Improved element-free Galerkin method for three-dimensional transient heat conduction problems[J]. Science China Physics， Mechanics and Astronomy， 2013， 56： 1568-1580. DOI： 10.1007/s11433-013-5135-0.

[7] CHENG Y M， PENG M J. Boundary element-free method for elastodynamics[J]. Science in China Series G： Physics and Astronomy， 2005， 48： 641-657. DOI： 10.1360/142004-25.

[8] SUN Y Z， LIEW K M. Analyzing interaction between coplanar square cracks using an efficient boundary element-free method[J]. International Journal for Numerical Methods in Engineering， 2012， 91（11）： 1184-1198. DOI： 10.1002/nme.4309.

[9] PENG M J， LI R X， CHENG Y M. Analyzing three-dimensional viscoelasticity problems via improved element-free Galerkin（IEFG） method[J]. Engineering Analysis with Boundary Elements， 2014， 40： 104-113. DOI： 10.1016/j.enganabound.2013.11.018.

[10] ZHANG Z， LIEW K M， CHENG Y M， et al. Analyzing 2D fracture problems with improved element-free Galerkin method[J]. Engineering Analysis with Boundary Elements， 2008， 32（3）： 241-250. DOI： 10.1016/j.enganabound.2007.08.012.

[11] 魏慶節， 彭妙娟. 三維彈性動力學的改進的無單元Galerkin方法[J]. 應用力學學報， 2019， 36（1）： 22-27.

[12] 蔡小杰， 彭妙娟， 程玉民. 彈塑性大變形問題的改進的無單元Galerkin方法[J]. 中國科學， 2018， 48（2）： 53-62. DOI： 10.1360/SSPMA2017-00231.

[13] 鄒詩瑩， 席偉成， 彭妙娟， 等. 運用改進的無單元Galerkin方法分析機場道面斷裂力學問題[J]. 物理學報， 2017， 66（12）： 36-44. DOI： 10.7498/aps.66.120204.

[14] 劉學聰， 章青， 夏曉舟. 一種新型裂尖加強函數的顯式動態擴展有限元法[J]. 工程力學， 2017， 34（10）： 10-18. DOI： 10.6052/j.issn.1000-4750.2016.05.0399.

[15] 石路楊， 余天堂. 多裂紋擴展的擴展有限元法分析[J]. 巖土力學， 2014， 35（1）： 263-272. DOI： 10.16285/j.rsm.2014.01.001.

[16] 楊永濤， 徐棟棟， 鄭宏. 動載下裂紋應力強度因子計算的數值流形元法[J]. 力學學報， 2014， 46（5）： 730-738. DOI： 10.6052/0459-1879-14-024.

[17] 張洪武， 盧夢凱， 鄭勇剛. 非均質飽和多孔介質彈塑性動力分析的廣義耦合擴展多尺度有限元法[J].計算力學學報， 2016， 33（4）： 454-461. DOI： 10.7511/jslx201604005.

（編輯 武曉英）