郭 晴1,，劉東海，陳 輝

(1.河北工程大學(xué) 水利水電學(xué)院，河北 邯鄲 056021； 2.天津大學(xué) 水利工程仿真與安全國家重點(diǎn)實(shí)驗(yàn)室，天津 300350 )

1 研究背景

瀝青混凝土心墻壩是以堆石為主體，瀝青混凝土作為壩體防滲材料的一種土石壩型[1]。瀝青混凝土自身良好的防滲性、變形性和黏結(jié)性，以及瀝青混凝土心墻可與壩殼材料同步施工以縮短施工工期的性能，使得瀝青混凝土心墻壩在世界范圍內(nèi)的建設(shè)規(guī)模不斷擴(kuò)大，已經(jīng)成為當(dāng)前具有極大發(fā)展?jié)摿Φ男滦蛪涡汀Ｒ虼?，有必要通過有限元計(jì)算，模擬大壩的運(yùn)行情況，找出壩體可能存在的安全隱患，為大壩正常運(yùn)行提供分析依據(jù)。涂幸等[2]、代凌輝[3]、Coroller[4]基于有限元分析軟件對瀝青混凝土心墻壩進(jìn)行了三維有限元分析，探討了在不同工況下壩體的應(yīng)力、變形特性；張蕓蕓等[5]、Valstad等[6]對考慮堆石料濕化效應(yīng)或受地震動(dòng)荷載作用下的瀝青混凝土心墻壩進(jìn)行數(shù)值模擬，提出了改善壩體受力特性的有效方案。

但是，目前國內(nèi)外的研究在采用有限元分析時(shí)多是假定同一分區(qū)采用相同的壩料力學(xué)參數(shù)。大壩實(shí)際施工質(zhì)量在空間上存在差異性，這種差異性勢必會(huì)導(dǎo)致壩體材料物理力學(xué)參數(shù)在空間分布上的不確定性，如果按壩體材料的設(shè)計(jì)參數(shù)計(jì)算分析大壩的應(yīng)力和變形將會(huì)與實(shí)際情況存在差異。如楊鴿等[7]采用局部平均細(xì)分法模擬了二維的堆石料物理力學(xué)性質(zhì)隨機(jī)場，表明忽略材料的不確定性可能導(dǎo)致大壩的地震反應(yīng)被低估。

本文采用考慮空間差異的隨機(jī)有限元分析方法，研究瀝青混凝土心墻壩的應(yīng)力、變形情況；統(tǒng)計(jì)壩體應(yīng)力、應(yīng)變的空間分布規(guī)律，進(jìn)而為大壩安全分析與設(shè)計(jì)優(yōu)化提供參考依據(jù)。

2 瀝青心墻壩隨機(jī)有限元分析方法

2.1 壩體隨機(jī)有限元計(jì)算流程

隨機(jī)有限元法是將蒙特卡羅(Monte-Carlo)技術(shù)與有限元分析方法相結(jié)合，通過有限次的循環(huán)計(jì)算參數(shù)統(tǒng)計(jì)量，進(jìn)而研究參數(shù)的分布特征，并且能夠保證概率結(jié)果具有較高的置信度。在模擬的過程中，每一個(gè)仿真循環(huán)是完全獨(dú)立的，而且任何一組循環(huán)與其他仿真循環(huán)結(jié)果無關(guān)，因此非常適合并行計(jì)算[8]。本文采用蒙特卡羅直接抽樣法模擬壩體在施工過程中任何材料參數(shù)下的響應(yīng)，一個(gè)仿真循環(huán)得到壩體在某組材料參數(shù)下的應(yīng)力、變形等響應(yīng)。隨機(jī)有限元計(jì)算流程如圖1所示。

(1)

2.2 材料物理力學(xué)參數(shù)的不確定性

在結(jié)構(gòu)的隨機(jī)有限元計(jì)算中，涉及的不確定性因素有：幾何不確定性、初始條件不確定性、材料參數(shù)不確定性、邊界條件不確定性，以及載荷不確定性等[9]。對于瀝青混凝土心墻壩，它的幾何模型大小、邊界條件和受力情況都基本不變化。天然條件下存在地基覆蓋層材料的非均勻性與復(fù)雜性，因此，本文將研究考慮壩體與覆蓋層材料參數(shù)的不確定性對瀝青混凝土心墻壩應(yīng)力、變形的影響。

