冉 蠡 程 琳，2 楊 杰 仝 飛 陳詩怡

(1. 西安理工大學(xué) 西北旱區(qū)生態(tài)水利工程國家重點(diǎn)實(shí)驗(yàn)室培育基地，陜西 西安 710048；2. 河海大學(xué) 水文水資源與水利工程科學(xué)國家重點(diǎn)實(shí)驗(yàn)室，南京 210098；3. 福州大學(xué)至誠學(xué)院，福建 福州 350002)

流變在混凝土中一般被稱為徐變，在巖石力學(xué)中一般被稱為蠕變，均表示變形隨時間發(fā)生變化[1]．對流變的分析，目前主要是采用流變模型．由于混凝土的力學(xué)性質(zhì)(包括徐變性質(zhì)在內(nèi))與材料的齡期有關(guān)，當(dāng)混凝土的齡期較小時，材料受力后徐變發(fā)展得很快，只有當(dāng)它的齡期達(dá)到相當(dāng)大的時候(如1年以上)，徐變的發(fā)展才可以近似地認(rèn)為與齡期無關(guān)，所以在建立混凝土的徐變理論時，應(yīng)該考慮材料齡期的影響．混凝土的流變，可用專門的理論來分析，也可采用元件流變模型來分析．對于巖石的蠕變，則一般用元件流變模型[2]來模擬．圖1是典型的巖石蠕變過程線，由圖可知，巖石的蠕變可分為3個階段.第1階段(a-b)，減速蠕變階段：應(yīng)變速率隨時間增加而減小；第2階段(b-c)，等速蠕變階段：應(yīng)變速率保持不變；第3階段(c-d)，加速蠕變階段：應(yīng)變速率隨時間增加而增加．

圖1 巖石蠕變曲線

流變計(jì)算較為復(fù)雜，一般需借助相關(guān)的計(jì)算軟件．在結(jié)構(gòu)有限元分析的各種商業(yè)軟件中，由于MSC. Marc軟件[3-4]擁有多種高級非線性有限元求解器、充足的單元類型、豐富的材料庫以及計(jì)算精度高等優(yōu)點(diǎn)，本文以MSC. Marc軟件為基礎(chǔ)，進(jìn)行二次開發(fā)，編制相關(guān)的程序，來實(shí)現(xiàn)流變計(jì)算．然后采用某抽水蓄能電站的面板堆石壩為算例，對開發(fā)的程序進(jìn)行了驗(yàn)證．

1 流變問題的基本原理

設(shè)流變過程中結(jié)構(gòu)的靜力平衡方程為[5-6]

(1)

式中，B為應(yīng)變矩陣；σ為應(yīng)力；V為體積；f為外力．

采用小應(yīng)變和小位移假設(shè)，即矩陣BT和坐標(biāo)與時間無關(guān)．將上式對時間t求導(dǎo)，并根據(jù)應(yīng)力應(yīng)變轉(zhuǎn)換關(guān)系進(jìn)行推導(dǎo)，可以得到平衡方程的微分形式為

(2)

(3)

處理流變問題一般有顯示方法和隱式方法兩種．所謂的顯式方法(初應(yīng)力法)就是如式(3)所示蠕變項(xiàng)被完全處理成偽荷載，這樣不需要重新組裝求解剛度矩陣．所謂的隱式方法就是將式(3)右端第二項(xiàng)移到方程的左邊，使剛度矩陣變成與流變應(yīng)變相關(guān)的量，這樣每個時間步都要重新組裝剛度矩陣．

2 流變分析的伯格斯模型

在解決流變問題時，通常采用元件模型進(jìn)行有限元分析[7]．元件模型是一種可模擬材料粘彈性性質(zhì)的力學(xué)模型，其中伯格斯模型較為常用．

伯格斯模型是由一個Kelvin模型和一個Maxwell模型串聯(lián)得到的[8]，如圖2所示．

圖2 伯格斯流變模型

則由Kelvin模型和Maxwell模型可推導(dǎo)出，在全量形式下，t0時刻到t=t0+Δt時段總流變應(yīng)變?yōu)?/p>

(4)

式中，EK為彈性模量；ηK為偏應(yīng)力下粘性系數(shù)；ηM為靜水應(yīng)力下粘性系數(shù)；{εV，K}t0、{εV，Mt0}分別為t0時刻Kelvin模型與Maxwell模型的粘性應(yīng)變；[C]為泊松比矩陣．

