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

王哲學(xué)1, 2,王建華1, 2,程星磊1, 2

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

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

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

中圖分類號：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—)， 男, 江蘇泰興人, 高級工程師, 主要從事水泵、水閘設(shè)計與研究工作。

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

1本構(gòu)模型簡述

采用Mises屈服準則構(gòu)建邊界面約束方程：

(1)

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

(2)

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

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

(3)

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

(4)

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

(5)

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

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

(6)

式中：σz為試樣軸向應(yīng)力；σr為試樣徑向應(yīng)力。

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

(7)

(8)

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

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

(9)

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

(10)

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

(11)

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

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

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

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

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

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

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

2.1平衡迭代過程

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

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

2.2應(yīng)力更新算法

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

(12)

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

(13)

由此得到應(yīng)力增量大?。?dσij=dσmδij+dsij

(14)

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

(15)

得到增量形式的彈塑性應(yīng)力應(yīng)變關(guān)系，整理成矩陣形式：dσ=Depdε

(16)

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

因此，本構(gòu)關(guān)系的積分式為:Δσ=∫Depdε

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

(17)

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

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

(18)

圖4　有限元預(yù)測結(jié)果與單元預(yù)測結(jié)果比較 Fig.4 Comparison between FEM results and unit predicted results

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

3有限元程序驗證

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

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

4三軸試驗?zāi)M

表2　三軸試驗?zāi)Ｐ蛥?shù)

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

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

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

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

圖6　三軸應(yīng)力-應(yīng)變曲線 Fig.6 Triaxial stress-strain curves

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

5結(jié)語

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

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

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

(4)通過有限元方法進行二次開發(fā)，成功地將此本構(gòu)模型應(yīng)用到了數(shù)值計算中，對三軸壓縮和拉伸試驗進行預(yù)測，結(jié)果證明此本構(gòu)模型數(shù)值實現(xiàn)方法合理。對于更復(fù)雜的應(yīng)力狀態(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]．寧夏工程技術(shù)， 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)應(yīng)力應(yīng)變響應(yīng)的彈塑性模擬[J].巖土工程學(xué)報，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]張學(xué)言，閆澍旺. 巖土塑性力學(xué)基礎(chǔ)[M]. 天津：天津大學(xué)出版社，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]．工程力學(xué)，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]朱伯芳. 有限單元法原理與應(yīng)用[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]. 北京：清華大學(xué)出版社，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]. 水利水運工程學(xué)報, 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