本構模型對面板堆石壩應力變形的影響

張 飛， 費文平

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

本構模型對面板堆石壩應力變形的影響

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

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

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

本構模型；面板堆石壩；應力變形特性

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 引言

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

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

1 基本理論

1.1 理想彈塑性模型

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

其中參數α和k可以用土的粘聚力C和內摩擦角φ來表示：

1.2 鄧肯-張彈性靜力本構模型

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

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

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

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

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

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

2 ANSYS二次開發

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

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

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

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

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

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

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

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

*create,E-B.mac !創建E-B計算的宏命令

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

*set,Pa,1e5

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

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

sigm3=max(sigm3,0.1*Pa)

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

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

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

*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)!保存最大應力水平

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

*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 !結束創建宏文件

3 工程實例分析

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

3.1 工程概況

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

3.2 有限元模型

采用二維實體單元，進行堆石壩的二維有限元靜力分析。圖1為堆石壩結構的二維有限元計算網格圖，共797個單元，834個節點。

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

3.3 計算參數與工況

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

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

表1 E-B模型參數表

表2 D-P模型參數表

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

3.4 計算成果及對比

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

表3 壩體及面板應力、位移最大值

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

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

圖4 E-B模型堆石體主應力圖

圖5 D-P模型堆石體主應力圖

3.4.1 本構模型對應力變形結果影響分析

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

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

3.4.2 E-B模型參數及荷載子步對堆石體應力變形結果影響分析

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

表4 參數K對應力和位移的影響

圖6 荷載子步與位移關系對比曲線

圖7 荷載子步與應力關系對比曲線

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

表5 參數n對應力和位移的影響

表6 荷載子步對應力和位移的影響

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

4 結語

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

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

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

[1]孫陶，高希章，楊建.紫坪鋪混凝土面板堆石壩應力-應變分析[J].巖土力學，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]孫明權，陳姣姣，劉運紅.鄧肯-張E-B模型的ANSYS二次開發及應用[J].華北水利水電學院學報，2013，（2）：30-34.

[3]顧凎臣，黃金明.混凝土面板堆石壩的堆石本構模型與應力變形分析 [J].水力發電學報，1991，（1）：13-24.

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

[5]吳長彬，張巖，許小東.APDL二次開發功能在土石壩非線性分析中的應用研究[J].災害與防治工程，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

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

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

費文平（1972-），男，湖北天門人，副教授，主要從事水工結構數值模擬、現場監測和室內試驗研究。