軟黏土不排水循環(huán)應力應變關系的數(shù)值模擬

軟黏土不排水循環(huán)應力應變關系的數(shù)值模擬

王哲學1, 2,王建華1, 2,程星磊1, 2

(1. 天津大學水利工程仿真與安全國家重點實驗室, 天津300072； 2. 天津大學 巖土工程研究所，天津300072)

摘要：對程星磊等提出的總應力形式增量彈塑性本構模型進行二次開發(fā)，以模擬復雜應力狀態(tài)下軟黏土的響應。通過Newton-Raphson算法，對材料非線性問題進行迭代求解；針對本構模型中應力反向等一系列關鍵性問題，應用歐拉切線算法編寫了有限元程序，并結合子增量方法提高了計算精度。預測了軟黏土在軸對稱應力狀態(tài)下的響應，得到了應力應變關系曲線，將其與單元預測結果進行比較，二者趨于一致，從而驗證了有限元程序編寫的合理性。利用該程序模擬三軸不固結不排水試驗，模擬結果與試驗吻合良好，表明該本構模型有限元程序可以反映軸對稱應力狀態(tài)下軟黏土的不排水應力應變特性，可應用于更加復雜邊值問題的模擬計算。

關鍵詞：循環(huán)荷載； 飽和軟黏土； 應力應變； 本構模型； 歐拉切線算法； 三軸試驗

中圖分類號：TU43

文獻標志碼：A

文章編號：1009-640X(2015)03-0081-07

Abstract：An incremental elasto-plastic constitutive model in the form of total stress developed by CHENG Xing-lei et al can effectively reflect the undrained cyclic stress-strain response of the saturated soft clay. In order to simulate the response of the soft clay under complex stress states，a secondary development for the constitutive model is conducted by use of the finite element method. The Newton-Raphson method is applied to solve iteratively the material nonlinear problem; considering a series of key problems such as stress reversing, in programming the model, Euler tangent algorithm is adopted，with a sub-incremental method improving calculation accuracy. The response of the soft clay under axisymmetric stress states is predicted, and the stress-strain relationship curves are received, which is compared with the results given by the unit predicting. The analysis comparison results show that it is in good agreement, verifying the rationality of the finite element program. The finite element program is applied to simulate unconsolidated undrained triaxial tests， and the results are in good agreement with the experiments, which indicates that the finite element program of the constitutive model can describe the undrained stress-strain behavior of the soft clay under triaxial stress states. Furthermore, it can be applied to solve more complex boundary value problems.

DOI:10.16198/j.cnki.1009-640X.2015.03.013

收稿日期：2014-09-05

基金項目：江蘇省水利科技項目(2013057)

作者簡介：丁軍(1967—)， 男, 江蘇泰興人, 高級工程師, 主要從事水泵、水閘設計與研究工作。

嵌固在海洋軟黏土中的錨固基礎在風、浪等荷載作用下，為維持上部結構工作性能，必須保持自身穩(wěn)定，因此需要合理分析軟黏土中錨固基礎在循環(huán)荷載作用下的變形失穩(wěn)過程。錨固基礎的變形不僅與基礎本身有關，海洋中軟土地基對其影響更為重要。因此研究循環(huán)荷載下軟黏土的不固結不排水應力應變關系具有重要意義。J. H. Prévost[1-2]提出了基于運動硬化塑性理論的多面模型，但模型的變形分析很復雜，需要記憶所有假設的多個次加載面形狀，并分析它們的移動，導致在數(shù)值應用中計算過程復雜。Y. F. Dafalias等[3-4]建立了包括1個正常屈服面和1個次加載面的兩面模型，大大簡化了多面模型；后來他們又將兩面模型中的次加載面簡化為一點，提出了單面模型[5]。兩面和單面模型都叫邊界面模型[6-8]。由于改進初期的邊界面模型映射中心固定在原點，無法合理描述土單元的應力應變滯回曲線[9]，且硬化模量場的演化過程表述較為復雜，不便于應用到數(shù)值分析中。程星磊等[10]采用總應力形式在偏應力空間中構建硬化模量場，硬化模量插值函數(shù)數(shù)學表達簡潔，引入模型參數(shù)較少，硬化模量場的演化規(guī)律表述簡單，跟蹤循環(huán)應力路徑所需記憶的參數(shù)較少，便于計算應用。通過上述分析，針對程星磊等建立的軟黏土循環(huán)彈塑性本構模型進行二次開發(fā)，并對循環(huán)三軸試驗進行數(shù)值模擬，通過計算結果與試驗結果的比較，驗證該本構模型數(shù)值分析方法的可行性。

