空間差異性下瀝青混凝土心墻壩應(yīng)力與變形的隨機有限元分析

郭 晴1,，劉東海，陳 輝

(1.河北工程大學 水利水電學院，河北 邯鄲 056021； 2.天津大學 水利工程仿真與安全國家重點實驗室，天津 300350 )

1 研究背景

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

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

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

2 瀝青心墻壩隨機有限元分析方法

2.1 壩體隨機有限元計算流程

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

(1)

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

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

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

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

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

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

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

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

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

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

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

3 空間差異性下的壩體應(yīng)力、變形達標判斷準則

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

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

(2)

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

P=m/M；

(3)

P′=1-P。

(4)

4 工程實例分析

某瀝青混凝土心墻壩最大壩高82.5 m，上下游壩坡比均為1∶1.8，采用漸變式瀝青混凝土心墻，在心墻上下游側(cè)設(shè)置有厚度為4.5 m的過渡層。圖2中標明了大壩材料分區(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 計算模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖7 設(shè)計工況下心墻應(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 心墻特征點位置Fig.9 Typical locations on core wall

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

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

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

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

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

5 結(jié) 論

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

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

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

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