基于粒子群優(yōu)化算法的軟土場地Davidenkov模型參數(shù)擬合與應用

王沿朝， 陳清軍

(同濟大學 土木工程防災國家重點實驗室，上海 200092)

土的動力非線性本構模型是土動力學研究的重要課題之一，也是土層及土-地下結構相互作用非線性地震反應分析中必須要考慮的內容。土體是一種非線性很強的材料，在循環(huán)荷載作用下，其應力-應變關系主要表現(xiàn)出非線性、滯后性和變形累積的特性。研究表明，在地震動的作用下，土體幾乎不存在線彈性階段，且隨著地震動強度的增大，其非線性程度亦越強[1]。因此，如何對土體在動荷載作用下的應力-應變關系進行正確描述一直是研究者們所關注的問題。

目前，對于土體動力非線性本構模的研究主要集中在以下三個方面[2-7]：①等效線性化方法；②真非線性動力本構模型；③時域滯回非線性分析方法。在這三種方法中，等效線性化方法形式簡單直觀，計算量少，但由于不能計算永久變形，在強震動時的計算結果誤差大，且容易產(chǎn)生“虛共振”現(xiàn)象等缺點，因而也制約了其進一步的發(fā)展。真非線性動力本構模型雖然能夠準確的對土體的動應力-應變關系進行描述，然而由于其模型參數(shù)較多，實現(xiàn)起來較為復雜，且通常很難通過常規(guī)試驗獲得，因此其應用范圍也不是很廣。時域滯回非線性分析方法兼有上述二者的優(yōu)點[8]，比如Davidenkov模型，既能較好的對土體的非線性行為進行描述，同時其表達形式相對于真非線性動力本構模型更為簡單，便于推廣應用。

為探究Davidenkov模型在軟土場地地震反應分析中的應用，本文將推導Davidenkov模型的非線性動力本構模型在初始加載、卸載和再加載階段的增量剪切模量表達式，并將土體動應力-應變關系推廣到三維?；诹Ｗ尤簝?yōu)化算法對Davidenkov模型的參數(shù)進行擬合，利用MATLAB編制擬合程序，通過典型砂土和黏土驗證擬合效果。在此基礎上，將利用ABAQUS軟件提供的二次開發(fā)平臺UMAT，編制基于Davidenkov模型的土體非線性動力本構模型計算子程序，利用復雜加載路徑驗證該子程序的正確性，實現(xiàn)在ABAQUS軟件中軟土非線性動力本構模型的二次開發(fā)，并通過對某典型上海軟土場地的地震反應計算，對比分析上海軟土場地基于非線性動力本構模型的計算結果和等效線性化結果間的差異。

1 Davidenkov模型表達式及其參數(shù)擬合

1.1 Davidenkov骨架曲線理論公式

Hardin等[9]提出，土體的應力-應變關系可以通過一種雙曲線形式來描述，其動剪切模量比的表達式如下

G/Gmax=1-H(γ)

(1)

(2)

式中：G為剪切模量；Gmax為最大剪切模量；γ為剪應變幅；γref為參考剪應變，一般可取動剪切模量比為0.5時所對應的剪應變幅。

Martin等[10]采用Davidenkov骨架曲線來描述上述關系，將式(2)改寫如下

(3)

式中：A，B和γ0均為和土性有關的擬合參數(shù)。此時，γ0已不再是具有明確物理意義的參考剪應變，而僅為一個擬合參數(shù)。由式(3)可知，當A=1.0，B=0.5，γ0=γref時，Davidenkov模型即退化成Mashing雙曲線模型。

因此，Davidenkov模型的應力-應變關系的骨架曲線可用下式表示

τ(γ)=Gγ=Gmaxγ[1-H(γ)]

(4)

將式(3)代入到式(4)中，即可得到

(5)

由式(5)可以發(fā)現(xiàn)，當B<0.5時，若γ→ ∞，則τ(γ) → ∞，而這與土體實際的應力應變關系是不相符的。為此，參考文獻[11]的研究方法，通過分段函數(shù)的形式來對Davidenkov模型的骨架曲線進行改進，即引入一個剪應變上限值γult，當土的剪應變γ超過γult時，便可認為該土體處于破壞狀態(tài)，此時若γ繼續(xù)增加，土體的剪應力τ也不再增加，甚至可能會減小。若不考慮軟化效應，則骨架曲線可表示為