1本構模型簡述

采用Mises屈服準則構建邊界面約束方程：

(1)

依據(jù)經(jīng)典彈塑性理論[11]，有總偏應變增量計算式：

(2)

式中：deij為總偏應變增量；G為彈性剪切模量；H′為塑性硬化模量；dskl為當前應力點的偏應力增量。

引入彈塑性硬化模量H，令：

(3)

在初始加載或應力反向時，塑性變形模量為無窮大：

(4)

圖1　彈塑性模量的映射規(guī)則圖解 Fig.1 Elasto-plastic mapping rule for loading and unloading

(5)

式中：Hmax為彈塑性模量最大值，取2G；μ為反映剪應變累積速率和大小的模型參數(shù)。

在不固結不排水條件下，邊界面半徑A0和彈性模量G可通過不固結不排水循環(huán)三軸試驗得到。在三軸應力狀態(tài)下，式(1)簡化為下式：

(6)

式中：σz為試樣軸向應力；σr為試樣徑向應力。

由于初始應力應變響應均遵循同一初始加荷曲線，且可用式(7)擬合，則擬合曲線的漸近線所對應的強度就是邊界面半徑A0，其值等于1/b；曲線的初始切線模量記為E，其值等于1/a。另外，彈性剪切模量G可由式(8)求得：

(7)

(8)

式中：a和b為擬合函數(shù)中的待定系數(shù)；ν為泊松比，不排水時取0.5。

參數(shù)μ可以用下式表示：

(9)

式中：μ0為反映累積偏應變大小的待定系數(shù)；dγ8p為八面體累積剪應變增量，以此描述偏應變累積變化的速率和趨勢。對試驗數(shù)據(jù)采用式(10)擬合：

(10)

式中：c，d為待定系數(shù)；N為循環(huán)次數(shù)。則有dγ8p=cdNd-1

(11)

下面結合某一具體三軸試驗，給出待定系數(shù)c，d，μ0的確定方法。圖2給出了該三軸試驗得到的靜應力比τ8,a/τ8,f=0.3和0.5時八面體累積剪應變隨振次的變化關系。 用式(10)擬合試驗數(shù)據(jù)，如圖中虛線所示，得到不同靜應力比τ8,a/τ8,f，不同循環(huán)應力比τ8,cy/τ8, f下的待定系數(shù)c和d值。通過式(11)可得應力循環(huán)一周八面體累積剪應變增量dγ8p。

圖2　三軸壓縮時八面體累積剪應變隨循環(huán)次數(shù)的變化關系曲線 Fig.2 Relation curves of octahedral cumulative shearing strain and cycle times in triaxial compression

進一步，可通過試算方法確定μ0，即利用式(9)將dγ8p代入式(5)，給定一個μ0，通過式(2)計算得到破壞振次下的八面體累積剪應變。這里規(guī)定按照靜偏應力作用下試樣的八面體靜剪應變與靜偏應力與循環(huán)應力共同作用下八面體循環(huán)累積剪應變之和達到14%的標準確定循環(huán)破壞次數(shù)。若模型計算所得破壞振次下的八面體累積剪應變與試驗結果一致，那么此時的μ0值即為該試驗應力條件下待定系數(shù)μ0的取值。

這樣，不同靜應力比τ8,a/τ8,f，不同循環(huán)應力比τ8,cy/τ8,f下的參數(shù)c，d，μ0都已給出(見表1)。則對任意循環(huán)應力比下的參數(shù)值，可通過線性插值方法確定。

表1　不同靜應力比、循環(huán)應力比下擬合函數(shù)中的待定系數(shù)c，d及μ 0取值

至此，軟黏土不固結不排水條件下循環(huán)本構關系及模型參數(shù)的確定方法已完全給出。

2本構模型二次開發(fā)

2.1平衡迭代過程

圖3　增量內(nèi)迭代過程 Fig.3 Iterative process in increment

