基于Abaqus的土工袋數(shù)值模擬算例分析

章致遠(yuǎn)

（江西省高安市水利局，江西 高安 330800）

近年來，有關(guān)土工袋技術(shù)在系統(tǒng)性理論研究、定量化分析計(jì)算、數(shù)值建模與預(yù)測(cè)以及大規(guī)?，F(xiàn)場(chǎng)試驗(yàn)等方面取得長足進(jìn)展。Recio J和OumeraciH[1]研究了采用土工充砂袋鋪筑的海岸結(jié)構(gòu)在穩(wěn)定水流作用下的有效變形問題，并基于實(shí)驗(yàn)室研究成果提出了土工袋技術(shù)砌筑結(jié)構(gòu)的水力特性及相應(yīng)概念模型。Hornsey W P等[2]人詳細(xì)研究了土工袋技術(shù)在海岸防護(hù)系統(tǒng)中穩(wěn)定性、耐久性和土工袋的壽命周期等問題及其在澳大利亞的實(shí)際應(yīng)用。Corbella S等[3]介紹了土工織物充砂袋在南非的應(yīng)用以及研究成果。高存貴等[4]提出了土工袋技術(shù)施工過程中的控制、檢測(cè)要求，并總結(jié)了土工袋技術(shù)的優(yōu)點(diǎn)；陶同康等[5]通過理論分析和試驗(yàn)研制了低價(jià)高摩擦性能的土工袋；劉斯宏[6～12]等人在土工袋力學(xué)性能研究的基礎(chǔ)上，采用土工袋技術(shù)治理南水北調(diào)工程輸水總干渠的膨脹土邊坡，深入分析了土工袋處理膨脹土的機(jī)理，為膨脹土處理開辟了新途徑。

目前針對(duì)土工袋的有限元分析方法主要有整體式、分離式兩種型式。前者是將袋子的張力作用等效為附加應(yīng)力作用在土骨架上，并將土體單元視為連續(xù)均勻的材料，這種方法無法揭示袋子與土之間相互作用的機(jī)理[13]；后者主要將土體單元與袋子分開采用不同的單元模型，各自劃分為單元，引進(jìn)接觸面單元實(shí)現(xiàn)模擬，采用這種模擬方法模擬多個(gè)土工袋時(shí)計(jì)算量巨大[13，14]。

采用有限元方法模擬土工袋涉及到土體和袋子單元的選取、土—袋子、袋子—袋子、袋子—結(jié)構(gòu)物之間的接觸設(shè)置、土體本構(gòu)的選擇等問題。下面結(jié)合前人研究成果與有限單元法元理論從模擬軟件介紹、土體本構(gòu)、袋子單元、接觸分析和算例計(jì)算5個(gè)方面概述土工袋有限元數(shù)值模擬。

1 Abaqus軟件概述

Abaqus是一款功能強(qiáng)大，操作便捷，適用性廣的有限元軟件。目前該軟件已經(jīng)被廣泛應(yīng)用于巖土、水利工程的數(shù)值計(jì)算中，有效性得到了實(shí)際工程的驗(yàn)證。Abaqus有限元軟件擁有眾多單元庫與模型庫，可以解決巖石力學(xué)、固體力學(xué)等不同領(lǐng)域的問題，特別是在非線性分析領(lǐng)域，可以解決復(fù)雜的工程實(shí)際問題，融結(jié)構(gòu)、傳熱學(xué)、流體以及熱固耦合、流固耦合于一體，具有駕馭龐大求解規(guī)模的能力。

2 土工袋模擬方法

在進(jìn)行土工袋有限元計(jì)算分析時(shí)，涉及袋內(nèi)材料的本構(gòu)模型、土工編織袋的模擬及土工袋層間接觸模擬等問題。受到袋內(nèi)材料和袋子材料的影響，土工袋呈現(xiàn)出的變形規(guī)律也各不相同。因此，在進(jìn)行土工袋有限元計(jì)算時(shí)，應(yīng)根據(jù)袋內(nèi)材料特性和袋子特性分別選取合適的計(jì)算參數(shù)。

2.1 本構(gòu)模型

土的本構(gòu)模型是影響有限元計(jì)算結(jié)果的一個(gè)重要因素[15～17]。土體變形影響因素很多，規(guī)律非常復(fù)雜。對(duì)于細(xì)砂材料，工程中常用鄧肯-張模型進(jìn)行有限元計(jì)算分析。對(duì)于本項(xiàng)目使用的灰渣，在進(jìn)行有限元計(jì)算分析時(shí)，工程中則常用Mohr-Coulomb彈塑性模型。

