本構(gòu)模型對面板堆石壩應(yīng)力變形的影響

張 飛， 費文平

(1.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室，四川 成都 610065；2.四川大學(xué)水利水電學(xué)院，四川 成都610065)

本構(gòu)模型對面板堆石壩應(yīng)力變形的影響

張 飛1，2， 費文平1，2

(1.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室，四川 成都 610065；2.四川大學(xué)水利水電學(xué)院，四川 成都610065)

鄧肯-張E-B模型，能考慮土體的非線性特征以及應(yīng)力歷史對變形的影響，比較符合土體的受力特征。理想彈塑性D-P模型是在摩爾-庫倫模型的基礎(chǔ)上考慮了靜力壓力的影響，廣泛應(yīng)用于巖體材料的模擬。堆石體是介于土體與巖體之間的一種土石混合體，其本構(gòu)模型的選擇一直存在較大的爭議?；贏PDL二次開發(fā)，將E-B本構(gòu)模型成功引入到ANSYS靜力分析中，與理想彈塑性D-P模型的計算成果進行對比分析，并對E-B模型的主要參數(shù)進行敏感性分析。研究成果表明：本構(gòu)模型對堆石體豎向位移影響較大，但對堆石體的應(yīng)力、面板應(yīng)力、及面板位移影響較??；E-B模型參數(shù)中K，n以及荷載子步對面板堆石壩應(yīng)力、位移的影響較大。

本構(gòu)模型；面板堆石壩；應(yīng)力變形特性

Abstract:Duncan-Zhang E-B model can consider the nonlinear characteristics of soil and the influence of stress history on deformation,which is in line with the characteristics of soil.The ideal elastoplastic D-P model is based on the Mohr-Coulomb model,which takes into account the effects of static pressure and is widely used in the simulation of rock mass materials.Rockfill is a mixture of soil and rock between the soil and rock,the constitutive model of choice has been a big controversy.Based on the secondary development of APDL,E-B constitutive model was successfully introduced into ANSYS static analysis,compared with the calculated results of ideal elastic-plastic D-P model,and the sensitivity of E-B model was analyzed.The results showthat the constitutive model has a great influence on the vertical displacement of the rockfill body,but it has little effect on the stress ofthe rockfill body,the stress ofthe panel and the displacement ofthe panel.The influence ofK,n and load sub-step on E-Bmodel parameters has great influence on the stress and displacement ofconcrete face rockfill dam.

Keywords:The constitutive model；panel rockfill damand stress and deformation characteristics

0 引言

混凝土面板堆石壩是以堆石為主要材料，以混凝土面板為防滲體及墊層和過渡層所組成的，在計算壩體應(yīng)力變形等方面時，構(gòu)成混凝土面板堆石壩的各部分材料的工程特性及所選用的本構(gòu)模型對數(shù)值分析的計算成果均有重要影響。目前用于面板堆石壩靜力計算的本構(gòu)模型大致分為兩類：其一是彈塑性模型，其二則是非線性彈性模型。其中非線性彈性模型的發(fā)展經(jīng)歷了漫長的歷程，主要包括早期的鄧肯-張的E-μ模型，隨著堆石料本構(gòu)性質(zhì)探索內(nèi)容的不斷深化，該模型也越來越滿足不了計算的要求；1970年鄧肯等人在E-μ模型的基礎(chǔ)上提出了E-B模型，該模型最核心的思想是把材料各種形式的變形均假定為彈性變形，經(jīng)過不斷改動材料的彈性參數(shù)來大致地考慮塑性變形以減小差異。由于計算中采用的回彈模量與加載模量不同，所得出的變形也不同。它的優(yōu)點是使用了增量算法，局部反映了變形應(yīng)力路徑的影響水平。彈塑性模型主要包括早期的Drucker-Prager模型，它是1952年由Drucker和Prager共同提出的，接近莫爾-庫倫屈服準(zhǔn)則，但卻是以最簡略的數(shù)學(xué)模式，在土力學(xué)塑性極限設(shè)計中得到很好的應(yīng)用，并被人們所熟知，這種彈塑性本構(gòu)關(guān)系在巖土力學(xué)與工程中應(yīng)用最廣。