非線性平衡方程采用Newton-Raphson迭代算法[12]求解。在迭代過程中，把整個加載過程分成若干分段，荷載以增量形式施加在結構上，在增量內(nèi)部進行若干迭代，直至達到所需精度為止。整個流程融合了增量與迭代的過程[13]。在增量內(nèi)的迭代過程(見圖3)作如下說明[14]：①根據(jù)當前的應力狀態(tài)nσ ，確定彈塑性矩陣Dep，進而確定剛度矩陣K 。②由Δu=K-1ΔR ，得到位移增量。其中，ΔR 為外加荷載增量，Δu 為位移增量。③由Δε=BΔu ，得到應變增量。④采用應力更新算法，通過當前應力狀態(tài)以及上一步得到的應變增量計算應力增量Δσ，進而得到新的應力狀態(tài)n+1σ。⑤判斷內(nèi)力、外力是否平衡，位移是否滿足收斂標準：若不滿足，返回①，繼續(xù)進行迭代；若滿足，進入下一增量計算。

2.2應力更新算法

對第④步的應力更新算法做出詳細闡述，由應力應變關系式(2)，可得：

(12)

由上述關系式可得到偏應力增量。靜水壓力增量可由彈性虎克定律得到：dσm=3Kdεm

(13)

由此得到應力增量大?。?dσij=dσmδij+dsij

(14)

將式(8)和(9)代入式(10)，有：

(15)

得到增量形式的彈塑性應力應變關系，整理成矩陣形式：dσ=Depdε

(16)

式中：dσ為應力張量；dε為應變張量；Dep為彈塑性矩陣，由于在應力應變關系曲線上Dep為切線方向，因此也叫切線矩陣。

因此，本構關系的積分式為:Δσ=∫Depdε

隨著機動車保有量的迅猛增長，人們就醫(yī)出行選擇開私家車前往醫(yī)院的現(xiàn)象越來越普遍，大大增加了醫(yī)院周邊交通流。隨著生活水平的改善，人們對于更好的醫(yī)療服務需求顯得更加重視，大型醫(yī)院有著更好的醫(yī)療硬件設施和專家門診服務，所以前往大型醫(yī)院看病的人員越來越多。在這樣的背景之下，涌向大型醫(yī)院的人流與車流已經(jīng)遠遠超出了醫(yī)院設計之初的承載能力，而醫(yī)院只是對建筑進行擴建來增加醫(yī)院的就醫(yī)能力，醫(yī)院容量的擴大則會帶來更大的交通需求，一味地對醫(yī)院進行擴建不但不能從根本上解決問題，反而會增加醫(yī)院周邊的交通流，導致供需失衡越來越嚴重。

(17)

在計算過程中，采用式(16)將式(17)線性化。但是這種計算方法必然帶來偏差[15]。可采取減小增量大小的方法減小這種偏差。

當應力狀態(tài)趨于邊界面附近時，由式(5)得硬化模量趨于零，由式(17)得偏應力增量趨于零，因此應力狀態(tài)點將始終保持在邊界面以內(nèi)。但實際計算中，由于采用式(16)所示的歐拉切線方法對式(17)進行線性化，而歐拉方法只有在偏應變增量趨于無限小的時候才能足夠準確，因此在邊界面附近時，應力增量dσij有可能大于零，從而使應力狀態(tài)超出邊界面。但邊界面是破壞面，土體的應力狀態(tài)不允許超過邊界面。這將導致計算出錯。采用子增量方法雖然大大減小了偏應變增量的大小，但是也不足以保證應力狀態(tài)不超過邊界面。為此，將超出邊界面的點沿邊界面的半徑方向徑向拉回至邊界面，并認為這一點的應力狀態(tài)即為當前的應力狀態(tài)。具體做法是，計算應力空間原點與當前應力點之間距離，按照邊界面半徑與該距離比例對當前應力點各個應力分量進行線性縮放，則縮放后的應力點必然落在邊界面上。進一步利用修正后的應力狀態(tài)建立與之相適應的彈塑性矩陣Dep。這樣就保證了任意時刻的應力狀態(tài)只能在邊界面以內(nèi)或者在邊界面之上；而到達邊界面的點在持續(xù)加載的條件下，應力狀態(tài)始終保持在邊界面上，類似于理想彈塑性模型中的塑性加載狀態(tài)。與之相適應的彈塑性矩陣，則保證了應力與應變狀態(tài)的一致性，有利于后續(xù)分析的收斂。

(18)

圖4　有限元預測結果與單元預測結果比較 Fig.4 Comparison between FEM results and unit predicted results

如果這一余弦值為正或零，說明應力增量張量與當前邊界面外法線張量呈銳角或直角，則此時沒有應力狀態(tài)改變；如果這一余弦值為負，說明應力增量張量與當前邊界面外法線張量呈鈍角，則認為此時應力狀態(tài)發(fā)生改變。三軸應力狀態(tài)只有夾角為0°(應力狀態(tài)不變)和180°(應力狀態(tài)反向)兩種應力狀態(tài)的改變形式。