由于瀝青混凝土與堆石料都是典型的散粒體材料，具有材料非線性特點(diǎn)，因此材料的本構(gòu)模型選取鄧肯張E-B模型[10]。將密度ρ與鄧肯張E-B模型材料參數(shù)(黏聚力c、內(nèi)摩擦角φ、破壞比Rf、彈性模量數(shù)K、體積模量指數(shù)m、彈性模量指數(shù)n、體積模量數(shù)Kb)作為隨機(jī)變量進(jìn)行抽樣，Kur的取值是根據(jù)K的概率分布某次抽樣后得到的數(shù)值，將2K賦值給Kur，即Kur=2K。然后再將這8個(gè)材料參數(shù)與密度作為該次抽樣的堆石料材料參數(shù)賦值給堆石料的單元，從而進(jìn)行該次有限元分析。

由于瀝青混凝土心墻厚度較薄，一般為0.5～1.2 m，并且作為大壩防滲結(jié)構(gòu)，其施工質(zhì)量控制嚴(yán)格，心墻料物理與力學(xué)參數(shù)差異不大。根據(jù)實(shí)際工程中的現(xiàn)場檢測，瀝青混凝土心墻的技術(shù)指標(biāo)大都符合設(shè)計(jì)控制指標(biāo)[11]。因此，本文只考慮壩體除心墻外的過渡層區(qū)、堆石區(qū)和覆蓋層區(qū)的密度與鄧肯張E-B模型材料參數(shù)的不確定性，而心墻區(qū)的材料物理力學(xué)參數(shù)仍舊使用設(shè)計(jì)值。

2.3 模型參數(shù)批量賦值方法

為了使得有限元分析充分反映大壩施工壓實(shí)質(zhì)量的空間差異性，本文采用每個(gè)有限元單元對應(yīng)不同的鄧肯張E-B模型參數(shù)的方法。有限元分析軟件ABAQUS可以通過調(diào)用全命令流式的inp文件，直接進(jìn)行結(jié)構(gòu)的有限元分析，因此可以通過改變inp文件中的命令來實(shí)現(xiàn)對大壩所有網(wǎng)格材料的重新賦值。批量賦值的過程采用開發(fā)的Fortran程序來實(shí)現(xiàn)，具體流程如下：

(1)確定賦值區(qū)域單元號(hào)E。通過ABAQUS軟件建立壩體不同分區(qū)的集合，從inp文件中采集堆石區(qū)、過渡層以及覆蓋層的單元號(hào)E。

(2)定義單元集合FE。在Fortran程序中讀入單元號(hào)E，按照軟件固定的格式循環(huán)，實(shí)現(xiàn)每個(gè)單元定義為一個(gè)單元集合的模式，生成任意單元集合FE，將結(jié)果輸出到文本文件1.txt。

(3)定義單元材料截面屬性。在Fortran程序中將鄧肯張E-B模型材料屬性的名稱ME分配給相應(yīng)的單元集合FE，按照軟件固定的格式進(jìn)行循環(huán)賦值，實(shí)現(xiàn)每個(gè)單元集合對應(yīng)一個(gè)截面屬性，將結(jié)果輸出到文本文件2.txt。

(4)定義單元材料特性。讀取2.2節(jié)每個(gè)單元隨機(jī)抽樣的密度ρi與鄧肯張E-B模型參數(shù)ci，φ0i，Rfi，Ki，mi，ni，Kbi，Δφi，Kuri，Pa(Pa為大氣壓，Pa=100 kPa)。對材料名稱ME、密度ρ和鄧肯張E-B模型參數(shù)進(jìn)行三重嵌套循環(huán)，實(shí)現(xiàn)對有限元模型任意單元的力學(xué)模型參數(shù)賦值，將結(jié)果寫入文本文件3.txt。

