王桔

（湘潭市水利水電勘測設(shè)計(jì)院　湘潭市　411100）

土石壩有限元分析中的應(yīng)用接觸單元法在混凝土心墻

王桔

（湘潭市水利水電勘測設(shè)計(jì)院湘潭市411100）

為研究混凝土心墻土石壩中混凝土心墻與土石壩壩體之間的摩擦接觸對心墻、壩體應(yīng)力應(yīng)變的影響，文章以某混凝土心墻土石壩為例，采用ANSYS有限元分析軟件建立二維模型，引入點(diǎn)點(diǎn)接觸單元模擬混凝土心墻與其周圍土體間的接觸問題，研究了心墻與壩體的受力特性，并將其與不考慮心墻與壩體之間摩擦接觸的情況比較。結(jié)果顯示，點(diǎn)點(diǎn)接觸單元可以模擬心墻及其周圍土體間的接觸面，且引入接觸單元后，心墻、壩體的受力均減小，運(yùn)行期心墻的受力減小最明顯。

接觸單元土石壩混凝土防滲墻應(yīng)力

引言

許多學(xué)者和專家對土石壩中防滲墻的應(yīng)用作了大量的研究，但是在研究過程中，他們大多沒有深入考慮防滲墻與壩體之間的接觸問題。實(shí)際上，混凝土防滲墻與其周圍土體之間存在著起固壁作用的泥漿，它具有一定的塑性和徐變性能，不考慮土體材料與混凝土材料間的接觸，也就忽略了泥漿固壁對混凝土防滲墻墻體受力的影響。為研究混凝土心墻土石壩的實(shí)際應(yīng)力情況，有必要引入一種方法考慮心墻與其周圍土石壩壩體之間的摩擦接觸問題，據(jù)此，本文根據(jù)接觸單元的力學(xué)機(jī)理，引入點(diǎn)-點(diǎn)接觸單元，基于某心墻土石壩，對該心墻土石壩進(jìn)行有限元分析。

1 接觸單元的力學(xué)機(jī)理

1.1接觸單元的選擇

在有限元分析中，通常在兩種材料接觸面上設(shè)置接觸單元［1～2］來模擬接觸面上應(yīng)力狀況。有限元分析中常用的接觸單元類型主要有以Goodman單元［3］為代表的無厚度單元和以Desai單元［4］為代表的有厚度單元。薄層單元的厚度難以取到合適的值，且?guī)Ш穸鹊膯卧谟邢拊治鲕浖幸膊蝗菀讘?yīng)用，而Goodman單元類型能夠較好模擬接觸面的錯(cuò)動(dòng)與張開，且單元參數(shù)相對容易確定，本次選用無厚度單元進(jìn)行分析。ANSYS支持的無厚度度單元包括面-面、點(diǎn)-面、點(diǎn)-點(diǎn)接觸三類［5］。面-面接觸單元只支持一般的靜態(tài)或瞬態(tài)、屈曲、模態(tài)、子結(jié)構(gòu)或譜分析問題；點(diǎn)-面接觸單元不適用于不光滑連續(xù)的接觸。本文選取點(diǎn)-點(diǎn)接觸單元CONTAC12。

1.2模擬接觸面的模型選擇

ANSYS軟件處理接觸面模擬問題時(shí)，需要采用合理的模型進(jìn)行模擬。有限元分析中通用的模型包括等效連續(xù)模型、接觸單元模型及接觸邊界模型。等效連續(xù)模型采用的連續(xù)性假設(shè)將導(dǎo)致參數(shù)定義困難和幾何模型過于簡化；接觸單元模型適用范圍有所限制，對精度也有一定影響［6］；接觸邊界模型往往需要進(jìn)行多次迭代求解，對計(jì)算能力的要求比較高。當(dāng)今電子計(jì)算機(jī)的快速發(fā)展，這種計(jì)算能力的要求是不難實(shí)現(xiàn)的，本次選用接觸邊界模型。

1.3接觸單元的力學(xué)機(jī)理