2.2 土工編織袋模擬方法

土工編織袋可以有兩種模擬方式：一種是將土工袋組合體簡化為一個(gè)等效的加筋土體；另一種是用桿單元模擬土工編織袋，袋子與袋中材料單獨(dú)模擬。

2.2.1 用組合體模擬土工袋

將土工袋組合體簡化為一個(gè)等效的加筋土體，加筋土體的抗剪強(qiáng)度中考慮袋子張力引起的附加凝聚力。由式（1）計(jì)算得到：

式中：c為粘聚力，kPa；T 為溫度，℃；Kp為被動(dòng)土壓力系數(shù)；H為水頭，m；

2.2.2 用桿單元模擬土工袋

土工編織袋具有能抗拉，但不能抗壓、彎、剪的普遍特性，可采用線彈性的桿單元進(jìn)行模擬，桿單元的勁度矩陣?yán)奂拥秸麄€(gè)土體結(jié)構(gòu)的整體勁度矩陣中。當(dāng)桿單元處于受壓狀態(tài)，則不計(jì)該桿單元的勁度貢獻(xiàn)。

2.3 土工袋層間接觸模擬

在結(jié)構(gòu)物中，如果兩種材料性質(zhì)相差很遠(yuǎn)，在一定的受力條件下有可能在其接觸面上產(chǎn)生錯(cuò)動(dòng)或開裂[7]。在土工袋組合體中就存在這樣3種形式的接觸面：袋子與袋子之間、袋子與結(jié)構(gòu)物之間以及袋內(nèi)材料與袋子之間的接觸面，這時(shí)可設(shè)置接觸單元進(jìn)行模擬。本項(xiàng)目擬采用Goodman無厚度接觸單元進(jìn)行模擬。

3 土工袋有限元分析算例

3.1 計(jì)算參數(shù)

為模擬該實(shí)際工程，選擇大土工袋尺寸為1 000mm×800mm×400mm（長×寬×高），小土工袋尺寸為500mm×800mm×200mm（長×寬×高），與設(shè)計(jì)采用的土工袋一致。在計(jì)算中土工袋的力學(xué)參數(shù)取用見表1，土工袋內(nèi)裝砂的物理力學(xué)參數(shù)取自拉姆火電廠地區(qū)的地質(zhì)工程勘探表（見表2）。

根據(jù)表2中土體的物理力學(xué)參數(shù)，結(jié)合以往工程經(jīng)驗(yàn)，選取鄧肯張模型的參數(shù)見表3。

計(jì)算過程中忽略編織袋與土體的接觸，只考慮袋子與袋子的接觸。為模擬編織袋之間的相互作用，在兩者之間設(shè)置無厚度的Goodman接觸單元，接觸面的參數(shù)如表4。

表1 黑色PP編織袋的性能指標(biāo)

表2 拉姆火電廠區(qū)域砂的物理力學(xué)參數(shù)

表3 拉姆火電廠區(qū)域砂鄧肯張模型參數(shù)

表4 接觸面單元參數(shù)

3.2 土工袋在ABAQUS中的模擬

算例中土工袋的模擬是根據(jù)實(shí)際工程中建議的大土工袋尺寸建立的模型，長度為100cm，高度為40cm。本次計(jì)算中用平面應(yīng)變單元模擬砂，用僅受張拉的桿單元和等效加筋兩種方法分別模擬土工袋。假設(shè)砂土和編織袋之間結(jié)合完好，砂土—編織袋之間不設(shè)置接觸單元，而在編織袋交接處設(shè)置接觸面單元。

3.3 土工袋無側(cè)限壓縮試驗(yàn)?zāi)M

土工袋無側(cè)限壓縮是對(duì)兩個(gè)上下堆疊的土工袋組合體施加豎向荷載，計(jì)算過程中在土工袋上、下方設(shè)置剛性加載板且約束下方加載板的位移和轉(zhuǎn)動(dòng)。加載板表面積遠(yuǎn)大于土工袋上、下表面，故忽略土工袋自重對(duì)壓縮的影響，將上層土工袋自重對(duì)下層土工袋的影響轉(zhuǎn)化為均布力施加在下土工袋頂部。為了模擬土工袋上部荷載的施加過程，計(jì)算中共分為10級(jí)加載，最后豎向荷載達(dá)到100kPa。

3.3.1 桿單元法有限元模型

土工袋模型由1 136個(gè)單元、1 034個(gè)節(jié)點(diǎn)組成，其中桿單元100個(gè)，接觸面單元120個(gè)，如圖1所示。圖中點(diǎn)V1～V6和H1～H6為位移監(jiān)測(cè)點(diǎn)。

