王沿朝， 陳清軍

(同濟(jì)大學(xué) 土木工程防災(zāi)國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海 200092)

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

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

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

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

1.1 Davidenkov骨架曲線理論公式

Hardin等[9]提出，土體的應(yīng)力-應(yīng)變關(guān)系可以通過(guò)一種雙曲線形式來(lái)描述，其動(dòng)剪切模量比的表達(dá)式如下

G/Gmax=1-H(γ)

(1)

(2)

式中：G為剪切模量；Gmax為最大剪切模量；γ為剪應(yīng)變幅；γref為參考剪應(yīng)變，一般可取動(dòng)剪切模量比為0.5時(shí)所對(duì)應(yīng)的剪應(yīng)變幅。

Martin等[10]采用Davidenkov骨架曲線來(lái)描述上述關(guān)系，將式(2)改寫(xiě)如下

(3)

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

因此，Davidenkov模型的應(yīng)力-應(yīng)變關(guān)系的骨架曲線可用下式表示

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

(4)

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

(5)

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

(6)

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

(7)

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

(8)

式中：τc和γc分別表示的是應(yīng)力-應(yīng)變滯回曲線中加卸載轉(zhuǎn)折點(diǎn)處的剪應(yīng)力與剪應(yīng)變幅值。

1.2 剪切模量的增量表達(dá)式

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

(9)

對(duì)上式求γ的導(dǎo)數(shù)，即可得到

(10)

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

G=GmaxMt

(11)

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

(|τ|≤τult)

(12)

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

(|τ|≤τult)

(13)

而當(dāng)|τ|>τult時(shí)，G=0。

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

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

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

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

(14)

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

(15)

式中:w為慣性權(quán)值；c1，c2為學(xué)習(xí)因子；c3為另一個(gè)慣性權(quán)值；rand1，rand2為兩個(gè)[0，1]之間的隨機(jī)數(shù)。

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

(16)

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

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

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

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

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

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

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

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

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

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

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

土類剪應(yīng)變×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所示，從圖中可以看出，對(duì)試驗(yàn)數(shù)據(jù)擬合的曲線走勢(shì)很好，Davidenkov模型參數(shù)擬合的結(jié)果及NCC值見(jiàn)表2所示，可以看出，兩種土體擬合結(jié)果的NCC值均十分接近于1，說(shuō)明了擬合效果很好。

圖2 砂土擬合曲線

圖3 黏土擬合曲線

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

在此基礎(chǔ)上，本文通過(guò)粒子群優(yōu)化算法對(duì)典型上海軟土進(jìn)行Davidenkov模型參數(shù)擬合，表3給出了典型上海軟土各層土的剪切模量比及阻尼比與剪應(yīng)變幅的試驗(yàn)數(shù)據(jù)。對(duì)模型參數(shù)進(jìn)行擬合時(shí)，其他擬合參數(shù)不變，粒子維數(shù)取為3。典型上海軟土各層土的擬合曲線如圖4所示，Davidenkov模型參數(shù)擬合的結(jié)果及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 典型上海軟土擬合結(jié)果

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

2.1 實(shí)現(xiàn)過(guò)程

ABAQUS為用戶提供了功能強(qiáng)大的用戶子程序接口，以幫助用戶開(kāi)發(fā)基于ABAQUS內(nèi)核的程序，而這其中應(yīng)用最為廣泛的子程序即為UMAT[18]，它主要用于用戶開(kāi)發(fā)自己的本構(gòu)模型，以彌補(bǔ)ABAQUS自帶材料模型的不足。ABAQUS子程序通過(guò)FORTRAN語(yǔ)言開(kāi)發(fā)。按照FORTRAN的語(yǔ)法，用戶可以自由的編輯代碼，并形成一個(gè)獨(dú)立的程序單元，該單元能被獨(dú)立的編譯和存儲(chǔ)，也可以被其他單元引用，利用UMAT能夠?qū)⒋罅康臄?shù)據(jù)帶回供引用程序使用。

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

(1) 子程序定義語(yǔ)句；

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

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

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

(5) 子程序返回和結(jié)束語(yǔ)句。

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

2.2 加卸載轉(zhuǎn)折點(diǎn)的實(shí)現(xiàn)思路

對(duì)類似于巖土體的顆粒破碎材料，其一維與三維條件下的力學(xué)行為的差別是很大的。因此，在進(jìn)行地下結(jié)構(gòu)三維的非線性動(dòng)力分析時(shí)，應(yīng)當(dāng)將前述的應(yīng)力-應(yīng)變關(guān)系從一維擴(kuò)展至三維。

