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

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