四節(jié)點矩形彈簧元及其特性研究

張青波，李世海，馮 春

（中國科學(xué)院力學(xué)研究所，北京 100190）

1 引 言

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

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

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

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

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

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

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

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

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

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

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

式中：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構(gòu)成，其中彈簧 s1、s3、s5的法向沿X軸正向，s2、s4、s6的法向沿Y軸正向。圖中A、B、C、D是4個插值點，其坐標(biāo)值和位移值均由4個節(jié)點的坐標(biāo)值和節(jié)點位移值進(jìn)行線性插值得到，4個插值點位于各邊的中點處。各基本彈簧的編號及彈簧首端編號和末端編號的對應(yīng)關(guān)系見表1。

表1 基本彈簧編號及首末端對應(yīng)關(guān)系Table 1 Relationships between the numbers ofbasic spring and its nodes

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

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

k =5時

k =6時

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

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

3 彈簧剛度系數(shù)的確定

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

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

令

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

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

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

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

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

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

式中：i=1，2。

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

式中：i =x，y。

5 算 例

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

5.1 重力作用下塊體位移場

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

圖3 計算模型Fig.3 Calculation model

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

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

5.2 側(cè)向線性壓力下位移場

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

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

圖5 計算模型Fig.5 Calculation model

5.3 懸臂梁算例

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

圖7 計算模型Fig.7 Calculation model

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

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

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

6 結(jié) 論

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

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

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

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

[1]方建瑞, 許志雄, 莊曉營. 三維邊坡穩(wěn)定彈塑性有限元分析與評價[J]. 巖土力學(xué), 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)定數(shù)值分析[J]. 巖石力學(xué)與工程學(xué)報, 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]周健, 崔積弘, 賈敏才, 等. 靜力觸探試驗的離散元數(shù)值模擬研究[J]. 巖土工程學(xué)報, 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]王桂萱, 秦建敏. 顆粒材料力學(xué)性質(zhì)的離散元數(shù)值模擬[J]. 大連大學(xué)學(xué)報, 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]劉曉宇, 梁乃剛, 李敏. 三維鏈網(wǎng)模型及其參數(shù)標(biāo)定[J].中國科學(xué)(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ù)介質(zhì)模型的離散元方法中彈簧性質(zhì)研究[J]. 巖石力學(xué)與工程學(xué)報,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]田振農(nóng), 李世海, 劉曉宇, 等. 三維塊體離散元可變形計算方法研究[J]. 巖石力學(xué)與工程學(xué)報, 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].巖土工程學(xué)報, 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), 張沖, 等. 考慮彈塑性本構(gòu)的三維模態(tài)變形體離散元方法斷裂模擬[J]. 工程力學(xué), 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]方明霽, 李國強(qiáng), 孫飛飛, 等. 基于多彈簧模型的空間梁柱單元(I): 理論模型[J]. 計算力學(xué)學(xué)報, 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]李國強(qiáng), 方明霽, 孫飛飛, 等. 基于多彈簧模型的空間梁柱單元(II): 參數(shù)分析[J]. 計算力學(xué)學(xué)報, 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ù)介質(zhì)力學(xué)的塊體單元離散彈簧法研究[J]. 巖石力學(xué)與工程學(xué)報, 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]王勖成, 邵敏. 有限單元法基本原理和數(shù)值方法[M].北京: 清華大學(xué)出版社, 1996.

[17]任鈞國, 熊龍飛. 低階高精度單元的發(fā)展[J]. 國防科技大學(xué)學(xué)報, 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.