考慮蠕變的任意多邊形雜交應(yīng)力單元

中圖分類號：0344.6 文獻標志碼：A 文章編號：1000-582X（2025）08-028-12

Arbitrary polygonal hybrid stress element for creep analysis in roadbed engineering

LI Zijun'，TIAN Jingjie2，HU Changhao2，TAN Dawen2，XIA Fan2，HOU Chunyao2，LIU Luyao2， NING Shaoqing2，GUO Ran1 （1.Faculty of Civil Engineering and Mechanics，Kunming University ofScience and Technology，Kunming 650500， P. R. China; 2. Xiluodu Hydropower Plant， China Yangtze Power Co.， Ltd.， Zhaotong， Yunnan 657300，P.R.China）

Abstract：Post-construction settlement of roadbeds isa critical factor influencing the long-term safetyand performance of highways. Accurate prediction of creep behavior is of great significance for ensuring the structural integrity of roadbed systems. However，conventional finite element simulations of creep behavior require extensive meshing for complex structural models，leading to high computational costs and time consumption.To addressthese challnges，this paper proposes a novel arbitrary polygonal hybrid stress element method （PHSEM） that incorporates creep effects for roadbed settlement analysis.Based onthe hybrid stress finite element method andthetime-dependent deformation characteristics of roadbed soil，the theoretical formulationof the PHSEM is derived.This element introduces a higher-order stress field，improving computational accuracy.Anumerical creep model is further established and validated against results obtained from MARC software.The simulation results show that the PHSEM achieves good agreement with benchmark solutions， with deviations within acceptable limits.Furthermore，the element accommodates polygons with variable edge counts，offering flexible meshing for complex，heterogeneous roadbed models and enabling realistic stress distribution analysis.The proposed PHSEM provides a promising approach for creep modeling in practical geotechnical engineering applications.

Keywords： hybrid stress element;arbitrary polygonal hybrid stress element;creep; higher-order stressfield; finite element; roadbed

在公路建設(shè)過程中，路基沉降會對所建道路的質(zhì)量產(chǎn)生重要影響。由于地理環(huán)境、氣候條件和地質(zhì)成因等原因，不同地區(qū)的土壤質(zhì)地特征差異明顯。在路基施工過程中及施工后，這些土壤特征可能導(dǎo)致嚴重問題，如不均勻差異沉降邊坡坍塌和山體滑坡等，影響道路安全、運輸效率和道路使用壽命。因此，路基沉降量的研究成了公路建設(shè)中的一個重要問題[2]。在沉降分析時，必須考慮路基填土的蠕變特性[34]，因為會產(chǎn)生重大影響[5]。

路基的沉降可以劃分為3個階段，初始瞬時沉降、固結(jié)沉降和次固結(jié)沉降。次固結(jié)沉降也稱蠕變沉降，這部分沉降是指地基在外載荷作用下，土體中超孔隙水壓力已完全消散，有效應(yīng)力不變的情況下，由土的固體骨架長時間緩慢蠕變所產(chǎn)生的沉降，沉降效果是隨時間長期緩慢體現(xiàn)的。土在外荷載作用下的長期變形行為一直是地質(zhì)工程和巖土工程的研究重點之一[6-7]。1943年Terzaghi基于線彈性的假說，首先提出了飽和土的一維固結(jié)理論。研究表明，由于土壤骨架的黏性特征，土壤的應(yīng)力-應(yīng)變行為具有高度的時間依賴性[9-10]，這種時間依賴現(xiàn)象主要來源于飽和土固結(jié)和土骨架蠕變。固結(jié)是由超孔隙壓力的消散引起的，而土壤骨架的黏度會產(chǎn)生蠕變行為[]。在土壤骨架具有黏彈性的假設(shè)下，，則可以描述為持續(xù)很長時間的蠕變行為。Biot和Bjerrum在早期的研究中采用黏彈性的本構(gòu)模型（如Kelvin，Maxwell和Merchant模型）來描述土壤的固結(jié)和蠕變過程。

