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

第一頁，共三十二頁，編輯于2023年，星期三本章要求教學(xué)目的

講解：偏微分方程離散格式及求解的一般過程教學(xué)要求熟記一階及二階偏微分方程的離散格式；精通用EXCEL迭代對偏微分方程求解；探索用兩數(shù)組交替更新的辦法進(jìn)行編程求解；延伸對化學(xué)反應(yīng)工程中物理場的模擬進(jìn)行嘗試。教學(xué)重點各種偏微分方程的離散與求解EXCEL循環(huán)迭代問題教學(xué)難點特殊邊界條件的引入與應(yīng)用第二頁，共三十二頁，編輯于2023年，星期三5.1偏微分方程簡介偏微分方程如果一個微分方程中出現(xiàn)多元函數(shù)的偏導(dǎo)數(shù)，或者說如果未知函數(shù)和幾個變量有關(guān)，而且方程中出現(xiàn)未知函數(shù)對幾個變量的導(dǎo)數(shù)，那么這種微分方程就是偏微分方程。在化工或化學(xué)動態(tài)模擬方程中，常常有一個自變量是時間，其它的自變量為空間位置。如果只考慮一維空間，則只有兩個自變量；如果考慮兩維空間，則有3個自變量。許多化工過程均是通過對偏微分方程的求解進(jìn)行工藝參數(shù)的確定或數(shù)值模擬。第三頁，共三十二頁，編輯于2023年，星期三5.1偏微分方程簡介偏微分方程的分類線性微分方程Linearpartialdifferencialequation擬線性微分方程Quasilinearpartialdifferencialequation非線性微分方程Nonlinearpartialdifferencialequation第四頁，共三十二頁，編輯于2023年，星期三5.1偏微分方程簡介數(shù)學(xué)上的分類：橢圓方程Elliptic拋物線方程Parabolic雙曲線方程Hyperbolic物理實際問題的歸類：波動方程(雙曲型）一維弦振動模型：熱傳導(dǎo)方程（拋物線型）一維線性熱傳導(dǎo)方程拉普拉斯方程（橢圓型）穩(wěn)態(tài)靜電場或穩(wěn)態(tài)溫度分布場）第五頁，共三十二頁，編輯于2023年，星期三5.1微分方程的求解思路求微分方程數(shù)值解的一般步驟：Step1區(qū)域剖分：首先按一定規(guī)則將整個定義域分成若干小塊Step2微分方程離散：構(gòu)造離散點或片的函數(shù)值遞推公式或方程Step3初始、邊界條件離散：根據(jù)遞推公式，將初值或邊界值離散化，補充方程，啟動遞推運算Step4數(shù)值解計算：求解離散系統(tǒng)問題微分方程的定解問題離散系統(tǒng)的求解問題第六頁，共三十二頁，編輯于2023年，星期三5.2離散化公式將自變量在時間和空間上以一定的間隔進(jìn)行離散化，則應(yīng)變量就變成了這些離散變量的函數(shù)。一階偏導(dǎo)的離散化公式一般采用歐拉公式表示有時為了保證系統(tǒng)和穩(wěn)定性，對時間的差分往往采用向后公式第七頁，共三十二頁，編輯于2023年，星期三5.2離散化公式對于二階偏導(dǎo)，我們可以通過對泰勒展開式處理技術(shù)得到下面離散化計算公式：第八頁，共三十二頁，編輯于2023年，星期三5.2離散化公式推導(dǎo)將uk+1在uk處按二階泰勒式展開：將uk-1在uk處按二階泰勒式展開：二式相加得：第九頁，共三十二頁，編輯于2023年，星期三5.3幾種常見偏微分方程的離散化計算1、

波動方程

其中：為初值條件為邊值條件

當(dāng)該波動方程只提供初值條件時，稱此方程為波動方程的初值問題，二者均提供時稱為波動方程的混合問題。第十頁，共三十二頁，編輯于2023年，星期三5.3.1波動方程求解對于初值問題，是已知t=0時，u與依賴于x的函數(shù)形式，求解不同位置，不同時刻的u值。而

u是定義在的二元函數(shù)，即上半平面的函數(shù)。

對于混合問題除初值外，還有邊值。是已知初值及x=0及x=l時u依賴于t的函數(shù)，求解不同位置x，不同時刻的u值。此時u是定義在的帶形區(qū)域上的二元函數(shù)。xt0a)初值問題tx0lb)混合問題第十一頁，共三十二頁，編輯于2023年，星期三5.3.1波動方程求解方程離散化整理可得：邊界條件初始條件離散化xxiτn第十二頁，共三十二頁，編輯于2023年，星期三5.3.1波動方程求解例5.1:

用數(shù)值法求解下面偏微分方程。

