張青波，李世海，馮 春

（中國科學院力學研究所，北京 100190）

1 引 言

傳統(tǒng)的數值計算方法將巖土體看作連續(xù)介質或者非連續(xù)介質，而忽略了巖土體宏觀連續(xù)微觀不連續(xù)的特點。有限元及有限差分方法的基礎是連續(xù)介質力學，適用于模擬巖土體的連續(xù)變形及塑性破壞[1－2]。顆粒離散元法可方便地實現(xiàn)材料內部破壞，適合模擬微觀尺度下材料的力學行為，但模擬連續(xù)介質時需對顆粒之間的彈簧進行復雜的標定[3－4]。劉曉宇等[5－6]構造了三維鏈網模型，給出了模型的幾何與物理參數標定公式。

發(fā)展新型離散元模擬巖土體漸進破壞過程的研究取得了很大的進展。李世海等[7－9]提出的基于連續(xù)介質力學的離散單元法(CDEM)和張沖等[10－11]提出的三維簡單變形體離散元方法(3SDEM)都是可以描述介質由連續(xù)到非連續(xù)的漸進破壞過程的數值計算方法。兩種方法都可將塊體看作簡單變形體，在接觸面上設置法向和切向彈簧，通過彈簧判斷實現(xiàn)非連續(xù)計算。3SDEM 中接觸彈簧剛度系數的選取具有人為性，CDEM法給出了接觸彈簧剛度的計算公式。方明霽等[12－13]將塊體離散為法向、剪切和扭轉彈簧，構造了一種多彈簧梁柱單元模型。馮春等[14]將長方體單元離散為 12根彈簧，通過離散彈簧系統(tǒng)和傳統(tǒng)有限元單元節(jié)點力的等效關系給出了彈簧剛度的表達式。Li Shihai等[15]構造了三角形和四面體彈簧元，對連續(xù)問題得到了與傳統(tǒng)有限元等價的一組彈簧系統(tǒng)，該方法通過有限元的單元剛度矩陣標定各彈簧的剛度。

新型離散元方法在計算連續(xù)介質問題時的計算能力多等同于傳統(tǒng)有限元中的常應變單元或者雙線性單元，對提高單元計算精度的研究較少。但對某些常見問題（如梁的彎曲）這些單元的計算精度是很差的[16]。在有限元中（如Wilson非協(xié)調元等）構造各種形式的非協(xié)調元就是一種構造低階高精度單元的有效手段[17]。

Li Shihai等[15]給出了一種對連續(xù)問題完全等價于相應有限元的彈簧系統(tǒng)，通過對比單元剛度矩陣得到各彈簧的剛度，其彈簧剛度的計算公式較其他幾種方法更加合理。本文在此基礎上發(fā)展了一種四節(jié)點矩形彈簧元，并求得了彈簧剛度的表達形式。通過改變彈簧剛度表達式中的系數可使該單元精度分別等同于常應變、雙線性和Wilson非協(xié)調單元，極大地提高了彈簧元法求解連續(xù)問題時的計算精度。本文主要包含5個部分，下面分別進行闡述。

2 四節(jié)點矩形彈簧元的構成

考慮圖1所示的四節(jié)點矩形單元，節(jié)點編號為1、2、3、4，節(jié)點坐標為單元長為 a，寬為 b。用 u表示X方向的位移，v表示Y方向的位移。

根據彈性力學理論[16]，四節(jié)點矩形單元的應變能表達式可寫為

圖1 四節(jié)點矩形單元Fig.1 Four-node rectangular element

式中：εx、εy、σx、σy為單元正應變和正應力；γxy、τxy為單元剪應變和剪應力；E為彈性模量；μ為泊松比；t為單元厚度。

彈簧元法是一種新型離散元，用離散系統(tǒng)描述連續(xù)介質力學行為是該方法的重要內容之一。與其他新型離散元方法不同，其核心思想是將單元看作一種結構，將單元應變能離散到一系列的彈簧上，用彈簧系統(tǒng)的剛度代替有限元的單元剛度描述單元節(jié)點力與節(jié)點位移之間的關系，使之在計算連續(xù)介質問題時與傳統(tǒng)有限元完全等價。同時可以利用動態(tài)松弛方法等進行求解，避免傳統(tǒng)有限元必須形成總體剛度矩陣帶來的問題。

按照彈簧元的基本思想，將式（1）的應變能表達式寫為求和形式，用彈簧的彈性勢能替代單元的應變能，可得如下表達式：

式中：G為剪切模量；ai、bi為彈簧長度；ui、vi為彈簧彈性變形。

