管涌過程中的流固耦合模型

管涌過程中的流固耦合模型

0管涌數(shù)學(xué)模型的研究滲透破壞是主要原因之一，比如大壩和地基等項目的穩(wěn)定性和破壞。滲透破壞發(fā)生后,土體內(nèi)部顆粒大量流失,滲透流量急劇增加,從而形成集中滲透通道,使土體內(nèi)部滲流場及應(yīng)力場發(fā)生變化,可能使堤防在很短的時間內(nèi)發(fā)生破壞,并可能造成的邊坡失穩(wěn)或洪水災(zāi)害等會對人民的生命財產(chǎn)安全造成極大的損失。對于滲透破壞機理以及防治方法的研究,一直是學(xué)者研究的重要課題。管涌是滲透破壞的主要形式,在工程中,也常將其它的滲透破壞形式統(tǒng)稱為管涌。管涌是在滲流作用下,土體內(nèi)部細(xì)小顆粒沿大顆粒通道移動并被帶出土體的過程。由于土體構(gòu)成具有很強的隨機性與復(fù)雜性,管涌發(fā)生后,由于顆粒的流失,土體的滲透性一直在發(fā)生變化,滲透性的變化會影響滲流場的分布,滲流場分布的變化又會對顆粒的流失產(chǎn)生影響,使得顆粒的流失與滲流場的變化相互影響,是一個耦合過程,從而很難利用數(shù)學(xué)語言對管涌的發(fā)生發(fā)展過程進行描述。因此,在以前的研究中,一般采用室內(nèi)試驗的方法對管涌的過程進行模擬研究。眾多學(xué)者在管涌的室內(nèi)試驗與工程應(yīng)用方面取得了豐碩的成果。近年來,隨著研究手段的發(fā)展以及研究思想的創(chuàng)新,新的研究方法及研究手段被利用到管涌機理研究中來,如有限元法、離散元方法、神經(jīng)網(wǎng)絡(luò)方法、隨機方法等,極大的推動了管涌機理的研究與應(yīng)用。但如何建立并求解管涌發(fā)展的數(shù)學(xué)模型依然是困擾研究人員的主要問題。本文從運動顆粒的守恒定律以及運動流體的守恒定律出發(fā),建立起顆粒流失與滲透性變化的相互聯(lián)系,建立了管涌發(fā)展過程的數(shù)學(xué)模型,并以一維情況為例,利用有限差分法對模型進行了求解,得到管涌發(fā)展過程中顆粒流失、滲流性變化、以及流速變化等的規(guī)律。1移動顆粒連續(xù)方程的建立鑒于管涌的定義,管涌的發(fā)生與發(fā)展要滿足以下基本條件:①土體為管涌型土體,即土體中的顆粒可以分為骨架顆粒與可以移動的顆粒,且骨架顆粒的孔隙可以提供可動顆粒被帶出的通道;②可動顆粒的靜力平衡達(dá)到被破壞的條件。在滲流場中,顆粒受到流體流動的拖曳力的作用,當(dāng)拖曳力的大小達(dá)到破壞顆粒靜力平衡的臨界值時,顆粒開始發(fā)生移動。顆粒受到的拖曳力與顆粒的大小以及滲流的孔隙流速有關(guān),對于粒徑一定的顆粒,使其發(fā)生移動的孔隙流速稱為顆粒的起動流速。在本模型的推導(dǎo)過程中,假定符合以上基本條件,并假設(shè):①可動顆粒流失后土體體積不發(fā)生變化,即流失顆粒的體積變?yōu)榭紫扼w積;②只要孔隙流速大于可動顆粒的起動流速,則可動顆?？梢酝耆魇?即顆粒的淤堵作用。并令可動顆粒的密度為ρf,其物理意義為單位體積土體內(nèi)可動顆粒的質(zhì)量;骨架顆粒的密度為ρs,物理意義為單位體積土體內(nèi)骨架顆粒的質(zhì)量。對于顆粒運動的質(zhì)量守恒定律,Sterpi與Cividini等曾分別利用室內(nèi)試驗法與有限單元法進行過分析。在圖1所示的微元體中,利用與文獻(xiàn)類似的方法,根據(jù)質(zhì)量守恒定律,建立移動顆粒的質(zhì)量守恒方程,圖1中,q為移動顆粒通量。根據(jù)質(zhì)量守恒,可得到在Δt時間內(nèi)單元體內(nèi)移動顆粒的質(zhì)量守恒定律:式中ρm為移動顆粒的密度,其物理意義為單位時間內(nèi)流經(jīng)單位土體的顆粒質(zhì)量與流體體積的比值;I為移動顆粒的源匯項,其物理意義為管涌過程中的單位時間內(nèi),單位體積土體內(nèi)可動顆粒轉(zhuǎn)化為移動顆粒的質(zhì)量,此時源匯項為正;如移動顆粒轉(zhuǎn)化為靜態(tài)的可動顆粒,則為顆粒的沉積過程,此時源匯項為負(fù)。根據(jù)q的定義得式中,vx,vy,vx為滲透流速矢量在3個方向上的分量。將式(2)代入式(1)整理得式(3)即為管涌過程中移動顆粒的連續(xù)性方程,方程中共有ρm,兩個未知數(shù),需聯(lián)立滲流的連續(xù)性方程才可以求解。假定流體不可壓縮,得到流體的連續(xù)性方程為假定管涌過程中滲流仍符合達(dá)西定律,則式(4)可以表述為式中,H為水頭,k為滲透系數(shù)。聯(lián)立式(3)~(5),結(jié)合邊界條件與初始條件,即為管涌發(fā)展的數(shù)學(xué)模型。在模型求解前,需要考慮兩式的耦合:①當(dāng)孔隙流速大于起動流速時,顆粒起動并被帶走,孔隙流速小于起動流速時移動的顆粒發(fā)生沉積,因此顆粒的源匯項與流體的流速是相關(guān)的;②無論顆粒的流失或沉積,都會引起土體中孔隙率的變化,而土體的滲透系數(shù)與孔隙率有關(guān),孔隙率越大滲透性越強,反之越弱。下面利用各參數(shù)的相互聯(lián)系,建立耦合過程的數(shù)學(xué)關(guān)系。對于移動顆粒連續(xù)性方程中的源匯項I,其實質(zhì)是可動顆粒的流失量,當(dāng)孔隙流速大于起動流速時,土體內(nèi)的可動顆粒起動,轉(zhuǎn)化為移動顆粒,源匯項為正,即顆粒流失過程;反之,孔隙流速減小,流體的拖曳力不足以維持移動顆粒的運動,則移動顆粒會發(fā)生沉積,此時源匯項為負(fù),為顆粒沉積過程。當(dāng)孔隙流速大于起動流速時,I隨孔隙流速的增加而增加,當(dāng)孔隙流速與起動流速相等時顆粒處于臨界狀態(tài);同時與單元體內(nèi)剩余可動顆粒的量有關(guān),剩余量越多,I越大,剩余量為0則I也為零?？紤]這兩種因素,構(gòu)建流失過程的線性模型為式中α為系數(shù),且α>0;rv為孔隙流速的大小;v*為顆粒的起動流速。根據(jù)I的物理概念,I為單位時間內(nèi),單位體積土體內(nèi)可動顆粒轉(zhuǎn)化為移動顆粒的質(zhì)量,因此I的量綱為M/(T·L3)。式(6)中ρf的量綱為M/L3,(vr-v*)的量綱為L/T,因此α也是有量綱的參數(shù),其量綱為1/L。當(dāng)實際流速小于等于臨界流速時,顆粒發(fā)生沉積,I為負(fù),且絕對值隨流速的減少而增加;同時與單元內(nèi)移動顆粒的量有關(guān),移動顆粒越多,I絕對值越大,濃度為0,則I也為0,考慮這兩個影響因素,構(gòu)建顆粒沉積的線性模型為式中,β為系數(shù),且β>0。綜合式(6)、(7)得流失量為剩余可動顆粒濃度與總的流失量有關(guān),即式中,ρf0(x,y,z)為初始可動顆粒濃度。孔隙流速vr(x,y,z,t)與達(dá)西流速和孔隙率有關(guān),即由于流失的顆粒轉(zhuǎn)變?yōu)榭紫扼w積,那么孔隙率n(x,y,z,t)也與流失量有關(guān),即式中,n0為初始孔隙率,γfs為可動顆粒的顆粒密度。在耦合模型中,滲透系數(shù)也是一個隨管涌發(fā)展過程變化的量,根據(jù)Kozeny–Carman方程得式中Ss為單位體積土體的顆粒表面積總量;Cs為形狀系數(shù);T為扭曲系數(shù)。如果忽略顆粒流失或沉積對顆粒表面積總量的影響,則CsSs2T2為常數(shù),即認(rèn)為滲透系數(shù)與n3/(1-n)2呈正比例關(guān)系。2初始模型材料及時間耦合模型在確定上述參數(shù)后,結(jié)合邊界條件與初始條件,即可對模型進行求解。由于模型的復(fù)雜性,難以解得解析解,因此利用數(shù)值方法對模型進行求解。在模型求解前,利用分時步法對模型進行解耦:即將時間過程按一個很小的時間間隔進行離散,在一個時步下,根據(jù)滲流場的邊界條件及初始條件,先對滲流場(式(5))進行求解,將求得的滲流場參數(shù)代入式(8)中得到源匯項,然后代入式(3)中對顆粒的運動過程進行求解,求解的結(jié)果代入式(9)~(12)中,對滲流場參數(shù)進行更新,然后將更新后的參數(shù)再代回式(5)中進行下一時步的求解,如此循環(huán),得到隨著時間的發(fā)展,管涌過程中可動顆粒的侵蝕過程。以一維管涌情況為例,利用有限差分法介紹模型的求解過程,模型高度1m,底部水頭H1,上部水頭H2,H1>H2,在穩(wěn)定水頭差作用下發(fā)生由下向上的一維滲流。模型初始狀態(tài)為均質(zhì),初始孔隙率0.3,初始滲透系數(shù)1×10-5m/s。土體初始干密度1800kg/m3,其中可動顆粒質(zhì)量占總固相的20%,其余為骨架顆粒,即初始可動顆粒密度為ρf0=360kg/m3。在初始條件下,流速、移動顆粒密度都為0(圖2)。在一維情況下,耦合模型退化為根據(jù)有限差分原理,將模型離散為間距為Δz的差分格點,同樣將時間離散為Δt的間隔,以差商代替微分,得顆粒運動方程的隱式差分格式為流體流動的隱式差分格式為將邊界條件和初始條件都離散到各差分格點上,對每個格點列式(15)的方程,如果有N個格點,就可以得到N個方程組成的線性方程組,利用簡單迭代法對方程組進行求解得到各格點的水頭值,然后利用達(dá)西定律計算各格點上的流速,再根據(jù)式(10)計算各格點的孔隙流速,代入式(8)中計算流失量,同時對滲透系數(shù)等參數(shù)進行更新,將得到的參數(shù)代入式(14)中,結(jié)合邊界條件與初如條件,同樣得到N維方程組,求解后即得的顆粒運動狀態(tài),從而完成一個時步的計算;計算完成后將相關(guān)參數(shù)進行更新,再代入式(15),從而進行下一時步的計算,如此循環(huán)從而實現(xiàn)對管涌過程的模擬。各參數(shù)在計算過程中第j時步的流失量值為將式(16)代入式(14)中解得移動顆粒密度,并按下列兩式對孔隙率及滲透系數(shù)進行更新:在本算例中,取Δz=0.01m,則共有101個差分格點;取Δt=0.5s,取總的模擬時間長度為65s,即共有130個時間差分時步。取α=β=1000/m,并假定顆粒起動流速為v*=1.25×10-5m/s。取水頭差ΔH=H1-H2=0.8m,則總體水力梯度i=ΔH/L=0.8。將以上各參數(shù)代入模型中,對于線性方程組的求解,利用SOR超松弛迭代法進行求解。求解的結(jié)果見圖3~5。圖3為模型中移動顆粒通量q隨時間的變化規(guī)律,通量方向與滲透流速方向相同,即從模型底部指向頂部,在同一時刻,移動顆粒通量隨z的增加而呈線性增長;隨著管涌的發(fā)展,在同一層位,q隨時間的增加而呈現(xiàn)先增加后減小的趨勢。圖4,5為在z=0.5m層位可動顆粒流失量、孔隙流速、孔隙率以及剩余可動顆粒密度隨時間的變化曲線。隨著時間的發(fā)展,顆粒流失量也呈現(xiàn)出先增加后減小的趨勢,伴隨著顆粒的流失,孔隙率是單調(diào)減小、孔隙流速是單調(diào)增加的,當(dāng)管涌發(fā)展到65s時,剩余顆粒密度接近0,即說明土體內(nèi)的可動顆粒大部分已被帶出土體,這也是造成流失量先增加而后減小的原因。當(dāng)顆粒接近流失時,各曲線都呈現(xiàn)出向恒定值的漸近性,即在流失后期,各參數(shù)的變化速度都在減小,而越來越接近穩(wěn)定狀態(tài)。3水力梯度對孔隙流失量的影響為了分析不同水力梯度條件下的顆粒流失規(guī)律的異同,又分別采用i=0.4、i=0.6、i=1.0對模型進行求解。求解的結(jié)果見圖6~9。圖6為在z=0.5層位不同水力梯度條件下流失量的變化規(guī)律,從圖中可以看出,水力梯度對流失量的峰值具有明顯的影響,i=1.0條件下的流失量峰值接近i=0.6條件下的兩倍;水力梯度同樣影響流失量到達(dá)的時間,i=1.0時15s左右即達(dá)到流失量峰值,而此時i=0.4時的流失量還不到i=1.0的1%。流失量曲線具有“拖長尾”的現(xiàn)象,即當(dāng)管涌后期,土體內(nèi)可動顆粒已大部分被帶出土體,剩余可動顆粒已很少,導(dǎo)致流失量下降。圖7,8為不同水力梯度條件下孔隙流速和孔隙率隨時間的變化曲線。在4次模擬過程中,孔隙流速和孔隙率都呈增長的趨勢,根據(jù)i=1.0條件下的孔隙流速變化看出,孔隙流速隨著時間的發(fā)展先加速增加,而后減速增加,并有保持穩(wěn)定的趨勢。雖然孔隙流速與孔隙率呈反比例關(guān)系,孔隙率的增加并沒有阻止孔隙流速的增加,這也間接說明管涌發(fā)展是加速發(fā)展的過程。由于孔隙流速一直呈增長趨勢,因此在一維管涌過程中沒有發(fā)生顆粒的沉積過程。土體中剩余的可動顆粒隨管涌的發(fā)展逐漸減少,在管涌初期,可動顆粒隨時間基本呈線性減小,而當(dāng)可動顆粒所剩很少時,其減少速率降低,并越來越趨近于0。從圖9中同樣可以看出,水力梯度值越高,顆粒流失總量越多,管涌發(fā)展越快。4流失量變化結(jié)果與分析本文根據(jù)管涌發(fā)生發(fā)展的基本規(guī)律,利用質(zhì)量守恒定律推得了移動顆粒的連續(xù)性方程,并得到了顆粒運動與滲流場中流體流動的耦合模型。根據(jù)影響顆粒流失的相關(guān)因素,提出了顆粒流失量及沉積量的計算模型,并針對一維管涌的算例,利用有限差分法對模型進行了求解,得到了管涌發(fā)生后管涌呈加速發(fā)展的過程。單位體積內(nèi)的顆粒流失量呈先增加后減小的趨勢,孔隙率與滲透流速不斷增加。通過對不同水力梯度下管涌發(fā)展