計(jì)算方法 偏微分方程數(shù)值解

1、Email: JPhone: 02985583997第五章 偏微分方程數(shù)值解 Numerical Methods for Partial Differential EquationsNumerical Methods for Partial Differential Equations&5.1 5.1 偏偏微分方程簡介微分方程簡介 &5.2 5.2 離散化公式離散化公式 &5.3 幾種常見偏微分方程的離散化計(jì)算幾種常見偏微分方程的離散化計(jì)算&5.4吸附床傳熱傳質(zhì)模型中偏微分方程求解吸附床傳熱傳質(zhì)模型中偏微分方程求解Email: JPhone: 029855839

2、97本章要求本章要求 教學(xué)目的教學(xué)目的 講解： 偏微分方程離散格式及求解的一般過程 教學(xué)要求教學(xué)要求E熟記 一階及二階偏微分方程的離散格式；C 精通 用EXCEL迭代對偏微分方程求解； 探索 用兩數(shù)組交替更新的辦法進(jìn)行編程求解；F 延伸 對化學(xué)反應(yīng)工程中物理場的模擬進(jìn)行嘗試。 教學(xué)重點(diǎn)教學(xué)重點(diǎn)各種偏微分方程的離散與求解EXCEL 循環(huán)迭代問題 教學(xué)難點(diǎn)教學(xué)難點(diǎn) 特殊邊界條件的引入與應(yīng)用Email: JPhone: 029855839975. 1 偏微分方程簡介&偏微分方程如果一個(gè)微分方程中出現(xiàn)多元函數(shù)的偏導(dǎo)數(shù)，或者說如果未知函數(shù)和幾個(gè)變量有關(guān)，而且方程中出現(xiàn)未知函數(shù)對幾個(gè)變量的導(dǎo)數(shù)，那

3、么這種微分方程就是偏微分方程。在化工或化學(xué)動態(tài)模擬方程中，常常有一個(gè)自變量是時(shí)間，其它的自變量為空間位置。如果只考慮一維空間，則只有兩個(gè)自變量；如果考慮兩維空間，則有3個(gè)自變量。 許多化工過程均是通過對偏微分方程的求解進(jìn)行工藝參數(shù)的確定或數(shù)值模擬。Email: JPhone: 029855839975.1 偏微分方程簡介&偏微分方程的分類22222( )( )( )( )( )( )( )0uuuuuabcdefugx yxyxy 線性微分方程 Linear partial differencial equation擬線性微分方程 Quasilinear partial differe

