專利名稱：地下水溶質(zhì)運(yùn)移數(shù)值模擬中減少誤差的方法
技術(shù)領(lǐng)域：
本發(fā)明涉及一種數(shù)值模擬中減少誤差的方法，具體是一種在地下水溶質(zhì)運(yùn)移數(shù)值模擬中減少誤差的方法。

背景技術(shù)：
解決環(huán)境污染控制問題的一個重要技術(shù)過程就是將所研究的環(huán)境系統(tǒng)行為抽象為數(shù)學(xué)模型，這是進(jìn)行定量研究工作的基礎(chǔ)。在水環(huán)境污染模擬問題中，對流-彌散方程是溶質(zhì)(污染物)遷移的基本方程，其求解結(jié)果的精度可能會影響對水環(huán)境系統(tǒng)的分析，其被廣泛用于對溶質(zhì)輸運(yùn)問題的刻畫。近年來，國內(nèi)外學(xué)者對對流-彌散方程的數(shù)值解法進(jìn)行了廣泛而深入的研究，提出了許多方法。這些方法各有利弊，當(dāng)對流-彌散方程中的彌散項占主導(dǎo)地位時，使用各種數(shù)值方法均能得到較為滿意的結(jié)果，反之，若對流項占主導(dǎo)地位，各種常見的數(shù)值解法都會遇到困難(孫訥正.地下水污染-數(shù)學(xué)模型和數(shù)值方法[M].北京地質(zhì)出版社，1989)。因此，對流-彌散方程的數(shù)值解通常被認(rèn)為是“令人困窘”的難題。
對流-彌散方程的數(shù)值法大體分三類Euler法、Lagrange法和Euler-Lagrange法(成建梅，胡進(jìn)武.飽和水流溶質(zhì)運(yùn)移問題數(shù)值解法綜述[J].水文地質(zhì)工程地質(zhì)，2003.99.)。各種數(shù)值方法都包含程度不同的數(shù)值彌散和振蕩兩類誤差，這兩類誤差有共同的形成原因。Pinder指出，當(dāng)-種差分格式使數(shù)值彌散最小時，就會出現(xiàn)解的振蕩；而當(dāng)控制了解的振蕩，就要以增大數(shù)值彌散作為代價(Pinder GF，&Gray，W.G.Finiteelement simulation in surface and subsurface hydrology[M].Academic Press，1977.)。常用的Euler法主要包括有限差分法，有限單元法(Hossain M.A.，&Barber，M.E.Optimized Petrov-Galerkin model for advective-dispersive transport[J].Applied Mathematics and Computation，2000.115(1)，1-10.)等。這類方法應(yīng)用比較方便，但易受數(shù)值振蕩和數(shù)值彌散的影響。本發(fā)明主要是在Euler法中的上游加權(quán)法基礎(chǔ)上進(jìn)行改進(jìn)。
為減少數(shù)值振蕩與彌散，提高計算精度，通常采用兩類方法。一類方法是在空間離散時加密網(wǎng)格，并減小時間步長，但這種方法有時在實踐上會有困難，特別是對于一些區(qū)域性的污染問題，過小的空間步長和時間步長有時難以做到(薛禹群，謝春紅.地下水?dāng)?shù)值模擬[M].北京科學(xué)出版社，2007.)；另外一類方法是使用更高精度的差分格式，如20世紀(jì)80年代發(fā)展起來的一類流體力學(xué)計算的差分格式TVD(Total VariationDiminishing)(水鴻壽.一維流體力學(xué)差分方法[M].北京國防工業(yè)出版社，1998.)，這類高階方法常常減少了數(shù)值彌散，但是卻引入了振蕩。又有學(xué)者提出了通量限制器來消除振蕩，保存陡的濃度鋒面。如由Leonard等人開發(fā)的ULTIMATE(對流運(yùn)移方程瞬時插值模擬的通用限制器)(Zheng C.，& Bennett G.D.Applied Contaminant TransportModeling[M].San FranciscoWiley Inter-science，2002.)，這種格式不會引入過多的數(shù)值彌散，且基本沒有振蕩。但在實際計算時耗時過多。
數(shù)值彌散一般認(rèn)為是由時間和空間導(dǎo)數(shù)差分后的二階截斷誤差引起的，其大小可以通過分析二階截斷誤差近似得到。數(shù)值彌散使一個以彌散系數(shù)為D的對流-彌散問題變成為一個有新的彌散系數(shù)D*的對流-彌散問題(Noorishad J.，Tsang C.F.，Perrochet，P.，&Musy，A.A perspective on the numerical solution of convection-dominatedtransport problemsA price to pay for the easy way out[J].Water ResourcesResearch，1992.28(2)，551-561.)。