對圖1所示的四節(jié)點矩形單元，將其離散為圖2所示的四節(jié)點矩形彈簧元。

圖2 四節(jié)點矩形彈簧元Fig.2 Four-node rectangular spring element

如圖2所示，四節(jié)點矩形彈簧元由6個基本彈簧 s1、s2、s3、s4、s5、s6構成，其中彈簧 s1、s3、s5的法向沿X軸正向，s2、s4、s6的法向沿Y軸正向。圖中A、B、C、D是4個插值點，其坐標值和位移值均由4個節(jié)點的坐標值和節(jié)點位移值進行線性插值得到，4個插值點位于各邊的中點處。各基本彈簧的編號及彈簧首端編號和末端編號的對應關系見表1。

表1 基本彈簧編號及首末端對應關系Table 1 Relationships between the numbers ofbasic spring and its nodes

若設4個節(jié)點的位移分別為(ui，vi)(i=1, 2, 3, 4)，則6個基本彈簧的變形可由各彈簧兩端點的相對位移求得，即第k個基本彈簧的兩個方向的變形可表示為

式中：上標k =1, 2, 3, 4為彈簧編號；下標i為彈簧首端編號，下標j為彈簧末端編號。

k =5時

k =6時

則矩形單元的應變能可由圖2所示的6個基本彈簧的彈性勢能表示為

式（6）中含 Ki的項等價于式（1）中含的項，含Gi的項等價于式(1)中含的項，故稱Ki、Gi分別為彈簧si的法向和切向彈簧剛度系數。在式（6）中含KK的項等價于式（1）中含的項，含GG的項等價于式（1）中含的項，故稱KK、GG分別為四節(jié)點矩形彈簧元的泊松彈簧和純剪彈簧剛度系數。

3 彈簧剛度系數的確定

按照能量變分原理，式（6）表示的彈簧元應變能對節(jié)點位移求二次偏導數可得到彈簧元對應的單元剛度矩陣，其中

將有限元的形函數代入式（1）中并對節(jié)點位移求二次偏導數可得到有限元的單元剛度矩陣，其中

令

可以得到與四節(jié)點矩形有限單元對應的四節(jié)點矩形彈簧元中各彈簧的剛度系數表達式。式（10）給出了含待定系數的彈簧元剛度系數表達式。

式中：α、β為待定系數，其物理含義為單元應變能在各彈簧上的分配比例。

對任何一種確定的有限元形函數，都可通過式（9）得到相應的彈簧元的剛度系數表達形式。在有限元中可以通過給定含高階項的形函數得到高精度的有限元解，同理，在彈簧元中也可以通過給定對應于含高階項的形函數的待定系數，從而得到高精度的彈簧元解。對于四節(jié)點矩形單元，表2給出了有限元形函數為常應變非協(xié)調單元、雙線性單元和Wilson非協(xié)調單元時，對應的四節(jié)點矩形彈簧元的剛度系數表達式中待定系數的取值。

表2 不同單元彈簧元剛度系數的取值Table 2 Coefficients in the stiffness expressions of springs with different elements

4 彈簧力、節(jié)點力的確定

由6個基本彈簧的變形及各彈簧的剛度系數可求得基本彈簧兩個方向的彈簧力，其計算公式為

式中：i=1，2。

4個節(jié)點的節(jié)點力的計算公式為

式中：i =x，y。

5 算 例

將四節(jié)點矩形彈簧元與基于連續(xù)介質力學的離散單元法(CDEM)結合，形成了含四節(jié)點矩形彈簧元的計算程序。應用該程序計算了簡單模型在簡單載荷作用下的位移場，通過給定不同剛度系數，比較不同精度的單元對不同問題的適用性。

5.1 重力作用下塊體位移場

重力作用是邊坡穩(wěn)定性分析的主要載荷?？紤]一個100 m×100 m的塊體在重力作用下的位移，計算模型如圖3所示，通過給定不同的待定系數，用不同精度的單元計算A點Y方向的位移。材料的彈性模量 E =1 500 MPa，泊松比μ=0.25，密度ρ=2 g/cm3，計算結果如圖4所示。

圖3 計算模型Fig.3 Calculation model

圖4 不同單元計算的A點位移與網格數(n×n)的關系Fig.4 Relationships between Y-displacement at point A and grid numbers (n×n)with different elements

從圖中可以看出，對于計算重力載荷問題，雙線性單元和 Willson非協(xié)調單元的計算精度相差不大，但常應變單元的計算精度相比較差，且塊體劃分為單數個網格時的結果比劃分為偶數個網格時的結果好。