取接觸面單元由兩片長度為L的接觸面12和34組成，假設(shè)兩片接觸面之間由無數(shù)微小的彈簧連接。在受力前兩接觸面完全重合，即單元沒有厚度只有長度，是一種一維單元，并令兩相鄰接觸單元之間，只在結(jié)點(diǎn)處有力的聯(lián)系。每片接觸面兩側(cè)有兩個(gè)節(jié)結(jié)點(diǎn)。一個(gè)單元共有4個(gè)結(jié)點(diǎn)。

在結(jié)點(diǎn)力｛F｝e作用下，兩接觸面間的彈簧受內(nèi)應(yīng)力為：

相應(yīng)的在兩片接觸面之間產(chǎn)生相對位移：

角標(biāo)s表示切向，n表示法向。

在線性彈性假定下，應(yīng)力｛σ｝與相對位移｛ω｝成正比，寫成關(guān)系式為：，ks、kn分別為切向及法向的單位長度剛度系數(shù)（N/m3）。對于彈性材料為常量，若材料具有非線性特征，視為變量。

取線性位移模式，把每一片接觸面上沿長度方向各點(diǎn)的位移表示為結(jié)點(diǎn)的位移，用u、v分別表示為x、y方向上的位移，則頂面及底面的位移分別為：

接觸面單元內(nèi)各點(diǎn)的相對位移為（取壓為正）：

上式可寫成

由虛位移原理推得：

式中

式中L——接觸單元的長度；

kn——法向模量。

以上剛度矩陣是在標(biāo)準(zhǔn)坐標(biāo)系下建立起來的，而實(shí)際上接觸面不一定是水平的，這就需要進(jìn)行坐標(biāo)系的轉(zhuǎn)化。取接觸面仰角β，則局部坐標(biāo)系與整體坐標(biāo)系存在如下關(guān)系：

則結(jié)點(diǎn)力和結(jié)點(diǎn)位移的局部坐標(biāo)系與整體坐標(biāo)系間有如下關(guān)系：

｛F｝e=［Q］｛Fˉ｝e，｛δ｝e=［Q］｛δˉ｝e（13）

Fˉ和δˉ分別表示整體坐標(biāo)的結(jié)點(diǎn)力和結(jié)點(diǎn)位移。將式（1）～式（13）代入式（1）～式（9），即可得到用整體坐標(biāo)系表示的勁度關(guān)系：

｛Fˉ｝e=［Q］-1［k］［Q］｛δˉ｝e=［kˉ］｛δˉ｝e（14）在由｛δˉ｝e求應(yīng)力時(shí)也應(yīng)作坐標(biāo)變化。

1.4ANSYS接觸迭代算法

在有限元分析計(jì)算中，進(jìn)行接觸迭代的算法比較多，主要包括求解帶約束的泛函極值問題的方法、數(shù)學(xué)規(guī)劃的線性補(bǔ)償方法、直接求解的方法等。本文問題涉及到要實(shí)現(xiàn)接觸面上接觸力的傳遞以及接觸面沒有穿透，這兩點(diǎn)必須在接觸面的法向方向上實(shí)現(xiàn)。在用有限元軟件分析中，這兩點(diǎn)通常用求解帶約束泛函極值問題的方法來完成，帶約束的泛函極值問題的求解方法主要分為罰函數(shù)法，拉格朗日乘子法及推廣拉格朗日乘子法等［5］。

拉格朗日乘子法的剛度矩陣中有可能存在零對角，選擇求解器就會(huì)受一定的限制，且隨著自由度的增加，剛度矩陣有所增大，計(jì)算也就變得更加復(fù)雜；推廣拉格朗日乘子法需人為定義ANSYS中的TOLN（最大穿透值）值來確定接觸面允許的最大穿透值［7］，與罰函數(shù)法相比，它的迭代次數(shù)可能會(huì)更多。本此采用罰函數(shù)方法迭代求解。

2 接觸單元在混凝土心墻土石壩中的應(yīng)用

2.1工程概況

本工程實(shí)例為某一加固病險(xiǎn)粘土心墻土石壩。大壩典型橫斷面及壩體各特征數(shù)值見圖1，為簡化計(jì)算，本文將典型斷面近似看成對稱面。該土石壩已經(jīng)運(yùn)行近40年，壩基沉降基本趨于穩(wěn)定。此次加固是通過在粘土心墻中修建混凝土防滲墻的措施來解決壩體嚴(yán)重的滲漏問題，防滲墻沿壩軸線布置，墻厚0.6m，深入基巖1.0m，加固后橫斷面見圖2。

