魯麗雪，遲世春

（大連理工大學(xué) 海岸與近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 116024）

1 引 言

土體的動(dòng)力模型即動(dòng)荷載下土體的應(yīng)力-應(yīng)變關(guān)系是表征土體動(dòng)力學(xué)特性的基本關(guān)系，是土動(dòng)力學(xué)研究的中心問題之一，也是分析土動(dòng)力反應(yīng)等一系列工程問題的重要基礎(chǔ)。土體的動(dòng)力特性主要表現(xiàn)出非線性、滯后性、變形累積以及強(qiáng)度與剛度退化特性。到目前為止，所提出土體動(dòng)力模型大多不能令人滿意。Bouc-Wen土體動(dòng)力微分模型由Bouc[1]1967 年提出，Wen[2]和Baber等學(xué)者[3－4]對(duì)其進(jìn)行了改進(jìn)，使之能夠模擬土體的強(qiáng)度和剛度退化特性。Bouc-Wen模型既包含非線性阻尼，又包含非線性剛度，可通過合理選擇參數(shù)得到不同形狀的滯回圈，用以描述土體在動(dòng)力荷載作用下非線性、滯后性、變形積累以及強(qiáng)度與剛度退化等特性。Bouc-Wen土體動(dòng)力模型也能很好地反映剛度和阻尼隨剪應(yīng)變的變化，適于描述多種類型土的循環(huán)動(dòng)力特性。該模型還無(wú)需處理加載拐點(diǎn)，方便編程，具有良好的通用性。

圖1為GDS動(dòng)三軸試驗(yàn)得到的摻礫土心墻料滯回圈與Bouc-Wen模型滯回曲線的比較，可以看出，試驗(yàn)點(diǎn)與理論模型吻合良好。

國(guó)外學(xué)者已就Bouc-Wen模型的熱力學(xué)可容許性開展大量研究工作，指出盡管Bouc-Wen模型獨(dú)立于內(nèi)時(shí)理論發(fā)展起來(lái)，但其仍屬于基于內(nèi)時(shí)概念的內(nèi)部變量理論類模型[5－6]。由于某些應(yīng)力路徑引起負(fù)的能量耗散而可能違反熱力學(xué)第二定律[7]。不過，單個(gè)滯回圈自動(dòng)閉合而不發(fā)生應(yīng)力或應(yīng)變滑移，因此服從Drucker（或 Ilyushin）公設(shè)[8－9]。

Collins等[10]基于熱力學(xué)基本原理建立了土體本構(gòu)關(guān)系模型的方法，采用量綱分析，給出耗散函數(shù)的可能形式，對(duì)受力變形過程中的能量耗散機(jī)制等問題也進(jìn)行了初步探討。郭曉霞[11]、遲世春[12]等從熱力學(xué)基本定律出發(fā)，對(duì)土體Hardin-Drnevich模型及 Ramberg-Osgood模型的增量耗散函數(shù)進(jìn)行了構(gòu)造，但由于采用了Masing準(zhǔn)則假設(shè)，并未考慮實(shí)際阻尼比隨剪應(yīng)變的變化等。因此，有必要結(jié)合壩料的實(shí)際阻尼研究土石料動(dòng)力變形機(jī)制。

本文以Bouc-Wen模型為研究對(duì)象，探討其適用條件以及參數(shù)影響，從熱力學(xué)基本定律出發(fā)，構(gòu)造壩料的自由能函數(shù)和耗散函數(shù)，研究其耗散特征及動(dòng)力變形機(jī)制，繪制其耗散應(yīng)力空間和真實(shí)應(yīng)力空間屈服曲線，探討其屈服面的變化規(guī)律。

圖1 理論滯回曲線與試驗(yàn)結(jié)果比較Fig.1 Comparison between theoretically predicted hysteretic loops and test results

2 Bouc-Wen模型耗散函數(shù)的構(gòu)造

2.1 Bouc-Wen土體動(dòng)力模型

Bouc-Wen土的動(dòng)力模型的應(yīng)力-應(yīng)變關(guān)系為

