中科大計算熱物理講義04擴(kuò)散方程的數(shù)值方法

4.1一維導(dǎo)熱[1,2,3]可變下，引入面積函數(shù)因子A=A(x)，一維非穩(wěn)態(tài)導(dǎo)熱的通用控制方其中S為源項，λ是導(dǎo)熱系數(shù)。面積函數(shù)因子A(x)對不同情況和不同坐標(biāo)，有不同表示形將(4-1)兩端乘以A(x)，在t至t+Δt時間間隔內(nèi)對圖3－2所示的控制體P作積分將源項線化處理為S=SC+SPT；對上式的非定常項和源項的空間積分，取函數(shù)沿空間為為書寫方便，去掉時間上角標(biāo)t+Δt，而上角標(biāo)t改寫為0。于是，上述積分式可以化成aP=faE+faW+a-fSPAP(Δx)P顯式aPTP=aET+aWT+baP=aE+aW+a-SPAP(Δx)P,b=aT+SCAP(Δx)P(4-9b)以上各格式中的系數(shù)aE,aW,a表達(dá)式均為(4-6)所示。非穩(wěn)態(tài)熱物理問題不同離散格式中aP=aE+aW-SPAP(Δx)P,b=SCAP(Δx)P(4-11b)與非穩(wěn)態(tài)問題的全隱離散格式（4-9）一樣，穩(wěn)態(tài)問題離散格式（4-11）也是包含三個未知無論穩(wěn)態(tài)或是非穩(wěn)態(tài)問題，離散格式系數(shù)anb(nb=W,P,E)之值都與其界面w或者e1.加權(quán)平均法：假定在節(jié)點P、E之間的λ隨x成線性分布，則有很大。但按加權(quán)平均法計算，λe≈λP(δx)e+(δx)e，如為均勻網(wǎng)格，λe≈λP/2，即此時2.調(diào)和平均法：該法在穩(wěn)態(tài)導(dǎo)熱、無源項、導(dǎo)熱系數(shù)為階梯分布情況下基于界面熱流按照界面當(dāng)量導(dǎo)熱系數(shù)的意義，在穩(wěn)態(tài)、無和導(dǎo)熱系數(shù)呈階梯變化的假定下導(dǎo)出的，但串聯(lián)熱阻迭加原則的適用性不應(yīng)受以上條件限按照調(diào)和平均法，一維導(dǎo)熱（穩(wěn)態(tài)和非穩(wěn)態(tài)）通用離散方程一種方法是把導(dǎo)熱系數(shù)的突變階躍面設(shè)為控制容積的界面，使用調(diào)和平均算法S=SC+SPTP形式的離散方程迭代求解收斂，要求源項線化斜率SP≤0。這樣就可使離散求解代數(shù)方程組的系數(shù)矩陣對角占優(yōu)，且SP越大，離散方程主對角系數(shù)aP越大，兩次迭代之間求解函數(shù)TP變化越小，收斂變慢，有利克服迭代過程可能的發(fā)散。反之SP小，可能引起發(fā)散。1)S=2-3T，取SC=2,SP=-32）S=4+7T，可取SC=4+7T*,SP=0；亦可取SC=4+10T*,SP=-33）S=3-5T4，通常選取SP為S=f(T)曲線在節(jié)點P的斜率，并由此定SC。依于是取SC=3+15T*4,SP=-20T*3離散方程組要能封閉適定，必須引入定解條件，即初始邊界條件。初始條件引入通常(a)點中心法（外節(jié)點法，A法）第一類邊界條件，給定邊界上的函數(shù)值TB。對點中心網(wǎng)格，將給定的邊界函線性插值，得到1、2兩節(jié)點的溫度與邊界溫度的一個關(guān)聯(lián)方程。對均分網(wǎng)格，有在處理邊界條件時，通常約定通過邊界流入系統(tǒng)的熱流值為正，流出為負(fù)。第二類邊界條件，給定邊界上的熱流值qB，即函數(shù)的導(dǎo)數(shù)值。按照熱流值約定符號，左邊界熱流表達(dá)式為右邊界為；而第三類邊界條件，給定的是邊界上熱流值qB與周圍溫度為Tf、換熱系數(shù)為h的流體介質(zhì)之間換熱關(guān)系式，即的補(bǔ)充節(jié)點法或是附加源項法，在引入邊界條件后，對于離散格式為（4-9）和（4-11）的綜合以上兩式消除T0，并利用(△x)B=δx/2，可以得到邊界節(jié)點二階精度的離散在邊界節(jié)點B對應(yīng)的半控制容積建立熱平衡關(guān)系，其右側(cè)界面用‘e’厚度(Δx)b=0。與點中心時一樣，建立該節(jié)點的代數(shù)方程，既可直接用有限差分法，也可對零控制體節(jié)點B建立熱平衡方程，仍以‘w’和‘e令邊界與節(jié)點1的距離為δxb，對第二類邊界條件，在溫度空間分布為線性假設(shè)下而(Δx)B=0，所以，添加的邊界節(jié)點離散方程為如果是第三類邊界條件，將qB=h(Tf-TB)代入上式，可得（4-22）TT/2(4-23)綜合以上兩式，消除T0，得（4-20如為第三類邊界，代qB=h(Tf-TB)入(4-20)，即得此法不需在邊界上設(shè)置節(jié)點并求出邊界上的待求函數(shù)值，而只是將邊界上的第二、第個零控制體節(jié)點‘B’以作為‘1’的左邊鄰節(jié)點。但要注意，由于所設(shè)虛擬邊界節(jié)點控制Ba1-aB=a2-SPδx綜合以上三式，可以得到在第二類邊界條件下，緊鄰邊____T2-SPδx,=SCδx+qB（4-25b）對第三類邊界，qB=h(Tf-TB)；而按Fourier導(dǎo)熱定律，于是-SPδx,b=SCδx線化了的代數(shù)方程組求解，通常有直接解法和迭代解法兩類，深入討論留到本書后面一維穩(wěn)態(tài)導(dǎo)熱問題或者一維非穩(wěn)態(tài)導(dǎo)熱的隱式格式，解域內(nèi)部節(jié)點離散都得到如下形程組，用的最多的是基于Gauss消元法的Thomas算法，又稱追趕法或三對角陣算法1,2,...N-1,N（4-28）TPTi將此式乘(?ai)，再與方程(4-28)相iQi?1P因此，求解過程歸結(jié)如下：取a1=0,cN=0，按照式(4-30)計算消元過程中的遞推系數(shù)P,Qi(i=1,2,...N)；再按照式(4-29),逐個回代計算Ti(N,N?1,...2,1)。4.2多維導(dǎo)熱[1,3,4]在ρc=常數(shù)情況下，直角坐標(biāo)系下二維非穩(wěn)態(tài)導(dǎo)熱方程為梯分布，線化源項中溫度隨時間和空間均為階梯分布，取全隱式格式，令線化源項S=SC+SPT，則積分得到aPTP=aETE+aWTW（4-33a）aP=aE+aW+aN+aS+a-SPΔxΔy,離散網(wǎng)格系統(tǒng)如圖4－4所示，東西為軸向x坐標(biāo)方導(dǎo)方法，亦可得到與(4-32)完全相同的離散方程，其相應(yīng)系數(shù)anb和非齊次項b的表達(dá)式則為λ(δx)e(δx)w(δr)n(δr)s取如圖4－5所示的離散網(wǎng)格系統(tǒng)，東西為周向θ坐標(biāo)方向，南北為徑向r坐標(biāo)方向。相應(yīng)系數(shù)anb和非齊次項b的表達(dá)式為以上三種坐標(biāo)的離散方程具有完全相同的形式，差別只在方程系數(shù)和非齊次項表達(dá)式改一下個別參數(shù)取值就行了的通用程序提供了基礎(chǔ)。我們只要把通用離散方程式(4-32)中的分析x-y,x-r,θ-r三種坐標(biāo)系數(shù)表達(dá)式的主要差別，一是東西方向，極坐標(biāo)中的角度坐標(biāo)θ是無量綱的，而其它坐標(biāo)相應(yīng)的x以在角度坐標(biāo)θ中引入一個有量綱的尺度因子SX=r，而規(guī)定其它坐標(biāo)相應(yīng)的尺度因子SX=1；二是在南北方向，直角坐標(biāo)中沒有半徑這個量，而其它坐標(biāo)系中卻有，為統(tǒng)一表xxθXyrrY1rrR11rrδθ(δX)(SX)δyR(ΔY)/SXrΔθRΔXRΔXΔY★aE λ(δx)e λeΔre λReΔY(SX)e2(δX)e▼aW λw(δx)wλ(δx)waN(δr)n(δr)nλRnΔX(δY)naSλsRsΔX(δY)s0aPρcRΔXΔY/Δt★bScRΔXΔY+aT★aPaE+aW+aN+aN+a-SPRΔXΔY★▼對軸對稱圓柱坐標(biāo)，Re,Rw代表的是rP取全隱式格式，則積分完成后可以整理成如λaE=(δx)e,aW=(δx)w,aNλaP=aE+aW+aN+aS+aT+aB+a-SPΔxΔyΔz(4-40c)P節(jié)點控制容積相應(yīng)的界面w-e,s-n,b-t分別對應(yīng)θ,r,x三個坐標(biāo)方向，控制體體采用類似直角坐標(biāo)系的積分離散推演方法，可以得到與（4-39）形式完相應(yīng)系數(shù)anb和非齊次項b的表達(dá)式為 就不規(guī)則計算區(qū)域邊界的處理問題作簡單類似推演內(nèi)部節(jié)點通用離散方程一樣，將上式移項合并aB=aE+aN+aS+a-SPΔxΔy/2(4-44b)對第三類邊界qB=h(Tf-TB)，將此代入以上兩式中，全一樣離散形式，其中，系數(shù)aE,aN,aS,a與第二類時的（4－44a）中形式一樣，但aB對應(yīng)控制體的大小為內(nèi)部節(jié)點控制體的---qB=qBxi+qByj，建立該控制體的熱量平衡關(guān)系，對第二類邊界，（4-47a）對第三類邊界，如設(shè)qBx=hx(Tfx-TB),qBy=hy(Tfy-TB)，其中參數(shù)新附加的下標(biāo)完全一樣離散形式，其中，系數(shù)aE,aS,a與第二類時的（4－47a）中形式一樣，但aB與baPTP=aETE+aBTB+(aP-aB)TP=aETE+aB(TB-TP)+aNTN+aSTS+b這里qB為從邊界進(jìn)入P控制體的熱流率，于是aPTP=aETE+aNTN+aSTS+b(4-49)aP=aP-aB=aE+aN+aS+a-SPΔxΔy,b=qBΔy+b(4-50)當(dāng)邊界條件為第三類時，qB=h(Tf-TB)；而由Fourier導(dǎo)熱定律，與非角點處不一樣的是緊鄰角點處的aPTP=aETE+aSTS+baPTP=aETE+aSTS+b以上討論的邊界條件處理方法均是針對規(guī)則的矩形求解區(qū)域問題。當(dāng)所研究的多維問對于曲線邊界，一般采用點中心網(wǎng)格，對于第二、第三類邊界條件，邊界節(jié)點的離散右兩面來自節(jié)點1、3的導(dǎo)熱Q1-2,直方向來自P節(jié)點的導(dǎo)熱QP-2，熱平衡方+Q1-2+Q3-2+QP-2=0（4-57）邊界條件q2=h(Tf-T2)、FourierP同一維問題一樣，多維問題離散后，方程也可能因為系數(shù)是求解函數(shù)的函數(shù)，而使方對以矢量形式寫出的線性代數(shù)方程組AT=b，所謂迭代求解，就是要構(gòu)造一個向量T(n)}，當(dāng)n→∞,它收斂于某T*=A?1b。為了要得到這個函數(shù)序列，一般認(rèn)為第n次迭代所得之值T(n)取決于A,b及前一次迭代之值T(n?1)，即T(n)=F(A,b,T(n?1))（4-60）所要求的收斂精度，迭代停止。如二維導(dǎo)熱問題，中或（4-61b）對計算節(jié)點總數(shù)為N的問題，ε1≈Nε。后一種方法從解域總體上限制誤差量，要求較前ε值一般取為10?3~10?6，選取主要依照具體問題所要求的精度，也要考慮所采用的迭代求解線性方程組的迭代解法類型很多，針對熱物理計算中用的較多方法，以下介紹點迭代、塊迭代和交替方向線迭代等幾類，而每類又分為簡單(Jacob)迭代此需要首先提及。為敘述方便，選擇直角坐標(biāo)下二維非穩(wěn)態(tài)導(dǎo)熱離散方程aPTP=aETE+aWTW+aN(ⅰ)簡單迭代：又稱Jacobi迭代。每個節(jié)點更新值由上一輪迭代所得到的與它相關(guān)TaP（4-62）凡已算出新值的鄰點，立即啟用新值。對離散方程（4-32如掃描方向從下至上，從左到T迭代的新值，而是將此值與上一輪迭代值做加權(quán)平均后的值作為TaP(4-64a)T (n+1)令TP (n+1)(ⅰ)簡單線迭代：對于離散方程(4-32)，如果逐列掃描，迭代式為 (n+1)之值與上一輪迭代之值作加權(quán)平均，所得值才算新一輪迭代的新值，令 (n+1)為一項(1-w)aPTP(n)就可以了。T不可約：如果n×n矩陣A不能通過行的次序調(diào)換和其相應(yīng)列的次序調(diào)換而成為其中A11,A22分別為r和n-r階的方陣（1<r<n則稱矩陣A不可約。顯然，若A為可約,則方程組Ax=b可以化為一些低階的方程組，表示某些求解函數(shù)不需要全解域關(guān)聯(lián)且其中至少有一個式子取不等號，則稱矩陣A對角占優(yōu)。對于導(dǎo)熱問題，無論(ⅱ)若方程組Ax=b的系數(shù)矩陣A正定，則Gauss-Seidel迭代必收斂。(ⅳ)若方程組Ax=b的系數(shù)矩陣A不可約，且具對角優(yōu)，松弛因子ω滿足4.3管道內(nèi)充分發(fā)展的對流換熱截面上流體的無量綱溫度分布也不一定滿足式（4-73這個稱之為復(fù)雜類型充分發(fā)展的流視為常數(shù)，試確定溫度均勻的流體沿管軸流入管道后在充分發(fā)展的管內(nèi)區(qū)域換熱的Nu數(shù)。物理假設(shè)：①流體不可壓，物性為常數(shù)；X=x/(R.Pe),η=r/R(4-79)其中Pe=2Rum/a，而a=λ/(ρcP)為導(dǎo)溫系數(shù)，于是有等式左端是X的函數(shù)，右端是η函數(shù)，因此只能是一個常數(shù)。又由于dTbdX與Tb-T∞的符號相反，該常數(shù)一定為負(fù)數(shù)，令其為-Λ,Λ>0，有對給定的定解條件，只有