(6)

τult=Gmaxγult[1-H(γult)]

(7)

而根據(jù)曼辛準則，以Davidenkov骨架曲線為基礎的土體應力-應變滯回曲線可表示為

(8)

式中：τc和γc分別表示的是應力-應變滯回曲線中加卸載轉折點處的剪應力與剪應變幅值。

1.2 剪切模量的增量表達式

由于土體在外荷載作用下幾乎不存在彈性階段，因此可以將土體假定為理想的黏彈性體，進而，可以從一般的模量衰減曲線得到歸一化的剪應力

(9)

對上式求γ的導數(shù)，即可得到

(10)

式中：Ms為歸一化的割線模量，γ為剪應變，Mt為歸一化的切線模量。則增量剪切模量G可以表示為

G=GmaxMt

(11)

對式(6)進行求導，即可得到Davidenkov骨架曲線中的增量剪切模量如下式

(|τ|≤τult)

(12)

對式(8)進行求導并令R=|γ-γc|/(2γ0)，則可得到卸載及再加載階段的增量剪切模量如下式

(|τ|≤τult)

(13)

而當|τ|>τult時，G=0。

由此，得到了基于Davidenkov骨架曲線及滯回曲線的完整的加載、卸載及再加載的增量剪切模量的表達形式。

1.3 基于粒子群優(yōu)化算法的土體參數(shù)的擬合

當采用基于Davidenkov骨架曲線的非線性動力本構模型時，必須根據(jù)土體的動剪切模量與動剪應變試驗曲線來對式(3)中的三個參數(shù)A、B、γ0進行擬合。參數(shù)擬合問題的本質是函數(shù)優(yōu)化，目前已有一些解決方法，如牛頓法、遺傳算法等。然而，當函數(shù)是超高維、多局部極值的復雜函數(shù)時，這些算法的效果并不理想。研究發(fā)現(xiàn)，對大多數(shù)的非線性函數(shù)，粒子群優(yōu)化算法在精度和優(yōu)化速度上均優(yōu)于上述算法[12]。本節(jié)基于粒子群優(yōu)化算法編寫了相應的優(yōu)化程序對Davidenkov模型的參數(shù)進行擬合。

粒子群優(yōu)化算法是一種群體智能算法，由Eberhart等[13-14]發(fā)明，在標準粒子群優(yōu)化算法中，每個優(yōu)化問題的解被看成是搜索空間的一個“粒子”，假設在一個D維目標搜索空間中，有M個粒子組成一個群體，其中在第t次迭代時，粒子pi的位置矢量為xi(t)=(xi1,…xiD)，速度矢量vi(t)=(vi1,…,viD)。在每次迭代中，粒子都通過追蹤兩個極值來對自己的速度和位置進行更新。一個是個體極值點(粒子自身找到的最好解)，表示為pbesti=(pbesti1,…,pbestiD)；一個是全局極值點(整個種群當前找到的最優(yōu)解)，表示為gbesti=(gbesti1,…,gbestiD)。在第t+1次迭代時，粒子pi根據(jù)式(14)和(15)來更新自己的速度和位置。

(14)

xid(t+1)=xid(t)+c3×vid(t+1)

(15)

式中:w為慣性權值；c1，c2為學習因子；c3為另一個慣性權值；rand1，rand2為兩個[0，1]之間的隨機數(shù)。

目前，常用來構建目標函數(shù)的方法包括互相關法、最小二乘法及基于信息熵的互信息法等。在本文中，目標函數(shù)衡量的是預測值與實際值差異的大小，目標值越大則擬合效果越好。因此，本文引入圖像配準[15]中經(jīng)常使用的歸一化互相關法來構建目標函數(shù)，該方法通過計算互相關度量值NCC來確定試驗數(shù)據(jù)與擬合結果的相似程度，其表達式如下

(16)

本文中粒子群優(yōu)化算法具體的計算流程如圖1所示。

圖1 粒子群優(yōu)化算法流程圖

基于以上理論，本文利用MATLAB編寫了擬合程序。為了驗證程序的可行性和正確性，首先對典型的砂土和黏土參數(shù)進行擬合，擬合中主要的參數(shù)變量為慣性權值w，加速因子c1，c2，種群數(shù)N，迭代次數(shù)M和粒子維數(shù)D。

