流體的質(zhì)量守恒方程

流體的質(zhì)量守恒方程

這是流量修正的恒定方程。若單元體積內(nèi)不發(fā)生源匯,即W=0,則:divJ=0(10)這是整個流體體系的質(zhì)量守恒方程,即連續(xù)性方程。各向同性介質(zhì)的滲透系數(shù)是一個表示多孔介質(zhì)輸送流體能力的標量。此時不可壓縮流體的三維運動的Darcy定律為:{qx=ΚΙx=-Κ?Η?xqy=ΚΙy=-Κ?Η?yqz=ΚΙz=-Κ?Η?z(11)對于各向異性介質(zhì),滲透系數(shù)是一個張量,三維情況下:[Κ]=[ΚxxΚxyΚxzΚyxΚyyΚyzΚzxΚzyΚzz](12)Darcy定律的表達為:q=[Κ]Ι(13)根據(jù)連續(xù)性方程式(10)及Darey定律,可以得出:divJ=div(gradH)=0(14)即:ΔΗ=?2Η=?2Η?x2+?2Η?y2+?2Η?z2=0(15)式(15)為Laplace方程,亦稱不可壓縮流體運動的微分方程,它描述均質(zhì)各向同性的不變形介質(zhì)中不可壓縮流體滲流場中水頭的分布H=H(x,y,z,t)。它對于穩(wěn)定流動和不可壓縮流體的非穩(wěn)定流動都適用。滿足Laplace方程的非穩(wěn)定流動是由于隨時間變化和邊界條件引起的。任何偏微分方程都有無窮多個解。要從眾多可能解中得到所研究問題的解,除微分方程外,還要提供求解條件,包括初始條件和邊界條件。3.1.1初始條件是指滲流過程開始的某一特定時刻研究區(qū)域上所有點必須滿足的初始條件。例如必須給出初始時刻區(qū)域內(nèi)的水頭分布。設(shè)水頭H=H(x,y,z,t),初始條件的表達式一般為H(x,y,z,t)=H0(x,y,z),(x,y,z∈Ω)(16)式(16)中,H0—已知函數(shù);Ω—滲流區(qū)域。3.1.2邊界條件①第一類邊界條件為已知水頭的邊界條件,在邊界的所有點上水頭是給定的,例如:H(x,y,z)=Hb(x,y,z),(x,y,z)∈s(17)H(x,y,z,t)=Hb(x,y,z,t),(x,y,z)∈s(18)式中,s—三維區(qū)域的邊界曲面。②第二類邊界條件為已知通量的邊界條件,即垂直于邊界面的流量給定的,一般有:qn=→q?→n=-qb(x,y,z)?(x,y,z)∈s(19)或qn=→q?→n=-qb(x,y,z,t)?(x,y,z)∈s(20)式中,qb—邊界面上沿法線方向的單位面積入量;n—邊界外法線的單位矢量。③第三類邊界條件為水頭和水頭的法向?qū)?shù)的組合在邊界上為已知,即:?Η?n+λ(x,y,z)Η=f(x,y,z)?(x,y,z)∈s(21)式(21)中,λ、f—均為已知函數(shù),λ稱為變換系數(shù)。3.2二維有限元方法首先,根據(jù)滲流運動微分方程與邊界條件約束建立二維數(shù)學(xué)模型:{??x(Τ?Η?x)+??y(Τ?Η?y)=0,在Ω內(nèi)Η=Ηb(x,y),在Γ1上Τ?Η?xcos(n,x)+Τ?Η?ycos(n,y)=q,在Γ2上(22)式(22)中,Ω—計算區(qū)域;Γ1—第一類邊界;Γ2—第二類邊界;Hb—第一類邊界上的已知水頭;n—第二類邊界的外法線方向;q—第二類邊界上法向單寬流量,流入為正,流出為負。采用伽遼金加權(quán)剩余法,設(shè)ˉΗ(x,y)為區(qū)域Ω上地下水頭的試探解?！ウ?x,y)=ne∑L=1ΝLΗL(23)式(23)中,NL—在單元e上分片定義的基函數(shù);HL—單元節(jié)點上的水頭值;ne—單元結(jié)點數(shù)目。令基函數(shù)為權(quán)函數(shù),代入試函數(shù)ˉΗ。令權(quán)剩余為零,則權(quán)剩余方程為:∫ΩΝL(L(ˉΗ-f)dΩ=Μ∑e=1∫ΩeΝL(L(ˉΗ)-f)dΩ=0(24)式(24)把權(quán)剩余在整個區(qū)域上的積分,化為在各個單元上的積分,然后求和,這樣得到一個線性方程組,求解此線性方程組,便可得到各節(jié)點的水頭值。3.2.1基函數(shù)三角形單元節(jié)點i,j,m按逆時針編號,則:ˉΗ(x,y)=12A[(ai+bix+ciy)Ηi+(aj+bjx+cjy)Ηj+(am+bmx+cmy)Ηm](ΝiΗi+ΝjΗj+ΝmΗm)(25)于是,三角形單元基函數(shù)表達式為:NL=(aL+bLx+cLy)/2A,(L=i,j,m)(26)3.2.2滲透矩陣由上面給出的基函數(shù),故試函數(shù)的表達式為:ˉΗ=ne∑i=1ΝLΗL=ΝiΗi+ΝjΗj+ΝmΗm(27)權(quán)剩余方程為:?ΩΝL[??x(Τ?ˉΗ?x)+??y(Τ?ˉΗ?y)]dxdy=0(L=1,2,?,Ν)(28)對式(28)進行分部積分,并應(yīng)用格林公式,得:∫Γ2ΝLqds-?DL(Τ?ˉΗ?x?ΝL?x+Τ?ˉΗ?y?ΝL?y)dxdy=0(29)對于每一個單元,有:QLe=?Ω(Τe?Ηˉ?x?ΝL?x+Τe?Ηˉ?y?ΝL?y)dxdy-∫Γ2eΝLds(30)根據(jù)式(25)、式(26)易于導(dǎo)出:{?Ηˉ?x=12A(biΗi+bjΗj+bmΗm)?Ηˉ?y=12A(ciΗi+ciΗj+cmΗm)?ΝL?x=bL2A,?ΝL?y=cL2A,(L=i,j,m)(31)將式(31)代入式(30),整理得如。式(32)的矩陣方程:[QieQjeQme]=Τe4AB-[FiFjFm],B=[bibi+cicibibi+cicjbibm+cicmbjbi+cjcibjbj+cjcjbjbm+cjcmbmbi+cmcibmbj+cmcjbmbm+cmcm]{ΗiΗjΗm}(32)式(32)可記為:[Qe]=[G]e{Η}e-{F}e(33)式(33)中,[G]e—單元滲透矩陣。由于區(qū)域Ω是由M個單元組成,現(xiàn)將各單元滲透矩陣集合起來,便得到總滲透矩陣。式(33)又可寫作:[G]{Η}={F}(34)3.2.3邊界條件在第一類邊界上的節(jié)點,水頭H是給定的,無須計算。在第二類邊界上的節(jié)點將與節(jié)點L有關(guān)的單元e的列向量{F}作如下處理:FLe=∫Γe2qNLds(35)對于二維滲流情況:q=Τ?Η?xcos(n,x)+Τ?Η?ycos(n,y)(36)4基于2D-FLOW軟件滲流場數(shù)值模擬計算4.1中線法尾礦壩滲流場數(shù)值計算4.1.1壩體計算的幾何模型與工況正確確定計算幾何模型,是保證計算結(jié)果真實有效的根本。根據(jù)中線式尾礦壩的設(shè)計資料,以及筑壩極限高度,來構(gòu)造計算幾何模型,如圖2所示,計算最終標高為+2687m,即壩高185m,分兩個材料區(qū)考慮。幾何參數(shù):ab=877m,bc=382.8m,ac坡度12%,ad坡度1:3,ad=897mde=ef=10m,df=20m,g點為水位點,fg坡度1%,oe=185m(壩高度)。計算范圍按照400m考慮。底部為不透水層。同時按照下列情況進行模擬計算:(a)堆積子壩按照20m寬度考慮;(b)考慮到尾礦庫內(nèi)水位的不斷變化,對干灘面的距離分別按70m(規(guī)范中規(guī)定最干灘面長度)、100m和300m三種情況計算;(c)初期壩按照導(dǎo)水情況考慮;對壩體按照設(shè)置排滲系統(tǒng)和不設(shè)排滲系統(tǒng)兩種情況考慮。4.1.2壩體材料的計算參數(shù)根據(jù)筑壩情況,按照兩種材料考慮,即1號材料區(qū)按照尾粉砂考慮;2號材料區(qū)按照尾粉土考慮。庫底為基巖,計算時按照不透水層考慮。不同資料中尾礦壩體各土層的計算參數(shù)見表1。初期壩為堆石透水壩,重度γ=21.0kN/m3、內(nèi)摩擦角φ=38.0°、黏結(jié)力為C=0.0、滲透系數(shù)為K=20×10-2cm/s。4.1.3中線式尾礦壩地下水滲流場的數(shù)值計算結(jié)果按照2D-FLOW程序使用說明,將上述計算模型和材料參數(shù)輸入計算機,進行數(shù)值計算。這次分2個大類總共6個小類進行計算,即:(A)第1大類—沒有排滲設(shè)施(a)先按照干灘面長度為300m計算,即dg=300m。計算結(jié)果見圖3和圖4。(b)然后,按照干灘面長度為100m計算,即dg=100m。計算結(jié)果見圖5、圖6和圖7。(c)最后按照干灘面70m計算,即dg=70m。計算結(jié)果見圖8和圖9。(B)第2大類——按照有排滲設(shè)施(井)考慮,排滲設(shè)計(井)按照直徑1.0m考慮。干灘面長度分別為300m、100m和70m三種情況考慮。其他條件不變。其結(jié)果分別見圖10、圖11、圖12、圖13和圖14。5結(jié)論從這些計算結(jié)果中可以得出:(1)不設(shè)排滲設(shè)施,則洪水情況下,壩體的浸潤線會與壩坡相交(圖8和圖9),表明壩坡上有水滲出,對壩體的穩(wěn)定性產(chǎn)生很不利的影響。(2)采用垂直排滲和水平排滲措施后,則壩體浸潤線大大改觀,對壩體的穩(wěn)定極為有利。為確保尾礦庫的生產(chǎn)安全,了解尾礦庫滲流場的變化規(guī)律非常重要。通過數(shù)值模擬可以預(yù)先知道滲流場的變化規(guī)律,為防止尾礦壩產(chǎn)生滲流變形和壩面溢出提供依據(jù),也為采取措施降低壩體浸潤線指明方向。壩坡浸潤線是尾礦壩的生命線,它是直接影響壩體安全的一個非常重要的因素。而影響地下水位的最重要因素是尾礦庫基礎(chǔ)相對于尾礦沉積層的滲透性、顆粒分級程度和尾礦沉積層滲透性的側(cè)向變化,以及庫內(nèi)水位相對于壩頂?shù)奈恢?。當庫?nèi)水浸入沉積灘時,靠近壩面則產(chǎn)生很高的地下水位,危及穩(wěn)定性。所以選擇合適的尾礦排放方法和控制溢水的高度,可以把庫內(nèi)水位從壩頂向后推移,加大干灘面的長度;提高泄水速度,從而降低庫內(nèi)水位高程,可以增大庫內(nèi)水與壩頂之間的相對距離;在庫內(nèi)水浸入沉積灘的關(guān)鍵地段增大尾礦排放量,通過形成較高和較寬的沉積灘也可向后推移池水位置。通過地下水滲流場的計算,可以模擬尾礦壩地下水位的變化情況,這樣有助于尾礦庫的安全生產(chǎn)與管理。1尾礦壩浸潤線的特點尾礦庫工程的最突出特點是,地下水滲流狀態(tài)成為控制壩坡穩(wěn)定性和污染物遷移的決定性因素因此,尾礦庫地下水滲流狀態(tài)的可靠分析與評價是尾礦庫工程研究的關(guān)鍵性問題之一。地下水滲流分析的主要目的有兩個,一是為了估計孔隙壓力,為壩體穩(wěn)定性分析提供地下水的數(shù)據(jù)。對這種情況一般是假設(shè)尾礦壩內(nèi)滲流是在重力流動及穩(wěn)態(tài)條件下發(fā)生的地下水的流動;二是為了確定尾礦庫滲漏損失,以預(yù)測污染潛勢,這種情況需要進行非穩(wěn)態(tài)、瞬態(tài)或非飽和滲流評價,非飽和滲流是指在毛細作用而不是在重力梯度下發(fā)生的地下水的流動。在穩(wěn)態(tài)重力滲流假定的條件下,確定壩內(nèi)地下水位(或頂面流線),得到穩(wěn)定性分析所必需的孔隙水壓力分布。一般確定地下水位的最普通方法是采用流網(wǎng)分析。當然根據(jù)Darcy定律和LaPlace方程推導(dǎo)的基本滲流原理通常是應(yīng)用于普通水壩,如果考慮到邊界條件的差異,也可以應(yīng)用于尾礦壩。2滲流對尾礦壩穩(wěn)定性影響分析當尾礦庫中的滲流場形成后,如果浸潤線在壩前坡面逸出,坡內(nèi)尾礦除了受自重力作用外,還受滲透力的作用,因而會降低尾礦壩的穩(wěn)定性。當滲流從坡面逸出時,如果水流的方向與水平面夾角為則沿水流方向的滲透力式中為水的重度為水力坡度在壩坡面上取土體V中的土骨架為隔離體,其有效重量為γ′V(γ′為土的浮重度)。分析壩坡土體土骨架的穩(wěn)定性,則作用在土骨架上的滲透力J為:沿壩坡面的全部下滑動力T為重力和滲透力式(5)中:i—滲透坡降;γ′—土的浮重度;γw—水的重度;φ—土的內(nèi)摩擦角。若滲流在溢出段順坡面流動,即θ=α,這時,流經(jīng)途徑ds的水頭損失為dh,故有:代入式(5)中,得:式(7)中,γsat—土的飽和重度;其它符號意義同前。降低γ′/γsat,通常γ′/γsat約為0.5,則安全系數(shù)降低一半。因此若要保持同樣安全度,有滲流溢出時的坡角