, ,,

(三峽大學(xué) a.水利與環(huán)境學(xué)院; b.土木與建筑學(xué)院，湖北 宜昌 443002)

1 研究背景

有限單元法在土石壩的應(yīng)力應(yīng)變分析中得到了廣泛的應(yīng)用[1]，同時(shí)也為國內(nèi)外重大工程的建設(shè)提供了理論指導(dǎo)[2-5]。由于土體的應(yīng)力應(yīng)變關(guān)系具有明顯的非線性特性[6]，在進(jìn)行土石壩的有限元分析中一般采用增量法來對施工全過程進(jìn)行仿真分析，這種方式不僅能反映出施工過程中各階段的應(yīng)力應(yīng)變情況，而且能模擬土石壩的施工填筑及蓄水運(yùn)行的全過程。

采用有限單元法對土石壩施工填筑過程進(jìn)行數(shù)值仿真時(shí)，整個(gè)分析對象(壩體、地基等)的計(jì)算網(wǎng)格是一次性建成的，根據(jù)實(shí)際工程中壩體分區(qū)施工次序來對施工單元和未施工單元進(jìn)行不同處理，從而達(dá)到施工過程的模擬。大型商用軟件往往采用生死單元的方式[7-8]，通過控制單元的施工(激活)并一次性施加自重荷載來進(jìn)行施工過程的模擬，這種處理方式的計(jì)算效率不高，而且在激活施工單元時(shí)需要對其位移進(jìn)行處理，應(yīng)該說不是一種最佳的施工填筑模擬方法。

同時(shí)，在土石壩施工填筑仿真模擬中，考慮到壩體高度、仿真精度及計(jì)算效率的因素，計(jì)算模擬的施工填筑不可能按照實(shí)際填筑層來模擬，往往將幾個(gè)實(shí)際填筑層合并成一個(gè)計(jì)算填筑層，整個(gè)壩體分為十幾個(gè)或幾十個(gè)填筑層來進(jìn)行模擬。而在采用有限單元法對施工填筑級進(jìn)行模擬時(shí)，施工填筑單元的自重荷載是一次性施加的，而實(shí)際施工時(shí)填筑級是一層層碾壓逐級加載的，這2種方式下土壩的變形機(jī)理是不同的。在不考慮土體固結(jié)等時(shí)間因素對變形影響的前提下(即假定新填筑層的變形在施工過程中瞬時(shí)完成)，已填筑層的自重不影響其上部未施工土層的變形。填筑層內(nèi)某一點(diǎn)土體的變形來自于該點(diǎn)上方土體的壓力，該點(diǎn)下面的土體重力對其變形沒有影響，則每層填筑層頂部位移應(yīng)該為零。但按照施工填筑單元的自重荷載一次性施加的情況下，填筑層頂部的位移最大，并由此最終得到的累積豎向位移等值線出現(xiàn)不合理的鋸齒狀分布的現(xiàn)象，這樣就導(dǎo)致填筑仿真得到的位移場與實(shí)際不一致。針對這一問題，國內(nèi)外學(xué)者[8-11]提出在數(shù)值仿真計(jì)算位移的基礎(chǔ)上乘以修正系數(shù)的方法將其改善，但是對修正公式適用性的研究不完善。

本文基于上述研究背景，對有限單元法在土石壩施工填筑數(shù)值仿真中的施工填筑過程模擬及填筑位移修正進(jìn)行了研究，提出了一種對施工仿真單元的處理方式，并對填筑位移修正方法進(jìn)行了完善。

2 施工填筑中施工單元的處理

基于有限單元法的施工填筑仿真是根據(jù)實(shí)際工程中壩體分區(qū)施工次序來對施工單元和未施工單元進(jìn)行不同處理，從而達(dá)到施工過程的模擬。