ANSYS軟件是一個功能強大的大型通用有限元分析軟件，它包括金屬、混凝土等多種材料的本構(gòu)模型，但尚無土工數(shù)值分析中廣泛應(yīng)用的鄧肯-張模型，這一遺憾無疑使ANSYS軟件在土工有限元分析中的應(yīng)用范圍和適用性受到限制。本文通過APDL參數(shù)化設(shè)計語言對ANSYS進行二次開發(fā)，在ANSYS中加入了目前廣泛應(yīng)用于土體的鄧肯-張E-B模型，結(jié)合具體工程，用E-B模型和D-P模型進行結(jié)果對比分析。為堆石壩選取方便合適的本構(gòu)模型提供參考和依據(jù)。

1 基本理論

1.1 理想彈塑性模型

D-P模型是一種理想彈塑性模型，其表達式如下：

其中參數(shù)α和k可以用土的粘聚力C和內(nèi)摩擦角φ來表示：

1.2 鄧肯-張彈性靜力本構(gòu)模型

基于ANSYS的二次開發(fā)實現(xiàn)面板堆石壩的鄧肯-張E-B模型，該模型中彈性模量E計算見下式

式中：c為材料凝聚力；φ為材料內(nèi)摩擦角；Pa為標(biāo)準(zhǔn)大氣壓；σ1，σ3分別為單元的大主應(yīng)力、小主應(yīng)力；Rf為破壞比；K為彈性模量系數(shù)；n為彈性模量指數(shù)。修正后的切線體積模量如下：

式中，Kb是體積模量系數(shù)，m是一無量綱的量，取值范圍0～1.0，泊松比為

由于堆石體的強度包線有彎曲，呈非線形，故對摩爾-庫倫準(zhǔn)則中的φ角進行修正為：

關(guān)于卸載情況，E-B模型使用回彈模量Eur繼續(xù)計算，回彈模量的計算表達式如下：

式中：kur為卸載模量基數(shù)；nur為卸載模量指數(shù)，且通常情況下nur≈n。當(dāng) σ1與 σ3的差值小于歷史上的最大值(σ1-σ3)0，且 s低于歷史上最大值s0時用Eur，否則用Et。

2 ANSYS二次開發(fā)

APDL即ANSYS參數(shù)化設(shè)計語言（ANSYS Parametric Design Language)是一種類似FORTRAN的解釋性語言，可用來自動完成一些通用性強的任務(wù)，也可以用于參數(shù)化建模，為結(jié)構(gòu)分析提供了很多便利。

本文根據(jù)土石壩有限元的計算原理和分析方法，在ANSYS中利用APDL技術(shù)成功的實現(xiàn)了土石壩非線性有限元分析。具體步驟：

（1）利用ANSYS的“生死單元”功能，殺死全部壩體單元，激活當(dāng)前填筑層單元。

（2）當(dāng)前填筑層單元施加重力加速度荷載9.81/nstep（nstep為荷載子步數(shù)）進行計算，調(diào)用宏命令來提取第一及第三主應(yīng)力，計算每一個單元的彈性常數(shù)Et、μ、Bt等，并修改單元的材料參數(shù)。

（3）依次施加當(dāng)前層的重力加速度荷載2*9.81/nstep、3*9.81/nstep，…，9.81，仿照第（2）步求解并逐步調(diào)整單元的材料參數(shù)。

（4）利用ANSYS重啟動技術(shù)，重新讀取參數(shù)，恢復(fù)上級荷載時的單元生死狀態(tài)，再激活新填層單元，重復(fù)（2）～（3）步，利用命令PARSAV保存計算參數(shù)，直至壩體填筑完成。

（5）編寫后處理的命令流，繪制應(yīng)力等值線和變形等值線等。

其中，采用APDL語言編寫的E-B 模型的宏命令如下：

*create,E-B.mac !創(chuàng)建E-B計算的宏命令

*afun,deg !角度單位用度，不是弧度

*set,Pa,1e5

*set,sigm1,- S1(i) !S1為第一主應(yīng)力，拉負(fù)壓正

*set,sigm3,-S3(i) !S3為第三主應(yīng)力，拉負(fù)壓正