(5)將文本文件1.txt、2.txt、3.txt合并為一個(gè)inp文件，通過ABAQUS軟件直接調(diào)用不同抽樣結(jié)果下的inp文件進(jìn)行有限元計(jì)算。

3 空間差異性下的壩體應(yīng)力、變形達(dá)標(biāo)判斷準(zhǔn)則

在進(jìn)行壩體應(yīng)力、變形有限元計(jì)算時(shí)，隨著壩體的填筑，荷載增量逐步施加，部分壩體單元的應(yīng)力結(jié)果可能會(huì)超過其材料的極限應(yīng)力而產(chǎn)生破壞?；谏鲜龅膲误w破壞方式和有限元計(jì)算方法，在進(jìn)行壩料參數(shù)不確定時(shí)的壩體應(yīng)力、變形有限元計(jì)算時(shí)，可以根據(jù)壩體材料參數(shù)在選取設(shè)計(jì)值時(shí)的確定性計(jì)算結(jié)果作為每一次循環(huán)模擬計(jì)算結(jié)果是否達(dá)標(biāo)的判斷標(biāo)準(zhǔn)。

根據(jù)隨機(jī)有限元每次的循環(huán)計(jì)算結(jié)果，得出每次計(jì)算的壩體最大沉降、上下游水平位移最大值、大小主應(yīng)力最大值,分別與材料參數(shù)設(shè)計(jì)值計(jì)算得出的相應(yīng)最值進(jìn)行比較，若小于等于設(shè)計(jì)階段的計(jì)算結(jié)果，則說明達(dá)標(biāo)，否則，說明不達(dá)標(biāo)。判斷式如式(2)所示。

(2)

根據(jù)前述判斷準(zhǔn)則確定壩體應(yīng)力、變形最大值中沒有超過設(shè)計(jì)情況計(jì)算結(jié)果的次數(shù)m占總有限元模擬次數(shù)M的比例，由此確定在材料參數(shù)不確定性下的壩體應(yīng)力、變形達(dá)標(biāo)概率P，如式(3)所示。進(jìn)一步可以得出超標(biāo)概率P′，如式(4)所示。

P=m/M；

(3)

P′=1-P。

(4)

4 工程實(shí)例分析

某瀝青混凝土心墻壩最大壩高82.5 m，上下游壩坡比均為1∶1.8，采用漸變式瀝青混凝土心墻，在心墻上下游側(cè)設(shè)置有厚度為4.5 m的過渡層。圖2中標(biāo)明了大壩材料分區(qū)，包括堆石料Ⅰ區(qū)、堆石料Ⅱ區(qū)、過渡層區(qū)、瀝青混凝土心墻區(qū)。壩址區(qū)河床存在厚度為42～50 m的覆蓋層。

圖2 大壩材料分區(qū)與筑壩順序Fig.2 Materials partition and building sequence of dam

圖3 大壩整體三維有限元網(wǎng)格剖分Fig.3 Three-dimensional finite element meshing of the whole dam

4.1 計(jì)算模型

利用ABAQUS有限元模擬軟件提供的C3D8單元對模型進(jìn)行網(wǎng)格剖分，共剖分27 142個(gè)單元，大壩三維整體網(wǎng)格剖分見圖3?；A(chǔ)邊界采用截?cái)噙x取，豎直方向向下截取50 m，并在其底部施加固定位移約束；水平向截?cái)嚅L度為50 m，并在其截?cái)嗝嫔鲜┘庸潭ㄎ灰萍s束。通過壩基地應(yīng)力平衡和9級(jí)加載步模擬壩體的填筑過程(見圖2)。計(jì)算中規(guī)定沿壩軸向從右岸到左岸為x坐標(biāo)正向，沿壩體高程方向規(guī)定為y坐標(biāo)正向，垂直于壩軸線從上游到下游為z坐標(biāo)正向。

表1 材料隨機(jī)參數(shù)統(tǒng)計(jì)特性Table 1 Random properties of dam material parameters

4.2 隨機(jī)有限元計(jì)算結(jié)果分析

4.2.1 材料隨機(jī)參數(shù)的統(tǒng)計(jì)特性