此微分方程，是在不考慮流體本身熱傳導(dǎo)時的套管傳熱微分方程.由計算結(jié)果可知，當(dāng)計算的時間序列進(jìn)行到72時，傳熱過程已達(dá)到穩(wěn)態(tài)，各點上的溫度已不隨時間的增加而改變。如果改變套管長度或傳熱系數(shù)，則達(dá)到穩(wěn)態(tài)的時間亦會改變。EXCEL第十三頁，共三十二頁，編輯于2023年，星期三5.3.2一維流動熱傳導(dǎo)方程與波動方程的情形類似，用差商近似代替偏商，可以得到一維流動傳熱傳導(dǎo)方程的混合問題的差分方程，以其解作為流動傳熱傳導(dǎo)方程的近似解。2、一維流動熱傳導(dǎo)方程的混合問題離散化第十四頁，共三十二頁，編輯于2023年，星期三將上式進(jìn)行處理得到：該式是顯式格式。只要保證式中各項系數(shù)大于零，一般情況下是穩(wěn)定的，可以獲得穩(wěn)定的解。分析上式可以發(fā)現(xiàn)，當(dāng)為了提高數(shù)值精度取適當(dāng)小的Δx

時，最有可能小于零的系數(shù)是uin的系數(shù)，若要保證此項系數(shù)大于零，此時Δt必須相應(yīng)地更小,會導(dǎo)致計算量將大大增加，這是顯式格式的缺點，為了克服此缺點，下面提出一種隱式格式：偏微分方程在點上進(jìn)行離散化，且對時間的偏微分采用向后歐拉公式得到原偏微分方程的離散化公式：5.3.2一維流動熱傳導(dǎo)方程第十五頁，共三十二頁，編輯于2023年，星期三從圖5-3中可見要由初值及邊界條件一排一排推上去是不行的，需解線性方程組,同時添上二邊界條件：

正好共有m+2個方程，同時有m+2個變量，就能解出n+1排上各點值。這樣，每解一個線性方程組，就可以往上推算一排點的u值，雖然引入了方程組的求解，有可能增加計算量，但由于隱式格式無條件穩(wěn)定，Δt的取法與Δx無關(guān)，可以少計算許多排節(jié)點上的u值，相應(yīng)于顯式格式來說，最終反而節(jié)省了計算量。

5.3.2一維流動熱傳導(dǎo)方程第十六頁，共三十二頁，編輯于2023年，星期三例5.2考慮縱向?qū)岬奶坠軗Q熱器內(nèi)管各點溫度分布微分方程：解：首先根據(jù)前面的知識，將所求的方程離散化：代入微分方程并化簡得：分析上式可知，如果知道了某一時刻的各點t，（j=0,1,2….10,11），就可以求下一時刻的各點溫度值t(j=1,2….10),現(xiàn)在已經(jīng)知道了零時刻管內(nèi)各點的溫度分布及入口處在任何時刻的溫度，如想求下一時刻的溫度值，根據(jù)上面的離散化計算公式，還需知道在j=11處的溫度，這個溫度可利用給定的邊界條件離散化求得：

有了以上各式，上面的微分方程就可以求解了。5.3.2一維流動熱傳導(dǎo)方程EXCEL第十七頁，共三十二頁，編輯于2023年，星期三5.3.3穩(wěn)態(tài)導(dǎo)熱/擴散方程3、穩(wěn)態(tài)導(dǎo)熱/擴散方程

在化工導(dǎo)熱及擴散過程中，沒有物流的流動，僅靠導(dǎo)熱及擴散進(jìn)行熱量及質(zhì)量的傳遞。如果此時系統(tǒng)達(dá)到穩(wěn)定狀態(tài)，也就是說系統(tǒng)中每一個控制單元的各項性質(zhì)如溫度、濃度等不再隨時間的改變而改變，系統(tǒng)中的各種性質(zhì)只與其所處的位置有關(guān)，利用化工知識，我們可以得到下面二維、三維的穩(wěn)態(tài)導(dǎo)熱或擴散偏微分方程：二維：三維：二維的穩(wěn)態(tài)導(dǎo)熱或擴散偏微分方程又稱調(diào)和方程。常見有三種邊界條件：第一類邊界條件：第二類邊界條件：第三類邊界條件：第十八頁，共三十二頁，編輯于2023年，星期三離散化公式：取，經(jīng)化簡得：外節(jié)點（邊界節(jié)點）和內(nèi)節(jié)點求解方法劃分網(wǎng)格建立節(jié)點離散方程迭代求解（或解稀疏方程組）xyΩ求解區(qū)域N節(jié)點邊界Γ五點格式示意圖5.3.3穩(wěn)態(tài)導(dǎo)熱/擴散方程求解第十九頁，共三十二頁，編輯于2023年，星期三5.3.3穩(wěn)態(tài)導(dǎo)熱/擴散方程求解常用的3種迭代格式：