式中：τ、γ分別為切應(yīng)力和切應(yīng)變；Gmax為最大剪切模量；α為屈服比，控制屈服后剪切剛度，取值范圍(0,1)，主要影響滯回圈的面積。α越大，滯回圈越瘦，α對(duì)曲線的影響見圖2。

圖2 參數(shù)α 對(duì)滯回曲線的影響Fig.2 Influence of parameter α on hysteretic loop

參數(shù)ζ=ζ(t)是遲滯變量，用來(lái)控制土的非線性，由下列微分方程確定：

式中：A、g、b、n均為模型形狀參數(shù)。

由式（1）、（2）可知，Bouc-Wen模型中加、卸載拐點(diǎn)由遲滯變量ζ和應(yīng)變速率d/dtγ控制，無(wú)需特殊處理，這樣可大大簡(jiǎn)化程序的編寫。

A控制0ζ=時(shí)滯回圈的斜率，影響滯回圈的高度和面積。描述土的滯回特性時(shí)，A取值(0,1]，其單獨(dú)變化時(shí)對(duì)滯回圈的影響見圖3。

圖3 參數(shù)A對(duì)滯回曲線的影響Fig.3 Influence of parameter A on hysteretic loop

g控制dγ/dt變號(hào)時(shí)曲線的剛度。為確保能量耗散值為正，應(yīng)滿足 g＞0，在此基礎(chǔ)上具有物理意義的5種滯回圈分別是：①g+b＞0和g?b＞0（見圖4）；②g+b＞0和g?b＜0（見圖5）；③ g+b＞0和g?b=0（見圖6）；④g+b = 0 和g?b＞0（見圖7）；⑤g+b＜0和g?b＞0（見圖8、9）。圖2～9中 n = 1，Gmax=600000 MPa。

由圖可知，當(dāng) g+b＜0 時(shí)滯回圈內(nèi)凹，g+b＞0時(shí)滯回圈外凸，g+b = 0 時(shí)滯回圈呈擬線性。

圖4 參數(shù)g對(duì)滯回曲線的影響Fig.4 Influence of parameter g on hysteretic loop

圖5 參數(shù)b對(duì)滯回曲線的影響Fig.5 Influence of parameter b on hysteretic loop

圖6 參數(shù)g+b對(duì)滯回曲線的影響(g-b=0,A=0.40)Fig.6 Influence of parameter g+b on hysteretic loop(g-b=0,A=0.40)

圖7 參數(shù)g-b對(duì)滯回曲線的影響(g+b=0,A=0.50)Fig.7 Influence of parameter g-b on hysteretic loop(g+b=0,A=0.50)

圖8 參數(shù)g+b對(duì)滯回曲線的影響(g-b=40000,A=0.15)Fig.8 Influence of parameter g+b on hysteretic loop(g-b=40000,A=0.15)

圖9 參數(shù)g-b對(duì)滯回曲線的影響(g+b=20000,A=0.10)Fig.9 Influence of parameter g-b on hysteretic loop(g+b=20000,A=0.10)

n為屈服指數(shù)，控制加載拐點(diǎn)初始加載階段線性至非線性轉(zhuǎn)換時(shí)的屈服前后滯回圈光滑程度，即屈服尖銳程度系數(shù)，n越小，過渡越平滑，其取值范圍為0至無(wú)窮，n＞10時(shí)，近似理想彈塑性；n趨向無(wú)窮時(shí)，應(yīng)力-應(yīng)變曲線趨向雙線性。n對(duì)單調(diào)加載曲線的影響見圖10。

圖10 參數(shù)n對(duì)滯回曲線的影響Fig.10 Influence of parameter n on hysteretic loop

2.2 模型耗散函數(shù)的構(gòu)造

為確保模型不違背熱力學(xué)定律，本文基于熱力學(xué)基本定律構(gòu)造了Bouc-Wen土體動(dòng)力模型的增量耗散函數(shù)表達(dá)式。

