中科大計(jì)算熱物理講義03離散方法基礎(chǔ)

3.1解域離散[1,2](a)點(diǎn)中心法（外節(jié)點(diǎn)法，A法）表示法，即研究的空間節(jié)點(diǎn)位置記作(i,j,k),時間或類時間坐標(biāo)離散節(jié)點(diǎn)記為n,兩相鄰節(jié)點(diǎn)間界面分別對應(yīng)i±1/2,j±1/2,k±1/2,n±1/2。用有限差分方法進(jìn)行離散和定性分析，都使用大寫字母P表示所研究的節(jié)點(diǎn),E,W,N,S,T,B分別表示節(jié)點(diǎn)P周圍在三個坐標(biāo)方向上距離稱為網(wǎng)格間距。對于均分網(wǎng)格系統(tǒng)，即等距網(wǎng)格系統(tǒng)，δx=⊿x，不強(qiáng)調(diào)兩者間差異，i+,j,k2(a)點(diǎn)中心法（外節(jié)點(diǎn)法，A法）(a)點(diǎn)中心法（外節(jié)點(diǎn)法，A法）3.2微分方程的有限差分法離散??（3-1）（3-2a）y=常值=u(x+Δx,y)?u(x,y)y=y=常值?y=常值yExu=u+1,E1u=u?1（3-3）xu=u+1?u,Δtu=u+1?u（3-4）xu=u?u?1,▽tu+1=u+1?u（3-5）（4）一步中心差算子δ,定義為（3-6）?1,δ=E1/2?E?1/2(3.10a)EΔ▽δE122+24E?1122241224E2?E2?11??24由一離散點(diǎn)（j＋1,n將該點(diǎn)的函數(shù)值u+1在離散點(diǎn)（j,n）展開，有稱為離散截?cái)嗾`差（T.E—Trunca如果選擇的離散點(diǎn)是(j?1,n)，將該點(diǎn)的函數(shù)值u?1在離散點(diǎn)(j,n)展開，有（3-16）如果在離散點(diǎn)（j,n）的左、右分別選?。╦－1/2,n）和（j＋1/（3-17）以上差分近似相當(dāng)一步中心差算子δ作用于函數(shù)u，稱為一步中心差分格式，它的截差是二階小量。類似，如果在離散點(diǎn)（j,n）的左、右分別選取(j?1,n)和（j＋1,n）兩點(diǎn)，作（3-18）?u)ΔxO(Δx))ΔxO(Δx))ΔxO(Δx2))2ΔxO(Δx2)2ΔxO(Δx2)2ΔxO(Δx2)O(Δx3)O(Δx3)O(Δx4)2un2,j)Δx2O(Δx2)+u)Δx2O(Δx)+u)Δx2O(Δx)j?2j?1jj?2j?1j/(12Δx2)O(Δx4)O(Δx,Δy2)O(Δx2,Δy)O(Δx2,Δy)O(Δx2,Δy2)解：①在離散參考點(diǎn)（j,n取時間向前、空間兩步中心的差分格式，即FTCS)（3-20）即ForwardSpace)格式。此格式截差對（3-23）u+1=u?c(u+1?u),T.E=O(Δt,Δx)（3-24）這是一個兩時層的三點(diǎn)顯式格式，全一階精度。在a正，c>0，節(jié)點(diǎn)j+1在節(jié)點(diǎn)j的下游，F(xiàn)TFS格式的離散方程（3-24）表明，節(jié)點(diǎn)j在下一時刻n+1的值是由該節(jié)點(diǎn)及它的下游節(jié)點(diǎn)j+1在前一時刻n之值來決定。這是不合雙曲型此種情況下的FTFS格式稱為逆風(fēng)格式，逆風(fēng)格式數(shù)值上是不穩(wěn)定的，因此不如果a<0，流動自右向左，流速為負(fù)，c<0，節(jié)點(diǎn)j+1在節(jié)點(diǎn)j的上游，是求解點(diǎn)的依賴andBackwardSpace)格式。此格u+1=u?c(u?u?1),T.E=O(Δt,Δx)（3-26）這也是一個兩時層的三點(diǎn)顯式格式，全一階精度。在a>0,即c>0時，F(xiàn)TBS格式是迎風(fēng)格T.E=O(Δt,Δx2)（3-28）即u+1=u+α(u+1?2u+u?1),T.E=O(Δt,Δx2)（3-29）u+1=u?1+2α(u+1?2u+u?1),T.E=O(Δt2,Δx2)（3-32）T.E=O(Δt,Δx2)（3-33）u+1=u+α(u?2u+1+u),T.E=O(Δt,Δx2)將求解層各節(jié)點(diǎn)待求函數(shù)移至左邊，合并整?αu+(2α+1)u+1?αu=u,T.E=O(Δt,Δx2)（3-34）C-N格式,又稱算術(shù)平均格式。此格式截差對時間、空間均為二階))點(diǎn)函數(shù)值的系數(shù)，用dj表示右端已知值，則方程變?yōu)閍ju+bju+1+cju=dj,T.E=O(Δt2,Δx2)（3-37）解：如果a>0時，迎風(fēng)格式要求在(j,n)的上游方向找依賴區(qū)，因此的空間導(dǎo)數(shù)離散點(diǎn)應(yīng)選為(j,n),(j?1,n),(j?2,n),...。現(xiàn)要求空間離散為二階精度，各離散點(diǎn)上點(diǎn)(j,n),(j-1,n),(j-2,n)進(jìn)行離散，令（3-38）將u-2,u-1,u在（j,n）點(diǎn)作Taylor展開（3-40）比較等式兩邊u和u的各階導(dǎo)數(shù)的系數(shù)，有聯(lián)立求解關(guān)于系數(shù)aj、bj、cj的上述方程，得（3-42）同理，如果a<0，迎風(fēng)格式要求在(j,n),(j+1,n),(j+2,n)三點(diǎn)離散成二階精度，可類對時間導(dǎo)數(shù)項(xiàng)，要求做前差，精度為一階，故取兩個離散點(diǎn)(j,n),(j,n+1)即可。u(x,tn)=a+bx(3-47)不失一般性，令節(jié)點(diǎn)(j,n)的空間坐標(biāo)x0=0，則節(jié)點(diǎn)(j,n)及鄰近節(jié)點(diǎn)(j-1,n),(j+1,n)u(x,tn)?a+bx+cx2j,j-1,u+1=a+b2故在(j,n)點(diǎn)有這與用基于Taylor展開方法，在節(jié)點(diǎn)(j,n)對x作兩步中心差分和二階中心差分所得到的離散形式完全一樣。其截差同樣通過將u+1、u?1對點(diǎn)(j,n)作Taylor展開分析得到，是T.E=O(Δx2)。事實(shí)上，函數(shù)u(x,tn)的精確解在(3-51)表明，函數(shù)在x=0附近所作的二次多項(xiàng)式分布假設(shè)等同于該函數(shù)精確解的Taylor展 中計(jì)算的組合式總截差仍是O(Δx3)，計(jì)算的組合式總截差變?yōu)镺(Δx4)，在分別除以導(dǎo)數(shù)項(xiàng)的乘積系數(shù)Δx、Δx2之后例3.4已知二維導(dǎo)熱區(qū)域固壁邊界y=0上的及區(qū)域內(nèi)部的溫度，確定邊界y=0上的溫度分布。解：如圖3－6所示，設(shè)在xi處壁面附近溫度分布近似為二次多項(xiàng)式T(xi,y)?a+by+cy2令壁面y=0，對應(yīng)的離散節(jié)點(diǎn)為(i,1)，依次向上而二維空間（x,y）解域邊界y=0上的壁面熱流qw為辦法來解。設(shè)在xi處壁面附近溫度分布近似為三次多項(xiàng)式T(xi,y)?a+by+cy2+dy3同樣方法分析，所得熱流表達(dá)式的截差為O(Δy3)。每步推進(jìn)半個時間間隔Δt/2，第一步，空間導(dǎo)數(shù)均取二階中心差，但對x取隱式、對y取顯；第二步，續(xù)接前一步計(jì)算結(jié)果，空間導(dǎo)數(shù)亦取二階中心差，但對y取隱式、對x取顯：上述二維導(dǎo)熱問題的ADI格式恒穩(wěn)定，對時間步長沒直接推廣至三維，每個方向計(jì)算前進(jìn)三分之一個時間間隔Δt/3，所得到的相應(yīng)格式只能有法?，F(xiàn)針對最簡單的雙曲型對流方程（3－19）說明它的基本思想。設(shè)a<0___u+1=u-c(u+1-u)(3-57)___階中心差分，對流方程取時間向前、空間單邊三算過程中插入了一個半時層（n+1/2每個方程計(jì)算時向前推進(jìn)Δt/2時間間隔，兩步程（3-19僅對空間導(dǎo)數(shù)項(xiàng)離散，在a>0時，取迎風(fēng)格式，可得對j點(diǎn)離散的如下半在原始方程對應(yīng)的(x,t)平面上，沿每條x=xj直線，構(gòu)成3.3微分方程的有限容積法離散[1,4,12,13]（1）將守恒性控制微分方程在任意選定的控制容熱物理問題通用形式的守恒型微分方程（2-（1）非穩(wěn)態(tài)項(xiàng)：取φ隨x作階梯變化，即同一控(ⅰ)若取(ρuφ)隨t作階梯顯式變化，即在積分時間(ⅱ)若取(ρuφ)隨t作階梯全隱式變化，即在積分時(ⅲ)若取(ρuφ)隨t作分段線性變化，在積分時間間隔Δt下限值的平均值替代積分函數(shù)值，則有算術(shù)平均格式（(ⅰ)階梯顯式(ⅱ)階梯全隱式（4）源項(xiàng)：選S隨x作階梯變化，隨t作階梯或分段線性變化，則類似有(ⅰ)階梯顯式(ⅱ)階梯全隱式 t+Δtt這里S和+S分別表示t+Δt和t時刻源項(xiàng)在積分控制體內(nèi)的平均值。將以上各項(xiàng)積分表達(dá)式按對應(yīng)的三種不同格式分別組合，并去掉t+Δt，而把積分下限t記為0，則可得由控制容(ⅰ)階梯顯式(ⅱ)階梯全隱式如果網(wǎng)格為均分網(wǎng)格，即(Δx)P=(δx)e=(δx)w=Δx，并假設(shè)這與在相同假設(shè)條件下，用有限差分法從(j,n)點(diǎn)出發(fā)作FTCS、BTCS及Crank-Nicolson(ⅲ)型線選擇的多樣性，可得離散方程的多樣性。（P控制體）里，在Δt時間間隔內(nèi)函數(shù)φ值的增加或減少，應(yīng)等于同一（常密度ρ)及擴(kuò)散系數(shù)Γφ也為常值的簡化條件下，以上論述的數(shù)學(xué)表示就是方案：t時刻、t+Δt時刻、t+Δt/2時刻，分別對應(yīng)顯式、全隱式和C-N三種不同格式。其中下標(biāo)nb代表離散中用到的與P關(guān)聯(lián)的相鄰節(jié)點(diǎn)。一般情況下，對于一維問題，nb=E,W；對于二維問題nb=E,W,N,S；對于三維問題nb=E,W,N,S,T,B。當(dāng)然3）與求解函數(shù)相關(guān)的源項(xiàng)線化時需取負(fù)斜率——只有這樣，只有SP<0，對aP才是正貢獻(xiàn)，否則，有可能使aP變負(fù)。3.4有限差分法離散和有限容積法離散比較有限容積法離散，從守恒型方程出發(fā)進(jìn)行積分或基于控制容3.5離散格式的定性分析3.5.1誤差與精度[4]非穩(wěn)態(tài)問題，其相應(yīng)的微分方程為L(φ(xj,tn))=0，L表示對連續(xù)函數(shù)φ(xj,tn)在點(diǎn)（j,n）作某些微分運(yùn)算的算子；微分方程對應(yīng)的某個差分方程（差分格式）為LΔ(φjn)微分算子和差分算子之間的差，即T.E=L(φ(xj,tn))?LΔ(φjn)(3-75)方程的FTCS格式，截差為O(Δt,Δx2)，則器的有效字長，字長越長，R.E越??；要求運(yùn)算的次數(shù)，次數(shù)越多，R.ED.E=T.E+B.E(3-76)即C.E=N.E=D.E+R.E(3-77)它的長度來度量。范數(shù)則是這些簡單度量指標(biāo)的推廣和抽象，用雙絕對值符號.表示，是3.5.2離散格式的相容性[1,4]根據(jù)離散方程的截?cái)嗾`差定義，顯然，如果格式相容，則網(wǎng)格步長趨于零時（h→0）離對于自變量為(x,t)的方程，當(dāng)離散格式的截差形式為T.E=O(Δtm,Δxn)時，易從?1解：將格式各離散點(diǎn)的函數(shù)值φ在(j,n)點(diǎn)作Taylor展開，并將各展開式代入離散格式物理性質(zhì)和數(shù)學(xué)特征有很大差異，方程類型變異，當(dāng)然不能得到原3.5.3離散格式的收斂性和穩(wěn)定性[1,3,4,14,15]解φ(xj,tn)，即長和衰減來定義，這就為實(shí)際進(jìn)行穩(wěn)定性分析提供了方便，可以采用決，又稱之為vonNeumann分析法。同時還假定，應(yīng)用vonNeumann分析方法證明是弱穩(wěn)3.5.4初值問題離散格式穩(wěn)定性分析方法[1,3,4]應(yīng)用廣泛的vonNeumann方法。但這種方法分析的只是離散節(jié)點(diǎn)采用時間向前、空間二階中心的顯式FTCS差入誤差）為ε,引入誤差后實(shí)際能夠計(jì)算的數(shù)值解為N，顯然有N=D+ε。實(shí)際數(shù)值解N是按照格式演化規(guī)律得到的，即N+1=N+α(N+1-2N＋N-1)將N=D+ε代入上式，有D是離散方程的精確解，自然應(yīng)該滿足離散方程，即D+1=D+α(D+1-2D＋D-1)(3-82)ε+1=ε+α(ε+1-2ε+ε-1)（3-83）vonNeumann穩(wěn)定性分析方法是將一個時層上在網(wǎng)格點(diǎn)的誤差函數(shù)分布展成為諧波組成的有限Fourier級數(shù)，再考察其在時間演化過程中這些令誤差函數(shù)ε(x,t)可以在函數(shù)定義區(qū)域[?l,l]展成以下有限Fourier級數(shù)ε(x,t)=Σbmeikmx（3-84）m波的數(shù)目。如果將解域均分為M個子域，并設(shè)M為偶數(shù)，對應(yīng)M+1個網(wǎng)格節(jié)點(diǎn)，步長為Δx=2lM，則波數(shù)km可以寫成代表級數(shù)展開的常數(shù)項(xiàng)；而m=1,M2時，分別對應(yīng)最小和最大的波數(shù)而相應(yīng)的波長為最大和最小波長，最大波長（最小波εm(x,t)=bm(t)eikmxεm(x,t)=gneikmx（3-85）g=eaΔt（3-86）:(x,0)=g0eikmx=eikmx=εm(x,0)(x,Δt)=geikmx=eaΔteikmx=gεm(x,0)t=nΔt(n=n)εm(x,nΔt)=gneikmx=eateikmx=gεm(x,(n?1)Δt)(εm)+1=gn+1eikmxj=g(εm)=eaΔt(εm)ε+1=gn+1eikxj=gε=eaΔtε（3-87）aΔt例3－6試用vonNeumann方法分Tj2=gTj1=g2Tj0氵Tjn+1=gTjn=gn+1Tj0此不等式右端恒成立；左端，對任意(kΔx)，sin2≥0，要使不等式成立，應(yīng)有?1≤1?2α即例3－7試用vonNeumann方法分析二i,j初始誤差：T0=ei(kxxi+kyyi,j要使任意kxΔx,kyΔy滿足，有例3－8試用vonNeumann方法分別對兩個空間方向的一維導(dǎo)熱問題作C-Nxxαy=σΔt/Δy2i,j初始誤差：T0=ei(kxxi+kyyi,jT又δx2Tij,=Ti+1,j-2Ti,j＋Ti?1,j=?2(1?cos(kxΔx))Ti,jδy2Ti,j=Ti,j+1-2Ti,j＋Ti,j-1=?2(1?cos(kyΔy))Ti,j例3－9試用vonNeumann方法分析無源g≤1，也就是即判斷左端部分極值點(diǎn)：令A(yù)=kΔx/2，左端函數(shù)為F,由dFdA=0，得極值點(diǎn)處為sinA=0orcosA=0易見sinA極大，則cosA極??；cosA極大，則sinA極小。于是有令考慮邊界條件后，需求解的全部離散節(jié)點(diǎn)個數(shù)為J?1，對應(yīng)的求解函數(shù)在時層n的BUn+1=AUn+ForUn+1=GUn+D(3-88)其中A、B、G=B?1A均為(J?1)×(J?1)矩陣，F(xiàn)、D=B?1F為方程的非齊次項(xiàng)矢量。Bεn+1=Aεnorεn+1=Gεn(3-89)采用矩陣模的定義，要使矢量函數(shù)在推進(jìn)演化中不致于因?yàn)檎`差的引G≤1orGn≤K(K為常數(shù)）(3-90)如果G是正規(guī)矩陣，即滿足GG*=G*G，其中G*為G的轉(zhuǎn)置共軛陣；或者存在一個非例3-10用矩陣方法分析一維非穩(wěn)定導(dǎo)熱方程在第T(0,t)=f1(t),T(1,t)=f2(t)取點(diǎn)中心網(wǎng)格，并對x和t坐標(biāo)分別作等距分割起始左右邊界節(jié)點(diǎn)分別為(0,0),(J,0)。采用FTCS差分格式離散，離散方程為22對鄰接邊界的節(jié)點(diǎn)j=1和j=J-1,方程分別為TTn+1=GTn+Dnεn+1=Gεn?1?1j=1,2,...J?1解：按要求格式離散，并按vonNeuma令β=csin(kΔx)，上式合并，得G?λI離散格式在保證相容和穩(wěn)定性前提下，雖能收斂于原方之為耗散效應(yīng)；另一種是由于格式里引入了數(shù)值頻散項(xiàng)而使準(zhǔn)確解呈u+1=u?c(u?u?1)(c=aΔtΔx)相對原始方程，上述方程多出等式右端各項(xiàng)，它們是離散格式uxuxxutttutxxuxxx1a0aΔx?226002?2aΔt?20240aΔtΔx04aΔt220aΔt40?42aΔt002323023?423?41a00(aΔx2000(aΔx26)*因此，對一維線性對流方程，在a>0對微分方程L(φ(xj,tn))=0，MPDE一般可寫成（3-92）這里省略詳細(xì)推演，只列出相關(guān)結(jié)論：修正方程(3-92)右端的2M+1≠0，則在(?1)Mμ2M+1<0時，格式為由上結(jié)論可以看出，如耗散主項(xiàng)是2,6,10,....階導(dǎo)數(shù)時，則其系數(shù)為正時，格式為耗散，否則逆耗散；如耗散主項(xiàng)是4,8,12,....階導(dǎo)數(shù)時，則其系數(shù)為正時，格式為逆耗散，3.5.6離散格式的守恒性[1,13]對離散方程定義域的有限空間作求和運(yùn)算，所得表達(dá)式滿足該區(qū)足物理量守恒。對界面物理量的這種連續(xù)性要求實(shí)際上是物理量局部將其在求解域某個區(qū)段[J1,J2]節(jié)點(diǎn)上求和，有