有限體積法在二維穩(wěn)態(tài)對(duì)流擴(kuò)散方程中的穩(wěn)定性分析

趙文娟，李春光，楊 程(1.寧夏大學(xué)土木與水利工程學(xué)院，寧夏銀川750021;2.寧夏節(jié)水灌溉與水資源調(diào)控工程技術(shù)研究中心，寧夏銀川750021;.北方民族大學(xué)土木工程學(xué)院，寧夏銀川750021)

隨著計(jì)算機(jī)技術(shù)的提升，近幾年土壤水鹽運(yùn)移模型的研究工作已從一維逐漸過(guò)渡到多維。而求解這類模型的計(jì)算方法成為研究的重點(diǎn)。傳統(tǒng)有限差分法中常規(guī)離散格式會(huì)對(duì)土壤水鹽運(yùn)移方程中的擴(kuò)散項(xiàng)產(chǎn)生較大的影響，使得離散后方程組系數(shù)矩陣極易呈現(xiàn)病態(tài)，計(jì)算結(jié)果呈現(xiàn)出不穩(wěn)定的偽振蕩現(xiàn)象［1］。筆者利用有限體積法求解二維穩(wěn)態(tài)對(duì)流擴(kuò)散方程，通過(guò)數(shù)值模擬計(jì)算分析研究該數(shù)值方法的數(shù)學(xué)穩(wěn)定性，為模擬計(jì)算土壤水鹽運(yùn)移過(guò)程提供前提條件。

1 有限體積法

有限體積法［2］(Finite volume method)是將計(jì)算區(qū)域劃分為一系列不重復(fù)的控制體積，并且使每個(gè)網(wǎng)格點(diǎn)周圍有一個(gè)控制體積;將待解的微分方程對(duì)每一個(gè)控制體積積分，便得出一組離散方程。積分方程中每項(xiàng)都有明確的物理意義，從而使得方程離散時(shí)，對(duì)各離散項(xiàng)給出一定的物理解釋。現(xiàn)就二維非飽和土壤水分運(yùn)動(dòng)方程為例，對(duì)其采用有限體積方法進(jìn)行線性離散。

對(duì)式(1)在圖1所示的控制體內(nèi)積分，可得

圖1 網(wǎng)格節(jié)點(diǎn)控制體圖

由于全隱式格式可以在較大時(shí)間步長(zhǎng)保持計(jì)算結(jié)果的穩(wěn)定性［3］，采用全隱格式對(duì)擴(kuò)散項(xiàng)進(jìn)行時(shí)間上的積分，即

按節(jié)點(diǎn)場(chǎng)內(nèi)的變量重新整理，可得

劃分的網(wǎng)格節(jié)點(diǎn)構(gòu)成的方程組再利用相應(yīng)的迭代法進(jìn)行計(jì)算求解。方程的數(shù)值解可通過(guò)增加迭代次數(shù)逼近到真實(shí)值。

有限體積法在計(jì)算偏微分方程中不僅可以得到較高精度的數(shù)值解，而且穩(wěn)定性好。這些特點(diǎn)使其成為當(dāng)前求解流動(dòng)問(wèn)題數(shù)值計(jì)算方法中最成功的方法［4］。1979年皇甫貴真［3］利用有限體積法求解Navier-Stokes方程組，解決了平板激波－附面層干擾和跨音速翼型黏性繞流問(wèn)題。至今，有限體積法被分別應(yīng)用到跨聲速流場(chǎng)的數(shù)值計(jì)算［5－6］、地下水中污染質(zhì)運(yùn)移數(shù)值模擬［7］、河床變形數(shù)值模擬［8］以及結(jié)構(gòu)動(dòng)力學(xué)和結(jié)構(gòu)可靠性等固體力學(xué)問(wèn)題［9］中。

2 數(shù)值模擬

利用有限體積法分別模擬計(jì)算二維穩(wěn)態(tài)對(duì)流擴(kuò)散問(wèn)題，驗(yàn)證該方法的數(shù)值解的穩(wěn)定性。

選取網(wǎng)格節(jié)點(diǎn)數(shù)為8、20、50、100，將模擬計(jì)算得到的數(shù)值解與精確解進(jìn)行對(duì)比。由圖2可知，相對(duì)誤差隨網(wǎng)格節(jié)點(diǎn)數(shù)目的增加而降低，范圍控制在1.2% ～6.9%。網(wǎng)格加密可以獲得高精度的數(shù)值解。當(dāng)網(wǎng)格節(jié)點(diǎn)較稀疏時(shí)，計(jì)算結(jié)果與精確值仍可以保持較高的吻合性。這說(shuō)明有限體積法下大步長(zhǎng)條件下是可以保障一定計(jì)算精度的。