目前對于該問題，采用較多的是“生死單元”[7-8]法。在進(jìn)行施工仿真前，分析模型的有限元網(wǎng)格在一次性準(zhǔn)備好的前提下，先對所有單元進(jìn)行“殺死”，然后根據(jù)填筑次序，逐步地對填筑區(qū)內(nèi)的單元進(jìn)行“激活”，從而達(dá)到對施工過程的真實(shí)模擬。對于“殺死”的單元，計(jì)算程序并不是將這些單元從計(jì)算模型中刪除，而是將其剛度矩陣乘以一個(gè)很小的因子，同時(shí)單元的荷載全部變?yōu)榱悖蝗绻麊卧_始被施工，則單元被重新激活，程序?qū)⑵鋭偠然謴?fù)，并施加作用于單元上的所有荷載。有些商用軟件在重新激活單元時(shí)，被激活的單元由于其邊界上有非零位移，這時(shí)軟件會把非零位移作為已知位移條件施加在新激活單元的邊界上，這樣就產(chǎn)生了剛激活而未施加任何荷載就產(chǎn)生了位移和應(yīng)力的現(xiàn)象，這顯然是不符合實(shí)際情況。

國內(nèi)學(xué)者在對大體積混凝土進(jìn)行施工期溫度場仿真時(shí)[12]，為了對施工過程進(jìn)行模擬，曾引入“空氣單元”的方法來模擬未修建的部分。對于未施工單元，“空氣單元”的彈性模量取相鄰單元材料彈性模量的1/100～1/1 000，泊松比取0.49，但是對其計(jì)算得到的位移不予以累計(jì)。單元施工時(shí)，將其材料參數(shù)恢復(fù)，此后開始累積其位移和應(yīng)力作為計(jì)算成果。

“生死單元”法和“空氣單元”法作為對施工過程中處理單元的2種模擬方法，其處理方式基本是一樣的。從上述2種方法的處理方式可以看出，未施工的單元只是單元剛度矩陣和荷載被賦予很小的值，但仍然參與整個(gè)計(jì)算過程，施工填筑過程每個(gè)荷載步都需要進(jìn)行包含未施工單元在內(nèi)的所有單元剛度矩陣的計(jì)算和組裝，以及求解包含未施工單元自由度在內(nèi)的所有自由度的平衡方程，乃至隨后所有單元的應(yīng)力計(jì)算。這種處理方法在整個(gè)施工過程中方程的求解上，方程的總階數(shù)是相同的。但隨著計(jì)算規(guī)模的增大，計(jì)算效率將會大大降低。

在某一荷載步內(nèi)，對于未施工單元，其單元自由度和荷載是不存在的，可以不參與到計(jì)算過程中。因此，本文建議的處理方式是：在土石壩施工填筑過程的有限元模擬中，對每個(gè)單元賦予1個(gè)施工信息或者施工序號，來規(guī)劃施工過程，并判斷單元是否被施工。對于已施工單元，對其自由度進(jìn)行考慮，對其荷載進(jìn)行處理，并參與有限元計(jì)算；對未施工單元，不考慮其自由度，對其荷載也不用處理，方程求解的階數(shù)只是施工單元的自由度。從計(jì)算工作量上講，這種處理方式不用計(jì)算未施工單元的剛度矩陣和應(yīng)力，不用構(gòu)建未施工單元的自由度方程；從計(jì)算結(jié)果來看，由于未施工單元的自由度沒有被納入計(jì)算，不管是位移還是應(yīng)力，均為0，符合實(shí)際情況。與“生死單元”法相比，這種處理方式只是在每個(gè)荷載步內(nèi)需要根據(jù)單元的施工信息對自由度進(jìn)行統(tǒng)計(jì)，得到自由度指示數(shù)組，從計(jì)算時(shí)間上來看，這個(gè)統(tǒng)計(jì)工作與未施工單元參與計(jì)算過程相比不值一提。因此，從計(jì)算效率上來講，該方法比“生死單元”法有顯著的提升。

可粗略地估計(jì)這種處理方式對計(jì)算效率的提高效果，假定對于線彈性問題，初步假定每個(gè)單元的計(jì)算分析所占的時(shí)間一致，設(shè)為t，計(jì)算模型共有NE個(gè)單元，分NB個(gè)荷載步進(jìn)行計(jì)算分析，每個(gè)荷載步的施工單元數(shù)相同，即NE/NB。對于采用“生死單元”法來處理施工過程，其需要的計(jì)算時(shí)間為

T1=NE×t×NB。

(1)

采用本文提出的方法來模擬施工過程，其需要的計(jì)算時(shí)間為

(2)

從上式可以看出，采用本文建議的方法，可節(jié)省計(jì)算時(shí)間接近50%。

