基于有限差分法的凍結(jié)過程數(shù)值分析

基于有限差分法的凍結(jié)過程數(shù)值分析

在中國(guó)東北部和西北部的寒冷地區(qū)，土木工程方興不斷，人工土壤冷凍技術(shù)在土木工程施工中發(fā)揮著越來越重要的作用。要分析自然凍土和人工凍土的冰凍工程，必須正確規(guī)定土壤溫度變化規(guī)律。例如，konrad分凝冰模型和daiyuwa的凍脹模型在計(jì)算凍脹前必須正確計(jì)算凍面的溫度場(chǎng)。然而，在使用有限差分法求解溫度場(chǎng)的過程中，由于時(shí)間間隔不正確，很容易導(dǎo)致“異常日期”。在這項(xiàng)工作中，我們提出了一種基于自適應(yīng)時(shí)間間隔、半隱藏和全隱藏格式的交替使用方法。這種方法可以有效。防止傳統(tǒng)計(jì)算過程中出現(xiàn)的相變遺漏,并且計(jì)算速度較之固定時(shí)間步長(zhǎng)方法有了很大的提高.1顯熱容法求解常見的一維土壤凍結(jié)過程熱傳導(dǎo)方程為,(1)式中:CV,λ為土壤體積熱容量(J/(m3·K))和熱導(dǎo)率(J/(m·s·K));L,T為融化潛熱(J/m3)和溫度(℃);ρi,θi為冰密度(kg/m3)和冰含量(m3/m3);z,t為空間坐標(biāo)(m)和時(shí)間(s).求解上述方程的一種有效方法是顯熱容法,即在發(fā)生相變的一個(gè)小的溫度范圍之內(nèi)構(gòu)造比熱函數(shù),以代替以上方程的右端第二項(xiàng),只以溫度為待求函數(shù),(2)其中:式中:C*V為等效熱容;λ*為等效導(dǎo)熱系數(shù);Cf,λf分別為已凍土的體積熱容和導(dǎo)熱系數(shù);Cu,λu分別為未凍土的體積熱容和導(dǎo)熱系數(shù);Tf-ΔT≤T≤Tf+ΔT為假設(shè)小的相變溫度范圍.2相變影響下的半隱格式方程(2)是與時(shí)間和空間都有關(guān)系的發(fā)展方程,采用有限差分法能對(duì)其有效求解.在常用的有限差分格式中,顯格式精度最差,隱格式次之,而Crank-Nicolson格式精度較高;但是從解的穩(wěn)定性角度分析,則隱格式的穩(wěn)定性最好,因此本文采用隱格式進(jìn)行求解.根據(jù)方程系數(shù)的取法不同,隱格式又有全隱格式和半隱格式之分.全隱格式公式如下式中j,n分別為空間和時(shí)間坐標(biāo)標(biāo)記.所謂半隱格式是將式(3)中系數(shù)C*v和λ*的時(shí)間步n+1換成n,通過半隱格式只需求解線性方程組即可.事實(shí)上,半隱格式和全隱格式對(duì)本方程的求解精度影響不大可忽略不計(jì).作者也分別編制了兩種不同格式的求解程序進(jìn)行了驗(yàn)證.由式(2),(3)可見,相變的影響在模擬中是通過在假設(shè)相變范圍的等效熱容來實(shí)現(xiàn)的,因此要全面考慮相變的影響,必須保證在每一計(jì)算時(shí)間步內(nèi)所有發(fā)生相變的空間單元的等效熱容能夠包含相變潛熱項(xiàng).而土體空間單元的等效熱容是否應(yīng)該包含相變潛熱項(xiàng),則是通過該單元的溫度值來判斷的.由式(2)知,只有其溫度落在Tf-ΔT≤T≤Tf+ΔT范圍之內(nèi),該單元的相變影響才能被考慮.這樣,就不免會(huì)出現(xiàn)‘相變遺漏’的情況.比如空間某單元在某一時(shí)間步之內(nèi)溫度由T1>Tf+ΔT變?yōu)門2<Tf-ΔT,顯然該單元在這一時(shí)間步內(nèi)發(fā)生了相變,但是根據(jù)T1或T2都無(wú)法判斷出該單元發(fā)生相變這一事實(shí),因此單元等效熱容也就沒能包含本應(yīng)包含大相變潛熱項(xiàng).由于半隱格式是根據(jù)時(shí)間步初的溫度值、全隱格式是根據(jù)時(shí)間步末的溫度值來判斷等效熱容是否包含相變潛熱項(xiàng)的,因此對(duì)于半隱和全隱格式,一旦相應(yīng)的溫度判據(jù)即時(shí)間步初的溫度值和時(shí)間步末的溫度值沒能落在Tf-ΔT≤T≤Tf+ΔT范圍之內(nèi),該單元的相變影響才能被考慮,否則就會(huì)出現(xiàn)下述的‘相變遺漏’現(xiàn)象.例如在任意一個(gè)計(jì)算時(shí)間步內(nèi),假設(shè)計(jì)算空間上某一點(diǎn)在時(shí)間步初的溫度值為T1,在時(shí)間步末的溫度值為T2.如果時(shí)間步長(zhǎng)控制不好,不同的差分格式就會(huì)有不同的‘相變遺漏’.采用半隱格式,若T1>Tf+ΔT而Tf-ΔT≤T2≤Tf+ΔT(時(shí)段初未凍而時(shí)段末處于過渡期)則相變作用就會(huì)遺漏;采用全隱格式,若Tf-ΔT<T1<Tf+ΔT而T2<Tf-ΔT相變作用也會(huì)遺漏.但是如果T1>Tf+ΔT且T2<Tf-ΔT則不論采用半隱或是全隱格式相變作用都會(huì)被遺漏掉.可見如果僅僅采用一種差分格式,那么無(wú)論是半隱格式還是全隱格式,對(duì)于本文所述問題都會(huì)不可避免地存在‘相變遺漏’這一事實(shí),因此是不正確和不可取的.由于凍結(jié)過程伴隨著大量相變潛熱的釋放,所以原來采用固定格式算法得出溫度將會(huì)偏低,而凍結(jié)鋒面推進(jìn)速度則會(huì)偏快.3速度速的數(shù)值分析步驟為了避免數(shù)值分析中相變作用的遺漏,必須有效地控制時(shí)間步長(zhǎng).本文根據(jù)半隱格式和全隱格式的不同對(duì)方程的求解精度影響不大這一特點(diǎn),提出了自調(diào)節(jié)時(shí)間步長(zhǎng)、半隱和全隱格式交替使用的方法(ModifiedFDM以下簡(jiǎn)稱MFDM),其中半隱格式直接求解線性方程組即可,全隱格式則要采用迭代法求解.由于本文提出的方法采用的是穩(wěn)定的計(jì)算格式,而且采用了自調(diào)節(jié)時(shí)間步長(zhǎng)的技術(shù),所以該方法具有穩(wěn)定和計(jì)算速度快的特點(diǎn).采用該方法進(jìn)行數(shù)值分析的步驟如下:1)對(duì)第一時(shí)間步的計(jì)算進(jìn)行特殊處理:第一時(shí)間步采用較小的步長(zhǎng),由于第一時(shí)間步內(nèi)溫度變化比較劇烈,而且土體從初始的正溫開始變化,所以先采用半隱格式進(jìn)行計(jì)算,然后令所有溫度落在T<Tf+ΔT的空間點(diǎn)的溫度為Tf,再對(duì)第一時(shí)間步進(jìn)行全隱格式計(jì)算;2)將上一時(shí)間步末的空間點(diǎn)按其溫度值分為兩類:一類是溫度落在T≤Tf+ΔT的空間點(diǎn),另一類是溫度大于Tf+ΔT的點(diǎn);同時(shí)采用與上一時(shí)間步相同步長(zhǎng)的時(shí)間步,用半隱格式對(duì)本時(shí)間步進(jìn)行試算;3)對(duì)第2)步中劃分出的溫度大于Tf+ΔT的點(diǎn)的試算結(jié)果進(jìn)行判斷:如果這些點(diǎn)中的溫度最小值(在自上而下凍結(jié)的情況下,此點(diǎn)應(yīng)為這些點(diǎn)中最上端的點(diǎn))小于Tf-ΔT,就減小時(shí)間步長(zhǎng),再用半隱格式對(duì)此時(shí)間步進(jìn)行試算;若此最小值大于Tf+ΔT,則增大時(shí)間步長(zhǎng),再用半隱格式對(duì)該時(shí)間步進(jìn)行試算,直至該最小值落在[Tf-ΔT,Tf+ΔT]之中;4)用第3)步中確定下的最終時(shí)間步長(zhǎng)對(duì)2)中劃分出的溫度大于Tf+ΔT的點(diǎn)進(jìn)行全隱格式計(jì)算;5)重復(fù)2)～4)步,直至計(jì)算時(shí)間結(jié)束.可以設(shè)想,采用了上述改進(jìn)辦法,相變遺漏避免了,對(duì)于凍結(jié)問題也就是說釋放出來的相變潛熱被全面地考慮了.采用固定格式算法得出的溫度值和凍結(jié)鋒面都會(huì)得到相應(yīng)的糾正,溫度會(huì)相應(yīng)升高,鋒面推進(jìn)會(huì)相應(yīng)減緩.4時(shí)間步初的溫度值采用文獻(xiàn)提供的基于FEM的經(jīng)典算例.空間總長(zhǎng)為1m,計(jì)算總時(shí)間為2.88×105s,采用第一類邊界條件,TC=-20℃,TW=10℃;初始溫度T0=10℃.參數(shù)如下:Tf=0℃,ΔT=1℃,L=338000000(J/m3),λf=2.22(J/(m·s·℃)),λu=0.556(J/(m·s·℃)),Cf=1762000(J/(m3·℃)),Cu=4226000(J/(m3·℃)).1)采用固定時(shí)間步長(zhǎng)半隱差分格式(FDM)求解.將空間等步劃分為400個(gè)單元,時(shí)間等步劃分為192個(gè)時(shí)間步長(zhǎng).計(jì)算結(jié)果與文獻(xiàn)提供的結(jié)果幾乎相吻合.通過比較可以看出多處存在相變遺漏,其中計(jì)算時(shí)間初始階段出現(xiàn)較多,隨著計(jì)算時(shí)間的增加相變遺漏現(xiàn)象逐漸變少,見表1和表2.表1中在1.5×103～3×103s時(shí)間步之內(nèi),0.02～0.0275m處4個(gè)空間點(diǎn)顯然發(fā)生了相變,且它們?cè)跁r(shí)間步初的溫度值都大于,也即超出了考慮相變的溫度范圍,由于采用半隱格式,根據(jù)時(shí)間步初的溫度值判斷的結(jié)論是這4個(gè)空間點(diǎn)都沒有發(fā)生相變(表中黑體數(shù)字),因此它們的等效熱容將不會(huì)包含相變潛熱項(xiàng),也就是說,本時(shí)間步之內(nèi)這4個(gè)點(diǎn)的相變都被遺漏掉了.同樣在3×103～4.5×103s時(shí)間步之內(nèi),0.03m和0.0325m兩點(diǎn)的相變被遺漏掉了;在4.5×103～6×103s時(shí)間步之內(nèi),0.035m和0.0375m兩點(diǎn)的相變被遺漏掉了;在6×103～7.5×103s時(shí)間步之內(nèi),點(diǎn)0.04m處的相變被遺漏掉了.在表2中因?yàn)橄嘧儾幌癖?所示的時(shí)間初那樣劇烈,因此相變遺漏相應(yīng)減少,只有黑體數(shù)字表示的兩點(diǎn)發(fā)生相變遺漏.2)用MFDM求解將空間變步劃分為113個(gè)單元,初始時(shí)間步長(zhǎng)設(shè)為0.25s,采用此種方法只計(jì)算了68個(gè)時(shí)間步計(jì)算就結(jié)束了,計(jì)算結(jié)果見表3.由于對(duì)每一時(shí)間步內(nèi)黑體數(shù)據(jù)以上的點(diǎn)采用半隱格式,而黑體數(shù)據(jù)及其以下的點(diǎn)則采用全隱格式,所以‘相變遺漏’被有效避免了.比如在表3中的3.17×103～3.68×103s時(shí)間步之內(nèi),0.025～0.028m處3個(gè)空間點(diǎn)采用了半隱格式,點(diǎn)0.03m處則采用了全隱格式,它們相應(yīng)的溫度判據(jù)都落在了[Tf-ΔT,Tf+ΔT],因此實(shí)際發(fā)生的相變都被考慮到了.在3.17×103～3.68×103s時(shí)間步之內(nèi),點(diǎn)0.025m處已不存在相變,點(diǎn)0.026m和0.028m處均存在相變,對(duì)這兩個(gè)點(diǎn)必須用半隱格式才能考慮相變的影響,而點(diǎn)0.03m和0.031m處因?yàn)闀r(shí)間步初末的溫度值都在[Tf-ΔT,Tf+ΔT],采用半隱或全隱格式都能考慮到相變影響,點(diǎn)0.033m處也存在相變但只有時(shí)間步末的溫度值在相變范圍內(nèi),因此只能采用全隱格式.根據(jù)設(shè)計(jì)的算法在0.033m點(diǎn)以上各點(diǎn)采用半隱格式,在0.033m點(diǎn)(包括該點(diǎn))以下各點(diǎn)采用全隱格式,是可以滿足全部考慮相變影響要求的.圖1(橫軸表示距土體低溫端的空間距離,縱軸表示溫度值)比較了用FEM,FDM和本文方法計(jì)算所得結(jié)果中的兩個(gè)時(shí)刻的溫度場(chǎng).可以看出,FEM和FDM計(jì)算結(jié)果幾乎相同,而用本文方法計(jì)算的土體中同一點(diǎn)的溫度值則比FEM和FDM的均要高.這符合本文前述的分析和設(shè)想:傳統(tǒng)FEM和傳統(tǒng)的固定格式的FDM均存在相變遺漏,致使沒能全面考慮釋放的相變潛熱,因此計(jì)算出的溫度偏低.從圖1還可看出傳統(tǒng)算法和修正算法之間最大的誤差出現(xiàn)在凍結(jié)鋒面附近,這也是符合實(shí)際情況的.因?yàn)閷?shí)際凍結(jié)過程中相變發(fā)生最劇烈的地方就是凍結(jié)鋒面附近,固定格式算法遺漏相變最多的地方自然也在此處.由于本文提出的方法考慮了原來固定時(shí)間步長(zhǎng)法沒有考慮的‘相變遺漏’的影響,所以比傳統(tǒng)方法更加全面地考慮了水相變成冰釋放出來的潛熱,這必然延緩凍結(jié)鋒面的推進(jìn)速度.表4比較了采用半隱格式固定時(shí)間步長(zhǎng)方法和采用自調(diào)節(jié)時(shí)間步長(zhǎng)方法求出的凍結(jié)鋒面的變化情況.從圖2(橫軸表示凍結(jié)時(shí)間,