全量形式下，總流變應(yīng)變增量，即減去初值({εV，K}t0+{εV，M}t0)后的值為

(5)

式(5)右邊第一項(xiàng)是Kelvin模型對應(yīng)的流變增量，第二項(xiàng)為Maxwell模型對應(yīng)的流變增量．用偏應(yīng)力和靜水應(yīng)力可分別表達(dá)為

(6)

(7)

式中，GK為剪切模量；BK為體積模量；Sij為偏應(yīng)力；pm為平均應(yīng)力．

3 MSC. Marc的流變分析與二次開發(fā)

3.1 MSC. Marc軟件的對流變問題的處理

對于流變問題，當(dāng)結(jié)構(gòu)不存在塑性，即處于圖1所示的蠕變的第一、二階段時，會有以下3種情況：①僅采用Maxwell形式的模型；②僅采用Kelvin形式的模型(VISOELASTIIC)；③采用兩種模型的組合．

MSC. Marc軟件處理這種元件流變模型的常用思路見表1[9-11]．

表1 MSC. Marc流變分析二次開發(fā)快速參考

1)如果流變模型是Kelvin模型和Maxwell模型的組合(如伯格斯模型)，可采用以下方法：將總流變分為Kelvin模型部分和Maxwell模型部分，Kelvin模型部分的采用子程序CRIVIPS定義(只能定義偏應(yīng)力產(chǎn)生的)，Maxwell部分采用CRPLAW定義或其它方法．如果還要考慮體積流變則采用VSWELL定義．(先將應(yīng)力分為偏應(yīng)力和平均應(yīng)力，平均應(yīng)力產(chǎn)生的體積流變用VSWELL定義，偏應(yīng)力產(chǎn)生的對應(yīng)于Maxwell部分的流變用②中的方法定義，偏應(yīng)力產(chǎn)生的Kelvin形式的流變采用CRIVIPS定義．當(dāng)流變應(yīng)變的方程是采用應(yīng)力全量的形式給出時，由于隱含著體積流變和偏應(yīng)力流變形式相同，需將全量形式轉(zhuǎn)化為偏應(yīng)力和靜水應(yīng)力的分開形式進(jìn)行開發(fā))．

2)如果僅是Maxwell形式的模型，或是經(jīng)驗(yàn)?zāi)Ｐ?，指?shù)模型等，(如三參數(shù)模型，經(jīng)驗(yàn)流變模型，Maxwell模型)可以采用方案一中的顯式或隱式方法．

3)如果是僅是Kelvin形式的模型(典型的如廣義Kelvin模型)這時就直接采用CRIVIPS進(jìn)行定義偏應(yīng)力流變，如果還要考慮體積流變，則采用VSVELL定義．

4)當(dāng)計(jì)算較為粗略時，可以直接采用CRPLAW子程序進(jìn)行開發(fā)，計(jì)算中所需要的中間量(如上一步的Kelvin流變，Maxwell流變)均由用戶自己定義．這時對于Kelvin模型和Maxwell模型的組合形式的模型而言，采用CRPLAW直接定義應(yīng)變增量相當(dāng)于將兩種類型的模型產(chǎn)生的流變直接相加統(tǒng)一由用戶定義，不再進(jìn)行區(qū)分．

3.2 MSC.MARC軟件的二次開發(fā)實(shí)現(xiàn)

根據(jù)上述分析，采用MSC. Marc軟件實(shí)現(xiàn)伯格斯模型主要有兩種方式：

1)采用全量形式下的表達(dá)式(5)直接利用CRPLAW子程序來實(shí)現(xiàn)；

2)根據(jù)式(6)和式(7),分別用CRPLAW子程序?qū)崿F(xiàn)Maxwell偏應(yīng)變流變分量，即兩式的后半部分；用CRIVIPS子程序?qū)崿F(xiàn)Kelvin偏應(yīng)變流變分量，即兩式的后前部分．

為保證模型精度，本文采用第二種方法進(jìn)行二次開發(fā)．采用CRPLAW子程序定義顯式粘塑性模型的非彈性應(yīng)變率，CRIVIPS子程定義顯式的廣義Kelvin模型．考慮到結(jié)構(gòu)同時產(chǎn)生的體積流變，則采用VSWELL子程序?qū)崿F(xiàn)體積流變分量，并且認(rèn)為體應(yīng)力引起粘性變形的規(guī)律與應(yīng)力偏量引起粘性變形的規(guī)律都和一般應(yīng)力張量引起粘性變形的規(guī)律相同．