圖1　大壩加固前典型橫斷面圖

圖2　大壩加固后典型橫斷面圖

2.2模型及參數(shù)選擇

（1）模型及網(wǎng)格劃分。

根據(jù)典型斷面建立二維模型（見圖3），此次計(jì)算中，土料及混凝土心墻均采用plane42單元，混凝土與土體間接觸問題采用CONTAC12接觸單元模擬。以水平向?yàn)閄軸方向，指向壩體下游面為正；豎直向?yàn)閅軸向，向上為正；Z軸垂直紙面（Z≈0）。用四節(jié)點(diǎn)四邊形及少數(shù)的三節(jié)點(diǎn)三角形單元?jiǎng)澐帜Ｐ停Ｐ凸矂澐譃?個(gè)單元組，1377個(gè)平面單元，1 554個(gè)節(jié)點(diǎn)及60個(gè)接觸單元。考慮到防滲墻所受周圍土體之間的摩擦約束主要在壩基面以上，故接觸單元設(shè)置在壩基面以上防滲墻與土體接觸處。為了能更好的應(yīng)用點(diǎn)點(diǎn)接觸單元，在墻體與周圍土體相接處的網(wǎng)格應(yīng)盡量趨于一致。

圖3　大壩網(wǎng)格劃分圖

（2）模型材料參數(shù)及計(jì)算方案。

模型材料分為壩殼、粘土心墻、覆蓋層、基巖及混凝土防滲墻五種材料，材料參數(shù)根據(jù)相關(guān)實(shí)例確定［8～10］，如表1所示。

表1　模型材料參數(shù)表

本次按工況①：防滲墻竣工期；工況②：防滲墻正常運(yùn)行期兩種工況分別對混凝土心墻壩考慮接觸問題和不考慮接觸問題，即有、無引入接觸單元（分別為方法一、方法二）進(jìn)行平面應(yīng)力分析。后一工況的計(jì)算將以前一工況的應(yīng)力場作為初始條件，而將模型僅在壩基（壩基沉降近似趨于穩(wěn)定）土體重力荷載下并分六次加載壩體重力計(jì)算得出的應(yīng)力場作為工況①的初始應(yīng)力場，水荷載一次加載。此外，計(jì)算模型中，粘土心墻、壩殼材料采用Duncan-ChangEB模型，壩基與防滲墻看作線彈性模型，并借助ansys軟件進(jìn)行計(jì)算。

2.3計(jì)算結(jié)果

2.3.1初始應(yīng)力場

在壩基（壩基沉降近似趨于穩(wěn)定）土體重力荷載下并分六次加載壩體重力計(jì)算得壩體水平向及豎直向位移變形云圖如圖4。

圖4　水平及豎向位移變形云圖

由圖4可知，壩體水平位移和垂直位移沿壩軸線基本成對稱分布。水平位移在土體自重的作用下，上游部分水平位移指向上游，下游部分水平位移指向下游，最大水平位移為2.14cm。垂直位移自壩體兩坡面向壩體心墻逐步增大，壩體最大沉降為20.89cm，最大沉降發(fā)生壩體中部，約1/2壩高處。壩體變形符合一般土石壩受力變形規(guī)律［11］，故ANSYS能實(shí)現(xiàn)有限元中土體非線性鄧肯張E-B模型的加載。

2.3.2防滲墻竣工期應(yīng)力情況

有、無引入接觸單元求得壩體及墻體水平及豎直向應(yīng)力情況見表2。

表2　防滲墻完工期壩體及墻體水平及豎向應(yīng)力情況MPa

從Ansys計(jì)算云圖可以看出兩種方法壩體的最大應(yīng)力分布情況基本一致，最大豎向壓應(yīng)力均發(fā)生在壩基中央附近。兩種方法計(jì)算得出的防滲墻應(yīng)力沿高程分布情況見圖5、圖6。

2.3.3防滲墻正常運(yùn)行期應(yīng)力情況

在水荷載、圍土擠壓及基巖約束下，兩種方法求得壩體及墻體水平及豎直向應(yīng)力情況見表3。