5.2 側向線性壓力下位移場

側向線性壓力載荷是土石壩、擋土墻穩(wěn)定性分析中須考慮的重要載荷形式。如圖 5所示一個100 m×100 m的塊體，底部固定，左側承受靜水壓力。通過給定不同的待定系數，用不同精度的單元計算 A點 X方向的位移，材料的彈性模量 E =1 500 MPa，泊松比μ=0.25，密度ρ=2 g/cm3，計算結果如圖6所示。

從圖中可以看出，對于計算靜水壓力問題，非線性單元和 Wilson非協(xié)調單元的計算精度相差不大，使用常應變單元計算時塊體劃分為偶數個網格時的結果比劃分為奇數個網格時的結果好。

圖5 計算模型Fig.5 Calculation model

5.3 懸臂梁算例

在工程結構中經常將問題簡化為懸臂梁在自重作用下的撓度問題，考慮如圖7所示懸臂梁：平面尺寸為10 m×2 m，梁左端固定，關注A點在重力作用下的撓度，計算結果如圖 8、9所示。材料的彈性模量E=1 500 MPa，泊松比μ=0.25，密度ρ=2 g/cm3。

圖7 計算模型Fig.7 Calculation model

從圖8可以看出，懸臂梁的變形與理論解一致。從圖9可以看出，對懸臂梁問題傳統(tǒng)雙線性單元的精度與常應變單元的精度都較差，而Wilson非協(xié)調元的精度較好，這與文獻[16]中的結論一致，也從側面驗證了程序的可靠性。

圖8 Wilson非協(xié)調元計算的位移 (網格數5×1)Fig.8 Displacement calculated by Wilson incompatible elements (with grid number of 5×1)

圖9 不同單元計算的A點位移與網格數(5n×n)的關系Fig.9 Relationships between Y-displacement at point A and grid numbers(5n×n)with different elements

6 結 論

（1）將四節(jié)點矩形單元離散為6個基本彈簧，給出了一種四節(jié)點矩形彈簧元的構造形式。該單元的彈簧剛度表達形式、節(jié)點力公式簡單，易于程序化，同傳統(tǒng)有限元方法相比，可提高計算效率。

（2）通過改變彈簧剛度表達式中的待定系數可使該單元分別等價于常應變、雙線性和Wilson非協(xié)調單元，表明該單元具有良好的擴展性，當待定系數取其他值時，可得到更高精度或者更低精度的結果。與其他新型離散元法相比，在計算連續(xù)介質問題時該單元具有更高的計算精度。

（3）不同算例表明：對重力作用下一般邊坡的穩(wěn)定性及側向線性荷載作用下一般土石壩、擋土墻的穩(wěn)定性問題可采用雙線性及非協(xié)調元進行計算；對重力作用下一般懸臂梁的彎曲問題，應采用非協(xié)調元進行計算。

（4）本文研究了四節(jié)點矩形彈簧元的特性，對四節(jié)點任意四邊形單元，可按類似的方法得到對應于雙線性單元的離散彈簧剛度表達式，但對其特性的研究還有待深入。

[1]方建瑞, 許志雄, 莊曉營. 三維邊坡穩(wěn)定彈塑性有限元分析與評價[J]. 巖土力學, 2008, 29(10): 2667－2672.FANG Jian-rui, XU Zhi-xiong, ZHUANG Xiao-ying.Appraisal and analysis of three-dimensional slope stability based on elastoplastic FEM[J]. Rock and Soil Mechanics, 2008, 29(10): 2667－2672.

[2]戴自航, 劉志偉, 劉成禹, 等. 考慮張拉與剪切破壞的土坡穩(wěn)定數值分析[J]. 巖石力學與工程學報, 2008,27(2): 375－382.DAI Zi-hang, LIU Zhi-wei, LIU Cheng-yu, et al.Numerical anslysis of soil slope stability considering tension and shear failures[J]. Chinese Journal of Rock Mechanics and Engineering, 2008, 27(2): 375－382.

