ABAQUS滲流應(yīng)力耦合分析中滲透荷載施加問(wèn)題探討

，，

(1.三峽大學(xué) 水利與環(huán)境學(xué)院，湖北 宜昌 443002；2.西北農(nóng)林科技大學(xué) 水利與建筑工程學(xué)院，陜西 楊凌 712100)

1 研究背景

目前，關(guān)于巖體或混凝土滲流應(yīng)力耦合問(wèn)題的研究較多[1-4]，鑒于理論分析的繁瑣和試驗(yàn)?zāi)M的局限，數(shù)值分析成為重要的研究手段，然而不同的軟件對(duì)于水力邊界處理方法不同。ABAQUS作為一款計(jì)算功能強(qiáng)大的有限元分析軟件，提供了滲流應(yīng)力耦合分析功能，能夠求解多孔介質(zhì)的飽和滲流、非飽和滲流，以及二者的混合問(wèn)題[5-7]。該軟件的滲流應(yīng)力耦合分析理論是在Biot固結(jié)理論的基礎(chǔ)上建立起來(lái)的，可以直接考慮骨架顆粒變形與孔隙水之間的相互作用，即“直接耦合”；用戶還可以通過(guò)編寫(xiě)子程序的方式定義滲透系數(shù)與應(yīng)力、應(yīng)變等力學(xué)參數(shù)之間的關(guān)系，進(jìn)而實(shí)現(xiàn)直接耦合過(guò)程中引起的水力參數(shù)和力學(xué)參數(shù)的變化，即“間接耦合”；因此在利用該軟件進(jìn)行完全滲流應(yīng)力耦合分析時(shí)包含了“直接耦合”和“間接耦合”2種情況。ABAQUS中專門(mén)設(shè)有孔壓?jiǎn)卧?，較常規(guī)應(yīng)力位移單元多一個(gè)孔隙水壓自由度，能夠方便進(jìn)行滲流應(yīng)力耦合分析。然而在模型上下游邊界指定隨高程線性變化的孔壓來(lái)滿足水頭條件的同時(shí)，是否還需要再次對(duì)該邊界定義相應(yīng)的靜水壓力荷載，一直是學(xué)者不斷爭(zhēng)論的問(wèn)題之一。

一些學(xué)者認(rèn)為只要定義了孔壓邊界，軟件將自行根據(jù)水頭計(jì)算出模型內(nèi)滲透力，然后按照體力考慮水荷載的作用[8]；而有的學(xué)者則認(rèn)為由于材料存在孔隙，故水荷載不可能百分之百作用在邊界面上，因此從安全角度出發(fā)，認(rèn)為應(yīng)該再次施加靜水壓力荷載[9-10]。另外，利用該軟件對(duì)混凝土重力壩進(jìn)行滲流應(yīng)力耦合分析時(shí)，如果考慮了壩基的滲透作用，是否需要再單獨(dú)定義壩體揚(yáng)壓力荷載也成為學(xué)者爭(zhēng)論的問(wèn)題。在不考慮壩體透水性時(shí)，林凱生等[11]認(rèn)為由于考慮壩基滲透性，壩基面必存在孔隙水壓力，因此不需要人為定義壩體揚(yáng)壓力；而常曉林等[8]認(rèn)為建基面的孔隙水壓力并不能傳遞到壩體，需借助ABAQUS的Dload子程序在壩底面定義揚(yáng)壓力荷載?？梢?jiàn)，對(duì)上述2個(gè)爭(zhēng)論問(wèn)題尚需進(jìn)一步的探討。

本文從理論和數(shù)值模擬2個(gè)方面出發(fā)，并結(jié)合實(shí)例分析來(lái)探討應(yīng)用ABAQUS進(jìn)行混凝土重力壩滲流應(yīng)力耦合分析時(shí)滲流荷載和水邊界條件的施加問(wèn)題，旨在為研究此類問(wèn)題提供參考。

2 ABAQUS滲流應(yīng)力耦合原理

2.1 有效應(yīng)力方程

基于有效應(yīng)力原理，總應(yīng)力σ由有效應(yīng)力σ′、孔隙水壓力uw和孔隙氣壓力ua共同組成，其表達(dá)式為

σ=σ′+[χuw+(1-χ)ua]I。

(1)

式中：χ為飽和度s的函數(shù)；I為單位矩陣。當(dāng)材料完全飽和時(shí)，χ=1；當(dāng)材料非飽和時(shí)，χ=χ(s)。

2.2 應(yīng)力平衡方程

固相材料的應(yīng)力平衡方程由虛功原理表示，即在某一時(shí)刻t, 固相材料的虛功與作用在該材料上的作用力(面力和體力)產(chǎn)生的虛功相等[6]，即