3.3.2 等效加筋法有限元模型

采用等效加筋法進(jìn)行土工袋模擬時(shí)，袋子張力的約束作用增強(qiáng)了內(nèi)部土顆粒間的接觸作用，相當(dāng)于在內(nèi)部顆粒間增加了一個(gè)附加粘聚力。根據(jù)前述的土工袋強(qiáng)度原理，單個(gè)土工袋的凝聚力可以按式（1）計(jì)算：

圖1 土工袋網(wǎng)格模型圖

等效加筋法模擬土工袋的模型仍由1 136個(gè)單元組成，如圖2所示。圖中點(diǎn)V1～V6和H1～H6為位移監(jiān)測(cè)點(diǎn)。

3.3.3 計(jì)算結(jié)果對(duì)比分析

圖2是典型節(jié)點(diǎn)在加載過程中豎向位移與荷載的關(guān)系圖。從圖2（a）中可見，典型節(jié)點(diǎn)的豎向位移隨荷載增加均呈現(xiàn)變大的規(guī)律，且基本呈線性分布。由于底部無位移，靠近底部的典型節(jié)點(diǎn)豎向位移幾乎為0；越靠近頂部，位移越大。最大位移出現(xiàn)在靠近頂部的V1節(jié)點(diǎn)，為-2.48cm。圖2（b）中各典型節(jié)點(diǎn)位移變化與圖3（a）基本一致，最大位移同樣出現(xiàn)在V1節(jié)點(diǎn)，為-2.79cm，略大于桿單元法。

圖3是典型節(jié)點(diǎn)的水平位移與荷載的關(guān)系圖。從圖中可見，由于模型對(duì)稱加載，水平位移呈兩兩反對(duì)稱的形式。在荷載作用下，土工袋兩端水平移動(dòng)較為明顯，說明土工袋向兩端伸展，最大水平位移分布在每層土工袋中部，分別向上下兩側(cè)減小，方向均朝土工袋外側(cè)。兩種方法計(jì)算中等效加筋法的位移差值較小且略大于桿單元法，其中桿單元法典型節(jié)點(diǎn)最大水平位移出現(xiàn)在H2與H4，為-3.99cm；等效加筋法中典型節(jié)點(diǎn)最大水平位移也是H2與H4，為-4.55cm。

圖2 典型節(jié)點(diǎn)的豎向位移和加載關(guān)系

圖3 典型節(jié)點(diǎn)的水平位移和加載關(guān)系

圖4 加載后的位移矢量圖

圖4繪制了加載后袋內(nèi)土體的位移矢量圖，由于模型具有對(duì)稱性，故計(jì)算結(jié)果亦呈現(xiàn)出良好的對(duì)稱性。從圖4中可以看出，土工袋在外荷載作用下，整體發(fā)生壓縮變形，袋內(nèi)土體在土工袋受壓過程中向四周移動(dòng)，具有明顯的向兩側(cè)鼓出的趨勢(shì)，土工袋逐漸變成扁平狀，袋子周長變長，該趨勢(shì)必然受到編織袋的約束，限制其變形，從而在編織袋內(nèi)產(chǎn)生張力，這與實(shí)際中的土工袋受力變形情況完全一致。相比之下，等效加筋法的結(jié)果偏大，這是由于桿單元模擬的袋子對(duì)袋內(nèi)土起到約束作用，限制了其位移。

3.4 土工袋水平剪切試驗(yàn)?zāi)M

土工袋水平剪切試驗(yàn)是對(duì)兩個(gè)上下堆疊的土工袋組合體施加水平剪切力，計(jì)算過程采用剛性板加載法。在下方土工袋底部和左部設(shè)置剛性板并約束位移和轉(zhuǎn)動(dòng)，忽略自重對(duì)壓縮的影響，將上土工袋自重對(duì)下土工袋的影響轉(zhuǎn)化為均布力施加在下土工袋頂部。為了模擬土工袋荷載的施加過程，模型頂部施加100 kPa的豎向均布力，對(duì)上土工袋施加從右向左的水平推力，計(jì)算中共分為10級(jí)加載，最后水平荷載達(dá)到50kPa。

（1）桿單元法有限元模型。土工袋模型由1 136個(gè)單元組成，如圖5所示。圖5中點(diǎn)V11～V62和H1～H6為位移監(jiān)測(cè)點(diǎn)。

（2）等效加筋法有限元模型。土工袋模型由1 136個(gè)單元組成，如圖5所示。圖中點(diǎn)V11～V62和H1～H6為位移監(jiān)測(cè)點(diǎn)。

圖5 土工袋網(wǎng)格模型圖