假定塑性中心移動(dòng)沿骨架曲線建立耗散函數(shù)，采用Ziegler正交條件，得到耗散應(yīng)力空間中屈服函數(shù)的表達(dá)式。引入耗散應(yīng)力及真實(shí)應(yīng)力間的差別項(xiàng)，即遷移應(yīng)力，得到真實(shí)空間中的屈服函數(shù)。由于Bouc-Wen模型的表達(dá)式是微分形式，為簡(jiǎn)化表達(dá)式，取n=1構(gòu)造耗散函數(shù)。

滯回曲線統(tǒng)一表達(dá)為

骨架曲線為

式（3）與第1象限中骨架曲線走向相同的再加荷曲線g、γ1取負(fù)值，與第3象限中骨架曲線走向相同的卸荷曲線b取負(fù)值，與第3象限中骨架曲線走向相同的再加荷曲線g、γ1、b取負(fù)值；式（4）對(duì)于第3象限中骨架曲線g、b取負(fù)值，γ0為剪應(yīng)變幅值，γ1為滯回曲線與橫軸的交點(diǎn)，可由反對(duì)稱性得出，下同。

綜合公式（3）、（4）可得到塑性中心移動(dòng)為即時(shí)骨架曲線時(shí)的耗散函數(shù)：

真實(shí)應(yīng)力表達(dá)式為

式中：xbl、ρbl分別為耗散應(yīng)力和遷移應(yīng)力，遷移應(yīng)力，其公式為

3 屈服曲線的繪制

結(jié)合某心墻壩工程選擇粗堆石料、細(xì)堆石料及反濾料Ⅰ的動(dòng)力試驗(yàn)成果，使用上述方法探討其屈服面的變化規(guī)律。粗堆石料為角礫巖與花崗巖的混合物，細(xì)堆石料為弱風(fēng)化以下花崗巖開挖料加工而成，二者采用大三軸動(dòng)力試驗(yàn)機(jī)進(jìn)行動(dòng)力試驗(yàn)；反濾料I則為弱風(fēng)化以下花崗巖加工而成，采用共振柱試驗(yàn)機(jī)及中型三軸儀進(jìn)行動(dòng)力試驗(yàn)。

3.1 最大動(dòng)剪切模量

最大動(dòng)剪切模量[13]由下式確定：

式中：Gmax為最大動(dòng)剪切模量，σm為土體所受的初始平均靜應(yīng)力；pa為大氣壓力；k（與Gmax同量綱）、n1（無(wú)量綱）為試驗(yàn)常數(shù)，均由動(dòng)力試驗(yàn)確定，筑壩粗堆石料、細(xì)堆石料及反濾料Ⅰ的模量系數(shù)k分別為2455.6、1303.9、976.0，指數(shù) n1分別為0.6004、0.6109、0.4980。

壩料動(dòng)剪模量比隨動(dòng)剪應(yīng)變和阻尼比隨動(dòng)剪應(yīng)變的變化見圖11～13。

圖12 細(xì)堆石料對(duì)應(yīng)的模量衰減和阻尼增長(zhǎng)曲線Fig.12 Corresponding modulus decline and damping growth curves of fine rockfill materials

圖13 反濾料Ⅰ對(duì)應(yīng)的模量衰減和阻尼增長(zhǎng)曲線Fig.13 Corresponding modulus decline and damping growth curves of the filter materialⅠ

3.2 Bouc-Wen模型參數(shù)辨識(shí)

根據(jù)動(dòng)力試驗(yàn)提供的動(dòng)剪切模量 Gs、等效阻尼比ξeq與動(dòng)剪應(yīng)變[14]的關(guān)系（見圖11～13），用遺傳算法對(duì)Bouc-Wen模型進(jìn)行參數(shù)辨識(shí)。

遺傳算法中選擇算法采用比例選擇方法，交叉采用非均勻算數(shù)交叉的方法，變異采用均勻變異，種群數(shù)量為80，迭代次數(shù)為1000，雜交概率為0.8，變異概率為0.1。由于參數(shù)n僅影響遲滯曲線的光滑度，取 n = 1，因此需辨識(shí)參數(shù)有α、A、g、b。