(1) 種群規(guī)模。種群規(guī)模不能太大也不能太小，綜合考慮算法的精度及效率，本次擬合選擇種群規(guī)模為20。

(2) 最大迭代次數(shù)。迭代次數(shù)與解的收斂性成正比，但過大會影響運算速度，本文中選1 000次。

(3) 慣性權值。慣性權值w表征了粒子保留原來的速度的程度。w較大，全局收斂能力強；w較小，局部收斂能力強。本文選0.6。

(4) 加速因子。加速因子c1，c2分別用于控制粒子朝自身或鄰域最佳位置運動。文獻[16]建議φ=c1+c2≤4.0，并通常取c1=c2=2.0。本文也取c1=c2=2.0。

(5) 粒子維數(shù)。粒子維數(shù)即待優(yōu)化函數(shù)的維數(shù)。

下面，通過粒子群優(yōu)化算法對典型的砂土和黏土進行Davidenkov模型參數(shù)擬合，表1給出了Seed等[17]所給出的典型砂土和黏土的剪應變與剪切模量及阻尼比與剪應變幅的試驗數(shù)據(jù)。在對模型參數(shù)進行擬合時，其他擬合參數(shù)不變，粒子維數(shù)取為3(3個參數(shù)需要擬合)。

表1 典型砂土和黏土剪切模量比和阻尼比試驗數(shù)據(jù)

Tab.1 Shear modulus ratio and damping ratio of sand and clay

土類剪應變×10-5砂土黏土G/GmaxλG/Gmaxλ0.11.02.501.00.500.3160.9132.500.9840.8010.7612.500.9341.703.160.5653.500.8263.20100.4004.750.6565.6031.60.2616.500.44310.01000.1529.250.24615.53160.07613.80.11521.01 0000.03720.00.04924.63 1600.01326.00.04924.610 0000.00429.00.04924.6

典型砂土和黏土的擬合曲線如圖2和圖3所示，從圖中可以看出，對試驗數(shù)據(jù)擬合的曲線走勢很好，Davidenkov模型參數(shù)擬合的結果及NCC值見表2所示，可以看出，兩種土體擬合結果的NCC值均十分接近于1，說明了擬合效果很好。

圖2 砂土擬合曲線

圖3 黏土擬合曲線

土類ABγ×10-4NCC砂土1.100.403.250.999 6黏土1.010.812.000.999 7

在此基礎上，本文通過粒子群優(yōu)化算法對典型上海軟土進行Davidenkov模型參數(shù)擬合，表3給出了典型上海軟土各層土的剪切模量比及阻尼比與剪應變幅的試驗數(shù)據(jù)。對模型參數(shù)進行擬合時，其他擬合參數(shù)不變，粒子維數(shù)取為3。典型上海軟土各層土的擬合曲線如圖4所示，Davidenkov模型參數(shù)擬合的結果及NCC值如表4所示。

表3 典型上海軟土各土層剪切模量比和阻尼比

(a) 土1

(b) 土2

(c) 土3

(d) 土4

(e) 土5

(f) 土6

圖4 典型上海軟土各層土擬合曲線

Fig.4 Fitting curves of typical Shanghai soft soils

表4 典型上海軟土擬合結果

2 Davidenkov模型在ABAQUS中的實現(xiàn)

2.1 實現(xiàn)過程

ABAQUS為用戶提供了功能強大的用戶子程序接口，以幫助用戶開發(fā)基于ABAQUS內核的程序，而這其中應用最為廣泛的子程序即為UMAT[18]，它主要用于用戶開發(fā)自己的本構模型，以彌補ABAQUS自帶材料模型的不足。ABAQUS子程序通過FORTRAN語言開發(fā)。按照FORTRAN的語法，用戶可以自由的編輯代碼，并形成一個獨立的程序單元，該單元能被獨立的編譯和存儲，也可以被其他單元引用，利用UMAT能夠將大量的數(shù)據(jù)帶回供引用程序使用。

UMAT的主要任務是根據(jù)ABAQUS主程序傳入的應變增量更新應力增量和狀態(tài)變量(如有必要)，并給出雅克比矩陣供ABAQUS求解使用，主程序結合當前荷載增量求解位移增量，繼而進行平衡校核，如果不滿足指定的誤差，ABAQUS將進行迭代直到認為收斂，然后進行下一增量步的求解。因此，UMAT文件主要包含以下幾個部分：