發(fā)明內(nèi)容
本發(fā)明所要解決的技術(shù)問題是提供一種運(yùn)用對流-彌散方程對地下水溶質(zhì)運(yùn)移進(jìn)行數(shù)值模擬過程中，在保證無數(shù)值振蕩的前提下，可以較大程度地減少數(shù)值彌散的方法。
本發(fā)明所述地下水中溶質(zhì)運(yùn)移數(shù)值模擬中減少誤差的方法，包括以下步驟 1)從水文地質(zhì)概念模型出發(fā)，用對流-彌散方程描述地下水流運(yùn)動和溶質(zhì)運(yùn)移的基本特征，對流-彌散方程的一維溶質(zhì)運(yùn)移數(shù)學(xué)模型為 0≤x≤L，0≤t≤Tsum (1) C(x，0)＝Cx(0) C(0，t)＝C0，C(l，t)＝Cl 式中參數(shù)，C為溶質(zhì)(污染物)濃度，t為時間，DL為彌散系數(shù)，u為實際平均流速，Cx(0)為初始時刻x處的濃度，C0為邊界x＝0上的濃度，Cl為邊界x＝l上的濃度； 2)為了得到(1)式左端濃度的時間導(dǎo)數(shù)項的差分格式，對j時刻和j+1時刻在i點(diǎn)的濃度進(jìn)行泰勒展開，并引入時間權(quán)因子θ得 (2) 將(2)式右端第二項濃度對時間的二階導(dǎo)數(shù)項可以轉(zhuǎn)化成對空間的導(dǎo)數(shù) (3) (3)式是個加權(quán)六點(diǎn)格式，當(dāng)θ＝0時，為顯式格式，當(dāng)θ＝1/2時，為Crank-Nicolson格式，當(dāng)θ＝1時，為隱式格式； 3)對(3)式的對流項使用泰勒展開得到第i和i+1節(jié)點(diǎn)的表達(dá)式，并引入空間加權(quán)因子α得 (4) 上式右端第三項為二階截斷誤差； 4)將(3)、(4)式代回(1)式，且將彌散項用傳統(tǒng)差分格式代替，得 (5) 其中截斷誤差εL為 (6) 其中Cr為科朗數(shù)，(6)式中二階截斷誤差大小與時間加權(quán)因子θ，空間加權(quán)因子α，空間步長Δx，地下水實際平均流速u，以及時間步長t有關(guān)，且時間步長越長，空間步長以及地下水實際平均流速越大，則二階截斷誤差越大，亦即數(shù)值彌散越大； 令D′＝uΔx[(θ-0.5)Cr+(α-0.5)]，則 (6)式右端高于二階的截斷誤差，以及空間與時間的離散對計算結(jié)果的影響可以通過在D’上乘以一個因子體現(xiàn)出來，定義這個因子為NDF，且令Dnum＝NDF×D′ (7) (8) 將(8)代回上游加權(quán)差分格式，令修正彌散系數(shù)D*＝D-Dnum，得到新計算格式如下 (9) 本發(fā)明提出了一個新的參數(shù)數(shù)值彌散因子(NDF)。其取值范圍在0～1之間，根據(jù)不同的問題其取值不同。數(shù)值計算結(jié)果顯示，使用修正彌散系數(shù)對傳統(tǒng)的隱式格式有較大的改進(jìn)。在保持無數(shù)值振蕩的前提下，可以較大程度地減少數(shù)值彌散。本發(fā)明既不需要以加密網(wǎng)格為代價，又不需要引入復(fù)雜的差分格式，可以比較簡單、方便地提高模擬計算精度。