雜交應(yīng)力元法是Pian（卞學(xué)璜）2在1964年基于修正余能原理首次提出的，克服了傳統(tǒng)位移有限元法無法解決單元邊數(shù)大于4時位移場插值的困難。為了避免位移型多邊形單元的位移插值函數(shù)構(gòu)造困難問題，雜交元方法被引入到多邊形單元的構(gòu)建中。

1994年Ghosh等[13]提出了基于最小余能原理的 Voronoi有限元法（Voronoi cellfinite element method，VCFEM），該方法利用域內(nèi)的應(yīng)力插值，不需要構(gòu)造域內(nèi)位移插值函數(shù)，使單元構(gòu)造變得簡單。1995年Zhang等[14-15提出了基于Hellinger-Reissner兩場變分原理的雜交多邊形有限元法，用于含夾雜和孔洞材料的分析。2010年Fu等[基于變分原理開發(fā)了新型8節(jié)點平面單元，在極端網(wǎng)格畸變下也能保持精度不變。2011年Cen等首次提出形狀自由（shape-free）有限元法，用以描述有限元性能不受單元形狀影響的特性。同年Cen等8針對各向同性問題成功構(gòu)造了2個平面8結(jié)點單元和2個平面12結(jié)點單元，針對各向同性和各向異性問題，構(gòu)造了平面4結(jié)點帶旋轉(zhuǎn)自由度膜元[19]。2015年Zhou等[20]成功構(gòu)造了一種新的二次多邊形雜交應(yīng)力函數(shù)單元，構(gòu)造方法簡單且精度高于許多其他方法。2019年楊鋒等2提出了一種基于多邊形雜交應(yīng)力單元的優(yōu)化接觸算法（polygonal hybrid stress element method，PHSEM）。由于單元邊數(shù)可變，PHSEM能夠靈活適應(yīng)復(fù)雜的模型邊界與材料邊界，更容易劃分網(wǎng)格。

筆者基于PHSEM提出了一種考慮蠕變的任意多邊形雜交應(yīng)力新單元，用于分析不同巖土種類的路基產(chǎn)生的蠕變變形行為。同時，討論了任意多邊形雜交應(yīng)力單元的高效性，以及此單元對宏觀復(fù)雜的復(fù)合介質(zhì)材料研究的可行性。

1單元基本公式推導(dǎo)

多邊形單元的形狀具有任意性，在處理復(fù)雜結(jié)構(gòu)時具有較大優(yōu)勢。雜交應(yīng)力法不需要構(gòu)造域內(nèi)的位移插值函數(shù)，利用這種方法構(gòu)造的多邊形單元不僅精度高，而且具有很強的抗畸變敏感性能[20]。利用雜交應(yīng)力方法的優(yōu)勢，將其應(yīng)用到任意多邊形單元的構(gòu)造中，以解決目前多邊形單元構(gòu)造困難等問題。

任意多邊形雜交應(yīng)力單元的示意圖如圖1所示，網(wǎng)格劃分由邊數(shù)不定的多邊形構(gòu)成，可根據(jù)材料的空間分布及結(jié)構(gòu)的特點，高效便捷地劃分網(wǎng)格，可有效模擬復(fù)雜不規(guī)則的模型。典型任意多邊形雜交應(yīng)力單元的構(gòu)成如圖2所示。

圖1任意多邊形雜交應(yīng)力單元的示意圖 Fig.1Schematic diagram of arbitrary polygonal hybridstresselements

圖2任意多邊形雜交應(yīng)力單元

Fig.2Themodelofarbitrarypolygonalhybridstresselement

1.1 提出泛函

公式推導(dǎo)的理論框架基于Tong和Pian[2-24以及Atluri[25]建立的雜交應(yīng)力元增量模型。

式中： ΔB 為余能增量； n 為單元外邊界上外法線向量;增量步開始前已知的初始場變量： σ 為單元域 內(nèi)滿足平衡的應(yīng)力場； ε 為單元域內(nèi)滿足協(xié)調(diào)關(guān)系的應(yīng)變場； 為單元外邊界 的位移 和 分別為力邊界 ?Ω_t 給定外力和給定外力增量。未知的場變量： Δσ 為單元域內(nèi)的應(yīng)力增量； Δu 為單元外邊界上的位移增量。

對于式（1），滿足給定面力邊界條件式（2）和單元間的面力連續(xù)條件式（3）。

n^-?（σ_e+Δσ_e）+n⁺?（Δdownσ_e+Δσ_e）=0

在該模型的基礎(chǔ)上，考慮滿足單元域內(nèi)平衡關(guān)系的蠕變應(yīng)變，得到了修正后的增量形式的余能泛函。

式中， Δε^c 為單元域內(nèi)的蠕變應(yīng)變增量。

余能增量 ΔB 滿足

式中： Δε^e 為單元域內(nèi)彈性應(yīng)變增量； Δε^p 為單元域內(nèi)塑性應(yīng)變增量。

1.2 單元的應(yīng)力場/邊界位移假設(shè)

在無外力作用的靜力條件下，不考慮體力時應(yīng)力場滿足平衡方程

abla?σ=0

對于二維問題，引人Airy應(yīng)力函數(shù)，得到滿足平衡方程的應(yīng)力場以及應(yīng)力增量的矩陣形式為

σ_e=Pβ，

Δσ_e=PΔβ，

式中： β 為單元內(nèi)部假定獨立應(yīng)力場的應(yīng)力系數(shù)矩陣； P 為應(yīng)力插值矩陣。

單元的邊界位移增量可以表示為各個邊界段節(jié)點的廣義位移增量的插值

Δu=LΔq，

式中： L 為單元邊界上的形函數(shù)矩陣； Δq 為單元邊界節(jié)點位移。

將式 （7）_～（9） 帶入式（4）得

式中， s 為彈塑性柔度張量。其中，

1.3 單剛推導(dǎo)

根據(jù)修正的余能駐值條件

可以得到每個單元內(nèi)運動關(guān)系的弱表達式為

在上一增量步中滿足

則式（16）可以表示為

HΔβ+G^c-GΔq=0

在迭代求解過程中，第 i 迭代步的應(yīng)力參數(shù)增量表示為

在迭代求解過程中，第 i 迭代步的位移增量表示為

由式（18）～（20）可得到應(yīng)力參數(shù)迭代增量的位移迭代增量表達式。第1迭代步：

dβ¹=H^-1Gdq¹-H^-1G^c，

第2～i迭代步：

根據(jù)修正余能駐值條件

可以得到相應(yīng)的面力邊界條件的弱表達式為

將式（19）帶入式（24）可得

將式（21）帶入式（25）可整理得第1迭代步的求解方程

將式（22）帶入式（25）可整理得第 2～i 迭代步的求解方程

其中，剛度矩陣為

K=G^TH^-1G

2 土壤蠕變方程

式（4）中的蠕變應(yīng)變增量，采用以下的土壤蠕變本構(gòu)方程進行計算。對于一維蠕變問題，等效蠕變應(yīng)變率與任何相關(guān)參數(shù)的相關(guān)性可以直接以冪定律形式表示，采用合適的指數(shù)，參考MSCMARC手冊。

式中：A為蠕變材料常數(shù)； 為等效蠕變應(yīng)變率； σ，ε^c，T， t分別為等效應(yīng)力、等效蠕變應(yīng)變、溫度和時間。

m、n、 p 、 q 分別為與等效蠕變應(yīng)變、等效應(yīng)力、溫度以及時間相關(guān)的指數(shù)。

由式（29）可知，土壤蠕變過程是一個非常復(fù)雜的非線性力學(xué)過程，其應(yīng)力與應(yīng)變都隨時間而變化，是時間的函數(shù)。由于蠕變屬于不可逆變形且遵循塑性理論發(fā)展，對于三維蠕變問題，通常將塑性理論推廣到蠕變情況，可以建立如下蠕變應(yīng)變增量方程。

式中： E_ij^c 為蠕變應(yīng)變張量； S_ij 為應(yīng)力偏張量； 為等效應(yīng)力； 為等效蠕變應(yīng)變。

3 數(shù)值算例

建立了2個數(shù)值算例，第1個算例是簡單的泥板試件模型，運用一個任意多邊形雜交應(yīng)力單元模擬泥板試件，主要為了驗證考慮蠕變的PHSEM有限元模型的有效性。第2個算例是包含2種材料的路基模型，用于分析多種材料路基的蠕變沉降問題，以及此單元對復(fù)雜材料研究和在路基工程應(yīng)用中的可行性。

3.1 數(shù)值算例一

模型結(jié)構(gòu)如圖3（a）所示，是一塊長 2m ，寬 2m 的正方形泥板。泥板左下角和右下角采用完全固定的約束，對泥板施加豎直向下 1kN 的載荷。PHSEM模型采用4個節(jié)點，1個任意多邊形雜交應(yīng)力單元，250項的應(yīng)力函數(shù)多項式（21階的Airy應(yīng)力函數(shù)），網(wǎng)格劃分如圖3（b）所示;MARC模型采用10201個節(jié)點，10000個平面應(yīng)變?nèi)e分四邊形單元，網(wǎng)格劃分如圖3（c）所示。

模型所用材料為凍土[2]，，其材料力學(xué)參數(shù)如表1所示。計算模型考慮在施加載荷下產(chǎn)生的蠕變，采用的蠕變參數(shù)如表2所示[2]，模擬泥板試件 ^100h 的蠕變過程。

Fig.3 Model structuralconditionsand meshingdiagrams

表1凍土力學(xué)參數(shù)表

Table1 Mechanical parametersoffrozen soil

表2凍土蠕變參數(shù)表

圖3模型結(jié)構(gòu)工況圖和網(wǎng)格劃分圖

Table2 Creepparametersoffrozen soil

注：蠕變參數(shù)擬合時， σ 為等效應(yīng)力， Pa;T 為溫度， ^°C;t 為時間，h。

為了驗證任意多邊形雜交應(yīng)力單元的有效性，建立結(jié)構(gòu)相同和邊界條件一致的PHSEM模型和MARC模型。圖4顯示了PHSEM模型和MARC模型考慮蠕變時的水平方向應(yīng)力云圖 σ_x ，豎直方向應(yīng)力云圖 σ_y 和等效馮米塞斯應(yīng)力云圖 。圖5顯示了PHSEM模型和MARC模型考慮蠕變時的水平方向蠕變應(yīng)變云圖 E_x ，豎直方向蠕變應(yīng)變云圖 ε_y 和等效蠕變應(yīng)變云圖ε。由應(yīng)力云圖和蠕變應(yīng)變云圖可以看出，PHSEM模型的應(yīng)力和蠕變應(yīng)變分布與MARC模型得到的結(jié)果吻合，即計算結(jié)果與理論吻合，驗證了PHSEM模型考慮蠕變的有效性。

通過觀察路徑曲線上應(yīng)力值和蠕變應(yīng)變值的擬合情況可以看出，PHSEM模型的變化趨勢與MARC模型的結(jié)果一致，這證明了任意多邊形雜交應(yīng)力單元考慮蠕變的準確性。

圖6擬合了PHSEM模型和MARC模型在路徑曲線上的應(yīng)力值。可以注意到2個模型在水平方向應(yīng)力值、豎直方向應(yīng)力值和等效馮米塞斯應(yīng)力值的擬合趨勢都是高度一致的。由此可見，PHSEM的計算結(jié)果和MARC的計算結(jié)果較吻合，由PHSEM計算的整體應(yīng)力場分布是正確的。

圖7擬合了PHSEM模型和MARC模型在路徑曲線上的蠕變應(yīng)變值?？梢宰⒁獾?個模型在水平方向蠕變應(yīng)變值、豎直方向蠕變應(yīng)變值和等效蠕變應(yīng)變值的擬合趨勢基本都是一致的，誤差也在合理范圍內(nèi)。由此可見，PHSEM的計算結(jié)果和MARC的計算結(jié)果比較吻合，由PHSEM計算的整體蠕變應(yīng)變場分布是符合客觀規(guī)律的。

圖6應(yīng)力路徑圖

Fig.6Diagram of stress on thepath

圖7蠕變路徑圖

Fig.7Diagram ofcreep strain on the path

3.2 數(shù)值算例二

模型結(jié)構(gòu)如圖8所示，地基為長 70m 寬 20m 的長方形，路堤為上底 18m 、下底 26m 高 6m 的梯形。為了方便對比分析，PHSEM模型和MARC模型所受邊界條件和載荷完全相同，地基底部采用完全固定的約束，限制地基頂部的水平位移，對路堤施加豎直向下 400kN 的載荷。PHSEM模型采用7個節(jié)點，2個任意多邊形雜交應(yīng)力單元，250項的應(yīng)力函數(shù)多項式（21階的Airy應(yīng)力函數(shù)），網(wǎng)格劃分如圖9（a）所示;MARC模型采用1653個節(jié)點，1556個平面應(yīng)變?nèi)e分四邊形單元，網(wǎng)格劃分如圖9（b）所示。

圖8模型結(jié)構(gòu)工況圖

Fig.8Model structuralconditions

圖9模型網(wǎng)格劃分圖

Fig.9Model meshingdiagrams

模型所用材料為雜填土和暗綠草黃色粉質(zhì)黏土[27]，其力學(xué)參數(shù)如表3所示。計算模型考慮在施加載荷下產(chǎn)生的蠕變，采用的蠕變參數(shù)如表4所示[27]，模擬路基 1000h 的蠕變過程。

表3各土層力學(xué)參數(shù)表

Table3Mechanical parametersofeachsoillayer

表4各土層蠕變參數(shù)表

Table4Creepparametersofeachsoil layer

注：蠕變參數(shù)擬合時， σ 為等效應(yīng)力， Pa;T 為溫度， ^°C;t 為時間， ?h

同樣為了驗證任意多邊形雜交應(yīng)力單元的有效性，建立結(jié)構(gòu)相同和邊界條件一致的PHSEM模型和MARC模型。圖10顯示了PHSEM模型和MARC模型考慮蠕變時的水平方向應(yīng)力云圖 σ_x ，豎直方向應(yīng)力云圖 σ_y 和等效馮米塞斯應(yīng)力云圖 。圖11顯示了PHSEM模型和MARC模型考慮蠕變時的水平方向蠕變應(yīng)變云圖 ε_x ，豎直方向蠕變應(yīng)變云圖 ε_y 和等效蠕變應(yīng)變云圖ε。由應(yīng)力云圖和蠕變應(yīng)變云圖可以看出，PHSEM模型的應(yīng)力和蠕變應(yīng)變分布與MARC模型得到的結(jié)果吻合，即計算結(jié)果與理論吻合，驗證了PHSEM模型考慮蠕變的有效性。

圖10 應(yīng)力云圖

Fig. 10 Contoursof stress

圖11 蠕變云圖

Fig.11 Contoursof creep strain

從圖10可以注意到，2個模型的水平應(yīng)力分布圖是相似的，并且在路基的棱角和邊界處可觀察到明顯的應(yīng)力集中情況，說明PHSEM采用的高階應(yīng)力場能準確捕捉到應(yīng)力集中現(xiàn)象，雖然MARC也得到了明顯的應(yīng)力集中情況，但是在地基和路堤接觸的地方應(yīng)力等值線不夠光滑，應(yīng)力界限不規(guī)則。任意多邊形雜交應(yīng)力單元獲得高階應(yīng)力場以及應(yīng)力集中現(xiàn)象的優(yōu)勢由此可以看出。2個模型的豎直應(yīng)力分布圖和等效馮米塞斯應(yīng)力分布圖都是基本一致的，同時可以看出材料分層對應(yīng)力的分布趨勢有一定影響，PHSEM模型和MARC模型在地基和路堤接觸的地方應(yīng)力分布比較吻合。圖11顯示了蠕變應(yīng)變云圖，可以注意到2個模型的蠕變應(yīng)變云圖基本一致。PHSEM模型得到的整體應(yīng)力分布和蠕變應(yīng)變分布明確并符合客觀規(guī)律。

通過上述2個算例的驗證可以看出，任意多邊形雜交應(yīng)力單元（PHSEM）采用250項的應(yīng)力函數(shù)完備多項式（21階的Airy應(yīng)力函數(shù)），使用較少單元就可以模擬出較為準確的結(jié)果;普通位移有限元采用的為低階單元，若使用少量單元則捕捉不住應(yīng)力集中的現(xiàn)象，要想達到同樣的計算精度，需要劃分細密的網(wǎng)格來實現(xiàn)，從網(wǎng)格劃分圖中可以明顯看出在MARC模型中的單元數(shù)量要遠遠多余PHSEM模型中的單元數(shù)量。

4結(jié)論

在PHSEM的基礎(chǔ)上，提出一種考慮蠕變的任意多邊形雜交應(yīng)力單元，通過泥板試件模型和2種材料的路基模型對該單元進行了驗證，將計算結(jié)果與商用軟件MARC的計算結(jié)果進行比較，驗證了任意多邊形雜交應(yīng)力單元的有效性和高效性。得出以下結(jié)論：

1）PHSEM采用的是高階應(yīng)力場，單元數(shù)量較少，計算精度更高。在大規(guī)模工程計算中的優(yōu)勢將是顯而易見的。

2）PHSEM單元的邊數(shù)是任意的，模擬復(fù)雜結(jié)構(gòu)和多種材料時只需要幾個甚至一個單元，普通的位移有限元在復(fù)雜結(jié)構(gòu)處需劃分細密的網(wǎng)格。

3）考慮蠕變的任意多邊形雜交應(yīng)力單元在分析路基蠕變問題時，能夠準確反應(yīng)應(yīng)力場和應(yīng)變場，在邊界和材料界面能夠捕捉到應(yīng)力集中的現(xiàn)象。

任意多邊形雜交應(yīng)力單元法在處理復(fù)雜結(jié)構(gòu)的模型時更方便快捷，彌補了目前傳統(tǒng)有限元方法在路基沉降工程應(yīng)用中的一些不足，為實際工程問題的數(shù)值模擬提供了一種可參考的便利高效新單元。

參考文獻

[1]Lin QQ，ZhangYH，YangCJ，etal.Evaluationofdiferential setlementof subgrade forhighway-idening projects[J]. Sustainability，2023，15（4）：2950.

[2]Liu D W，ZhangWM，Tang Y，etal.Evolving support vectorregresson basedonimproved grey wolfoptimization for predicting setlementduring constructionofhigh-filled roadbed[J].Transportation Geotechnics，2024，45：101233.

[3]Yao Y P，F(xiàn)ang Y F.Negative creep of soils[J]. Canadian Geotechnical Journal，2020，57（1）： 1-16.

[4]Wang SM，Zhan QH，Wang L，et al.Unsaturatedcrepbehaviors and creep modelofslip-surfacesoilofalandslide nTree Gorges Reservoir Area，China[J].Buletin ofEngineering Geologyand the Environment，2021，80（7）： 5423-5435.

[5］冷伍明，鄧志龍，徐方，等.基于路基土蠕變效應(yīng)的路基預(yù)應(yīng)力損失模型研究[J].巖土力學(xué)，2022，43（6）：1671-1682. Leng WM，DengZL，XuF，etal.Aprestresslossmodelforsubgradecosidering crepeectofsubgradesoil[J]Rockand Soil Mechanics，2022，43（6）： 1671-1682. （in Chinese）

[6]Biot MA.Theoryof deformationofa porous viscoelasticanisotropic solid[J].Journal ofApliedPhysics，1956，27（5）： 459-467.

[7]MitchellJK，Soga K.Fundamentals of soil behavior[M].3rd ed.Hoboken，NJ： John Wileyamp; Sons，2005.

[8]Terzaghi K. Theoretical Soil Mechanics[M].New York： Wiley，1943.

[9]Bjerrum L.Engineering geologyof Norwegian normally-consolidated marine clays asrelated tosetlementsof buildings[J]. Géotechnique，1967，17（2）： 83-118.

[10]CabalarAF，HasanRA.Compresionalbehaviourof various size/shapesand-claymixtureswithdifferentporefluids[J]. EngineeringGeology，2013，164：36-49.

[11]YinJH，F(xiàn)engWQ.Anewsimplified methodanditsverificationforcalculationofconsolidationsetlementofaclayeysoil with creep[J]. Canadian Geotechnical Journal，2017， 54（3）： 333-347.

[12]PianTHHissbsdssrtis[J]Aal6，7）：3.

[13]Ghosh S，MallettRL.Voronoi cellfinite elements[J].Computersamp; Structures，1994，50（1）：33-46.

[14]ZhangJ，KatsubeN.Ahybridfinite element methodforheterogeneous materials withrandomlydispersedelastic inclusions[J]. Finite Elements inAnalysis and Design，1995，19（1/2）： 45-55.

[15]Zhng J，KatsubeN.Ahybridfinite elementmethodforheterogeneous materials withrandomlydispersed rigidinclusions[J]. International Journal for Numerical Methods in Engineering，1995，38（10）：1635-1653.

[16]FuXR，CenS，LiCF，etal.Analyticaltrialfunction method fordevelopmentofnew8-node planeelementbasedonthe variational principle containing Airy stress function[J].Engineering Computations，2010，27（4）： 442-463.

[17]CenS，F(xiàn)uXR，ZhouGH，etal.Shape-fre finiteelement method：theplanehybrid stres-function（HS-Felementmethod for anisotropic materials[J]. Science China Physics，Mechanics and Astronomy，2011，54（4）： 653-665.

[18]Cen S，F(xiàn)uXR，ZhouMJ.8-and2-node planehbridstres-functionelements immunetoseverelydistortedmeshcontaining elements withconcaveshapes[J].ComputerMethodsinApplied MechanicsandEnginering，20120（29/30/31/32）：2321- 2336.

[19]Cen S，ZhuMJ，F(xiàn)uXR.A4-odehbrdtre-function（HS-F）planeelementwithdrillngdegresoffreeomlesssiive to severe mesh distortions[J].Computers amp; Structures，2011，89（5/6）：517-528.

[20]ZhouPL，CenS.Anovelshape-freeplanequadratic plygonal hybridstres-functionelement[J].MathematicalProblemsin Engineering，2015，2015：491325.

[21]楊鋒，郭然.多邊形應(yīng)力雜交單元的接觸算法研究[J].應(yīng)用數(shù)學(xué)和力學(xué)，2019，40（10）：1059-1070. Yang F，GuoR.Studyoncontact algorithms for the polygonal hybrid stresselement method[J].Applied Mathematics and Mechanics，2019，40（10）：1059-1070.（in Chinese）

[22]TongP，Pian THH.Avariational principleandtheconvergenceofafinite-element methodbasedonasumed stress distribution[J]. International Journal of Solids and Structures，1969，5（5）： 463-472.

[23]PianTHH，TongP.Basisoffiniteelementmethodsforsolidcontinua[J].InternationalJournalforNumericalMethodsin Engineering，1969，1（1）： 3-28.

[24]Pian THH.Nonlinearcreepanalysis byassumed stressfinite elementmethods[J].AIAAJoural，1974，2（12）：1756-1758.

[25]AtluriS.Onthehybridstressfiniteelementmodelforincrementalanalysisoflargedeflectionproblems[J].Interational Journal of Solids and Structures，1973，9（10）：1177-1191.

[26]李方政，樓根達，余志松，等.基于蠕變效應(yīng)的穿越隧道凍結(jié)帷幕開挖與支持三維數(shù)值模擬[J].巖土力學(xué)，2005，26（S1）： 121-132. LiFZ，Lou GD，YuZS，etal.3Dcreep numericalsimulationof excavationand supportinfrozen soil curatin wallofrossing tunnel[J]. Rock and Soil Mechanics，2005，26（S1）：121-132. （in Chinese）

[27］郭海柱，張慶賀，郭健.時間硬化與 Druker-Prager耦合蠕變模型的應(yīng)用研究[J].結(jié)構(gòu)工程師，2008，24（3）：117-121. Guo H Z， Zhang Q H，Guo J.Application studyon time hardening and druker-prager yield criterioncoupled creep model of soil[J]. Structural Engineers，2008，24（3）： 117-121. （in Chinese）

（編輯 鄭潔）