兩種方法計(jì)算得出的防滲墻應(yīng)力沿高程分布情況見圖7、圖8。

圖5　防滲墻竣工時(shí)水平應(yīng)力Sx沿高程變化

圖6　防滲墻竣工期墻體豎直應(yīng)力Sy沿高程變化

表　3　防滲墻運(yùn)行期壩體水平及豎向應(yīng)力情況　MPa

圖7　防滲墻運(yùn)行期墻體水平應(yīng)力Sx隨高程變化情況

圖8　防滲墻運(yùn)行期墻體豎直應(yīng)力Sy隨高程變化情況

此外，兩種方法兩種工況下第一主應(yīng)力與第三主應(yīng)力的最大值情況見表4。

由表2～表4及圖2～圖8，點(diǎn)-點(diǎn)接觸單元可以來模擬混凝土防滲墻與其周圍土體之間的接觸問題。引入接觸單元后，兩種工況下計(jì)算得出防滲墻、壩體所受的應(yīng)力均小于未引入接觸單元時(shí)求得的應(yīng)力，防滲墻表現(xiàn)得更為明顯。運(yùn)行期較竣工期減少量更明顯，水平向應(yīng)力較豎直向應(yīng)力減少量更大，最大可以減少38.0%。

表4　防滲墻竣工期及運(yùn)行期墻體最大主應(yīng)力情況　MPa

3 結(jié)　論

本文借助有限元軟件ANSYS，引入接觸單元對混凝土心墻土石壩進(jìn)行有限元分析，主要結(jié)論如下：

（1）接觸單元可以實(shí)現(xiàn)力在不連續(xù)位移間的傳遞，將它引入到有限元分析中，是比較符合墻體實(shí)際受力情況的。

（2）當(dāng)引入接觸單元考慮防滲墻與壩體之間的接觸面問題時(shí)，心墻及壩體計(jì)算得出的應(yīng)力均小于未引入接觸單元時(shí)計(jì)算得出的應(yīng)力，防滲墻表現(xiàn)更為明顯，且運(yùn)行期受水荷載作用，應(yīng)力減少更大，此外，運(yùn)行時(shí)期的防滲墻水平向拉應(yīng)力表現(xiàn)最突出，減少量達(dá)38%。

［1］Clough G W， Duncan J M. Finite element analysis of retaining wall behavior ［J］. Journal of the Soil Mechanics and Foundations Division，1971， 97（12）： 1657-1673.

［2］安關(guān)峰，高大釗.接觸面彈粘塑性本構(gòu)關(guān)系研究［J］.土木工程學(xué)報(bào)，2001（1）：88-91.

［3］Goodman R F， Tayor R L， Brekle T L.A model for the mechanics of jointed rock ［J］. Journal of Soil Mechanics and Foundation Engineering Division， 1968， 94 （3）： 637-660.

［4］Desai C S， Zaman M M， Lightner J G， et al. Thin layer element for interfaces and joints ［J］. International Journal for Numerical and Analytical Methods in Geomechanics， 1984， 8 （1）： 19-43.

［5］劉濤，楊鳳鵬.精通ANSYS［M］.北京：清華大學(xué)出版社，2002.

［6］張楚漢，金峰.巖石和混凝土離散-接觸-斷裂分析［M］.北京：清華大學(xué)出版社，2008.

［7］王偉，盧廷浩，宰金珉，等.土與混凝土接觸面反向剪切單剪試驗(yàn)［J］.巖土力學(xué)，2009，（5）：1303-1306.

［8］路陽.土石壩加固中混凝土防滲墻應(yīng)用研究［D］.武漢：武漢大學(xué)，2010.

［9］劉娜.土石壩三維非線性有限元分析及防滲墻應(yīng)力狀態(tài)研究［D］.西安：西安理工大學(xué)，2007.

［10］華東水利學(xué)院力學(xué)教研室主編.土工原理與計(jì)算.上冊［M］.北京：水利電力出版社，1982.

王桔（1985-），女，碩士研究生，工程師，主要從事水利工程設(shè)計(jì)工作，手機(jī)：18273223164，E-mail：wj_whu@yeah.net。

［11］陳慧遠(yuǎn).土石壩有限元分析［M］.南京：河海大學(xué)出版社，1988.（2016-07-11）