從模型參數(shù)識(shí)別的結(jié)果來(lái)看，隨著動(dòng)力水平的提高，阻尼增加，滯回圈變高、變胖，α、A、g?b、g+b 均逐漸降低。圖14給出半對(duì)數(shù)坐標(biāo)下粗堆石料對(duì)應(yīng)的模型參數(shù)與動(dòng)剪應(yīng)變的變化規(guī)律。

圖14 參數(shù)α、A、g-b、g+b隨剪應(yīng)變的變化規(guī)律Fig.14 Variation trend of parameters α，A，g-b，g+b with the change of shear strain

在半對(duì)數(shù)坐標(biāo)下擬合α-γ、A-γ、gb?-γ及g b+-γ曲線，得

3.3 屈服曲線的繪制

以粗堆石料為例，結(jié)合模量衰減和阻尼比增長(zhǎng)建立Bouc-Wen模型，討論p-q（單位為MPa）的應(yīng)力空間中塑性中心的移動(dòng)為骨架曲線時(shí)的屈服曲線的變化規(guī)律。隨著γm不斷增大，繪制屈服曲線如圖15～17所示，其中（γm，τm）為滯回曲線的頂點(diǎn)。

圖15 γm= 5×10-6時(shí)的屈服曲線Fig.15 Yield curves when γm= 5×10-6

圖16 γm= 4×10-5時(shí)的屈服曲線Fig.16 Yield curves when γm= 4×10-5

圖17 γm= 8×10-4時(shí)的屈服曲線Fig.17 Yield curves when γm= 8×10-4

圖中分別繪制了γm=5× 10?6、4 × 10?5和8×10?4時(shí)試樣的屈服曲線。由圖可見，應(yīng)變水平為5× 10?6時(shí)，不同γ/γm的屈服函數(shù)幾乎為一條直線，滯回圈應(yīng)變達(dá)到 4 × 10?5時(shí)，屈服面為斜率不同的直線。不同γ/γm的屈服面在p-q應(yīng)力空間為同一直線說明土體屈服時(shí)土體顆粒之間的常摩擦系數(shù)起控制作用。因此，試樣滯回圈應(yīng)變?cè)?5 × 10?6～4 × 10?5之間存在一個(gè)剪切屈服曲線的直線斜率，不再保持一致的界限應(yīng)變，稱之為第1閾值應(yīng)變[15]，記為γt1。隨著應(yīng)變的增大，且達(dá)到 8 × 10?4時(shí)，屈服面已明顯變彎。剪切屈服面彎曲代表控制堆石體試樣塑性變形的機(jī)制除摩擦外，還有結(jié)構(gòu)變化等機(jī)制。因此，試樣滯回圈應(yīng)變?cè)?4 × 10?5～8 × 10?4之間存在剪切屈服曲線由直變彎的界限應(yīng)變，稱之為第2閾值應(yīng)變，記為γt2[15]。

求出 3種壩料的閾值應(yīng)變，并將其繪制于圖11～13，得到第1、第2壩料閾值應(yīng)變?cè)贐ouc-Wen模型對(duì)應(yīng)的模量衰減及阻尼衰減情況。從粗堆石料及細(xì)堆石料的第1閾值應(yīng)變所對(duì)應(yīng)的歸一化動(dòng)剪切模量看，第1閾值應(yīng)變對(duì)應(yīng)的 G/Gmax在0.93～0.99之間，阻尼比為1%～4%。第2閾值應(yīng)變對(duì)應(yīng)的割線模量與最大動(dòng)剪切模量之比范圍為 0.54～0.79，與傳統(tǒng)意義上以孔壓開始升高或體積開始變化為標(biāo)準(zhǔn)定義的門檻應(yīng)變相當(dāng)。

4 結(jié) 語(yǔ)