4 算 例

4.1 工程概況

某水利樞紐主要由上水庫、輸水系統(tǒng)、發(fā)電廠房及下水庫等4部分組成．上水庫布置于牢山寨北坡位置，擋水建筑物為主壩和副壩各1座，均為混凝土面板堆石壩，水庫采用全庫盆表面水平防滲型式，即庫周采用鋼筋混凝土面板防滲、庫底采用土工膜防滲．上水庫的主壩壩頂高程351 m，壩頂寬度10 m，壩頂長度1 110 m，最大壩高131 m(壩軸線處)，上游面坡比1∶1.4，下游面每20 m高設(shè)一寬2 m的馬道，馬道間坡比1∶1.35，其綜合坡比約1∶1.45．

4.2 計(jì)算模型和參數(shù)

建立該抽水蓄能電站上水庫面板堆石壩最大斷面的二維平面有限元流變計(jì)算分析模型，如圖3所示，有限元模型單元數(shù)278個，節(jié)點(diǎn)數(shù)307個，設(shè)置接觸單元6個．不考慮流變時，堆石體的本構(gòu)采用E-B模型來模擬，混凝土材料采用彈性本構(gòu)關(guān)系．對巖石和混凝土采用伯格斯模型來計(jì)算分析結(jié)構(gòu)的流變特性．

圖3 壩體最大斷面的有限元模型

國內(nèi)外的試驗(yàn)和觀測表明，堆壩的流變變形可以占到壩體瞬時變形的30%～50%．考慮大壩10年后，流變引起的變形占總沉降量的30%左右的情形，并參考其它類似工程的試驗(yàn)和計(jì)算結(jié)果擬定了壩體、回填區(qū)和面板材料的流變計(jì)算參數(shù)見表2．

表2 伯格斯流變模型參數(shù)

4.3 計(jì)算成果

壩體部分結(jié)點(diǎn)在大壩運(yùn)行10年后考慮和不考慮流變時的水平向位移和垂直向位移見表3．相應(yīng)結(jié)點(diǎn)的位移歷時曲線如圖4～5所示．大壩蓄水期末和運(yùn)行10年后考慮流變時的垂直位移分布如圖6～7所示．水平位移向上游為正，向下游為負(fù)；垂直位移向上為正．由圖表可知：考慮流變時大壩位移的分布規(guī)律與不考慮流變時基本相同．大壩運(yùn)行10年后，垂直向位移最大值出現(xiàn)在壩體中部1/3～1/2壩高處．流變的效應(yīng)使壩體水平向和垂直向位移的量都有所增加，但大部分結(jié)點(diǎn)后期的位移趨向平穩(wěn)．

表3 考慮和不考慮流變時的位移比較

圖4 壩體部分結(jié)點(diǎn)的水平位移歷時曲線(單位：m)

圖5 壩體部分結(jié)點(diǎn)的垂直位移歷時曲線(單位：m)

圖6 蓄水期末垂直向位移(單位：m)

圖7 大壩運(yùn)行10年后考慮流變時的垂直向位移(單位：m)

壩體部分結(jié)點(diǎn)的最大主應(yīng)力和最小主應(yīng)力的歷時曲線如圖8所示．大應(yīng)力以拉為正，壓為負(fù)．由圖可知：考慮流變效應(yīng)時壩體主應(yīng)力的分布規(guī)律與不考慮流變時基本相同．根據(jù)計(jì)算第一主應(yīng)力的最小值為-2.73 MPa，第二主應(yīng)力的最小值為-1.41 MPa．各向應(yīng)力在后期都趨向平穩(wěn)．

圖8 壩體部分結(jié)點(diǎn)的最小主應(yīng)力歷時曲線(單位：MPa)

圖9是主壩面板考慮流變和不考慮流變時的撓曲線的對比．考慮流變前后面板撓曲線的分布規(guī)律基本相同，都呈現(xiàn)出兩端大中間小的特點(diǎn)．流變使面板法向位移有所增加，面板最大撓度由297 mm變?yōu)?55 mm．根據(jù)面板應(yīng)力的計(jì)算結(jié)果，不計(jì)流變時，主壩面板的最大法向壓應(yīng)力為1 805 kPa，考慮流變時大壩運(yùn)行10年后面板最大的法向壓應(yīng)力為1 820 kPa，最大切向應(yīng)力為1 810 kPa，較未考慮流變時的計(jì)算結(jié)果均有所增加．