圖6 典型節(jié)點(diǎn)的豎向位移和加載關(guān)系

（3）計(jì)算結(jié)果對(duì)比分析。圖6是典型節(jié)點(diǎn)在加載過程中豎向位移與荷載的關(guān)系圖。從圖6（a）中可見，土工袋右側(cè)的典型節(jié)點(diǎn)豎向位移隨著水平荷載的增加逐漸減少。這是由于在加載過程中，隨著水平荷載的增加，土工袋受層間摩擦作用，上層土工袋右側(cè)豎向位移總體呈現(xiàn)隨荷載減小的趨勢(shì)。土工袋左側(cè)典型節(jié)點(diǎn)豎向位移變化與右側(cè)不同，下層袋子會(huì)向左側(cè)剛性板移動(dòng)。上方土工袋變化隨水平推力小幅增加，下方土工袋由于設(shè)置了左側(cè)擋板，袋內(nèi)土體向上下發(fā)散。總體而言，土工袋的變形基本上與加載的大小成正比，當(dāng)水平荷載達(dá)到50kPa時(shí)，土工袋右側(cè)豎向最大位移出現(xiàn)在頂部的V11節(jié)點(diǎn)，為-0.42cm，左側(cè)豎向最大位移出現(xiàn)在上方土工袋底部的V12節(jié)點(diǎn)（由于下方土工袋袋內(nèi)土體發(fā)散，不考慮其規(guī)律性），為-0.41cm。圖6（b）中各典型節(jié)點(diǎn)位移變化基本與（a）一致，右側(cè)最大位移同樣出現(xiàn)在V11節(jié)點(diǎn)，-0.40cm，左側(cè)豎向最大位移出現(xiàn)在頂部的V12節(jié)點(diǎn)，為-0.48cm。

圖7 典型節(jié)點(diǎn)的水平位移和加載關(guān)系

圖7是典型節(jié)點(diǎn)在加載過程中水平位移與荷載的關(guān)系圖。從圖7中可見，隨著水平推力的增大，土工袋典型節(jié)點(diǎn)的水平位移先向右逐漸減小后向左逐漸增大，下層土工袋水平位移變化速度慢于上層土工袋，說明上下層土工袋之間產(chǎn)生相對(duì)滑動(dòng)。加載初期，水平荷載相對(duì)豎向荷載作用小，土工袋仍表現(xiàn)為向外伸展，隨著水平力的增加，土工袋開始沿剪力方向移動(dòng)。由于底部及左側(cè)設(shè)置剛板，故加載初期土工袋底部水平位移較大，且變化速度最慢。從圖7（a）和圖7（b）計(jì)算的結(jié)果對(duì)比來看，桿單元法和等效加筋法計(jì)算結(jié)果也保持高度一致。其中，桿單元法典型節(jié)點(diǎn)最大水平位移出現(xiàn)在H2，為-2.85cm；等效加筋法中典型節(jié)點(diǎn)最大水平位移也是H2，為-3.08cm。

圖8 加載后的位移矢量圖

圖8繪制了加載后袋內(nèi)土體的位移矢量圖，土工袋受豎向荷載和水平力作用，上層土工袋受力端首先受到擠壓，局部產(chǎn)生滑動(dòng)。隨著水平剪力的增大，袋體其余部位沿受力方向逐漸發(fā)生滑動(dòng)。加載過程中，上層土工袋由于受到豎向荷載、水平剪力和層間摩擦力作用，袋內(nèi)土體發(fā)生側(cè)向壓縮和平動(dòng)，袋內(nèi)土體的位移從上到下逐漸減小。由于Abaqus中在模型頂部施加的剛性板約束袋子的滾動(dòng)，導(dǎo)致矢量圖中上層土工袋右側(cè)并沒有出現(xiàn)向上的鼓起，且從左到右逐漸減小。

4 結(jié)論

本文通過有限元數(shù)值模擬對(duì)土工袋無側(cè)限壓縮試驗(yàn)和剪切試驗(yàn)進(jìn)行分析，結(jié)果表明：

（1）用桿單元法和等效加筋法模擬土工袋，計(jì)算結(jié)果合理并基本一致。使用桿單元法時(shí)，僅受張拉的桿單元對(duì)袋內(nèi)土體起到約束作用，限制其位移，因此計(jì)算結(jié)果略小于等效加筋法。

（2）桿單元法中計(jì)算出的土工袋袋子張力充分闡明了土工袋的作用機(jī)理，通過與等效加筋法進(jìn)行對(duì)比計(jì)算分析，也表明桿單元法的計(jì)算結(jié)果要更為理想。