3有限元程序驗證

把有限元程序和單元預測循環(huán)三軸不固結不排水試驗得到的應力應變曲線進行對比(圖4)。實線為有限元計算結果，虛線為單元預測結果。圍壓為100 kPa，靜強度大小為τ8,f=16 kPa，靜剪應力比τ8,a/τ8,f=0.5，循環(huán)剪應力比τ8,a/τ8,f=0.392，其他參數(shù)為：A0=31 kPa,G=745.76 kPa,c=0.861,d=0.606,μ0=8。

圖4表明數(shù)值算法結果與基于此本構模型進行的單元預測結果吻合良好，說明數(shù)值算法程序能夠正確反映本構模型的應力應變關系，此程序編寫是合理的。

4三軸試驗模擬

表2　三軸試驗模型參數(shù)

針對前文所述三軸試驗，靜應力比τ8,a/τ8,f=0.5的三軸壓縮試驗和三軸拉伸試驗分別進行數(shù)值計算。根據(jù)表1，通過插值的方法確定模型所采用的參數(shù)(表2)。

將循環(huán)八面體累積剪應變隨循環(huán)次數(shù)的變化關系進行比較(圖5)。對于圖5(a)中τcy/τf=0.334的試驗模擬，當循環(huán)應變大于10%時，計算得到的曲線與試驗曲線略有差異，但循環(huán)累積應變小于10%時，二者吻合良好。如上所述，當靜應變與循環(huán)累積應變之和達到14%時，即視為土體破壞。這里靜應變約為4%，因此循環(huán)累積應變達到10%時，試樣已經(jīng)破壞，即在土樣破壞之前可以較好地模擬試驗結果。圖5(b)中τcy/τf=0.180的試驗模擬，循環(huán)次數(shù)達900次，模擬結果與試驗有較好的一致性，且曲線趨于平穩(wěn)，因此該模型可以描述循環(huán)次數(shù)較大時土體的應力應變響應。綜上，相同振次下的累積變形預測結果與試驗結果較為一致，此本構模型的有限元程序可以較好地反應循環(huán)變形的累積特性。

圖5　循環(huán)累積剪應變隨振次的變化 Fig.5 Cyclic cumulative shear strain vs. the number of cycles

將靜剪應力比τ8,a/τ8,f=0.5、循環(huán)剪應力比τ8,cy/τ8,f=0.392的應力應變曲線進行比較(圖6)。結果表明，模型可較好地描述循環(huán)加載時應力應變曲線的滯回特性及應變累積特性。

圖6　三軸應力-應變曲線 Fig.6 Triaxial stress-strain curves

上述分析表明，該本構模型有限元程序可以描述循環(huán)三軸應力狀態(tài)下軟黏土的應力應變關系，可進一步推廣至更加復雜的應力狀態(tài)。

5結語

(1)在數(shù)值實現(xiàn)過程中，采用Euler切線方法對本構關系積分進行線性化，并采用子增量法提高了計算精度，加快了后續(xù)迭代的收斂速度。

(2)由于歐拉方法不可避免會使應力狀態(tài)超出邊界面，因此采用徑向返回的方法將其拉至邊界面上，以保證任何單元的應力狀態(tài)均不會超出具有破壞意義的邊界面。

(3)在偏應力空間中，采用應力增量與像應力點在邊界面上外法線方向的夾角判斷應力狀態(tài)改變，確定加載卸載狀態(tài)。

(4)通過有限元方法進行二次開發(fā)，成功地將此本構模型應用到了數(shù)值計算中，對三軸壓縮和拉伸試驗進行預測，結果證明此本構模型數(shù)值實現(xiàn)方法合理。對于更復雜的應力狀態(tài)還有待于進一步研究。

參考文獻：

[1]PRéVOST J H. Mathematical modelling of monotonic and cyclic undrained clay behaviour[J]. International Journal for Numerical and Analytical Methods in Geomechanics, 1977, 1(2): 195- 216.

[2]PRéVOST J H. Anisotropic undrained stress-strain behavior of clays[J]. Journal of the Geotechnical Engineering Division, ASCE, 1978, 104(GT8): 1075- 1090.