圖9 考慮和不考慮流變情況下面板撓曲線的對比(單位：m)

由以上分析可知，在考慮堆石料與混凝土的流變情況時，其壩體與混凝土面板的應(yīng)力較不考慮材料流變時均有所增加．表明堆石料與混凝土的流變對壩體與混凝土面板的變形特性具有較大的影響，堆石料與混凝土的流變使壩體與面板的變形增加，應(yīng)力增大；而壩體變形的增加又會影響到混凝土面板應(yīng)力變形情況，進(jìn)而加大混凝土面板的應(yīng)力值．

5 結(jié) 論

本文基于MSC. Marc軟件進(jìn)行二次開發(fā)，編制相關(guān)程序，采用伯格斯模型進(jìn)行流變分析，再以某面板堆石壩為例，計(jì)算其流變情況下的應(yīng)力變形情況．根據(jù)以上理論分析和工程實(shí)例的計(jì)算結(jié)果，可得以下結(jié)論：

1)考慮流變作用后的壩體位移和應(yīng)力分布規(guī)律與未考慮流變時基本相同，位移和應(yīng)力的量值在考慮流變后都有一定增加，但后期趨向平穩(wěn)．

2)考慮流變作用后大壩各部位的位移和應(yīng)力計(jì)算結(jié)果處于允許水平，面板和連接板應(yīng)力和位移均有所增加，應(yīng)加強(qiáng)相關(guān)結(jié)構(gòu)部位的檢查和維護(hù)．

3)在分析面板堆石壩的變形時考慮堆石料與混凝土的流變更能真實(shí)的反映出面板堆石壩的應(yīng)力變形情況．對于高面板壩的變形分析，若不考慮流變，則會過小估計(jì)混凝土面板的應(yīng)力變形，從而增大面板開裂風(fēng)險．

4)本文開發(fā)程序的計(jì)算結(jié)果符合一般工程規(guī)律，證明二次開發(fā)是成功的，相關(guān)程序可以應(yīng)用于實(shí)際工程問題．

參考文獻(xiàn)：

[1] 劉 鈺，黃耀英，唐騰飛，等．大壩混凝土分?jǐn)?shù)階徐變模型探討[J]．長江科學(xué)院院報(bào)，2017，5(34)：1-5．

[2] 殷宗澤，等．土工原理[M]．北京：中國水利水電出版社，2007：368-141．

[3] 關(guān)天定，崔 杰，李亞東．基于MSC.Marc二次開發(fā)的土體等價非線性粘彈性模型[J]．土工基礎(chǔ)，2011，25(6)：57-59．

[4] 屠立峰，包騰飛，陳 波．基于MSC.Marc軟件的面板堆石壩加高的可行性研究[J]．三峽大學(xué)學(xué)報(bào)(自然科學(xué)版)，2015，37(5)：9-13．

[5] 佘亞鵬，王永明，盧繼旺，等．200m級面板堆石壩流變分析[J]．水電能源科學(xué)，2010，4(28)：90-92+160．

[6] 王德信．面板堆石壩流變研究[D]．南京：河海大學(xué)，2003．

[7] 王瑞駿，李 陽，丁占峰．堆石料流變模型參數(shù)敏感性分析的正交試驗(yàn)法[J]．水利學(xué)報(bào)，2016，47(2)：245-252．

[8] 楊振偉，金愛兵，周 喻，等．伯格斯模型參數(shù)調(diào)試與巖石蠕變特性顆粒流分析[J]．巖土力學(xué)，2015，36(1)：240-248．

[9] 黃耀英，沈振中，吳中如．基于MSC．Marc二次開發(fā)進(jìn)行粘彈性問題分析黑[J]．長江學(xué)院院報(bào)，2006，23(6)：95-98．

[10] 郭典塔，周翠英，石漢生，等．對Marc軟件進(jìn)行二次開發(fā)及其在地下結(jié)構(gòu)中的應(yīng)用[J]．地下空間與工程學(xué)報(bào)，2009，5(4)：786-791．

[11] 周 偉，常曉林．高混凝土面板堆石壩流變的三維有限元數(shù)值模擬[J]．巖土力學(xué)，2006，27(8)：1389-1393．