式中：δv為虛位移場(chǎng)；δε為虛應(yīng)變,δε=sym(?δv/?x)；T為單位面積的表面力；f為單位體積的體積力(不含流體重量)；n為孔隙率；V為體積；S為面積邊界；ρw為流體的密度；g為重力加速度。

2.3 滲流連續(xù)性方程

流經(jīng)dV的流體應(yīng)滿足連續(xù)性方程，使得在某時(shí)間增量?jī)?nèi)流入的流體流量等于流體體積的增加速率，其表達(dá)式[5]為

(3)

利用ABAQUS進(jìn)行滲流應(yīng)力耦合分析時(shí)，無(wú)需進(jìn)行滲流場(chǎng)與應(yīng)力場(chǎng)的反復(fù)迭代，只要按時(shí)間過(guò)程連續(xù)求解便可得到全部結(jié)果。即在上述應(yīng)力平衡方程和滲流連續(xù)方程的基礎(chǔ)上，將節(jié)點(diǎn)位移和孔隙水壓力作為節(jié)點(diǎn)自由度進(jìn)行空間離散，將應(yīng)力平衡方程和滲流連續(xù)方程寫(xiě)成矩陣形式，并在同時(shí)滿足位移邊界和滲流邊界的條件下，對(duì)每個(gè)時(shí)間步內(nèi)的方程進(jìn)行求解。

圖1 計(jì)算模型Fig.1 Computation model

3 孔壓邊界和靜水壓力荷載施加問(wèn)題

以圖1所述工程問(wèn)題為例來(lái)進(jìn)行理論推導(dǎo)，并采用數(shù)值分析，比較二者結(jié)果，驗(yàn)證不同水力邊界設(shè)置的影響。某巖基處于水下15 m，巖體的密度為2 000 kg/m3，彈性模量為27 GPa，泊松比為0.2，滲透系數(shù)為1×10-8m/s[3]，孔隙率為0.091，巖基初始飽和度為1。數(shù)值模型取計(jì)算范圍為100 m×100 m。

首先，根據(jù)豎向應(yīng)力平衡和孔隙靜水壓力平衡推導(dǎo)得豎向應(yīng)力的計(jì)算公式為

S22=ρg(z0-z)-γws(1-n)(z0-z)=

[ρg-γws(1-n)](z0-z)=γ浮(z0-z) 。(4)

式中：S22為豎向有效應(yīng)力；ρ為巖體的密度；g為重力加速度，取為10 m/s2；z0為巖基表面至計(jì)算面的深度，z0-z為巖基表面以下埋深；γw為水的重度；γ浮為巖體的浮重度。

分別采用以下2種方案計(jì)算地面以下不同深度處的豎向應(yīng)力，并與理論解進(jìn)行比較。

方案1：在巖基表面定義孔壓邊界即考慮滲透力作用，不再施加靜水壓力荷載。

方案2：在巖基表面同時(shí)定義孔壓邊界和靜水壓力荷載。

表1為按以上2種方案計(jì)算所得的豎向應(yīng)力，由表1可知，方案1在每一個(gè)深度處計(jì)算結(jié)果與理論值相差均為0.15 MPa，恰好是地基表面靜水壓力值。而方案2計(jì)算結(jié)果與理論值完全相同，即考慮了式(2)中的面力項(xiàng)(靜水壓力荷載)后是正確的，說(shuō)明應(yīng)用ABAQUS進(jìn)行滲流應(yīng)力耦合分析時(shí)，定義孔壓邊界同時(shí)必須施加相應(yīng)的靜水壓力荷載。另外，本算例巖石是沒(méi)有左右邊界的，不需要考慮靜水壓力；而對(duì)于有邊界的物體，在水下，物體的上邊界和左右邊界均需定義孔壓邊界同時(shí)必須施加相應(yīng)的水荷載。

表1 不同方案下豎向應(yīng)力Table 1 Vertical stresses in different load cases

注:正值表示拉應(yīng)力；負(fù)值表示壓應(yīng)力

關(guān)于考慮材料滲流作用后，為何還要施加表面水壓力荷載，在理論上可從以下2個(gè)方面進(jìn)行說(shuō)明：