[3]周健, 崔積弘, 賈敏才, 等. 靜力觸探試驗的離散元數值模擬研究[J]. 巖土工程學報, 2007, 29(11): 1604－1610.ZHOU Jian, CUI Ji-hong, JIA Min-cai, et al. Numerical simulation of cone penetration test by discrete element method[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(11): 1604－1610.

[4]王桂萱, 秦建敏. 顆粒材料力學性質的離散元數值模擬[J]. 大連大學學報, 2007, 28(6): 27－31.WANG Gui-xuan, QIN Jian-min. Numerical simulation of mechanical behaviors of ellipsoid granular materials by discrete element method[J]. Journal of Dalian University, 2007, 28(6): 27－31.

[5]劉曉宇, 梁乃剛, 李敏. 三維鏈網模型及其參數標定[J].中國科學(A輯), 2002, 32(10): 887－894.LIU Xiao-yu, LIANG Nai-gang, LI Min. 3D network model and its parameter calibration[J]. Science in China(Ser. A), 2002, 32(10): 887－894.

[6]GUSEV A A. Finite element mapping for spring network representations of the mechanics of solids[J]. Physical Review Letter, 2004, 93(3): 1－4.

[7]LI Shi-hai, ZHAO Man-hong, WANG Yuan-nian, et al. A continuum-based discrete element method for continuous deformation and failure process[C]//WCCM VI in Conjunction with APCOM’04. Beijing: [s. n.], 2004.

[8]魏懷鵬, 易大可, 李世海, 等. 基于連續(xù)介質模型的離散元方法中彈簧性質研究[J]. 巖石力學與工程學報,2006, 25(6): 1159－1169.WEI Huai-peng, YI Da-ke, LI Shi-hai, et al. Study of spring properties of continuum-based discrete element method[J]. Chinese Journal of Rock Mechanics and Engineering, 2006, 25(6): 1159－1169.

[9]田振農, 李世海, 劉曉宇, 等. 三維塊體離散元可變形計算方法研究[J]. 巖石力學與工程學報, 2008, 27(增刊 1): 2832－2840.TIAN Zhen-nong, LI Shi-hai, LIU Xiao-yu, et al.Research on deformable calculation method based on three-dimensional block discrete element[J]. Chinese Journal of Rock Mechanics and Engineering, 2008,27(Supp. 1): 2832－2840.

[10]張沖, 金峰, 侯艷麗. 三維簡單變形體離散元方法[J].巖土工程學報, 2007, 29(2): 159－163.ZHANG Chong, JIN Feng, HOU Yan-li. 3-D simple deformable distinct element method[J]. Chinese Journal of Geotechnical Engineering, 2007, 29(2): 159－163.

[11]金峰, 胡衛(wèi), 張沖, 等. 考慮彈塑性本構的三維模態(tài)變形體離散元方法斷裂模擬[J]. 工程力學, 2011, 28(5): 1－7.JIN Feng, HU Wei, ZHANG Chong, et al. A fracture simulation using 3-D mode distinct element method(3MDEM)with elastoplastic constitutive model[J].Engineering Mechanics, 2011, 28(5): 1－7.

[12]方明霽, 李國強, 孫飛飛, 等. 基于多彈簧模型的空間梁柱單元(I): 理論模型[J]. 計算力學學報, 2008, 25(1):129－133.FANG Ming-ji, LI Guo-qiang, SUN Fei-fei, et al.Multi-spring model beam element(I): Theoretical model[J]. Chinese Journal of Computational Mechanics, 2008, 25(1): 129－133.

[13]李國強, 方明霽, 孫飛飛, 等. 基于多彈簧模型的空間梁柱單元(II): 參數分析[J]. 計算力學學報, 2008, 25(1):134－138.LI Guo-qiang, FANG Ming-ji, SUN Fei-fei, et al.Multi-spring model beam element(II): Parameter analysisl[J]. Chinese Journal of Computational Mechanics, 2008, 25(1): 134－138.

[14]馮春, 李世海, 姚再興. 基于連續(xù)介質力學的塊體單元離散彈簧法研究[J]. 巖石力學與工程學報, 2010, 29(增刊 1): 2690－2705.FENG Chun, LI Shi-hai, YAO Zai-xing. Study of block-discrete-spring method based on continuum mechanics[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(Supp.1): 2690－2705.

[15]LI Shi-hai, ZHANG Ya-nan, FENG Chun. A spring system equivalent to continuum model[C]//Discrete Element Methods, Simulation of Discontinua：Theory and Applications. London: Queen Mary, University of London, 2010: 75－85.

[16]王勖成, 邵敏. 有限單元法基本原理和數值方法[M].北京: 清華大學出版社, 1996.

[17]任鈞國, 熊龍飛. 低階高精度單元的發(fā)展[J]. 國防科技大學學報, 1998, 20(4): 109－114.REN Jun-guo, XIONG Long-fei. Development of low order element with high accuracy[J]. Journal of National University of Defense Technology, 1998, 20(4): 109－114.