在程序編制過(guò)程中，如何確定加卸載轉(zhuǎn)折點(diǎn)是一個(gè)關(guān)鍵的問(wèn)題，對(duì)于三維問(wèn)題，至少存在6個(gè)應(yīng)變率張量分量[19]，即：

(17)

式中：Δeij為應(yīng)變?cè)隽繌埩俊?/p>

(18)

(19)

式中，vi表示各分量上的應(yīng)變?cè)隽?，下?biāo)i表示1～6個(gè)分量。

(20)

因此，可以判斷，當(dāng)d小于0時(shí)，應(yīng)變發(fā)生轉(zhuǎn)向。

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

按照彈性力學(xué)理論，應(yīng)力-應(yīng)變關(guān)系為

(21)

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

(22)

式中：ν為泊松比。

為了對(duì)土體的黏性效應(yīng)加以考慮，本文按Rayleigh阻尼的概念定義黏性阻尼矩陣為

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

(23)

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

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

(24)

式中：ξi為第i振型阻尼比，對(duì)于典型上海軟土，通?？扇?%～10%；ωi為第i振型的自振頻率。

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

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

(25)

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

(26)

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

綜上所述，最終的應(yīng)力-應(yīng)變關(guān)系為

(k=1,2,3)

(27)

反應(yīng)在雅克比矩陣?Δσ/?Δε中即為

(28)

式中：dt為積分時(shí)間步長(zhǎng)，[I]為6×6的單位矩陣。

2.3 模型算法與程序編制

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

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

(2) 調(diào)用主程序，讀取當(dāng)前步的應(yīng)力、應(yīng)變及應(yīng)變?cè)隽浚?/p>

(4) 計(jì)算d值，若d值小于0，則荷載轉(zhuǎn)向，此時(shí)記憶新的廣義剪應(yīng)變?chǔ)胏；

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

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

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

該模型的計(jì)算流程圖如圖5所示，在Davidenkov本構(gòu)模型子程序的編輯中，會(huì)用到大量的狀態(tài)變量數(shù)組，這些狀態(tài)變量數(shù)組可以用來(lái)存儲(chǔ)與求解過(guò)程有關(guān)的參變量，這些參變量會(huì)隨著求解過(guò)程而更新，如在判斷是否荷載轉(zhuǎn)向時(shí)，應(yīng)變?cè)隽考靶碌摩胏即可用狀態(tài)變量保存調(diào)用，并可在結(jié)果文件中輸出，用戶可以清楚的知道荷載在何時(shí)發(fā)生轉(zhuǎn)向，極大的方便了用戶的編輯及查錯(cuò)過(guò)程。

圖5 子程序計(jì)算流程圖

2.4 程序驗(yàn)證

為驗(yàn)證子程序的可靠性，本文首先利用一個(gè)單元進(jìn)行驗(yàn)證計(jì)算，因?yàn)橐粋€(gè)單元的應(yīng)力應(yīng)變狀態(tài)簡(jiǎn)單，能夠更直觀的反應(yīng)子程序是否能夠?qū)崿F(xiàn)剛度折減及荷載轉(zhuǎn)向等功能。在ABAQUS中建立一個(gè)1 m×1 m×1 m的單元，通過(guò)施加位移荷載的形式來(lái)研究復(fù)雜加載條件下的應(yīng)力-應(yīng)變變化情況，荷載路徑采用如下形式，路徑圖如圖6所示。

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

(29)

圖6 荷載路徑圖

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

圖7 剛度折減曲線

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

圖8 應(yīng)力-應(yīng)變關(guān)系曲線

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

3 軟土計(jì)算結(jié)果對(duì)比

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

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

本文利用作者文獻(xiàn)[20]中得到的基巖地震波作為輸入，進(jìn)行二維軟土場(chǎng)地地震反應(yīng)分析，并將非線性動(dòng)力本構(gòu)模型的計(jì)算結(jié)果和等效線性化結(jié)果進(jìn)行對(duì)比。選用基巖地震波的加速度時(shí)程曲線及傅里葉譜分別如圖10和圖11所示。

圖10 基巖地震波時(shí)程曲線

圖11 基巖地震波傅里葉譜

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

圖12 加速度響應(yīng)對(duì)比

圖13 相對(duì)位移響應(yīng)對(duì)比

Tab.6 Comparison between equivalent linearization results and subroutine results

對(duì)比類別地表加速度響應(yīng)/gal地表相對(duì)位移響應(yīng)/m等效線性化結(jié)果94.460.044本文子程序結(jié)果83.340.039相對(duì)差13.34%12.82%

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

4 結(jié) 論

通過(guò)本文的研究，可以得到以下結(jié)論：

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

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

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