2.2 算例2 已知二維穩(wěn)態(tài)對(duì)流擴(kuò)散方程為:

其中，ρu=5，ρv=2，Γ =1，Δx= Δy=0.2 m，計(jì)算區(qū)域東、南、西、北的邊界值為固定值，分別為50、300、100和200。采用有限體積法對(duì)計(jì)算區(qū)域(圖3)內(nèi)的節(jié)點(diǎn)進(jìn)行積分離散，最后得到二維穩(wěn)態(tài)對(duì)流擴(kuò)散方程計(jì)算節(jié)點(diǎn)的離散系數(shù)值和源匯項(xiàng)數(shù)值，迭代求解計(jì)算。

圖2 不同網(wǎng)格數(shù)條件下的數(shù)值解與精確解對(duì)比

圖3 計(jì)算區(qū)域圖

模擬結(jié)果見(jiàn)圖4。圖a為9個(gè)計(jì)算節(jié)點(diǎn)模擬曲面圖。8號(hào)節(jié)點(diǎn)溫度在模擬計(jì)算后為最高值，2號(hào)節(jié)點(diǎn)的溫度值為最小，與所給出的邊界條件相吻合;當(dāng)將網(wǎng)格節(jié)點(diǎn)加密，數(shù)目為25個(gè)時(shí)，模擬結(jié)果見(jiàn)圖b，東側(cè)接近邊界的計(jì)算節(jié)點(diǎn)模擬曲線高出西側(cè)的節(jié)點(diǎn)曲線值，計(jì)算斷面呈現(xiàn)東高西底的斷面形態(tài)。繼續(xù)加密網(wǎng)格節(jié)點(diǎn)到49和100個(gè)(圖c、d)。計(jì)算曲面隨節(jié)點(diǎn)數(shù)目的增加逐漸穩(wěn)定和光滑，且從計(jì)算結(jié)果呈現(xiàn)該區(qū)域溫度在位置坐標(biāo)(0.4，0.65)處為最低值為78，在位置坐標(biāo)(0.8，0.56)處為節(jié)點(diǎn)最高值為300。在不同網(wǎng)格節(jié)點(diǎn)條件下，有限體積法模擬計(jì)算二維穩(wěn)態(tài)對(duì)流擴(kuò)散方程的數(shù)值收斂好。模擬結(jié)果驗(yàn)證了該數(shù)值方法是可用于實(shí)際工程方面的。

3 結(jié)論

采用有限體積法求解二維穩(wěn)態(tài)對(duì)流擴(kuò)散問(wèn)題是可行的。通過(guò)算例驗(yàn)證，發(fā)現(xiàn)在加密網(wǎng)格條件下，數(shù)值解與精確解吻合度高且可較詳實(shí)地呈現(xiàn)模擬形態(tài);在稀疏網(wǎng)格條件下，數(shù)值解與精確值吻合較好且穩(wěn)定性強(qiáng)，表明有限體積法具有良好的計(jì)算效果，可為模擬二維非飽和土壤水鹽運(yùn)移提供數(shù)學(xué)參考。

圖4 不同網(wǎng)格數(shù)目條件下溫度場(chǎng)模擬值

［1］錢凌志，馮新龍.對(duì)流占優(yōu)擴(kuò)散問(wèn)題的特征AGE方法［J］.計(jì)算數(shù)學(xué)學(xué)報(bào)，2011，33(4):347 －357.

［2］VERSTEEG H K，MALALASEKERA W.An introduction to Computational Fluid Dynamics:The Finite V Method［M］.London:Longman Group Ltd，1995.

［3］姜啟源，謝金星，葉俊.數(shù)學(xué)模型［M］.北京:高等教育出版社，2011.

［4］李人憲.有限體積法基礎(chǔ)［M］.北京:國(guó)防科技出版社，2008.

［5］皇甫貴真.有限體積－時(shí)間分裂差分法及其對(duì)粘性流動(dòng)的應(yīng)用［J］.航空計(jì)算技術(shù)，1979(1):75－115.

［6］周新海，朱方元，劉松齡.跨聲速流場(chǎng)計(jì)算的時(shí)間分裂有限體積法［J］.空氣動(dòng)力學(xué)學(xué)報(bào)，1981(3):82－91.

［7］PUTTI M，YEH W，MULDER W A.用三角形有限體積法高密迎風(fēng)格式求解地下水運(yùn)移方程［J］.世界地質(zhì)，1992(2):157－172.

［8］陸永軍，張華慶.平面二維河床變形的數(shù)值模擬［J］.水動(dòng)力學(xué)研究與進(jìn)展(A輯)，1993(3):78 －79.

［9］陳浩.有限體積法在結(jié)構(gòu)動(dòng)力及可靠性分析中的應(yīng)用［D］.哈爾濱:哈爾濱工程大學(xué)，2012.