sigm3=max(sigm3,0.1*Pa)

St=2*(c*cos(Fai)+sigm3*sin(Fai))/(1-sin(Fai))

S=(sigm1-sigm3)/St!應(yīng)力水平

S=min(S,0.95)!應(yīng)力水平最大取值

*if,St_m(i),gt,sigm1-sigm3,and,S_m(i),gt,S,then

Et=Kur*Pa*(sigm3/Pa)**n !卸荷模量

*else

Ei=k*Pa*(sigm3/Pa)**n

Et=Ei*(1-Rf*S)**2 !加荷情況的切線模量

S_m (i)=max(S_m (i),S)!保存最大應(yīng)力水平

St_m(i)=max(St_m(i),sigm1-sigm3)!保存最大應(yīng)力

*endif

Bt=Kb*Pa*(sigm3/Pa)**m

Mu=(3*Bt-Et)/(6*Bt)

Fai=Fail0-dFail*log10(p/Pa)

Mu=min(Mu,0.49)

mp,ex,i,Et!修改單元i的彈性模量

mp,PRXY,i,Mu !修改單元i的泊松比

mp,dens,i,idens(emat(i))

!emat(i)用于儲層單元的初始材料號

mpchg,i,i!修改i單元的材料屬性，

*end !結(jié)束創(chuàng)建宏文件

3 工程實例分析

以某混凝土面板堆石壩為例，面板采用D-P模型，堆石體分別采用E-B模型以及D-P模型來模擬，進而探究不同的本構(gòu)模型對堆石體和面板應(yīng)力變形的影響，從而可以選擇相對更加合理的本構(gòu)模型，以便于對本工程數(shù)值模擬提供指導(dǎo)。

3.1 工程概況

該混凝土面板堆石壩最大壩高142 m，上游壩坡1∶1.4，上游蓋重坡度1∶2.5，下游壩坡1∶1.6，設(shè)三級馬道，馬道高程分別為3193.00 m，3223.00 m，3253 m。由蓋重區(qū)、墊層區(qū)、過渡區(qū)、主堆石區(qū)、下游堆石區(qū)、排水區(qū)、面板、下游護坡等部分組成。

3.2 有限元模型

采用二維實體單元，進行堆石壩的二維有限元靜力分析。圖1為堆石壩結(jié)構(gòu)的二維有限元計算網(wǎng)格圖，共797個單元，834個節(jié)點。

圖1 面板堆石壩的二維有限元計算圖

3.3 計算參數(shù)與工況

暫不考慮構(gòu)造地應(yīng)力的影響，模擬堆石壩在完建工況下的位移場及應(yīng)力場，分析其分布規(guī)律。正常蓄水位高程為3278 m。

堆石料分別用E-B構(gòu)模型和D-P模型來模擬。根據(jù)設(shè)計單位提供的實驗數(shù)據(jù)，E-B和D-P模型材料參數(shù)分別見表1和表2。

表1 E-B模型參數(shù)表

表2 D-P模型參數(shù)表

壩體計算工況僅計算完建期。計算采用生死單元模擬大壩實際填筑施工過程，壩體填筑荷載分15步施加，對比分析兩種本構(gòu)模型計算結(jié)果。

3.4 計算成果及對比

本次選擇大壩壩體最大橫斷面為典型剖面進行計算分析。兩種模型分析大壩壩體在完建期的應(yīng)力位移情況見圖2～圖5，壩體及面板應(yīng)力位移最大值如表3。以圖2～圖5中，應(yīng)力拉為正，壓為負(fù)，單位為Pa，堆石體水平位移指向下游為正，指向上游為負(fù)，鉛直沉降向上為正，向下為負(fù)，面板法向位移指向壩體為正，偏離壩體為負(fù)，位移單位均為米。

表3 壩體及面板應(yīng)力、位移最大值

圖2 E-B模型堆石體位移圖

圖3 D-P模型堆石體位移圖

圖4 E-B模型堆石體主應(yīng)力圖

圖5 D-P模型堆石體主應(yīng)力圖

3.4.1 本構(gòu)模型對應(yīng)力變形結(jié)果影響分析