[3]DAFALIAS Y F, POPOV E P. A model of nonlinearly hardening materials for complex loading[J]. Acta Mechanica, 1975, 21(3): 173- 192.

[4]DAFALIAS Y F, HEMNANN L R. A boundary surface soil plasticity model[C]∥ Proceedings of International Symposium on Soils under Cyclic Transient Load. Rotterdam: Balkema Publ, 1980: 335- 356.

[5]DAFALIAS Y F, POPOV E P. Cyclic loading for materials with a vanishing elastic domain[J].Nuclear Engineering and Design, 1977, 41(2): 293- 302.

[6]DAFALIAS Y F. Bounding surface plasticity Part I: Mathematical foundation and hypoplasticity[J]. Journal of Engineering Mechanics, 1986, 112(9):966- 987.

[7]DAFALIAS Y F, HERRMANN L R. Bounding surface plasticity. Part II: Application to isotropic cohesive soils[J]. Journal of Engineering Mechanics, 1986, 112(12): 1263- 1291.

[8]ANANDARAJAH A , DAFALIAS Y F. Bounding surface plasticity. Part III: Application to isotropic cohesive soils[J]. Journal of Engineering Mechanics, 1986, 112(12): 1292- 1318.

[9]花麗坤，孔亮. 幾種循環(huán)塑性模型的剖析和比較研究[J]．寧夏工程技術， 2003, 2(1)：19- 23. (HUA Li-kun, KONG Liang. Analysis and comparison of several cyclic plastic models[J]．Ningxia Engineering Technology, 2003, 2(1): 19- 23.(in Chinese))

[10]程星磊，王建華，李書兆. 軟黏土不排水循環(huán)應力應變響應的彈塑性模擬[J].巖土工程學報，2014, 36(5): 933- 941. (CHENG Xing-lei, WANG Jian-hua, LI Shu-zhao. Elastoplastic simulation for undrained cyclic stress-strain response of soft clay[J]. Chinese Journal of Geotechnical Engineering, 2014, 36(5): 933- 941. (in Chinese))

[11]張學言，閆澍旺. 巖土塑性力學基礎[M]. 天津：天津大學出版社，2003: 88- 115. (ZHANG Xue-yan，YAN Shu-wang. Fundamentals of geotechnics plasticity[M]. Tianjin: Tianjin University Press, 2003: 88- 115. (in Chinese))

[12]王軍祥，姜諳男. 完全隱式返回映射算法對巖土地基問題的求解[J]．工程力學，2013, 36(8): 83- 89. (WANG Jun-xiang, JIANG An-nan. Fully implicit return mapping algrithm for solving the problems of geotechnical foundation[J]. Engineering Machanics, 2013, 36(8): 83- 89. (in Chinese))

[13]朱伯芳. 有限單元法原理與應用[M]. 北京:中國水利水電出版社, 1998: 292- 302. (ZHU Bo-fang. The finite element method theory and applications[M]. Beijing: China WaterPower Press, 1998: 292- 302. (in Chinese))

[14]王勖成，邵敏. 有限單元法基本原理和數(shù)值方法[M]. 北京：清華大學出版社，1997: 510- 511.(WANG Xu-cheng, SHAO Min. The basic principles of the finite element method and numerical methods[M]. Beijing: Tsinghua University Press, 1997: 510- 511. (in Chinese))

[15]HU Cun, LIU Hai-xiao. Implicit and explicit integration schemes in the anisotropic bounding surface plasticity model for cyclic behaviours of saturated clay[J]．Computers and Geotechnics, 2014, 55: 27- 41.

Numerical simulation of undrained cyclic stress-strain response of soft clay

WANG Zhe-xue1,2, WANG Jian-hua1,2, CHENG Xing-lei1,2

(1.StateKeyLaboratoryofHydraulicEngineeringSimulationandSafety,TianjinUniversity,Tianjin300072,China； 2.GeotechnicalEngineeringInstitute,TianjinUniversity,Tianjin300072,China)

Key words： cyclic load; saturated soft clay; stress-strain curve; a constitutive model; Euler tangent algorithm; triaxial tests

丁軍， 丁慶朋. 胥浦活水泵站肘形進水流道流態(tài)分析及優(yōu)化[J]. 水利水運工程學報, 2015(3): 88-94. (DING Jun, DING Qing-peng. Flow pattern analysis and optimization of elbow inlet conduits of Xupu running water pumping station[J]. Hydro-Science and Engineering, 2015(3): 88-94.)

E-mail：djdj1031@163.com