根據(jù)陳輝等[12]先前對堆石料材料參數(shù)不確定性的研究，結(jié)合現(xiàn)有土石壩填筑標(biāo)準(zhǔn)和相似工程現(xiàn)場檢測資料，擬定不同分區(qū)鄧肯張E-B模型材料參數(shù)的變異系數(shù)(即標(biāo)準(zhǔn)差/均值)如表1所示，各參數(shù)的期望值根據(jù)設(shè)計(jì)參數(shù)確定。由于現(xiàn)有的資料只有針對某堆石壩主堆石區(qū)變異系數(shù)的研究，并且考慮到變異系數(shù)只是反映數(shù)據(jù)的離散程度，本文想通過變異系數(shù)與該工程的均值真值求出標(biāo)準(zhǔn)差進(jìn)行抽樣，而且雖是不同分區(qū)但同為堆石料，變異系數(shù)差異性較小，所以不同分區(qū)采用了同一個(gè)變異系數(shù)。

由不同材料分區(qū)的8個(gè)參數(shù)統(tǒng)計(jì)特性，假定其服從一定的分布規(guī)律，對材料物理力學(xué)參數(shù)進(jìn)行隨機(jī)抽樣。文獻(xiàn)[12]統(tǒng)計(jì)了實(shí)際壓實(shí)質(zhì)量下堆石壩主堆石區(qū)的鄧肯張E-B模型參數(shù)概率分布，認(rèn)為鄧肯張E-B模型參數(shù)服從正態(tài)或?qū)?shù)正態(tài)分布。

4.2.2 壩體隨機(jī)有限元計(jì)算結(jié)果分析

繪制出由蒙特卡羅模擬得到的壩體最大沉降的滑動(dòng)平均值隨蒙特卡羅模擬次數(shù)的變化曲線，如圖4所示。可知當(dāng)蒙特卡羅模擬次數(shù)達(dá)到300次后計(jì)算結(jié)果已經(jīng)較為平穩(wěn)，因此下文隨機(jī)模擬次數(shù)取足夠大的600次。

圖5是設(shè)計(jì)工況下壩體應(yīng)力、變形等值線。結(jié)合有限元計(jì)算結(jié)果可知，壩體竣工期最大沉降為1.130 m，占?jí)胃吲c覆蓋層厚度的0.85%，位于壩體的心墻中下部。壩體向上游水平位移最大值為0.412 m，向下游水平位移最大值為0.383 m，分布較為對稱，沒有呈現(xiàn)明顯的不均勻性。壩體大主應(yīng)力從壩頂向壩基呈現(xiàn)出逐漸增大的趨勢，其最大值為2 530 kPa，位于大壩底部。

注：頻數(shù)為變形或應(yīng)力在600次隨機(jī)有限元計(jì)算中出現(xiàn)的次數(shù)。圖6 壩體應(yīng)力、變形概率分布Fig.6 Probability distribution of dam stress and deformation

圖6為壩體應(yīng)力、變形概率分布。隨機(jī)有限元計(jì)算得到的壩體最大沉降概率分布如圖6(a)所示?？梢钥闯鲈诳紤]了堆石料與覆蓋層的空間差異性后，與設(shè)計(jì)情況相比壩體最大沉降超標(biāo)概率為50%，變異系數(shù)為9.6%(0.109/1.132×100%=9.6%)。隨機(jī)統(tǒng)計(jì)結(jié)果的最大沉降為1.43 m，占?jí)胃吲c覆蓋層厚度的1.0%(壩高與覆蓋層厚度為82.5 m+50 m=132.5 m，1.43/132.5×100%=1%)，雖然其發(fā)生概率僅有0.17%(通過600次隨機(jī)有限元計(jì)算，統(tǒng)計(jì)最大沉降1.43 m出現(xiàn)的次數(shù)與600的比值可得)，但仍然會(huì)對壩體結(jié)構(gòu)安全造成一定的威脅。對壩體沉降最大值進(jìn)行正態(tài)分布K-S檢驗(yàn)，sig值為0.744。由此可知壩體沉降最大值服從均值為1.132、標(biāo)準(zhǔn)差為0.109的正態(tài)分布。(sig值為統(tǒng)計(jì)顯著性，由SPSS數(shù)據(jù)分析軟件直接求得，在進(jìn)行正態(tài)分布K-S檢驗(yàn)時(shí)，若sig值>0.05，該樣本服從正態(tài)分布，否則不服從正態(tài)分布。)沉降樣本最大值1.43 m落在3σ區(qū)間內(nèi)，因此在設(shè)計(jì)階段很有必要考慮壩體沉降在3σ區(qū)間內(nèi)的情況，對于本工程即考慮壩體最大沉降范圍在0.803～1.457 m時(shí)壩體的應(yīng)力、變形情況。