圖1對隱式差分格式使用新方法計算結(jié)果對比，圖1a的NDF＝1，圖1b的NDF＝0.7～0.8 圖2Peclet數(shù)為50時新舊方法計算結(jié)果對比，圖2a是解析法、傳統(tǒng)上游加權(quán)算法以及新方法(修正彌散系數(shù)取NDF＝1)計算結(jié)果對比，圖2b是解析法、傳統(tǒng)上游加權(quán)算法以及新方法(修正彌散系數(shù)取分別NDF＝1和NDF＝0.9)計算結(jié)果對比， 圖3Peclet數(shù)為100與150時新舊方法對比，圖3a是Peclet數(shù)為100，圖3b是Peclet數(shù)為150， 圖4本方法與舊方法計算結(jié)果對比圖(lst year)(b，c)和研究區(qū)示意(a)， 圖5本方法與舊方法計算結(jié)果對比圖(2st year)(b，c)和流場圖(a)， 圖6觀測孔1(a)與觀測孔2(b)處穿透曲線處新舊方法結(jié)果對比。

具體實施例方式 隨著經(jīng)濟(jì)社會的發(fā)展和環(huán)境保護(hù)工作力度的持續(xù)加大，海水入侵、油罐滲漏、填埋場滲濾、工業(yè)廢水排放等引起的地下水污染問題得到越來越為廣泛的關(guān)注和重視，地下水污染預(yù)防和治理成為科學(xué)研究的熱點(diǎn)和政府管理的重點(diǎn)問題之一。建立地下水溶質(zhì)(污染物)運(yùn)移模型，對溶質(zhì)在不同地質(zhì)和水動力條件下的運(yùn)移過程進(jìn)行模擬和預(yù)報，是地下水污染科學(xué)預(yù)防和治理的重要途徑和方法手段。模型建立和模擬過程主要包括以下步驟， ①建立水文地質(zhì)概念模型 對研究區(qū)進(jìn)行水地質(zhì)測繪、鉆探、試驗和長期觀測，查明地質(zhì)、水文地質(zhì)條件，獲取相關(guān)資料。對復(fù)雜的天然地質(zhì)體忽略一些和當(dāng)前無關(guān)和關(guān)系不大的因素，得到天然地質(zhì)體的概念模型。
②建立數(shù)學(xué)模型 從水文地質(zhì)概念模型出發(fā)，用一組數(shù)學(xué)關(guān)系式來刻畫它的數(shù)量關(guān)系和空間形式，以反映所研究地質(zhì)體的地質(zhì)、水文地質(zhì)條件和所要描述的地下水流運(yùn)動和溶質(zhì)運(yùn)移的基本特征，達(dá)到再現(xiàn)或復(fù)制一個實際地下水系統(tǒng)基本狀態(tài)的目的。
數(shù)學(xué)模型包含描述水流和溶質(zhì)運(yùn)移規(guī)律的偏微分方程(控制方程)、研究區(qū)的范圍和所有相關(guān)參數(shù)、描述流體初始狀態(tài)的初始條件和流體與周圍環(huán)境相互關(guān)系的表達(dá)式，即邊界條件。如考慮對流、彌散的一維溶質(zhì)運(yùn)移數(shù)學(xué)模型如下， 控制方程 C(x，0)＝Cx (0) 初始條件 C(0，t)＝C0，C(l，t)＝Cl 邊界條件 式中參數(shù)，C為溶質(zhì)(污染物)濃度，t為時間，DL為彌散系數(shù)，u為實際平均流速，Cx(0)為初始時刻x處的濃度，C0為邊界x＝0上的濃度，Cl為邊界x＝l上的濃度。
③數(shù)學(xué)模型的識別和檢驗 第一步，對模型進(jìn)行求解。通常描述實際地下水溶質(zhì)運(yùn)移的模型都比較復(fù)雜，一般很難求得模型的解析解，需采用數(shù)值方法求解。數(shù)值方法有廣泛的適用性、經(jīng)濟(jì)的可行性、程序化的操作等優(yōu)點(diǎn)，但數(shù)值計算由于方法的計算誤差會引起解的振蕩和彌散現(xiàn)象，有效減小數(shù)值計算誤差，即充分抑制和弱化數(shù)值解的振蕩和彌散成為影響數(shù)值模擬結(jié)果和應(yīng)用的一個重要因素，這是本發(fā)明所要解決的問題。
第二步，模型識別。將模型通過數(shù)值模擬得到的結(jié)果和嘗試動態(tài)觀測資料或野外試驗(如抽水試驗、彌散試驗)對含水層施加某種影響所得的實際觀測結(jié)果進(jìn)行比較，看兩者是否一致，若不一致，就要對模型進(jìn)行修正，此即模型識別，也即通過數(shù)值方法不斷調(diào)整模型的輸入?yún)?shù)，直至模型輸出的模擬結(jié)果和野外觀測結(jié)果基本一致為止。
第三步，模型檢驗。為確保經(jīng)過識別后的模型能再現(xiàn)所研究的實際地質(zhì)體，還需在參數(shù)值不變的前提下根據(jù)另一時間段內(nèi)模型的計算結(jié)果與相應(yīng)時間段現(xiàn)場實際觀測資料的對比來進(jìn)一步檢驗、考核模型。
④模型運(yùn)轉(zhuǎn) 經(jīng)過識別檢驗后的模型說明它確實能代表所研究的地質(zhì)體中溶質(zhì)運(yùn)移發(fā)生的真實過程，表明它有能以足夠的精度預(yù)測未來的運(yùn)移狀況和特征，因而可以根據(jù)需要用這個模型來進(jìn)行相應(yīng)的預(yù)測(地下水污染發(fā)展趨勢)或計算，此即模型預(yù)報或模型運(yùn)轉(zhuǎn)。
本發(fā)明的步驟如下 在一維的情形下，笛卡爾直角坐標(biāo)系中，對流-彌散方程的理想形式(低濃度溶液，空隙率、地下水流速為常數(shù)、忽略吸附、離子交換、化學(xué)反應(yīng)、放射性衰變條件)為 0≤x≤L，0≤t≤Tsum(1) 為了得到上式左端濃度的時間導(dǎo)數(shù)項的差分格式，對j時刻和j+1時刻在i點(diǎn)的濃度進(jìn)行泰勒展開，并引入時間權(quán)因子θ得 (2) (2)式右端第二項濃度對時間的二階導(dǎo)數(shù)項可以轉(zhuǎn)化成對空間的導(dǎo)數(shù) (3) 這是個加權(quán)六點(diǎn)格式，當(dāng)θ＝0時，為顯式格式，當(dāng)θ＝1/2時，為Crank-Nicolson格式，當(dāng)θ＝1時，為隱式格式。對對流項使用泰勒展開得到第i和i+1節(jié)點(diǎn)的表達(dá)式，并引入空間加權(quán)因子α可得 (4) 上式右端第三項為二階截斷誤差。將(3)、(4)式代回(1)式，且將彌散項用傳統(tǒng)差分格式代替，得 (5) 其中截斷誤差εL為 (6) 其中Cr為科朗數(shù)。分析(6)式可知，二階截斷誤差大小與時間加權(quán)因子θ，空間加權(quán)因子α，空間步長Δx，地下水實際平均流速u，以及時間步長t有關(guān)。且時間步長越長，空間步長以及地下水實際平均流速越大，則二階截斷誤差越大，亦即數(shù)值彌散越大。
令D′＝uΔx[(θ-0.5)Cr+(α-0.5)]，則 (6)式右端高于二階的截斷誤差，以及空間與時間的離散對計算結(jié)果的影響可以通過在D’上乘以一個因子體現(xiàn)出來。定義這個因子為NDF(Numerical DispersionFactor)，且令Dnum＝NDF×D′ (7) (8) 將(8)代回上游加權(quán)差分格式，令修正彌散系數(shù)D*＝D-Dnum，得到新計算格式如下 (9) 1、一維算例研究 (1)一維濃度鋒面的推進(jìn)問題 設(shè)有一半無限長均質(zhì)砂柱，柱中原有溶質(zhì)濃度C0＝0，地下水在其中以實際平均流速u穩(wěn)定流動，在原點(diǎn)持續(xù)注入定濃度的溶質(zhì)，數(shù)學(xué)模型如下[10] (10) 對該模型進(jìn)行數(shù)值求解，參數(shù)選擇如下彌散系數(shù)DL＝0.04；實際平均流速u＝2m/d；研究區(qū)域為x＝120m的半無限空間；空間網(wǎng)格剖分Δx＝0.5m；時間步長Δt＝0.05d，時間0≤t≤25。本問題是對流占優(yōu)問題，會形成較陡的濃度鋒面?？臻g加權(quán)因子α取1。本例中，表征對流和彌散相對大小的Peclet數(shù)Pe＝(u*Δx)/DL＝25?？评蕯?shù)Cr＝(u*Δt)/Δx＝0.2。
將計算結(jié)果與解析解(Wang H.Q.& Lacroix M.Optimal weighting in the finitedifference solution of the convection-dispersion equation[J].Journal ofHydrology，1997.200(1-4)，228-242.)對比如圖1所示，對隱式格式使用修正的水動力彌散系數(shù)，且令參數(shù)NDF＝1進(jìn)行計算，得到的結(jié)果雖然數(shù)值彌散大幅減少，濃度鋒面和解析解非常接近，但卻同時產(chǎn)生了明顯的數(shù)值振蕩(圖1a)。為消除數(shù)值振蕩，利用試算法，對參數(shù)NDF進(jìn)行調(diào)整。將數(shù)值彌散因子(NDF)分別取0、0.7、1進(jìn)行計算。當(dāng)NDF＝0時，格式變?yōu)閭鹘y(tǒng)的上游加權(quán)格式，當(dāng)NDF＝1，即是圖1a中的新方法。當(dāng)NDF＝0.7時，數(shù)值振蕩被消除，數(shù)值彌散也比原來的上游加權(quán)方法減小了許多。對本例，NDF在0.7～0.8之間可以保證結(jié)果無數(shù)值振蕩且數(shù)值彌散盡量小(圖1b)。
(2)一維均勻流場中初始濃度丘的輸運(yùn)問題 許多資料顯示，當(dāng)Peclet數(shù)小于4時，標(biāo)準(zhǔn)差分法無論上游加權(quán)還是中心加權(quán)都比較精確，可以放心使用。但當(dāng)Peclet數(shù)大于4時，也就是對流占優(yōu)時，數(shù)值彌散與振蕩會比較明顯，為了討論新的參數(shù)NDF與Peclet數(shù)的關(guān)系，本發(fā)明模擬計算了一個經(jīng)典的初始濃度丘輸運(yùn)問題。
對于無限長均質(zhì)砂柱中的穩(wěn)定均勻流場，初始條件為形狀如高斯分布的單鋒波形，本例所用參數(shù)X0＝20，α0＝8，Δt＝0.05，u＝2，DL＝0.04，時間0＜t＜25，研究區(qū)長度120m。數(shù)學(xué)模型為 (11) 解析解為其中 應(yīng)用修正彌散系數(shù)代替?zhèn)鹘y(tǒng)的水動力彌散系數(shù)進(jìn)行計算，討論P(yáng)eclet數(shù)分別為50、100、150時的模擬結(jié)果。取Δx＝1，Pe＝50時，在濃度丘的后端產(chǎn)生不明顯的振蕩，將NDF調(diào)至0.9，振蕩現(xiàn)象被消除(圖2)。隨著Peclet數(shù)增加，上游加權(quán)法的誤差逐漸增大，當(dāng)Δx＝2，Pe＝100時，調(diào)節(jié)NDF的大小(圖3a)，通過試算法，NDF取0.8～0.9之間的值，可以得到較為滿意的結(jié)果。當(dāng)Δx＝3，Pe＝150時，NDF取在0.7～0.8之間，可以完全消除數(shù)值振蕩，并得到較好結(jié)果(圖3b)。
2、二維算例研究 對二維問題的研究，本發(fā)明結(jié)合美國Alabama大學(xué)的Zheng和Wang所開發(fā)的MT3DMS源代碼來完成。對MT3DMS中的DSP4CF子程序進(jìn)行了修改，在其中加入了新的數(shù)值彌散項。
本算例研究區(qū)為一矩形區(qū)域，東西方向長1400m，南北方向長1800m，垂直方向10m，在研區(qū)內(nèi)有注水井與抽水井各一口(圖4a)。模型的東西邊界為隔水邊界；北部為給定水頭邊界，值為250m；南部也為給定水頭邊界，取值如圖4a。承壓含水層的滲透系數(shù)為1.474×10-4m/s。把一定濃度的水通過完整井注入到含水層中，同時位于下游的抽水井從含水層中將溶質(zhì)抽出。在兩口井之間存在一個滲透系數(shù)為1.474×10-7m/s的低滲透系數(shù)區(qū)。地下水流場恒定，水流大致流向西南。兩口井抽注水率保持恒定，第一年持續(xù)注入溶質(zhì)濃度57.78ppm，在第二年停止注入。
為達(dá)到較好的模擬結(jié)果，本發(fā)明利用試算法對NDF參數(shù)進(jìn)行了分區(qū)(圖4a)，濃度陡鋒面處容易出現(xiàn)數(shù)值振蕩，故將NDF值取得較小，處理為低NDF區(qū)，保證振蕩被消除；其余部分濃度梯度較小，可處理為高NDF區(qū)，盡可能減小數(shù)值彌散。對本例，低NDF區(qū)取值0.3；高NDF區(qū)取值0.8。
兩個應(yīng)力期模擬結(jié)果表明，新方法在保證無數(shù)值振蕩的條件下，數(shù)值彌散比原上游加權(quán)方法有大幅減少，計算精度明顯提高。傳統(tǒng)上游加權(quán)法誤差較大，濃度鋒面由于誤差被大幅拉平。新方法保證了較陡的濃度鋒面(圖4，圖5)，并與TVD(ULTIMATE)方法擬合很好。計算表明，合理劃分NDF參數(shù)分區(qū)可以使新方法得到較好模擬結(jié)果。
本算例的觀測孔2位于注水井附近，靠近濃度鋒面，觀測孔1則遠(yuǎn)離濃度鋒面(圖4a)。對比觀測孔1和觀測孔2處不同方法所得的濃度穿透曲線可知，新方法對靠近濃度鋒面處的改進(jìn)效果非常明顯，可以得到較陡的濃度鋒面(圖6b)，而對遠(yuǎn)離濃度鋒面處的模擬結(jié)果也有一定改進(jìn)(圖6a)，但改進(jìn)效果略遜于鋒面較陡處。
本發(fā)明首先使用所提出的新方法對上游加權(quán)方法進(jìn)行改進(jìn)，并對兩個對流占優(yōu)的一維對流-彌散問題進(jìn)行研究。提出了一個新的參數(shù)數(shù)值彌散因子(NDF)。其取值范圍在0～1之間，根據(jù)不同的問題其取值不同。數(shù)值計算結(jié)果顯示，使用修正彌散系數(shù)對傳統(tǒng)的隱式格式有較大的改進(jìn)。在保持無數(shù)值振蕩的前提下，可以較大程度地減少數(shù)值彌散。
本發(fā)明對一個非均質(zhì)含水層中的二維溶質(zhì)運(yùn)移問題進(jìn)行了模擬計算。計算中對NDF參數(shù)進(jìn)行分區(qū)處理。計算結(jié)果表明，對二維溶質(zhì)運(yùn)移問題，新方法同樣比傳統(tǒng)方法有更少的數(shù)值彌散，精度更高，進(jìn)一步說明了新方法的可靠性。
權(quán)利要求
1、一種地下水溶質(zhì)運(yùn)移數(shù)值模擬中減少誤差的方法，其特征在于包括以下步驟
1)從水文地質(zhì)概念模型出發(fā)，用對流-彌散方程描述地下水流運(yùn)動和溶質(zhì)運(yùn)移的基本特征，對流-彌散方程的一維溶質(zhì)運(yùn)移數(shù)學(xué)模型為
0≤x≤L，0≤t≤Tsum(1)
C(x，0)＝Cx(0)
C(0，t)＝C0，C(l，t)＝Cl
式中參數(shù)，C為溶質(zhì)濃度，t為時間，DL為彌散系數(shù)，u為實際平均流速，Cx(0)為初始時刻x處的濃度，C0為邊界x＝0上的濃度，Cl為邊界x＝l上的濃度；
2)為了得到(1)式左端濃度的時間導(dǎo)數(shù)項的差分格式，對j時刻和j+1時刻在i點(diǎn)的濃度進(jìn)行泰勒展開，并引入時間權(quán)因子θ得
將(2)式右端第二項濃度對時間的二階導(dǎo)數(shù)項轉(zhuǎn)化成對空間的導(dǎo)數(shù)
(3)式是個加權(quán)六點(diǎn)格式，當(dāng)θ＝0時，為顯式格式，當(dāng)θ＝1/2時，為Crank-Nicolson格式，當(dāng)θ＝1時，為隱式格式；
3)對(3)式的對流項使用泰勒展開得到第i和i+1節(jié)點(diǎn)的表達(dá)式，并引入空間加權(quán)因子α得
上式右端第三項為二階截斷誤差；
4)將(3)、(4)式代回(1)式，且將彌散項用傳統(tǒng)差分格式代替，得
其中截斷誤差εL為
其中Cr為科朗數(shù)，(6)式中二階截斷誤差大小與時間加權(quán)因子θ，空間加權(quán)因子α，空間步長Δx，地下水實際平均流速u，以及時間步長t有關(guān)，且時間步長越長，空間步長以及地下水實際平均流速越大，則二階截斷誤差越大，亦即數(shù)值彌散越大；
5)令D′＝uΔx[(θ-0.5)Cr+(α-0.5)]，則
(6)式右端高于二階的截斷誤差，以及空間與時間的離散對計算結(jié)果的影響通過在D’上乘以一個因子體現(xiàn)出來，定義這個因子為NDF，且令Dnum＝NDF×D′，則
6)將(8)式代回上游加權(quán)差分格式，令修正彌散系數(shù)D*＝D-Dnum，得到新計算格式如下
。
全文摘要
本發(fā)明公開了一種地下水溶質(zhì)運(yùn)移數(shù)值模擬中減少誤差的方法，其是在水動力彌散系數(shù)上加上一個數(shù)值彌散估算值，得到一個修正彌散系數(shù)，用其替代方程中有明確物理意義的水動力彌散系數(shù)進(jìn)行計算。并提出了一個新的參數(shù)—數(shù)值彌散因子(NDF)，可以根據(jù)研究需要進(jìn)行參數(shù)分區(qū)并適當(dāng)調(diào)節(jié)該因子的大小，從而達(dá)到控制數(shù)值振蕩，減小數(shù)值彌散的目的。從一維到二維的多個數(shù)值算例的模擬計算結(jié)果表明，新方法能在消除數(shù)值振蕩的基礎(chǔ)上，較好地減少數(shù)值彌散，達(dá)到令人滿意的精度。
文檔編號G01N33/18GK101582096SQ20091003359
公開日2009年11月18日 申請日期2009年6月24日 優(yōu)先權(quán)日2009年6月24日
發(fā)明者吳吉春, 祝曉彬, 一 梅 申請人:南京大學(xué)