本文探討了Bouc-Wen土體動(dòng)力微分模型的適用條件以及參數(shù)影響。從熱力學(xué)基本定律出發(fā)，克服Hardin-Drnevich模型及Ramberg-Osgood模型未考慮實(shí)際阻尼比隨剪應(yīng)變變化的局限，結(jié)合模量衰減和阻尼增長(zhǎng)給出基于Bouc-Wen微分模型的增量耗散函數(shù)，研究了筑壩土石料的動(dòng)力耗散特征及動(dòng)力變形機(jī)制，得到筑壩土石料的閾值應(yīng)變。其中第二閾值應(yīng)變與傳統(tǒng)試驗(yàn)方法得到的門檻應(yīng)變相當(dāng)。

[1]BOUC R.Forced vibration of mechanical systems with hysteresis[C]//Proceedings of the Conference on Nonlinear Oscillators.Prague: [s.n.],1967: 315.

[2]WEN Y K.Method for random vibration of hysteretic systems[J].Journal of Engineering Mechanics Division,ASCE,1976.102: 249－263.

[3]BABER T T,WEN Y K.Random vibration of hysteretic systems[J].Journal of Engineering Mechanics Division,ASCE,1981,107(6): 1069－1086.

[4]BABER T T,NOORI M N.Random vibration of degrading,pinching systems[J].Journal of Engineering Mechanics,ASCE,1985,111(8): 1011－1027.

[5]CASCIATI F.Stochastic dynamics of hysteretic media[J].Structural Safety,1989,6(2－4): 259－269.

[6]SIVASELVAN M V,REINHORN A M.Hysteretic models for deteriorating inelastic structures[J].Journal of Engineering Mechanics,ASCE,2000,126(6): 633－640.

[7]CHARALAMPAKIS A E,KOUMOUSIS V K.A Bouc-Wen model compatible with plasticity postulates[J].Journal of Sound and Vibration,2009,322(4－5): 954－968.

[8]CAPECCHI D,DE FELICE G.Hysteretic systems with internal variables[J].Journal of Engineering Mechanics,2001,127(9): 891－898.

[9]ERLICHER S,BURSI O S.Bouc-Wen-Type models with stiffness degradation: Thermodynamic analysis and applications[J].Journal of Engineering Mechanics,2008,134(10): 843－855.

[10]COLLINS I F,HOULSBY G T.Application of thermomechanical principles to the modeling of geomaterials[C]//Proceedings of the Royal Society of London,Series A: Mathematical,Physical and Engineering Sciences.London: [s.n.],1997: 1975－2001.

[11]遲世春,郭曉霞,楊峻,等.土的動(dòng)力 Hardin-Drnevich模型小應(yīng)變特性及其閾值應(yīng)變研究[J].巖土工程學(xué)報(bào),2008,30(2): 243－249.CHI Shi-chun,GUO Xiao-xia,YANG jun,et al.Small strain characteristics and threshold strain of dynamic Hardin-Drnevich model for soils[J].Rock and Soil Mechanics,2008,30(2): 243－249.

[12]郭曉霞,遲世春,林皋.基于熱力學(xué)定律的土體動(dòng)力Hardin-Drnevich模型再認(rèn)識(shí)[J].巖土力學(xué),2008,29(9):2335－2340.GUO Xiao-xia,CHI Shi-chun,LIN Gao.Recognition of dynamic Hardin-Draevich model for soils based on generalized thermodynamics[J].Rock and Soil Mechanics,2008,29(9): 2335－2340.

[13]謝定義.土動(dòng)力學(xué)[M].西安: 西安交通大學(xué)出版社,1988.

[14]王志良,韓清宇.黏彈塑性土層地震反應(yīng)的波動(dòng)分析法[J].地震工程與工程振動(dòng),1981,1(1): 117－136.WANG Zhi-liang,HAN Qing-yu.Sticky elastoplastic seismic response of soil layer wave analysis method[J].Earthquake Engineering and Engineering Vibration,1981,1(1): 117－136.

[15]郭曉霞.熱力學(xué)方法在土體本構(gòu)模型中的應(yīng)用研究[D].大連: 大連理工大學(xué),2009.