為了驗(yàn)證本文提出的施工單元處理方式與“生死單元”在計(jì)算效率上的差別，選用某面板堆石壩工程填筑期仿真計(jì)算進(jìn)行算例比較，該工程計(jì)算模型節(jié)點(diǎn)總數(shù)為50 294，單元總數(shù)為48 250，填筑期采用52級進(jìn)行模擬，計(jì)算平臺為HP Z400工作站 (處理器為至強(qiáng)W3670@3.20G六核)。采用本文提出方法計(jì)算時(shí)間為28 min，采用“生死單元”的處理方式計(jì)算所需時(shí)間為50 min，由此可見本文提出的處理方式與“生死單元”法相比，有較高的計(jì)算效率，大約可節(jié)省50%的計(jì)算時(shí)間，其與每個(gè)荷載步內(nèi)新施工單元數(shù)及每個(gè)單元的計(jì)算所需時(shí)間相關(guān)。

3 施工填筑位移修正

在土石壩施工填筑仿真模擬中，往往將整個(gè)壩體分若干個(gè)填筑層來模擬，每級單元施工時(shí)，施工填筑單元的自重荷載是一次性施加的，而實(shí)際施工時(shí)填筑級是一層層碾壓逐級加載的，這2種方式下土壩的變形機(jī)理是不同的：自重荷載一次性施加時(shí)，豎向某點(diǎn)的沉降為該點(diǎn)下部土層的壓縮量，這時(shí)頂部位移最大；而實(shí)際施工的逐級加載時(shí)，豎向某點(diǎn)的沉降為該點(diǎn)以上土重所引起的該點(diǎn)以下土層的壓縮量。Clough和Woodward[13]依據(jù)一次性加載和逐步加載的機(jī)理不同，從理論上推導(dǎo)了2種方式的沉降量，首次提出了在有限元模擬時(shí)一次性加載得到計(jì)算結(jié)果的基礎(chǔ)上乘以修正系數(shù)的方式來得到逐步多級加載的變形效果。

3.1 單填筑層的位移修正

假定土體為彈性體，考慮體積壓縮系數(shù)mv為常量情況下，分析一次性加載和逐步加載2種不同作用下豎向位移的理論值。

如圖1所示，某一填筑層的高度為h，A點(diǎn)位于壩頂下z深度，A點(diǎn)的沉降是A點(diǎn)以下厚度為h-z的土層AB的壓縮量。

圖1 單填筑層壩體豎向沉降受力示意圖

一次性加載條件下，引起壓縮的應(yīng)力為自重應(yīng)力，在壓縮層AB內(nèi)自重應(yīng)力呈梯形分布，如圖1中ABCD所示。根據(jù)土體單向壓縮沉降理論，一次性加載條件下A點(diǎn)的沉降為

(3)

對逐級加載來說，當(dāng)填土到A點(diǎn)的高度時(shí)，A點(diǎn)以下土體自重所引起的沉降已在填土過程中完成，A點(diǎn)的沉降為零；當(dāng)填土超過A點(diǎn)的高度時(shí)，A點(diǎn)才有沉降。此時(shí)，A點(diǎn)的沉降是A點(diǎn)以上土層OA的土重所引起的A以下土層AB的壓縮，引起A點(diǎn)沉降的應(yīng)力在壓縮層AB內(nèi)呈矩形分布，如圖1中ABED所示。則逐級加載機(jī)理下A點(diǎn)的沉降為

(4)

(5)

3.2 多填筑層的位移修正

上述修正公式只適用于單填筑層，并不適用于多填筑層情況，而在土石壩施工填筑仿真模擬中，往往將壩體分為十幾個(gè)或幾十個(gè)填筑層來進(jìn)行模擬，如果按照上述公式來修正是會存在一定誤差。

在多填筑層下，引起壩頂以下深度為z的某點(diǎn)產(chǎn)生沉降位移不僅僅是在本填筑層會發(fā)生沉降，已填筑層均會發(fā)生沉降。如圖2所示，A點(diǎn)的沉降為新填土層AB的沉降和已填土層BB′的沉降,即可表示為

(6)

式中：z為A點(diǎn)距新填土層頂部O點(diǎn)的深度；h為新填土層的厚度;H為已填筑土層的厚度;σz為引起A點(diǎn)土體沉降的應(yīng)力。

圖2 多填筑層壩體豎向沉降受力示意圖