(1)同步迭代：

(2)異步迭代：

(3)超松弛迭代：當(dāng)計算范圍R為矩陣區(qū)域，x方向m等分，y方向n等分，最佳松弛因子為：由數(shù)學(xué)知識可知，用這些迭代法求解上面的偏微分方程均收斂。

緊湊迭代第二十頁，共三十二頁，編輯于2023年，星期三5.3.3穩(wěn)態(tài)導(dǎo)熱/擴散方程求解例5.3

：處于傳熱平衡狀態(tài)的某保溫，假設(shè)其形狀為長方體，在x，y兩個方向上存在熱傳導(dǎo)，且導(dǎo)熱系數(shù)相等，已知邊界溫度分布如下圖所示：

解：取某一微元進(jìn)行能量衡算，由于已達(dá)傳熱平衡狀態(tài)，故可得：傳導(dǎo)入熱量-傳導(dǎo)出熱量=0

1xy10(1,1)溫度分布xyz第二十一頁，共三十二頁，編輯于2023年，星期三5.3.3穩(wěn)態(tài)導(dǎo)熱/擴散方程求解MicrosoftExcel迭代計算公式中的循環(huán)引用在“工具”菜單上，單擊“選項”，再單擊“重新計算”選項卡。選中“迭代”復(fù)選框。若要設(shè)置MicrosoftExcel進(jìn)行重新計算的最大次數(shù)，請在“最多迭代次數(shù)”框中鍵入迭代次數(shù)。迭代次數(shù)越高，Excel用于計算工作表的時間越多。若要設(shè)置兩次迭代結(jié)果之間可以接受的最大誤差，請在“最大誤差”框中鍵入所需的數(shù)值。數(shù)值越小，結(jié)果越精確，Excel用于計算工作表的時間也越多。第二十二頁，共三十二頁，編輯于2023年，星期三5.4吸附床傳熱傳質(zhì)模型中

偏微分方程求解實例5.4.1基本設(shè)定及假設(shè)

1.吸附器結(jié)構(gòu)參數(shù)的設(shè)定

上圖所示的是套筒式吸附器，該吸附器的有效長度為L，其有效內(nèi)徑為D,環(huán)隙寬度為δ，吸附器壁厚為δb。導(dǎo)熱流體通過環(huán)隙將熱量傳入或傳出吸附器，吸附質(zhì)通過吸附器上端的小管進(jìn)入或離開吸附器。吸附器結(jié)構(gòu)示意圖

δDL熱流體第二十三頁，共三十二頁，編輯于2023年，星期三5.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)各計算微元內(nèi)達(dá)到吸附平衡。吸附量可利用回歸方程計算；

4.

吸附熱利用微分吸附熱，隨吸附量和吸附溫度的改變而改變；比熱采用有效比熱，亦隨溫度改變，但在計算微元內(nèi)，可認(rèn)為是常數(shù)；5.床層活性炭導(dǎo)熱系數(shù)采用當(dāng)量導(dǎo)熱系數(shù)，可由實驗測量得到。

第二十四頁，共三十二頁，編輯于2023年，星期三5.4.2流體傳熱模型的建立在軸方向上取一環(huán)隙微元，作能量分析如下：1.流體通過流動流入環(huán)隙微元的能量為

2.流體通過流動流出環(huán)隙微元的能量3.流體熱傳導(dǎo)在x

處的熱量導(dǎo)入

7總能量平衡方程流體傳熱微元模型

其中：f——流體的密度uf——環(huán)隙的流體速度，

Sf——環(huán)隙的橫截面積，Cpf——流體的比熱。4.流體熱傳導(dǎo)在x+x處的熱量導(dǎo)入

5.微元體傳遞給吸附床的熱量

qt

6.微元體內(nèi)的能量變化率

為流體的橫截面積第二十五頁，共三十二頁，編輯于2023年，星期三5.4.3吸附床內(nèi)吸附劑傳熱傳質(zhì)模型的建立吸附床內(nèi)發(fā)生著熱量和質(zhì)量的傳遞，但質(zhì)量的傳遞是建立在熱量傳遞基礎(chǔ)上的，故只要建立熱量傳遞方程，就可以根據(jù)平衡吸附量方程求出各處的吸附量。吸附床內(nèi)的熱量傳遞主要以熱傳導(dǎo)為主，既有經(jīng)向的熱傳導(dǎo)，也有軸向的熱傳導(dǎo)，為了便于建模分析，選取如圖所示的吸附床微元體，進(jìn)行衡算：