4、ncial equation非線性微分方程 Nonlinear partial differencial equation,x y 1, ,/,nx yuxy , ,/,nx yuxy Email: JPhone: 029855839975.1 偏微分方程簡介&數(shù)學(xué)上的分類：橢圓方程 Elliptic拋物線方程 Parabolic雙曲線方程 Hyperbolic&物理實(shí)際問題的歸類：波動方程(雙曲型）一維弦振動模型：熱傳導(dǎo)方程（拋物線型）一維線性熱傳導(dǎo)方程拉普拉斯方程（橢圓型）穩(wěn)態(tài)靜電場或穩(wěn)態(tài)溫度分布場）240bac240bac240bac22222uutx 22uutx 22

5、220uuxyEmail: JPhone: 029855839975.1 微分方程的求解思路&求微分方程數(shù)值解的一般步驟：Step1區(qū)域剖分：首先按一定規(guī)則將整個(gè)定義域分成若干小塊Step2微分方程離散：構(gòu)造離散點(diǎn)或片的函數(shù)值遞推公式或方程Step3初始、邊界條件離散：根據(jù)遞推公式，將初值或邊界值離散化，補(bǔ)充方程，啟動遞推運(yùn)算Step4 數(shù)值解計(jì)算：求解離散系統(tǒng)問題&微分方程的定解問題 離散系統(tǒng)的求解問題Email: JPhone: 029855839975.2 離散化公式&將自變量在時(shí)間和空間上以一定的間隔進(jìn)行離散化，則應(yīng)變量就變成了這些離散變量的函數(shù)。&一階

6、偏導(dǎo)的離散化公式一般采用歐拉公式表示&有時(shí)為了保證系統(tǒng)和穩(wěn)定性，對時(shí)間的差分往往采用向后公式, ,( , , , )ni j kt n t x i x yj y z k zuu t x y z 1, , ,1, , ,1, , ,1, ,nni j ki j kt n t x i x yj y z k znnij ki j kt n t x i x yj y z k znni jki j kt n t x i x yj y z k znni j ki j kt n t x i x yj y z k zuuuttuuuxxuuuyyuuuzx 1, , ,(1),nni j ki j k

7、tnt x i x yj y z k zuuutt Email: JPhone: 029855839975.2 離散化公式P對于二階偏導(dǎo)，我們可以通過對泰勒展開式處理技術(shù)得到下面離散化計(jì)算公式：112, , , ,22,21, , ,1, ,22,2,1, ,1, ,22,222()2()2()nnni j ki j ki j kt n t x i x yj y z k znnnij ki j kij kt n t x i x yj y z k znnni jki j kij kt n t x i x yj y z k zt n tuuuuttuuuuxxuuuuyyuz , ,1, , ,1

8、2,2()nnni j ki j ki j kx i x yj y z k zuuuz Email: JPhone: 029855839975.2 離散化公式推導(dǎo)&將uk+1在uk處按二階泰勒式展開：&將uk-1在uk處按二階泰勒式展開：&二式相加得：22312()2!kkkkuuhuuhO hxx 22312()2!kkkkuuhuuhO hxx 211222()kkkuuuuxx Email: JPhone: 029855839975.3幾種常見偏微分方程的離散化計(jì)算&1 1、 波動方程波動方程 其中： 為初值條件 為邊值條件 當(dāng)該波動方程只提供初值條件時(shí)，

9、稱此方程為波動方程的初值問題，二者均提供時(shí)稱為波動方程的混合問題。2222200012( , )( ),( )( ),( )ttxx luuaf x ttxuuxxtut ut 00012( ),( )( ),( )ttxx tuuxxtutut Email: JPhone: 029855839975.3.1 波動方程求解對于初值問題，是已知t=0時(shí)，u與 依賴于x的函數(shù)形式，求解不同位置，不同時(shí)刻的u值。而 u是定義在 的二元函數(shù)，即上半平面的函數(shù)。 對于混合問題除初值外，還有邊值。是已知初值及x=0及x=l 時(shí)u依賴于t的函數(shù)，求解不同位置x，不同時(shí)刻的u值。此時(shí)u是定義在 的帶形區(qū)域上的

10、二元函數(shù)。ut 0,tx xt 0a)初值問題初值問題tx0lb)混合問題混合問題Email: JPhone: 029855839975.3.1 波動方程求解2222200012( , )( ),( )( ),( )ttxx luuaf x ttxuuxxtut ut 1121122(i1,2,m-1)22( , ) ()()(1,2,)nnnnnniiiiiiuuuuuuaf x ttxn 22212221211222()()()(22)()( , )()()()nnnnniiiiitttuauauauuxf x txxx 100012(),() (i1,2,m)(),() (n1,2,)i

11、iinnmuuujxixtuntunt 方程離散化整理可得：邊界條件初始條件離散化xxinniu1niu 1niu Email: JPhone: 029855839975.3.1 波動方程求解&例例5.1:5.1: 用數(shù)值法求解下面偏微分方程。 0012()3150,30,2501,0WnWjxttTtxTtttxx 110.01,0.1nnjjnnjjttttttxxx 110.020.680.3nnnjWjjtTtt 112()3nnnnjjjjjWntttttTx 此微分方程，是在不考慮流體本身熱傳導(dǎo)時(shí)的套管傳熱微分方程.由計(jì)算結(jié)果可知，當(dāng)計(jì)算的時(shí)間序列進(jìn)行到7272時(shí)，傳熱過程

12、已達(dá)到穩(wěn)態(tài)，各點(diǎn)上的溫度已不隨時(shí)間的增加而改變。如果改變套管長度或傳熱系數(shù)，則達(dá)到穩(wěn)態(tài)的時(shí)間亦會改變。Email: JPhone: 029855839975.3.2 一維流動熱傳導(dǎo)方程 與波動方程的情形類似，用差商近似代替偏商，可以得到一維流動傳熱傳導(dǎo)方程的混合問題的差分方程，以其解作為流動傳熱傳導(dǎo)方程的近似解。222001( , ) 00 00 0 tx lxuuuabf u t (xl , t)txxu(x), (xl)u(t)xu (t) 0 (t) 12111201012(,)()() (i1 2)0 (0,1,2,)(nnnnnnniiiiiiiinnmmnuuuuuuuabf i

13、x n ttxxui x, ,muunxun ) (0,1,2,)tn 1, , ,21, , ,1, ,22,2,1, ,1, ,22,22,2()2()nni j ki j kt n t x i x yj y z k znnnij ki j kij kt n t x i x yj y z k znnni jki j kij kt n t x i x yj y z k zt n t x i x yj y z kuuuttuuuuxxuuuuyyuz , ,1, , ,122()nnni j ki j ki j kzuuuz 2 2、一維流動熱傳導(dǎo)方程的混合問題、一維流動熱傳導(dǎo)方程的混合問題E