一次加載條件下，新填土層AB內(nèi)引起沉降的壓縮應(yīng)力呈梯形分布，如圖中ABCD所示；已填土層 內(nèi)的壓縮應(yīng)力呈矩形分布，如圖中BB′C′C所示。則A點(diǎn)的沉降可表示為

mvγ(h2-z2)/2+mvγhH。

(7)

對逐級加載來說，A點(diǎn)的沉降是A點(diǎn)以上土層OA的土重所引起的A以下土層AB′的壓縮，引起A點(diǎn)沉降的應(yīng)力在壓縮層AB′內(nèi)呈矩形分布，如圖2中AB′E′D所示。則逐級加載機(jī)理下A點(diǎn)的沉降為

mvγz(h-z)+mvγzH。

(8)

則其修正系數(shù)如下式所示：

(9)

從以上分析可以看出，單填筑層中新填土層中某點(diǎn)的沉降只有該點(diǎn)以下土體的沉降，而多填筑層新填土層中某點(diǎn)的沉降包括已填筑土層在該填筑層作用下的沉降和新填土層中該點(diǎn)以下土體的壓縮變形2部分，由此，其位移修正公式有所不同。因此，施工填筑位移的修正公式如下：

(10)

表1 壩體填筑材料計(jì)算參數(shù)

式中：i為填筑級；z為計(jì)算沉降點(diǎn)至填筑頂高;h為本級填筑高度;H為已填筑高度;N為總填筑層數(shù)。

Clough和Woodward[13]提出的修正公式(如式(5)),只考慮了單填筑層的情況，而未考慮多填筑層時(shí)已填筑土層在新填筑土層作用下會產(chǎn)生沉降；而文獻(xiàn)[10]在研究位移場修正時(shí)，直接針對多填筑層模型來分析，沒有區(qū)分單填筑層和多填筑層的不同所引起的修正系數(shù)的差異。

從上述分析可以看出當(dāng)已填筑層厚度H=0時(shí)，式(5)和式(9)一致，同時(shí)當(dāng)已填筑層較厚時(shí)(即h/H→0)，式(9)可簡化如下：

(11)

為了分析上式簡化及考慮多填筑層與否所引起的誤差，假定修正系數(shù)為β，將式(9)按h/H的不同，繪制修正系數(shù)沿填筑層厚度的分布如圖3所示。從圖中可以看出，當(dāng)新填土層與已填土層厚度之比越大時(shí)(即h/H→∞)，多填筑層的修正系數(shù)趨于單填筑層的修正系數(shù)(β→2z/(h+z))；而當(dāng)已填筑層較厚時(shí)(即h/H→0)，多填筑層的修正系數(shù)β→2z/(h+z)。

圖3 修正系數(shù)沿土層深度分布

如果采用式(11)來簡化代替式(9)，隨著填筑層的增加，其誤差越來越小，修正系數(shù)在第二填筑層時(shí)誤差最大，在土石壩施工填筑仿真中，考慮到施工強(qiáng)度的因素，填筑層厚較為均勻，即在第二填筑層時(shí)，h/H=1.0，此時(shí)的簡化所引起的最大誤差為5%，而且隨著填筑層的增加而逐漸減小。

如果不考慮施工填筑層的影響，如文獻(xiàn)[10]和[13]所述，采用式(5)或式(11)，其最大誤差可達(dá)到17%，因此在進(jìn)行位移修正時(shí)，需要考慮填筑層的差別，建議在使用有限單元法進(jìn)行土石壩施工填筑仿真時(shí)位移修正系數(shù)如下式所示：

(12)

4 算例應(yīng)用

選用文獻(xiàn)[8]和[10]的算例：一心墻堆石壩，壩高100 m，頂寬10 m，心墻頂寬6 m，壩體上下游壩坡坡比為1∶2.0，心墻坡比為1∶0.2；網(wǎng)格沿壩高分為20層，每層厚度為5 m，壩體和心墻的本構(gòu)模型采用鄧肯E-B模型，材料參數(shù)如表1所示，計(jì)算網(wǎng)格如圖4所示。

圖4 有限元網(wǎng)格圖

為了比較位移修正考慮單、多填筑層的差別以及式(12)簡化前后的誤差，擬定以下5個(gè)計(jì)算方案：

(1) 施工填筑位移不修正，稱為未修正方案。

(2) 根據(jù)單填筑層推導(dǎo)出來的修正公式[9,13](式(5))，稱為已有修正方案Ⅰ。