壩體的順河向最大水平位移及大主應(yīng)力最大值的概率分布如圖6(b)—圖6(d)所示，最終統(tǒng)計(jì)結(jié)果如表2?？煽闯隹紤]筑壩材料與覆蓋層料的不確定性確實(shí)會(huì)使大壩應(yīng)力、變形產(chǎn)生一定程度的離散。假如在設(shè)計(jì)階段忽略這種不確定性，而仍用確定性的分析方法對結(jié)構(gòu)進(jìn)行分析，則有50%的可能性會(huì)低估壩體最大沉降，46%～47%的可能性低估上下游最大水平位移，43%的可能性低估壩體大主應(yīng)力最大值，51%的可能性會(huì)低估壩體小主應(yīng)力最大值。

表2 壩體應(yīng)力、變形最大值隨機(jī)有限元分析結(jié)果統(tǒng)計(jì)Table 2 Statistical results of dam stress and deformation by stochastic finite element analysis

4.2.3 心墻隨機(jī)有限元計(jì)算結(jié)果分析

設(shè)計(jì)工況下心墻應(yīng)力、變形分布如圖7所示。結(jié)合有限元計(jì)算結(jié)果可知，心墻竣工期豎直沉降最大值為1.130 m，發(fā)生于心墻中下部。大主應(yīng)力分布呈現(xiàn)出從頂部到底部逐漸增加的趨勢，最大值為1 740 kPa，位于心墻底部位置，且為壓應(yīng)力，說明心墻受力狀態(tài)良好。圖8是心墻順河向水平位移沿相對高程分布情況，可知設(shè)計(jì)工況下向上游水平位移最大值為0.050 m，主要位于心墻下部，向下游水平位移最大值為0.025 m，主要位于心墻上部。而隨機(jī)有限元模擬向上游水平位移最大值均值為0.052 m，向下游水平位移最大值均值為0.023 m，模擬均值的分布情況與設(shè)計(jì)值相同。

圖7 設(shè)計(jì)工況下心墻應(yīng)力、變形等值線Fig.7 Stress and deformation of core wall in design condition

圖8 心墻順河向水平位移沿高度分布Fig.8 Horizontal displacement of core wall along the river against elevation

圖9 心墻特征點(diǎn)位置Fig.9 Typical locations on core wall

在心墻典型剖面的上中下部位各選擇一個(gè)特征點(diǎn)(有限元單元編號(hào)分別為536，5724，1844)對其應(yīng)力、變形分布規(guī)律進(jìn)行分析。特征點(diǎn)選取位置如圖9所示。由于心墻向上下游方向水平位移較小，對大壩的正常運(yùn)營不會(huì)產(chǎn)生影響，在此不再分析心墻順河向水平位移的隨機(jī)有限元分析結(jié)果。

圖10 編號(hào)5724特征點(diǎn)應(yīng)力、變形概率分布Fig.10 Probability distribution of stress and deformation at characteristic point No.5724