(1) 子程序定義語句；

(2) ABAQUS定義的參數(shù)說明；

(3) 用戶定義的局部變量說明；

(4) 用戶編制的程序主體；

(5) 子程序返回和結束語句。

對于本文來講，ABAQUS的本構模型開發(fā)工作主要是修改用戶定義的局部變量及程序主體的編制，通過定義材料參數(shù)數(shù)組來實現(xiàn)初始參數(shù)的輸入，通過定義狀態(tài)變量數(shù)組來存儲與本構模型有關的參變量。對于程序主體，則按照骨架曲線和滯回曲線的加載、卸載及再加載的路徑模式進行編制，來體現(xiàn)土體剪切模量的變化。

2.2 加卸載轉折點的實現(xiàn)思路

對類似于巖土體的顆粒破碎材料，其一維與三維條件下的力學行為的差別是很大的。因此，在進行地下結構三維的非線性動力分析時，應當將前述的應力-應變關系從一維擴展至三維。

在程序編制過程中，如何確定加卸載轉折點是一個關鍵的問題，對于三維問題，至少存在6個應變率張量分量[19]，即：

(17)

式中：Δeij為應變增量張量。

(18)

(19)

式中，vi表示各分量上的應變增量，下標i表示1～6個分量。

(20)

因此，可以判斷，當d小于0時，應變發(fā)生轉向。

在實際情況中，當土體單元的廣義剪應力的絕對值大于土體的剪應力上限值τult時，土體的抗剪強度并不會完全喪失，其仍然具有一定的殘余強度，因此，編程時可假定各類土破壞時的剪切模量Gmin=0.05Gmax。

按照彈性力學理論，應力-應變關系為

(21)

式中，λ是拉密常數(shù)，可通過下式求得

(22)

式中：ν為泊松比。

為了對土體的黏性效應加以考慮，本文按Rayleigh阻尼的概念定義黏性阻尼矩陣為

[C]=α0[M]+α1[K]

(23)

式中：α0，α1分別為與質量矩陣和剛度矩陣有關的Rayleigh阻尼系數(shù)。

由于振型對Rayleigh阻尼矩陣是正交的，因此可得

(24)

式中：ξi為第i振型阻尼比，對于典型上海軟土，通常可取8%～10%；ωi為第i振型的自振頻率。

本文參考文獻[11]的研究成果，取第一振型作為計算Rayleigh尼系數(shù)的振型并假設阻尼矩陣只與剛度矩陣有關，則：

α0=0,α1=2ξi/ω1

(25)

因此，總應變速率產(chǎn)生的阻尼力為

(26)

式中：Del為初始彈性矩陣。

綜上所述，最終的應力-應變關系為

(k=1,2,3)

(27)

反應在雅克比矩陣?Δσ/?Δε中即為

(28)

式中：dt為積分時間步長，[I]為6×6的單位矩陣。

2.3 模型算法與程序編制

在ABAQUS軟件中編輯UMAT子程序來實現(xiàn)Davidenkov模型，其具體的算法步驟如下

(1) 在ABAQUS中輸入Davidenkov模型所需的土體初始的動力參數(shù)，包括初始剪切模量Gmax，泊松比ν，擬合參數(shù)A，B，γ0，剪應變上限γult，Rayleigh阻尼系數(shù)α1；

(2) 調用主程序，讀取當前步的應力、應變及應變增量；

(4) 計算d值，若d值小于0，則荷載轉向，此時記憶新的廣義剪應變γc；

(5) 計算土體在初始加載、卸載和再加載條件下的剪切模量G，為了保證計算的穩(wěn)定性，可取一剪切模量小值Gmin，如果G

(6) 計算拉密常數(shù)，更新雅克比矩陣；

(7) 更新內能及狀態(tài)變量。

該模型的計算流程圖如圖5所示，在Davidenkov本構模型子程序的編輯中，會用到大量的狀態(tài)變量數(shù)組，這些狀態(tài)變量數(shù)組可以用來存儲與求解過程有關的參變量，這些參變量會隨著求解過程而更新，如在判斷是否荷載轉向時，應變增量及新的γc即可用狀態(tài)變量保存調用，并可在結果文件中輸出，用戶可以清楚的知道荷載在何時發(fā)生轉向，極大的方便了用戶的編輯及查錯過程。