從圖2和圖3中可以看出，兩種模型水平位移對稱分布，發(fā)生在壩體的中下部，上下游壩坡面附近。E-B模型得出的完建期最大沉降位移為437.553 mm，約在壩體2/3高度處；D-P模型得出的完建期最大沉降位移為272.926 mm，約在壩體的中間部位。用E-B模型計算出的值要比D-P模型的值大1倍以上，這是因為壩體分層澆筑，用E-B模型模擬時，中上部的堆石體仍處于較低的應(yīng)力水平，彈性模量較小，因而位移值較大。從圖4和圖5中可以看出，E-B模型得出壩體斷面內(nèi)第一主應(yīng)力為0.219 MPa，發(fā)生在壩底底部中央，并且隨著高程的增加，第一主應(yīng)力主應(yīng)力逐漸減??；第三主應(yīng)力也呈現(xiàn)同樣的分布規(guī)律，其最大值為3.49 MPa；D-P模型得出壩體斷面內(nèi)第一主應(yīng)力為0.24 MPa，發(fā)生在壩底底部中央，并且隨著高程的增加，第一主應(yīng)力主應(yīng)力逐漸減小，第三主應(yīng)力也呈現(xiàn)同樣的分布規(guī)律，第三主應(yīng)力最大值為2.28 MPa。兩種模型得出的應(yīng)力分布規(guī)律大致相同，E-B模型得出壩體斷面內(nèi)第三主應(yīng)力比D-P模型稍大，這是由于E-B模型計算過程中，其底部堆石體的彈性模量隨著澆筑的進行逐漸變大，因而在壩底的第三主應(yīng)力比D-P模型得到的應(yīng)力值稍大。

從表3中可以看出，E-B模型得出完建期面板，向壩體內(nèi)最大變形為29.736 mm，向壩體外最大變形為24.148 mm；DP模型得出完建期面板向壩體內(nèi)最大變形為38.581mm，向壩體外最大變形為30.626 mm。E-B模型得出完建期面板第一主應(yīng)力為0.219 MPa，面板第三主力為0.826 MPa；D-P模型得出面板第一主應(yīng)力為0.24 MPa。第三主應(yīng)力為0.907 MPa。在完建期由于面板的自重作用，面板澆筑后基本處于受壓狀態(tài)，可以得出兩種模型分析結(jié)果基本符合實際情況。

3.4.2 E-B模型參數(shù)及荷載子步對堆石體應(yīng)力變形結(jié)果影響分析

當(dāng)參數(shù)取基準(zhǔn)值時 (見以下各表中變化率為0相應(yīng)的數(shù)值)，壩體應(yīng)力位移最大值如下：堆石壩體沉降0.437 m，河流向向上游水平位移0.041 m，河流向向下游水平位移0.039 m，堆石體第一主拉應(yīng)力為0.219 MPa，第三主壓應(yīng)力為3.49 MPa。在研究堆石壩應(yīng)力變形的影響因素時發(fā)現(xiàn)，E-B模型參數(shù)中，K、n對堆石體應(yīng)力變形結(jié)果影響較明顯，如表3和表4所示；每層荷載子步數(shù)對應(yīng)力變形結(jié)果影響也較明顯，如圖6和圖7所示。

表4 參數(shù)K對應(yīng)力和位移的影響

圖6 荷載子步與位移關(guān)系對比曲線

圖7 荷載子步與應(yīng)力關(guān)系對比曲線

從表4可知，K值變化對堆石體水平位移影響較小。K增加20%，堆石體最大沉降減小10.6%，堆石體第一主應(yīng)力增大8.6%，第三主應(yīng)力增大8.3%；K減小20%，壩體沉降增大約10.2%。這是由于，K值變大，相應(yīng)的彈性模型在變大，則堆石體的位移變小，應(yīng)力變大，反之亦然。從表5可知，當(dāng)n增加20%，堆石體最大沉降減小10.8%，堆石體第一主應(yīng)力增大6.2%，第三主應(yīng)力增大5%；n減小20%，壩體沉降增約7.6%，堆石體第一主應(yīng)力減小約8.6%，第三主應(yīng)力減小約12.3%。這是由于n值變大，相應(yīng)的彈性模量在變大，則堆石體的位移變小，應(yīng)力變大，反之亦然。

