連續(xù)介質(zhì)滲流區(qū)域非線性問題數(shù)值求解分析

潘廣良

(方正縣河湖運(yùn)行保障中心，黑龍江 方正 150800)

0 前 言

土石壩的非線性滲流求解在于滲流自由面和溢出面的確定。文章基于移動(dòng)網(wǎng)格法、單元滲透矩陣調(diào)整法和復(fù)合單元高斯點(diǎn)法，運(yùn)用Microsoft Developer Studio可視化集成開發(fā)二維無壓滲流分析有限元程序，結(jié)合具體算例求解連續(xù)介質(zhì)滲流區(qū)域的非線性問題，成果對(duì)滲流區(qū)域非線性問題研究具有重要意義。

1 緒 論

中國目前的巖土體水力學(xué)研究，為實(shí)際應(yīng)用和方便處理考慮，基本選用線性數(shù)學(xué)模型。線性數(shù)學(xué)模型通常將滲流運(yùn)動(dòng)和物體狀態(tài)簡化處理，因此助推巖土體水力學(xué)近些年快速發(fā)展[1]。但實(shí)際滲流由多種因素復(fù)合影響而成，絕大多數(shù)的滲流分析均呈非線性，故材料的非線性問題在巖土體水力學(xué)研究中為基本問題，具有十分重要的研究價(jià)值。

2 滲流基本原理和數(shù)學(xué)模型

連續(xù)介質(zhì)滲流區(qū)域非線性問題以二維非穩(wěn)定流的飽和滲流模型為例，4種表達(dá)方程式[2]如下：

1)Boussinesq方程：在自由面進(jìn)行流量補(bǔ)給，并將緩變形式簡化為水平向滲流，在區(qū)域地下水活動(dòng)研究中應(yīng)用廣泛。

(1)

2)Laplace方程：自由面既可進(jìn)行流量補(bǔ)給，同時(shí)為下降流速的限值條件，表達(dá)式為：

(2)

3)擴(kuò)散方程：自由面不進(jìn)行流量補(bǔ)給，該方程僅針對(duì)自由面微弱變化的均質(zhì)土壩，表達(dá)式為：

(3)

4)固結(jié)方程：因考慮到土體具有壓縮性，常應(yīng)用于粘土筑壩的固結(jié)，表達(dá)式為：

(4)

3 算例分析

3.1 移動(dòng)網(wǎng)格法

連續(xù)介質(zhì)滲流的二維有限元模型分為2種不同單元進(jìn)行劃分，三角形單元和8節(jié)點(diǎn)四邊形等參單元[3]。算例為一均質(zhì)土壩，壩面寬6m，上、下游水位分別為4m和3m，計(jì)算模型示意圖，見圖1。

圖1 計(jì)算模型示意圖

邊界條件：滲流邊界BC無補(bǔ)給，AD為封閉邊界。

模型進(jìn)行三角形單元?jiǎng)澐謺r(shí)，共劃分18個(gè)單元，16個(gè)結(jié)點(diǎn)；八節(jié)點(diǎn)四邊形單元?jiǎng)澐謺r(shí)，共劃分為9個(gè)單元，40個(gè)結(jié)點(diǎn)，三角形單元網(wǎng)格剖分，見圖2；八結(jié)點(diǎn)四邊形等參元網(wǎng)格剖分，見圖3。

圖2 三角形單元網(wǎng)格剖分 圖3 八結(jié)點(diǎn)四邊形等參元網(wǎng)格剖分

通過程序模擬計(jì)算，得到2不同單元的滲流自由面結(jié)點(diǎn)水頭計(jì)算結(jié)果值，區(qū)域滲流自由面結(jié)點(diǎn)水頭對(duì)比值，見表1。

表1 區(qū)域滲流自由面結(jié)點(diǎn)水頭對(duì)比值

3.2 單元滲透矩陣調(diào)整法[4]

該法的無壓滲流分析單元?jiǎng)澐衷硗揪W(wǎng)格法。算例為一均質(zhì)土壩，高8m，寬4m，上、下游水位6.0m和1m，該均質(zhì)土壩模型，見圖4。

圖4 均質(zhì)土壩計(jì)算模型

邊界條件：BC邊界無流量補(bǔ)給，流量為零，CD邊界流量可能溢出，AE為封閉邊界。本算例先選用四邊形單元剖分模型，共剖分24個(gè)單元，93個(gè)結(jié)點(diǎn)。四邊形單元網(wǎng)格剖分，見圖5。