圖5 子程序計算流程圖

2.4 程序驗證

為驗證子程序的可靠性，本文首先利用一個單元進行驗證計算，因為一個單元的應力應變狀態(tài)簡單，能夠更直觀的反應子程序是否能夠實現(xiàn)剛度折減及荷載轉向等功能。在ABAQUS中建立一個1 m×1 m×1 m的單元，通過施加位移荷載的形式來研究復雜加載條件下的應力-應變變化情況，荷載路徑采用如下形式，路徑圖如圖6所示。

y=2e-3×(sin(2x)+sin(4x))

(29)

圖6 荷載路徑圖

圖7是單元在荷載作用下剛度折減情況，可以看到，隨著荷載的增加，土體單元的剛度會隨之衰減，當荷載轉向時，土體的剛度也隨之發(fā)生變化，當土體剛度小于初始剛度的5%時，土體剛度為一定值(Gmin)不再降低。

圖7 剛度折減曲線

圖8是單元的應力-應變關系曲線，可以看到，土體單元的應力-應變關系存在明顯的滯回現(xiàn)象，滯回圈飽滿。圖9是通過狀態(tài)變量記錄到的荷載轉向的次數(shù)，將其與輸入荷載對比，可以看出，本文的子程序能正確對土體單元的加卸載情況進行判斷。

圖8 應力-應變關系曲線

圖9 荷載轉向次數(shù)

3 軟土計算結果對比

本文選取了典型上海軟土場地進行二維地震反應有限元分析，并將分析結果與等效線性化結果對比。典型上海軟土場地各層土的物理力學參數(shù)見表5所示，各層土的基于Davidenkov模型的參數(shù)擬合結果詳見表4。

表5 典型上海軟土場地物理力學參數(shù)

本文利用作者文獻[20]中得到的基巖地震波作為輸入，進行二維軟土場地地震反應分析，并將非線性動力本構模型的計算結果和等效線性化結果進行對比。選用基巖地震波的加速度時程曲線及傅里葉譜分別如圖10和圖11所示。

圖10 基巖地震波時程曲線

圖11 基巖地震波傅里葉譜

圖12和圖13分別是經(jīng)過等效線性化算法和本文子程序計算得到的地表的加速度時程響應曲線和地表相對位移時程響應曲線。二者的響應峰值及相對差別列于表6中。

圖12 加速度響應對比

圖13 相對位移響應對比

Tab.6 Comparison between equivalent linearization results and subroutine results

對比類別地表加速度響應/gal地表相對位移響應/m等效線性化結果94.460.044本文子程序結果83.340.039相對差13.34%12.82%

由表6可知，采用等效線性化和本文子程序計算的地表加速度響應峰值分別是94.46 gal和81.34 gal，二者相對差為13.34%；采用等效線性化和本文子程序計算的地表相對位移響應峰值分別是0.044 m和0.039 m，二者相對差為12.82%。這主要是由于等效線性化方法會產(chǎn)生“虛(擬)共振效應”,過高地估計地震反應。

4 結 論

通過本文的研究，可以得到以下結論：

(1) 本文推導了Davidenkov模型的非線性動力本構模型在初始加載、卸載和再加載階段的增量剪切模量表達式，并將土體動應力-應變關系推廣到三維。基于粒子群優(yōu)化算法對Davidenkov模型的參數(shù)進行擬合，并通過典型砂土和黏土驗證了擬合效果。結果表明，兩種土體擬合結果的互相關度量值NCC均接近于1，粒子群優(yōu)化算法對Davidenkov模型參數(shù)具有很好的擬合效果。

(2) 利用ABAQUS軟件提供的二次開發(fā)平臺UMAT，本文編制了基于Davidenkov模型的土體非線性動力本構模型計算子程序，并通過復雜加載路徑驗證了該子程序的正確性，實現(xiàn)了在ABAQUS軟件中軟土非線性動力本構模型的二次開發(fā)。

(3) 以某典型上海軟土場地為例，文中對比分析了上海軟土場地基于非線性動力本構模型的計算結果和等效線性化的結果。結果表明，上海軟土場地基于非線性動力本構模型的加速度響應值和相對位移響應值要小于等效線性化的結果。