(3) 根據(jù)多層填筑推導(dǎo)出來的修正公式[10](式(11))，稱為已有修正方案Ⅱ。

(4) 考慮填筑層不同的修正公式(式(10))，稱為本文修正方案。

(5) 考慮了填筑層不同的簡化修正公式(式(12))，稱為簡化修正方案。

以上5種計(jì)算方案的填筑計(jì)劃均按10級填筑，每級兩層單元。

按照上述施工填筑方案得到的豎向位移分布沿壩軸線對稱，因而不同修正方案的豎向位移的等值線見圖5所示。不同修正方案的位移統(tǒng)計(jì)成果見表2所示。

圖5 不同修正方案的豎向位移

從未修正方案的豎向位移等值線圖可以看出，由于每填筑層的頂部位移不為0，致使累加位移產(chǎn)生如圖所示的鋸齒狀，采用修正方案I修正后，雖有所改善，但由于該方案沒有考慮已填筑土層在該填筑層作用下的沉降這一因素，致使修正后的位移圖仍有波動性。修正方案Ⅱ、本文修正方案及簡化修正方案的修正效果基本一致，如圖5(c)所示，這與實(shí)際情況較為相符，修正方案Ⅱ雖然沒有考慮第一個(gè)填筑層的修正系數(shù)的不同，但由于后續(xù)填筑層中考慮了占絕對權(quán)重的已填筑土層在該填筑層作用下的沉降(牽移位移)，因此修正效果與本文修正方案相同，只是第一個(gè)填筑層內(nèi)部點(diǎn)的修正沉降位移偏小。

表2 不同方案豎向位移計(jì)算結(jié)果

從不同方案的計(jì)算結(jié)果來看，未修正方案的豎向位移最大，達(dá)到0.744 m，但由于每填筑層的頂部位移不為0，致使最終壩頂豎向位移達(dá)到0.192 m的沉降，這顯然不符合實(shí)際。修正方案Ⅰ雖然對每個(gè)填筑層的頂部位移進(jìn)行了修正，但沒有考慮已填筑土層在該填筑層作用下的沉降，因此修正后豎向位移偏大。修正方案Ⅱ、本文修正方案及簡化修正方案的豎向位移基本相同，只是由于由式(9)簡化成式(11)后引起較小的差別。

從計(jì)算結(jié)果和修正后豎向位移等值線分布圖來看，由于修正方案II考慮了占絕對權(quán)重的已填筑土層在該填筑層作用下的沉降(牽移位移)，其計(jì)算結(jié)果和修正效果均與本文修正方案相同，但從理論上講，單填筑層和多填筑層的沉降土層不同，其位移修正公式應(yīng)有所區(qū)分，因此本文提出的修正方案及簡化修正方案是更為合理。

5 結(jié) 論

本文針對有限單元法在土石壩施工填筑數(shù)值仿真中的施工填筑仿真模擬及填筑位移修正進(jìn)行了研究，得到結(jié)論如下：

(1) 土石壩施工填筑數(shù)值仿真中，建議給每個(gè)單元賦予1個(gè)施工信息來進(jìn)行控制，并且未施工單元不參與計(jì)算，這種處理方式從計(jì)算效率上來講，比“生死單元”法有顯著的提升，可節(jié)省計(jì)算時(shí)間達(dá)到50%。

(2) 施工填筑位移修正時(shí)，應(yīng)考慮單填筑級和多填筑級沉降土層的不同，由此對于第一個(gè)填筑層與后續(xù)填筑層的位移修正公式應(yīng)有所區(qū)分。

(3) 在多填筑層位移修正時(shí)，采用β=z/h的簡化方式進(jìn)行處理時(shí)，所引起的最大誤差為5%，而且隨著填筑層的增加而逐漸減小，從修正效果和計(jì)算結(jié)果來看，這種簡化方式是可行的。

參考文獻(xiàn)：

[1] 陳慧遠(yuǎn).土石壩有限元分析[M].南京:河海大學(xué)出版社,1988.(CHEN Hui-yuan. Finite Element Analysis of Earth-Rock Dam[M]. Nanjing: Hohai University Press, 1988. (in Chinese))