(1)從材料內(nèi)部結(jié)構(gòu)組成方面來(lái)看，不同材料的內(nèi)部結(jié)構(gòu)組成形式有以下3種情況：顆粒與顆粒的點(diǎn)接觸、封閉的實(shí)體表面及二者的混合，其截面如圖2所示。從宏觀荷載角度上來(lái)說(shuō)，對(duì)于圖2(a)所示材料，水壓力將表現(xiàn)為100%體積力的形式；對(duì)于圖2(b)所示材料，水完全不能滲入，水壓力將表現(xiàn)為100%表面水壓力的形式；對(duì)于圖2(c)所示材料，不僅要考慮滲透力，而且有必要考慮作用于材料表面的水壓力，而巖體和混凝土恰恰屬于這一類材料。需要說(shuō)明的是，對(duì)于圖2中的(a)和(c)所示的可滲透性材料而言，利用ABAQUS進(jìn)行滲流應(yīng)力耦合分析時(shí)，均需在其邊界上施加靜水壓力荷載，并同時(shí)定義孔隙水壓邊界。

圖3 壩體斷面尺寸Fig.3 Dimensions of the cross section of dam

(2)從研究對(duì)象受力平衡方面來(lái)看，若以整個(gè)巖基(固相骨架和孔隙水)為研究對(duì)象，則外荷載包括自身重力和頂部表面水壓力荷載；若以固相顆粒骨架為研究對(duì)象，根據(jù)有效應(yīng)力原理可知，此時(shí)的應(yīng)力平衡方程還應(yīng)包含孔隙水壓力項(xiàng)。利用ABAQUS進(jìn)行滲流應(yīng)力耦合計(jì)算時(shí)，在邊界條件中定義的孔壓并非以荷載的形式作用于骨架上，僅是用以提供孔隙水壓力項(xiàng)，而表面水壓力則是以外荷載的形式單獨(dú)存在的。

4 混凝土重力壩壩底面揚(yáng)壓力施加問(wèn)題

某重力壩，橫斷面見(jiàn)圖3，其上游水位為175 m，下游無(wú)水，壩體混凝土的重度為24 kN/m3，彈性模量為25 GPa，泊松比為0.167;壩基巖體的重度為27 kN/m3，彈性模量為20 GPa，泊松比為0.25。壩體混凝土滲透系數(shù)為1×10-9m/s，初始孔隙比為0.1，假定初始飽和度為0；壩基巖體滲透系數(shù)為1×10-7m/s，初始孔隙比為0.1，假定初始飽和度為1。

圖4 壩內(nèi)孔壓分布Fig.4 Distribution of pore pressure in the dam

4.1 壩體浸潤(rùn)線的確定

壩體和壩基均按可滲

材料考慮，計(jì)算得到的壩體浸潤(rùn)線如圖4所示。從圖4可以看出，壩體內(nèi)浸潤(rùn)線明顯偏高，與實(shí)際工程存在差異，原因是模型中未考慮壩內(nèi)的排水設(shè)施等因素，而本模型旨在研究揚(yáng)壓力的施加問(wèn)題，因此對(duì)下文的分析沒(méi)有影響。

圖5 壩內(nèi)應(yīng)力分布Fig.5 Distribution of stress in the dam

4.2 考慮壩基和壩體均透水情況

認(rèn)為壩基和壩體均為

透水介質(zhì)，該工況計(jì)算得到的應(yīng)力屬于有效應(yīng)力(見(jiàn)圖5)，即考慮到孔隙水壓力的影響，而建基面揚(yáng)壓力荷載的實(shí)質(zhì)就是認(rèn)為壩體不透水時(shí)作用于壩體與壩基膠結(jié)面處孔隙水壓力，因此不必再單獨(dú)考慮壩基揚(yáng)壓力荷載，故以該情況作為參考基準(zhǔn)來(lái)研究揚(yáng)壓力的施加問(wèn)題。

4.3 考慮壩基透水，壩體不透水情況

由于壩體混凝土的滲透性較小，計(jì)算中可以將其視為不透水介質(zhì)來(lái)考慮，僅認(rèn)為壩基透水。為了避免建基面上共用結(jié)點(diǎn)的單元類型不同的問(wèn)題，數(shù)值分析時(shí)通過(guò)定義壩體滲透系數(shù)為0來(lái)反映壩體的不透水性，故該工況下壩體和壩基的單元類型均為CPE8P(CPE8P是ABAQUS中的一種單元類型，即考慮孔壓的8節(jié)點(diǎn)平面應(yīng)變單元)，此時(shí)整體模型的滲流場(chǎng)分布如圖6所示。

圖6 考慮壩基的滲流作用，不考慮壩體滲流時(shí) 整體模型的滲流場(chǎng)Fig.6 Seepage field of the whole model in consideration of foundation seepage effect in the absence of dam body seepage

