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

中圖分類號(hào)：0344.6 文獻(xiàn)標(biāo)志碼：A 文章編號(hào)：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è)過程中，路基沉降會(huì)對(duì)所建道路的質(zhì)量產(chǎn)生重要影響。由于地理環(huán)境、氣候條件和地質(zhì)成因等原因，不同地區(qū)的土壤質(zhì)地特征差異明顯。在路基施工過程中及施工后，這些土壤特征可能導(dǎo)致嚴(yán)重問題，如不均勻差異沉降邊坡坍塌和山體滑坡等，影響道路安全、運(yùn)輸效率和道路使用壽命。因此，路基沉降量的研究成了公路建設(shè)中的一個(gè)重要問題[2]。在沉降分析時(shí)，必須考慮路基填土的蠕變特性[34]，因?yàn)闀?huì)產(chǎn)生重大影響[5]。

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

雜交應(yīng)力元法是Pian（卞學(xué)璜）2在1964年基于修正余能原理首次提出的，克服了傳統(tǒng)位移有限元法無(wú)法解決單元邊數(shù)大于4時(shí)位移場(chǎng)插值的困難。為了避免位移型多邊形單元的位移插值函數(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)造變得簡(jiǎn)單。1995年Zhang等[14-15提出了基于Hellinger-Reissner兩場(chǎng)變分原理的雜交多邊形有限元法，用于含夾雜和孔洞材料的分析。2010年Fu等[基于變分原理開發(fā)了新型8節(jié)點(diǎn)平面單元，在極端網(wǎng)格畸變下也能保持精度不變。2011年Cen等首次提出形狀自由（shape-free）有限元法，用以描述有限元性能不受單元形狀影響的特性。同年Cen等8針對(duì)各向同性問題成功構(gòu)造了2個(gè)平面8結(jié)點(diǎn)單元和2個(gè)平面12結(jié)點(diǎn)單元，針對(duì)各向同性和各向異性問題，構(gòu)造了平面4結(jié)點(diǎn)帶旋轉(zhuǎn)自由度膜元[19]。2015年Zhou等[20]成功構(gòu)造了一種新的二次多邊形雜交應(yīng)力函數(shù)單元，構(gòu)造方法簡(jiǎn)單且精度高于許多其他方法。2019年楊鋒等2提出了一種基于多邊形雜交應(yīng)力單元的優(yōu)化接觸算法（polygonal hybrid stress element method，PHSEM）。由于單元邊數(shù)可變，PHSEM能夠靈活適應(yīng)復(fù)雜的模型邊界與材料邊界，更容易劃分網(wǎng)格。

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

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

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

任意多邊形雜交應(yīng)力單元的示意圖如圖1所示，網(wǎng)格劃分由邊數(shù)不定的多邊形構(gòu)成，可根據(jù)材料的空間分布及結(jié)構(gòu)的特點(diǎn)，高效便捷地劃分網(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 為單元外邊界上外法線向量;增量步開始前已知的初始場(chǎng)變量： σ 為單元域 內(nèi)滿足平衡的應(yīng)力場(chǎng)； ε 為單元域內(nèi)滿足協(xié)調(diào)關(guān)系的應(yīng)變場(chǎng)； 為單元外邊界 的位移 和 分別為力邊界 ?Ω_t 給定外力和給定外力增量。未知的場(chǎng)變量： Δσ 為單元域內(nèi)的應(yīng)力增量； Δu 為單元外邊界上的位移增量。

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

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

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

式中， Δε^c 為單元域內(nèi)的蠕變應(yīng)變?cè)隽俊?/p>

余能增量 ΔB 滿足

式中： Δε^e 為單元域內(nèi)彈性應(yīng)變?cè)隽浚?Δε^p 為單元域內(nèi)塑性應(yīng)變?cè)隽俊?/p>

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

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

abla?σ=0

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

σ_e=Pβ，

Δσ_e=PΔβ，

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

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

Δu=LΔq，

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

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

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

1.3 單剛推導(dǎo)

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

可以得到每個(gè)單元內(nèi)運(yùn)動(dòng)關(guān)系的弱表達(dá)式為

在上一增量步中滿足

則式（16）可以表示為

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

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

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

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

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

第2～i迭代步：

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

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

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

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

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

其中，剛度矩陣為

K=G^TH^-1G

2 土壤蠕變方程

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

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

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

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

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

3 數(shù)值算例

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

3.1 數(shù)值算例一

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

模型所用材料為凍土[2]，，其材料力學(xué)參數(shù)如表1所示。計(jì)算模型考慮在施加載荷下產(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ù)擬合時(shí)， σ 為等效應(yīng)力， Pa;T 為溫度， ^°C;t 為時(shí)間，h。

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

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

圖6擬合了PHSEM模型和MARC模型在路徑曲線上的應(yīng)力值?？梢宰⒁獾?個(gè)模型在水平方向應(yīng)力值、豎直方向應(yīng)力值和等效馮米塞斯應(yīng)力值的擬合趨勢(shì)都是高度一致的。由此可見，PHSEM的計(jì)算結(jié)果和MARC的計(jì)算結(jié)果較吻合，由PHSEM計(jì)算的整體應(yīng)力場(chǎng)分布是正確的。

圖7擬合了PHSEM模型和MARC模型在路徑曲線上的蠕變應(yīng)變值?？梢宰⒁獾?個(gè)模型在水平方向蠕變應(yīng)變值、豎直方向蠕變應(yīng)變值和等效蠕變應(yīng)變值的擬合趨勢(shì)基本都是一致的，誤差也在合理范圍內(nèi)。由此可見，PHSEM的計(jì)算結(jié)果和MARC的計(jì)算結(jié)果比較吻合，由PHSEM計(jì)算的整體蠕變應(yīng)變場(chǎ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所示，地基為長(zhǎng) 70m 寬 20m 的長(zhǎng)方形，路堤為上底 18m 、下底 26m 高 6m 的梯形。為了方便對(duì)比分析，PHSEM模型和MARC模型所受邊界條件和載荷完全相同，地基底部采用完全固定的約束，限制地基頂部的水平位移，對(duì)路堤施加豎直向下 400kN 的載荷。PHSEM模型采用7個(gè)節(jié)點(diǎn)，2個(gè)任意多邊形雜交應(yīng)力單元，250項(xiàng)的應(yīng)力函數(shù)多項(xiàng)式（21階的Airy應(yīng)力函數(shù)），網(wǎng)格劃分如圖9（a）所示;MARC模型采用1653個(gè)節(jié)點(diǎn)，1556個(gè)平面應(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所示。計(jì)算模型考慮在施加載荷下產(chǎn)生的蠕變，采用的蠕變參數(shù)如表4所示[27]，模擬路基 1000h 的蠕變過程。

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

Table3Mechanical parametersofeachsoillayer

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

Table4Creepparametersofeachsoil layer

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

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

圖10 應(yīng)力云圖

Fig. 10 Contoursof stress

圖11 蠕變?cè)茍D

Fig.11 Contoursof creep strain

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

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

4結(jié)論

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

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

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

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

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

參考文獻(xiàn)

[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]李方政，樓根達(dá)，余志松，等.基于蠕變效應(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］郭海柱，張慶賀，郭健.時(shí)間硬化與 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）

（編輯 鄭潔）