[2] 孔憲京,鄒德高,徐 斌,等.紫坪鋪面板堆石壩三維有限元彈塑性分析[J].水力發(fā)電學(xué)報(bào),2013,(2):213-222.(KONG Xian-jing, ZOU De-gao, XU Bin,etal. Three-dimensional Finite Element Elasto-plastic Analysis of Zipingpu Concrete Faced Rock-fill Dam [J]. Journal of Hydroelectric Engineering,2013,(2):213-222.(in Chinese))

[3] 潘家軍,王觀琪,江 凌,等.基于ABAQUS的高混凝土面板堆石壩地震反應(yīng)三維非線性分析[J].水力發(fā)電學(xué)報(bào),2011,(6):80-84.(PAN Jia-jun, WANG Guan-qi, JIANG Ling,etal. ABAQUS Three-dimensional Nonlinear Analysis of Seismic Responses of High CFRD [J]. Journal of Hydroelectric Engineering, 2011, (6):80-84. (in Chinese))

[4] 劉萌成,高玉峰,劉漢龍.混凝土面板堆石壩應(yīng)力變形長期性狀有限元模擬[J].巖土力學(xué),2010,31(增1):412-418.(LIU Meng-cheng, GAO Yu-feng, LIU Han-long. Finite Element Analysis of Long-term Stress-deformation Behavior for Concrete-faced Rockfill Dam[J].Rock and Soil Mechanics,2010,31(Sup.1):412-418.(in Chinese))

[5] 楊啟貴,常曉林,周創(chuàng)兵,等.水布埡超高面板堆石壩變形控制方法研究[J].巖土力學(xué),2010,31(增2):247-253.(YANG Qi-gui,CHANG Xiao-lin,ZHOU Chuang-bing,etal.Study of Dam Deformation Control Method for Shuibuya High Concrete Faced Rockfill Dam [J].Rock and Soil Mechanics,2010,31(Sup.2):247-253.(in Chinese))

[6] 李廣信.高等土力學(xué)[M].北京:清華大學(xué)出版社,2004.(LI Guang-xin. Advanced Soil Mechanics[M]. Beijing: Tsinghua University Press, 2004. (in Chinese))

[7] 郝文化.ANSYS土木工程應(yīng)用實(shí)例[M].北京:中國水利水電出版社,2005.(HAO Wen-hua. ANSYS Application Examples in Civil Engineering[M]. Beijing: China Water Power Press, 2005. (in Chinese))

[8] 費(fèi) 康,劉漢龍.ABAQUS 的二次開發(fā)及在土石壩靜、動力分析中的應(yīng)用[J].巖土力學(xué),2010,31(3):881-890.(FEI Kang, LIU Han-long. Secondary Development of ABAQUS and Its Application to Static and Dynamic Analysis of Earth-rockfill Dam[J]. Rock and Soil Mechanics,2010,31(3):881-890. (in Chinese))

[9] 殷宗澤.土工原理[M].北京:中國水利水電出版社,2007.(YIN Zong-ze. Earthwork Principle[M]. Beijing: China Water Power Press, 2007. (in Chinese))

[10] 杜 斌，吳夢喜.土石壩分層填筑位移場修正的一種新方法[J].力學(xué)與實(shí)踐,2011,33(4):23-28. (DU Bin, WU Meng-xi. A New Approach for Modification of Displacement in Incremental Construction Analysis of Earth-rock Dam[J]. Mechanics in Engineering, 2011,33(4):23-28. (in Chinese))

[11] 王敏強(qiáng)，聶向珍.土(堆)石壩施工過程位移場有限元仿真分析[J].武漢大學(xué)學(xué)報(bào)(工學(xué)版),2003,36(2):74-78. (WANG Min-qiang, NIE Xiang-zhen. Displacement Simulation of Earth-rock Dam in Construction Process by Finite Element Method[J]. Engineering Journal of Wuhan University,2003,36(2):74-78. (in Chinese))

[12] 劉德富,黃達(dá)海,田 斌.拱壩封拱溫度場及溫控優(yōu)化[M].北京:中國水利水電出版社,2008.(LIU De-fu, HUANG Da-hai, TIAN Bin. Optimization of Joint Closure Temperature Field and Temperature Control Measures for Arch Dam[M]. Beijing: China Water Power Press, 2008. (in Chinese))

[13] CLOUGH R W, WOODWARD R J. Analysis of Embankment Stress and Deformations[J]. Journal of the Soil Mechanics and Foundation Division, ASCE,1967,93(SM4):529-549.