按照目前已有的理論[12-13]，應(yīng)該考慮的荷載有壩基滲透力、上下游壩面靜水壓力及建基面揚(yáng)壓力，而在運(yùn)用ABAQUS進(jìn)行分析時(shí)，考慮壩基的滲流作用，在壩底面存在孔隙水壓力，如圖7所示。

圖7 建基面孔隙水壓力Fig.7 Pore water pressure in the dam base surface

針對(duì)圖7所示孔隙水壓力能否以荷載的形式傳遞到壩上，即此時(shí)是否應(yīng)在壩底面再定義揚(yáng)壓力荷載的問(wèn)題，分別采用以下4種方案進(jìn)行分析。

方案1：考慮壩基的滲流作用，不考慮壩體滲流，也不施加壩底面揚(yáng)壓力。壩體內(nèi)無(wú)孔隙水壓作用，探討該軟件是否可自動(dòng)施加壩底面的揚(yáng)壓力作用。

方案2：考慮壩基的滲流作用，不考慮壩體滲流，但在建基面施加揚(yáng)壓力。揚(yáng)壓力按照?qǐng)D7所示面力分布進(jìn)行施加。

方案3：考慮壩基的滲流作用，不考慮壩體滲流，但壩體自身浸潤(rùn)線以下按照浮重度考慮[14]。該情況雖考慮壩體內(nèi)有效應(yīng)力，但不能考慮滲透力作用。

方案4：不考慮壩體和壩基的滲流作用，但在建基面施加揚(yáng)壓力，壩基和壩體內(nèi)均無(wú)孔隙水壓，相當(dāng)于按總應(yīng)力計(jì)算。

在以上4種分析方案中，均考慮了重力，壩體上、下游面水壓力，壩基上下游表面庫(kù)水壓力荷載。

采用以上4種方案及考慮壩體滲流作用情況計(jì)算所得壩底面豎向應(yīng)力如圖8所示，壩體內(nèi)部的豎向應(yīng)力分布如圖9所示。方案1、2和4均未考慮壩體內(nèi)的孔隙水壓力，因此計(jì)算得到的壩體應(yīng)力為總應(yīng)力，其中，方案4是不考慮壩體壩基滲透作用時(shí)，進(jìn)行重力壩靜力分析時(shí)所普遍采用的方法，故以方案4作為判斷依據(jù)，來(lái)比較僅考慮壩基滲流作用后，ABAQUS是否可自動(dòng)施加壩建基面的揚(yáng)壓力荷載。

圖8 壩底面豎向應(yīng)力Fig.8 Vertical stress in the bottom of dam

圖9 壩內(nèi)豎向應(yīng)力分布等值線Fig.9 Distribution of vertical stress in the dam

由圖8、圖9知：方案1計(jì)算得到的壩內(nèi)和壩底面豎向應(yīng)力與方案4差異較大；方案2計(jì)算得到的壩內(nèi)和壩底面的豎向應(yīng)力與方案4基本相同。其原因在于，方案2不僅考慮了建基面上的孔隙水壓力，而且在壩底面施加了揚(yáng)壓力荷載。這說(shuō)明定義壩基的滲流作用后，雖然建基面上存在孔隙水壓力，但是其并不能以外荷載的形式傳遞到壩底面，所以在利用ABAQUS進(jìn)行重力壩滲流應(yīng)力耦合分析時(shí)，僅考慮壩基透水，而不考慮壩體透水時(shí)，應(yīng)單獨(dú)定義壩體揚(yáng)壓力作用。

下面將針對(duì)揚(yáng)壓力的2種不同施加方式進(jìn)行探討，比較圖9(b)和圖9(a)發(fā)現(xiàn)，通過(guò)在壩底面施加面力來(lái)定義揚(yáng)壓力時(shí)，僅對(duì)壩底部分區(qū)域的應(yīng)力有影響，且造成壩踵處應(yīng)力集中；比較圖9(c)和圖9(a)發(fā)現(xiàn)，通過(guò)對(duì)浸潤(rùn)線以下壩體按浮重度考慮來(lái)定義揚(yáng)壓力時(shí)，雖然壩踵處也會(huì)出現(xiàn)應(yīng)力集中，但該工況下壩內(nèi)和壩底面的豎向應(yīng)力與考慮壩體滲流作用時(shí)(圖5)所得結(jié)果基本一致。

因此，筆者認(rèn)為在利用該軟件對(duì)重力壩進(jìn)行滲流應(yīng)力耦合分析時(shí)，若考慮壩基滲透性而不考慮壩體滲透性時(shí)應(yīng)對(duì)壩基施加揚(yáng)壓力，但揚(yáng)壓力不能以面力的形式施加，宜采用浸潤(rùn)線以下壩體按浮重度進(jìn)行計(jì)算的方式；若同時(shí)考慮壩體和壩基的滲透性時(shí)則無(wú)需再對(duì)建基面施加揚(yáng)壓力。