心墻上3個(gè)不同部位的特征點(diǎn)豎直沉降在以設(shè)計(jì)工況為判斷標(biāo)準(zhǔn)時(shí)，都是有50%的超標(biāo)概率，規(guī)律比較明顯。相比之下，心墻的大、小主應(yīng)力也有不同程度的超標(biāo)，但無明顯規(guī)律。從3個(gè)點(diǎn)應(yīng)力、變形結(jié)果的整體來看， 3個(gè)特征點(diǎn)位置的不同，心墻上部的特征點(diǎn)超標(biāo)概率較小，而心墻下部的特征點(diǎn)超標(biāo)概率較大。從確定性有限元分析可知，心墻的豎直沉降與大小主應(yīng)力最大區(qū)域更靠近心墻中下部，因此分析這種情況可能是由于下部點(diǎn)距離心墻的豎直沉降和大小主應(yīng)力變化最大的區(qū)域更近導(dǎo)致的，所以在心墻施工時(shí)應(yīng)對其中下部進(jìn)行嚴(yán)格的質(zhì)量控制。圖10是編號(hào)5724特征點(diǎn)隨機(jī)有限元計(jì)算得到的部分應(yīng)力、變形概率分布。心墻特征點(diǎn)應(yīng)力、變形隨機(jī)有限元結(jié)果統(tǒng)計(jì)如表3所示。

表3 心墻特征點(diǎn)應(yīng)力、變形隨機(jī)有限元分析 結(jié)果統(tǒng)計(jì)Table 3 Statistical results of stress and deformation at characteristic points of core wall by stochastic finite element analysis

這一現(xiàn)象也可以從536與1844特征點(diǎn)的豎直沉降不服從正態(tài)分布看出。由于這2個(gè)特征點(diǎn)位距離心墻豎直沉降最大區(qū)域較遠(yuǎn)，所以其沉降樣本點(diǎn)取值多等于該點(diǎn)沉降均值，沒有產(chǎn)生過度的離散，因而不服從正態(tài)分布?？梢婋m然本構(gòu)模型參數(shù)服從正態(tài)分布，但其應(yīng)力、變形結(jié)果可能因研究位置的不同，概率分布也會(huì)不同。

5 結(jié) 論

本文針對大壩實(shí)際施工過程造成的材料參數(shù)空間差異對瀝青混凝土心墻壩應(yīng)力、變形存在一定影響的問題，給出了隨機(jī)有限元分析的流程步驟和實(shí)現(xiàn)空間差異性計(jì)算的方法，并且結(jié)合工程實(shí)例，利用蒙特卡羅概率設(shè)計(jì)方法探討了在考慮堆石料與覆蓋層的空間差異性后，壩體和心墻特征點(diǎn)的應(yīng)力、變形統(tǒng)計(jì)規(guī)律。具體結(jié)論如下：

(1)基于隨機(jī)有限元分析各仿真循環(huán)相對獨(dú)立的特點(diǎn)，提出了在有限元軟件中實(shí)現(xiàn)材料參數(shù)空間差異性的方法，并且結(jié)合壩體破壞方式給出了判斷壩體隨機(jī)有限元結(jié)果達(dá)標(biāo)的判斷準(zhǔn)則。

(2)考慮壩體材料的空間差異性確實(shí)會(huì)使大壩應(yīng)力、變形發(fā)生一定程度的離散。而如果忽略這種差異性，仍用確定性分析方法對結(jié)構(gòu)進(jìn)行分析，則有50%左右的可能性會(huì)低估壩體的應(yīng)力、變形。其中沉降最大值有可能會(huì)對壩體正常服役產(chǎn)生威脅，所以很有必要分析在隨機(jī)有限元計(jì)算結(jié)果的最大值下，壩體結(jié)構(gòu)是否仍然安全，從而指導(dǎo)工程設(shè)計(jì)。由于壩體上下游水平位移變異系數(shù)較小且大小主應(yīng)力全為壓應(yīng)力，因此二者對壩體的安全運(yùn)行威脅較小。

(3)空間差異性下瀝青混凝土心墻上部主應(yīng)力的超標(biāo)概率小于下部主應(yīng)力超標(biāo)概率，不同位置處沉降的超標(biāo)概率都為50%。且表征大壩或心墻的性能指標(biāo)不一定都服從正態(tài)分布，即使該指標(biāo)的最值服從正態(tài)分布，但有可能因?yàn)橹笜?biāo)考察的結(jié)構(gòu)位置不同而使其概率分布規(guī)律改變。