表5 參數(shù)n對應(yīng)力和位移的影響

表6 荷載子步對應(yīng)力和位移的影響

從表6可知，E-B模型中荷載子步變化對堆石體位移和應(yīng)力影響較大。荷載子步數(shù)增加66%，堆石體最大沉降增加約32.4%，堆石體第一主應(yīng)力增大14.3%，第三主應(yīng)力增大15.3%；荷載子步數(shù)減小66%，壩體沉降增大約21.9%。分析結(jié)果表明：對E-B模型而言，加載歷史對計算結(jié)果會產(chǎn)生較大影響。

4 結(jié)語

（1）由兩種模型的垂直位移圖看出，各個模型計算出的最大位移都在壩體的中問部位，用鄧肯E-B模型時，壩體產(chǎn)生的最大位移處要比D-P模型偏高，約在壩體2/3高度處，而D-P模型約在壩體的中間位置。用鄧肯E-B模型計算出的值要比D-P模型的值大1倍以上，這是因為壩體分層澆筑，用E-B模型模擬時，中上部的堆石體仍處于較低的應(yīng)力水平，彈性模量較小，因而位移值較大。

（2）兩種模型分析二維壩體完建期面板應(yīng)力變形情況，EB模型得出面板的最大向壩體內(nèi)變形為29.736 mm，最大向壩體外變形為24.148 mm，最大壓應(yīng)力為0.826 MPa。D-P模型得出面板的最大向壩體內(nèi)變形為38.581 mm，最大向壩體外變形為30.628 mm，最大壓應(yīng)力為0.907 MPa。完建期由于面板的自重作用，面板澆筑后基本處于受壓狀態(tài)且面板的變形量較小。兩種模型對面板應(yīng)力和位移影響較小，均在工程允許范圍內(nèi)。

（3）E-B模型參數(shù)K、n對面板壩的應(yīng)力位移影響最顯著。因此，對材料參數(shù)取值時，一方面要參考類似的工程經(jīng)驗，另一方面通過大量室內(nèi)及現(xiàn)場試驗，反復(fù)核對參數(shù)的準(zhǔn)確性。

[1]孫陶，高希章，楊建.紫坪鋪混凝土面板堆石壩應(yīng)力-應(yīng)變分析[J].巖土力學(xué)，2006，（2）：247-251.

[2]Drucker D C，Prager W.Soil mechanics and plastic analysis of limit design [J].Quarterly of Applied Mathematics，1952，(10)：157-165.[8]孫明權(quán)，陳姣姣，劉運紅.鄧肯-張E-B模型的ANSYS二次開發(fā)及應(yīng)用[J].華北水利水電學(xué)院學(xué)報，2013，（2）：30-34.

[3]顧凎臣，黃金明.混凝土面板堆石壩的堆石本構(gòu)模型與應(yīng)力變形分析 [J].水力發(fā)電學(xué)報，1991，（1）：13-24.

[4]吳業(yè)飛，馬海霞.基于ANSYS的土石壩應(yīng)力變形有限元分析[J].水利與建筑工程學(xué)報，2008，（4）：210-212.

[5]吳長彬，張巖，許小東.APDL二次開發(fā)功能在土石壩非線性分析中的應(yīng)用研究[J].災(zāi)害與防治工程，2010，（1）：6-10.

Influence of the Constitutive Model on Stress and Deformation of Concrete Face Rockfill Dam

Zhang Fei,Fei Wenping

（StateKeyLaboratoryofHydraulicsand Mountain RiversDevelopmentand Protection,Sichuan University,Chengdu 610065,SichuanCollege ofWater Resources and Hydropower,Sichuan University Chengdu 610065,Sichuan）

TV641.43

A

1673-9000（2017）05-0151-04

2017-05-25

國家自然科學(xué)基金項目（51409179，51379139，51609163）

張飛（1990-），男，湖北孝感人，碩士研究生，主要從事水工結(jié)構(gòu)工程研究。

費文平（1972-），男，湖北天門人，副教授，主要從事水工結(jié)構(gòu)數(shù)值模擬、現(xiàn)場監(jiān)測和室內(nèi)試驗研究。