5 結(jié) 論

本文旨在結(jié)合算例探討應(yīng)用ABAQUS進(jìn)行滲流應(yīng)力耦合分析時(shí)模型中水荷載及邊界條件施加的問(wèn)題，分析結(jié)果表明：

(1)利用ABAQUS進(jìn)行滲流應(yīng)力耦合分析時(shí)，應(yīng)在壩體和壩基表面同時(shí)定義水荷載和孔壓邊界條件，定義的表面水壓力將以外荷載的形式作用于固相顆粒骨架上，而定義孔壓邊界的作用則是用以定義孔隙水壓力項(xiàng)，以及滿足滲流連續(xù)方程。而并不像有些學(xué)者認(rèn)為的只要定義滲流孔壓邊界即考慮了滲透力的作用，從而無(wú)需再施加靜水壓力。

(2)不考慮壩體滲透性時(shí)，在建基面雖然存在孔隙水壓力，但分析發(fā)現(xiàn)其并不能以荷載的形式傳遞到壩體上，此時(shí)需重新定義揚(yáng)壓力荷載，但不能以面力的形式施加，宜采用浸潤(rùn)線以下壩體按浮重度進(jìn)行計(jì)算的方式考慮；而考慮壩體滲透性時(shí)，壩內(nèi)存在孔隙水壓力作用，故此時(shí)無(wú)需再單獨(dú)定義揚(yáng)壓力。

參考文獻(xiàn)：

[1] 楊天鴻, 唐春安, 梁正召, 等. 脆性巖石破裂過(guò)程損傷與滲流耦合數(shù)值模型研究[J]. 力學(xué)學(xué)報(bào), 2003, 35(5): 533-541.

[2] 梁 通, 金 峰. 基于廣義有效應(yīng)力原理的混凝土壩分析[J]. 水力發(fā)電學(xué)報(bào), 2007, 28(2): 47-51.

[3] 柴軍瑞, 仵彥卿. 碾壓混凝土壩滲流場(chǎng)與應(yīng)力場(chǎng)耦合分析的數(shù)學(xué)模型[J]. 水利學(xué)報(bào), 2000, (9): 33-35.

[4] 沈振中, 陳小虎, 徐力群. 重力壩應(yīng)力-滲流相互作用的無(wú)單元耦合分析[J]. 巖土力學(xué), 2008, 29(增): 74-79.

[5] HIBBETT D, KARLSSON B, SORENSEN P. ABAQUS/Standard: User’s Manual[K]. U.S.: Hibbitt, Karlsson & Sorenson, 1998.

[6] 陳衛(wèi)忠, 伍國(guó)軍, 賈善坡. ABAQUS在隧道及地下工程中的應(yīng)用[M]. 北京: 中國(guó)水利水電出版社, 2010.

[7] 莊 茁, 由小川, 廖劍暉, 等.基于ABAQUS的有限元分析和應(yīng)用[M]. 北京: 清華大學(xué)出版社, 2009.

[8] 常曉林, 袁 曦. 基于流-固耦合的強(qiáng)度折減法在重力壩抗滑穩(wěn)定分析中的應(yīng)用[J]. 武漢大學(xué)學(xué)報(bào)(工學(xué)版), 2012, 45(5): 545-550.

[9] 柴軍瑞. 混凝土壩水荷載討論[J]. 水電能源科學(xué), 2000, 18(2): 18-21.

[10] 張曉詠, 戴自航. 應(yīng)用ABAQUS程序進(jìn)行滲流作用下邊坡穩(wěn)定分析[J]. 巖石力學(xué)與工程學(xué)報(bào), 2010, 29(增): 292-294.

[11] 林凱生, 李宗利. 高滲透孔隙水壓作用下混凝土損傷破壞過(guò)程數(shù)值分析[D]. 楊凌: 西北農(nóng)林科技大學(xué), 2010.

[12] 王 媛, 速寶玉. 壩基應(yīng)力計(jì)算中的水荷載組合形式[J]. 勘察科學(xué)技術(shù), 1995, (3): 3-7.

[13] 黃耀英, 沈振中, 田 斌. 地基水荷載對(duì)混凝土壩位移影響研究[J]. 水利水運(yùn)工程學(xué)報(bào), 2010, (1): 42-49.

[14] 范書(shū)立, 陳建云, 林 皋. 滲透壓力對(duì)重力壩有限元分析的影響研究[J]. 巖土力學(xué), 2007, 28(增): 575-581.