14、mail: JPhone: 02985583997將上式進(jìn)行處理得到： 該式是顯式格式。只要保證式中各項(xiàng)系數(shù)大于零，一般情況下是穩(wěn)定的，可以獲得穩(wěn)定的解。 分析上式可以發(fā)現(xiàn)，當(dāng)為了提高數(shù)值精度取適當(dāng)小的x 時(shí)，最有可能小于零的系數(shù)是 uin的系數(shù)，若要保證此項(xiàng)系數(shù)大于零，此時(shí)t必須相應(yīng)地更小,會導(dǎo)致計(jì)算量將大大增加，這是顯式格式的缺點(diǎn)，為了克服此缺點(diǎn)，下面提出一種隱式格式： 偏微分方程在 點(diǎn)上進(jìn)行離散化，且對時(shí)間的偏微分采用向后歐拉公式得到原偏微分方程的離散化公式：111111211122(,(1)()nnnnnnniiiiiiiuuuuuuuabf i x nttxx 122211222(,

15、)()(12)()()()nnnniiiitttttutf i x n tabuabuauxxxxx (,(1)i x nt 5.3.2 一維流動熱傳導(dǎo)方程Email: JPhone: 02985583997 從圖5-3中可見要由初值及邊界條件一排一排推上去是不行的，需解線性方程組,同時(shí)添上二邊界條件： 正好共有m+2個(gè)方程，同時(shí)有m+2個(gè)變量，就能解出n+1排上各點(diǎn)值。這樣，每解一個(gè)線性方程組，就可以往上推算一排點(diǎn)的u值，雖然引入了方程組的求解，有可能增加計(jì)算量，但由于隱式格式無條件穩(wěn)定，t的取法與x 無關(guān)，可以少計(jì)算許多排節(jié)點(diǎn)上的u 值，相應(yīng)于顯式格式來說，最終反而節(jié)省了計(jì)算量。 1110

16、11(1), nnnmmuntuu 5.3.2 一維流動熱傳導(dǎo)方程Email: JPhone: 02985583997&例5.2 考慮縱向?qū)岬奶坠軗Q熱器內(nèi)管各點(diǎn)溫度分布微分方程： 解：首先根據(jù)前面的知識，將所求 的方程離散化：代入微分方程并化簡得：分析上式可知，如果知道了某一時(shí)刻的各點(diǎn)t，（，（j=0,1,2.10,11），），就可以求下一時(shí)刻的各點(diǎn)溫度值t(j=1,2.10),現(xiàn)在已經(jīng)知道了零時(shí)刻管內(nèi)各點(diǎn)的溫度分布及入口處在任何時(shí)刻的溫度，如想求下一時(shí)刻的溫度值，根據(jù)上面的離散化計(jì)算公式，還需知道在j=11處的溫度，這個(gè)溫度可利用給定的邊界條件離散化求得： 有了以上各式，上面的微分

17、方程就可以求解了。220012()0.0013150,30,3001,0WnWjxtttTtxxTtttxx 11211222()0.01,0.1nnjjnnjjnnnjjjttttttxxttttxxx 1110.020.0010.6780.301nnnnjWjjjtTttt 110nnjjxtttxx 5.3.2 一維流動熱傳導(dǎo)方程Email: JPhone: 029855839975.3.3 穩(wěn)態(tài)導(dǎo)熱/擴(kuò)散方程&3 3、穩(wěn)態(tài)導(dǎo)熱、穩(wěn)態(tài)導(dǎo)熱/ /擴(kuò)散方程擴(kuò)散方程 在化工導(dǎo)熱及擴(kuò)散過程中，沒有物流的流動，僅靠導(dǎo)熱及擴(kuò)散進(jìn)行熱量及質(zhì)量的傳遞。如果此時(shí)系統(tǒng)達(dá)到穩(wěn)定狀態(tài)，也就是說系統(tǒng)中每一

18、個(gè)控制單元的各項(xiàng)性質(zhì)如溫度、濃度等不再隨時(shí)間的改變而改變，系統(tǒng)中的各種性質(zhì)只與其所處的位置有關(guān)，利用化工知識，我們可以得到下面二維、三維的穩(wěn)態(tài)導(dǎo)熱或擴(kuò)散偏微分方程： 二維： 三維：二維的穩(wěn)態(tài)導(dǎo)熱或擴(kuò)散偏微分方程又稱調(diào)和方程。22220uuxy2222220uuuxyzuu un ()uun 常見有三種邊界條件： 第一類邊界條件： 第二類邊界條件： 第三類邊界條件：Email: JPhone: 02985583997& 離散化公式：& 取 ，經(jīng)化簡得：& 外節(jié)點(diǎn)（邊界節(jié)點(diǎn)）和內(nèi)節(jié)點(diǎn)& 求解方法劃分網(wǎng)格建立節(jié)點(diǎn)離散方程迭代求解（或解稀疏方程組）1,1,1,12222

19、0()()iji jiji ji ji juuuuuuxy xy ,1,11,11()4i jiji jiji juuuuu xy求解區(qū)域N 節(jié)點(diǎn)邊界五點(diǎn)格式示意圖5.3.3穩(wěn)態(tài)導(dǎo)熱/擴(kuò)散方程求解Email: JPhone: 029855839975.3.3穩(wěn)態(tài)導(dǎo)熱/擴(kuò)散方程求解 常用的3種迭代格式： (1)同步迭代： (2)異步迭代： (3)超松弛迭代：當(dāng)計(jì)算范圍R 為 矩陣區(qū)域，x方向m等分，y方向n等分，最佳松弛因子為：由數(shù)學(xué)知識可知，用這些迭代法求解上面的偏微分方程均收斂。 (1)( )( )( )( ),1,11,11()4kkkkki jiji jiji juuuuu (1)( )

20、( )(1)(1),1,11,11()4kkkkki jiji jiji juuuuu ( )( )(1)(1)1,11,1(1)( ),1()4(1)kkkkiji jiji jkki ji juuuuuuwuw u (,)axb cyd 22coscos112wmn 緊湊迭代Email: JPhone: 029855839975.3.3 穩(wěn)態(tài)導(dǎo)熱/擴(kuò)散方程求解例例5.35.3 ：處于傳熱平衡狀態(tài)的某保溫，假設(shè)其形狀為長方體，在x，y兩個(gè)方向上存在熱傳導(dǎo)，且導(dǎo)熱系數(shù)相等，已知邊界溫度分布如下圖所示： 解：取某一微元進(jìn)行能量衡算，由于已達(dá)傳熱平衡狀態(tài)，故可得： 傳導(dǎo)入熱量-傳導(dǎo)出熱量=0 1x

21、y10(1,1)24( , ) 30 27001,0t x yyyx 1( , ) 30 27001,0t x yxxy 2( , )30001,0tx yyx3( , ) 30001,1t x yxy ()yttyyx zy ()xttxxy zx ytx zy xty zx xyz22220 xyttxy Email: JPhone: 029855839975.3.3 穩(wěn)態(tài)導(dǎo)熱/擴(kuò)散方程求解&Microsoft Excel 迭代計(jì)算公式中的循環(huán)引用在“工具工具”菜單上，單擊“選項(xiàng)”，再單擊“重新計(jì)算”選項(xiàng)卡。 選中“迭代”復(fù)選框。 若要設(shè)置 Microsoft Excel 進(jìn)行重新

22、計(jì)算的最大次數(shù)，請?jiān)凇白疃嗟螖?shù)”框中鍵入迭代次數(shù)。迭代次數(shù)越高，Excel 用于計(jì)算工作表的時(shí)間越多。 若要設(shè)置兩次迭代結(jié)果之間可以接受的最大誤差，請?jiān)凇白畲笳`差”框中鍵入所需的數(shù)值。數(shù)值越小，結(jié)果越精確，Excel 用于計(jì)算工作表的時(shí)間也越多。Email: JPhone: 029855839975.4 吸附床傳熱傳質(zhì)模型中偏微分方程求解實(shí)例&5.4.1 5.4.1 基本設(shè)定及假設(shè)基本設(shè)定及假設(shè) 1.吸附器結(jié)構(gòu)參數(shù)的設(shè)定 上圖所示的是套筒式吸附器，該吸附器的有效長度為L，其有效內(nèi)徑為D,環(huán)隙寬度為，吸附器壁厚為b。導(dǎo)熱流體通過環(huán)隙將熱量傳入或傳出吸附器，吸附質(zhì)通過吸附器上端的小管進(jìn)

23、入或離開吸附器。吸附器結(jié)構(gòu)示意圖吸附器結(jié)構(gòu)示意圖 DL熱流體Email: JPhone: 029855839975.4.1 基本設(shè)定及假設(shè)2.吸附床外流體傳熱的一些基本假設(shè)： 1). 忽略流體在環(huán)隙寬度上的溫度梯度； 2). 忽略熱損失； 3). 忽略吸附器壁厚b上的溫度梯度，用集中參數(shù)法求取吸附器壁面溫度。.吸附床內(nèi)傳熱傳質(zhì)的一些基本假設(shè)： 1). 吸附床內(nèi)的吸附質(zhì)氣體處于氣滯狀態(tài)； 2). 忽略蒸發(fā)器、冷凝器和吸附床之間的壓力差； 3). 吸附床內(nèi)各計(jì)算微元內(nèi)達(dá)到吸附平衡。吸附量可利用回歸方程計(jì)算； 4. 吸附熱利用微分吸附熱，隨吸附量和吸附溫度的改變而改變；比熱采用有效比熱，亦隨溫度改變

24、，但在計(jì)算微元內(nèi)，可認(rèn)為是常數(shù)；5. 床層活性炭導(dǎo)熱系數(shù)采用當(dāng)量導(dǎo)熱系數(shù)，可由實(shí)驗(yàn)測量得到。 Email: JPhone: 029855839975.4.2流體傳熱模型的建立在軸方向上取一環(huán)隙微元，作能量分析如下：1.流體通過流動流入環(huán)隙微元的能量為 2.流體通過流動流出環(huán)隙微元的能量3.流體熱傳導(dǎo)在x 處的熱量導(dǎo)入 7 總能量平衡方程fffTTTxxxxx 流體傳熱微元模型流體傳熱微元模型 infffpffquSCT ()fou tfffpffTquSCTxx outqinqxqfxffTqSx 其中： f 流體的密度 uf 環(huán)隙的流體速度， Sf 環(huán)隙的橫截面積，Cpf流體的比熱。4. 流

25、體熱傳導(dǎo)在x+ x處的熱量導(dǎo)入 5. 微元體傳遞給吸附床的熱量 qt 6. 微元體內(nèi)的能量變化率 為流體的橫截面積xxq(/)()fFxxffTTxqxSxx ()tffwqhxDTT cq()fpfffcCSx Tqt 221()ffffffwffffTuThDTTTtxSx fffpfC fSEmail: JPhone: 029855839975.4.3 吸附床內(nèi)吸附劑傳熱傳質(zhì)模型的建立&吸附床內(nèi)發(fā)生著熱量和質(zhì)量的傳遞，但質(zhì)量的傳遞是建立在熱量傳遞基礎(chǔ)上的，故只要建立熱量傳遞方程，就可以根據(jù)平衡吸附量方程求出各處的吸附量。吸附床內(nèi)的熱量傳遞主要以熱傳導(dǎo)為主，既有經(jīng)向的熱傳導(dǎo)，也有軸

26、向的熱傳導(dǎo)，為了便于建模分析，選取如圖所示的吸附床微元體，進(jìn)行衡算： 2bxaTqrrx x+ x x xrr+ r吸附床內(nèi)傳熱傳質(zhì)微元體1.軸向?qū)霟崃?:2.軸向?qū)С鰺崃?.徑向?qū)霟崃?4.徑向?qū)С鰺崃?5.微元體內(nèi)的能量變化率 其中 為吸附床層內(nèi)的有效比熱。6.總能量平衡方程 222()bbxxaTTqrrxxx 2braTqrxr 222 ()()bbdraTTqrrxrrr 2bceffaTqrrx Ct ()effpapbPmCCmCHT 222211bbbbeffTTTTtrrxraeffaeffC Email: JPhone: 029855839975 .4 .4 吸附器內(nèi)/

27、外無量綱化方程吸附器內(nèi)/外無量綱化方程2,021()()fwwwfwwbswswshDThDTTTTTxSSx ,sswpsssC 221()ffffffwffffTuThDTTTtxSx fffpfC 222211bbbbeffTTTTtrrxraeffaeffC *2000000,ffwbininintxrxrtLRLTTTTTTTTTTTT 無量綱化處理無量綱化處理Email: JPhone: 029855839975 .4 .4 吸附器內(nèi)吸附器內(nèi)/ /外無量綱化方程外無量綱化方程2*2()PeKA Bitxx 20*2()sswsfswsfDrBiBiBiBitx 22222,fffffbbefffbbffffwsswsfsswswuC LhDLSDrPeBiKASShDLhDLLDrBiBiCSSR 22*2*2*bDrCCtxrrr 整理可得：其中：Email: JPhone: 029855839975 .4 .4 吸附器內(nèi)吸附器內(nèi)/ /外無量綱化方程外無量綱化方程*(,0)0(,0)0(,0)0 xtxtxrt *0*1(0,)1 0(1,)(1 ,) 0 xxxtxx