圖5 四邊形單元網(wǎng)格剖分

為保證迭代的穩(wěn)定性和結(jié)果的精確性，進(jìn)行5次迭代后計(jì)算得出溢出點(diǎn)水頭值5.68m，滲流自由面各節(jié)點(diǎn)水頭統(tǒng)計(jì)值，見表2。

表2 滲流自由面結(jié)點(diǎn)水頭統(tǒng)計(jì)值

由結(jié)果可得，在相同節(jié)點(diǎn)不同數(shù)解值有較大差距數(shù)，故在四邊形單元基礎(chǔ)上再次進(jìn)行三角形網(wǎng)格加密，共劃分單元48個(gè)，結(jié)點(diǎn)35個(gè)，三角形單元網(wǎng)格剖分，見圖6。

圖6 三角形單元網(wǎng)格剖分

本次求解進(jìn)行10次迭代后得到滲流面溢出點(diǎn)水頭值為3.23m，模型的滲流自由面和節(jié)點(diǎn)水頭結(jié)果見圖7和圖8，滲流自由面各節(jié)點(diǎn)水頭統(tǒng)計(jì)值，見表3。

圖7 計(jì)算模型滲流自由面

圖8 節(jié)點(diǎn)水頭等值線圖

表3 滲流自由面結(jié)點(diǎn)水頭統(tǒng)計(jì)值

由上述結(jié)果可得，網(wǎng)格加密的滲流求解更加貼合實(shí)際，為更好證明自編程序的合理性，重新運(yùn)用四邊形單元加密劃分模型，本次單元?jiǎng)澐謫卧?50個(gè)，節(jié)點(diǎn)501個(gè)，模型單元?jiǎng)澐纸Y(jié)果見圖9。迭代9次后得到自由面水頭值為3.78m，計(jì)算模型滲流自由面，見圖10；節(jié)點(diǎn)水頭等值線圖，見圖11；滲流自由面結(jié)點(diǎn)水頭值，見表4。

表4 滲流自由面結(jié)點(diǎn)水頭值

圖10 計(jì)算模型滲流自由面

圖11 節(jié)點(diǎn)水頭等值線圖

圖9 模型四邊形單元網(wǎng)格剖分

由表4可知，單元進(jìn)行加密后的2解相差減少，對(duì)比上述方法的計(jì)算結(jié)果，本法更加吻合實(shí)際。

3.3 復(fù)核單元高斯點(diǎn)法[5]

對(duì)上算例進(jìn)行復(fù)核單元高斯點(diǎn)法求解，得到計(jì)算結(jié)果如下：滲流自由面節(jié)點(diǎn)水頭統(tǒng)計(jì)值見表5；對(duì)比不同方法的代表水頭結(jié)果見表6，經(jīng)迭代9次后得到自由面水頭值為3.35m。不同方法得到的滲流自由面上典型點(diǎn)水頭值及誤差見，計(jì)算模型滲流自由面和節(jié)點(diǎn)水頭計(jì)算值見圖12和圖13。

表5 滲流自由面結(jié)點(diǎn)水頭值

表6 滲流自由面代表節(jié)點(diǎn)水頭計(jì)算統(tǒng)計(jì)表

圖12 計(jì)算模型滲流自由面

圖13 節(jié)點(diǎn)水頭等值線圖

通過上述結(jié)果表明，復(fù)合單元高斯點(diǎn)求得的滲流面節(jié)點(diǎn)位置更貼合實(shí)際，誤差僅為5.2%，遠(yuǎn)<單元滲透矩陣調(diào)整法，二者均能體現(xiàn)滲流區(qū)域震蕩，在實(shí)際應(yīng)用中均有自身的利用優(yōu)勢(shì)，但在同樣的前提條件下，復(fù)合單元高斯點(diǎn)具有更適宜的合理性應(yīng)用。

4 結(jié) 論

文章在連續(xù)介質(zhì)滲流區(qū)域非線性問題求解的基本原理基礎(chǔ)上，建立不同算例進(jìn)行建模，采用3種不同方法有限元計(jì)算對(duì)應(yīng)算例的滲流區(qū)非線性問題，得出不同模型的滲流自由面和代表節(jié)點(diǎn)水頭值，認(rèn)為在相同條件下，復(fù)合單元高斯點(diǎn)法更能貼合實(shí)際應(yīng)用，同時(shí)驗(yàn)證了自編程序的合理性和準(zhǔn)確性。