第5章偏微分方程數(shù)值解

第5章

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

5.1問題的提出

5.2基本離散化公式

5.3幾種常見方程的離散化計(jì)算

5.4吸附床傳熱傳質(zhì)模型中偏微分方程求解實(shí)例

目錄5.1問題的提出包含有偏導(dǎo)數(shù)的微分方程稱為偏微分方程。從實(shí)際問題中歸納出來的常用偏微分方程可分為三大類：波動(dòng)方程、熱傳導(dǎo)方程和調(diào)和方程。對(duì)于它們特殊的定解條件，有一些解決的解析方法，而且要求方程是線性的、常系數(shù)的。但是在實(shí)際中碰到的問題卻往往要復(fù)雜得多，尤其在化工和化學(xué)模擬計(jì)算中，不僅偏微分方程的形式無一定標(biāo)準(zhǔn)，且邊界條件五花八門，方程中的系數(shù)隨工況改變而改變，想利用解析求解是不可能的。另一方面實(shí)際問題的要求不一定需要嚴(yán)格的精確解，只要求達(dá)到一定精度，所以就可借助于差分方法來求偏微分方程的數(shù)值解。在第4章里，我們介紹了一個(gè)套管式換熱器穩(wěn)態(tài)的傳熱問題。如果我們考慮一個(gè)動(dòng)態(tài)的傳熱過程，且不忽略縱向的熱傳導(dǎo)，就可以得到以下的偏微分方程：5.15.45.35.2總目錄5.1問題的提出上面方程中變量的含義如下：通過求解上面的偏微分方程，就可以得到傳熱管各點(diǎn)溫度隨時(shí)間的變化值，從而確定達(dá)到傳熱平衡所需的時(shí)間，為實(shí)驗(yàn)測(cè)量提供依據(jù)。想求解上述方程，就必須首先學(xué)會(huì)偏微分方程的求解方法，下面我們首先介紹如何對(duì)偏微分方程進(jìn)行離散化的工作，然后再對(duì)各類不同的偏微分方程進(jìn)行求解，我們一般只給出離散化的基本公式及計(jì)算方法，對(duì)離散化公式的具體推導(dǎo)工作一般不作詳細(xì)介紹，對(duì)這方面感興趣的讀者可自行參考有關(guān)數(shù)值計(jì)算的書籍。

5.15.45.35.2總目錄5.2基本離散化公式

在偏微分方程中，自變量都在兩個(gè)或兩個(gè)以上，應(yīng)變量隨兩個(gè)或兩個(gè)以上的自變量變化而變化。在化工或化學(xué)動(dòng)態(tài)模擬方程中，常常有一個(gè)自變量是時(shí)間，其它的自變量為空間位置。如果只考慮一維空間，則只有兩個(gè)自變量；如果考慮兩維空間，則有3個(gè)自變量。一般我們將自變量在時(shí)間和空間以一定的間隔進(jìn)行離散化，則應(yīng)變量就變成了這些離散變量的函數(shù)，以3維空間為例，我們將離散化的應(yīng)變量表示成，它所表示的真正含義如下:有了以上的定義，對(duì)于一階偏導(dǎo)我們可以利用第四章的歐拉公式直接得出向前歐拉公式：對(duì)于時(shí)間偏導(dǎo)而言，有時(shí)我們常常采用向后歐拉公式，時(shí)間的向后歐拉公式如下：5.15.45.35.2總目錄5.2基本離散化公式

這樣在以后的計(jì)算中，得到的是隱式的計(jì)算公式，需通過求解線性方程組才能求解。具體的計(jì)算過程我們?cè)谙旅鏁?huì)針對(duì)具體的偏微分方程進(jìn)行講解。對(duì)于二階偏導(dǎo)，我們可以通過對(duì)泰勒展開式處理技術(shù)得到下面離散化計(jì)算公式：有了以上的離散化公式，就可以進(jìn)行偏微分方程的數(shù)值求解工作。當(dāng)然，在具體求解時(shí)，還會(huì)碰到不同的問題，需要區(qū)別對(duì)待，同時(shí)在利用計(jì)算機(jī)編程計(jì)算時(shí)也會(huì)碰到困難，這些問題我們會(huì)通過具體的例子加以說明。

5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

1、

波動(dòng)方程其中：為初值條件為邊值條件

當(dāng)該波動(dòng)方程只提初值條件時(shí)，稱此方程為波動(dòng)方程的初值問題，二者均提時(shí)，稱為波動(dòng)方程的混合問題。對(duì)于初值問題，是已知t=0時(shí)，u與依賴于x的函數(shù)形式，求解不同位置，不同時(shí)刻的u值。而

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

對(duì)于混合問題除初值外，還有邊值。是已知初值及x=0及x=l時(shí)u依賴于t的函數(shù)，求解不同位置x，不同時(shí)刻的u值。此時(shí)u是定義在的帶形區(qū)域上的二元函數(shù)。如圖可以看出初值問題和混合問題的定義域。

5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

根據(jù)5.2節(jié)提供的公式，將上面波動(dòng)方程離散化，得到：(5-1)將式（5-1）進(jìn)行處理，把（n+1）時(shí)刻的變量留在右邊，其余放在左邊得到：

（5-2）

同時(shí)將邊界條件和初始條件也離散化，得到：

(5-3)xt0a)初值問題tx0lb)混合問題5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

這樣，由式（5-2），并結(jié)合式（5-3），就可以從n時(shí)刻的各點(diǎn)u值，計(jì)算得到下一時(shí)刻的u值，這樣層層遞推，就可以計(jì)算出任意時(shí)刻，任意位置的u值。而圖5.2則表明了這種層層遞推的計(jì)算過程，在圖5.2中*表示需求u值的點(diǎn)，○表示為了求x點(diǎn)的u

值必須已知u值的點(diǎn)。需要說明的是，在應(yīng)用式（5-2）進(jìn)行計(jì)算時(shí)，初值與邊值應(yīng)當(dāng)滿足相容性條件。由初值得到，由邊值得到

，,但在利用式（5-2）進(jìn)行第一輪計(jì)算時(shí)，若取n=0，則發(fā)現(xiàn)等式右邊出現(xiàn)了，這是一個(gè)無法計(jì)算的值。這時(shí)可以利用另一個(gè)初值條件算得，這樣，可在第一輪計(jì)算的時(shí)候，取n=1，計(jì)算得到，由，遞推得到，這樣就可由式（5-2）一排一排往上推，計(jì)算得到所有希望得到的u值。對(duì)于式（5-2）取n=0計(jì)算中碰到的，也可利用另一種方法進(jìn)行計(jì)算，解決的辦法是將另一個(gè)初值條件利用向后歐拉離散化算得，這樣利用式（5-2），取n=0就可以得到，取n=1，，和前一種處理方法一樣一排一排往上推，計(jì)算得到所有希望得到的u值。象這樣可以用已知點(diǎn)上函數(shù)值直接推出所有點(diǎn)上函數(shù)值的格式，稱為顯式格式。當(dāng)方程非齊次時(shí)，，式（5-1）可寫為5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

當(dāng)方程是初值問題時(shí)，邊界條件沒有了，由于在t=0時(shí)，u與值是已知的，若需要求某的值，只要按“波及原則”多算一些初值，即可推得，所圖5-2所示。為了保證差分方程的解在時(shí)收斂于原來波動(dòng)方程的解，要求式（5-2）中等式右邊的各項(xiàng)系數(shù)均大于0，即：化簡(jiǎn)得：

而且，可以證明，只要初始條件，邊界條件滿足一定的光滑性要求，且滿足收斂關(guān)系式時(shí)，差分格式是穩(wěn)定的。

圖5-2x5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

例5.1:

用數(shù)值法求解下面偏微分方程，并寫出VB程序。解：首先根據(jù)前面的知識(shí)，將所求的方程離散化，先假設(shè)以下各式：代入微分方程并化簡(jiǎn)得：（5-3）分析式（5-3）可知，如果知道了某一時(shí)刻的各點(diǎn)t，（j=0,1,2….100），就可以求下一時(shí)刻的各點(diǎn)溫度值。有了以上各式，上面的微分方程就可以求解了。其實(shí)這個(gè)微分方程，是在不考慮流體本身熱傳導(dǎo)時(shí)的套管傳熱微分方程。由計(jì)算結(jié)果可知，當(dāng)計(jì)算的時(shí)間序列進(jìn)行到７２時(shí)，傳熱過程已達(dá)到穩(wěn)態(tài)，各點(diǎn)上的溫度已不隨時(shí)間的增加而改變。如果改變套管長度或傳熱系數(shù)，則達(dá)到穩(wěn)態(tài)的時(shí)間亦會(huì)改變。

5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

2、一維流動(dòng)傳熱傳導(dǎo)方程的混合問題與波動(dòng)方程的情形類似，用差商近似代替偏商，可以得到差分方程，以其解作為流動(dòng)傳熱傳導(dǎo)方程的近似解。一維流動(dòng)傳熱傳導(dǎo)方程的混合問題：

上面的偏微分方程其實(shí)就是在5.1節(jié)中提出的偏微分方程，利用5.2節(jié)中的離散化公式進(jìn)行離散化，得到其離散化公式：

5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

將上式進(jìn)行處理得到：

（5-4）利用初始條件和邊界條件，可以得到零時(shí)刻各點(diǎn)的（i=0,1,2,…m）及

，這樣就可以利用公式（5-4）計(jì)算得到，依次類推，可以得到其它時(shí)刻的各點(diǎn)值，所以式（5-4）也是顯式格式。只要保證式（5-4）中各項(xiàng)系數(shù)大于零，一般情況下，式（5-4）的計(jì)算公式是穩(wěn)定的，可以獲得穩(wěn)定的解。分析式（5-4）可以發(fā)現(xiàn)，當(dāng)為了提高數(shù)值精度取適當(dāng)小的時(shí)，最有可能小于零的系數(shù)是的系數(shù)，若要保證此項(xiàng)系數(shù)大于零，此時(shí)必須相應(yīng)地更小，這樣，計(jì)算量將大大增加，這是顯式格式的缺點(diǎn)，為了克服此缺點(diǎn)，下面提出一種隱式格式。偏微分方程在點(diǎn)上進(jìn)行離散化，且對(duì)時(shí)間的偏微分采用向后歐拉公式得到原偏微分方程的離散化公式：

（i=1,2….m）

5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

從圖5-3中可見要由初值及邊界條件一排一排推上去是不行的，需解線性方程組,同時(shí)添上二邊界條件：

正好共有m+2個(gè)方程，同時(shí)有m+2個(gè)變量，就能解出n+1排上各點(diǎn)值。（至于線性方程組的求解方法我們?cè)诘?章中已作過介紹，請(qǐng)讀者自行參照第3章的內(nèi)容）。這樣，每解一個(gè)線性方程組，就可以往上推算一排點(diǎn)的u值，雖然引入了方程組的求解，有可能增加計(jì)算量，但由于隱式格式無條件穩(wěn)定，

的取法與無關(guān)，可以少計(jì)算許多排節(jié)點(diǎn)上的u值，相應(yīng)于顯式格式來說，最終反而節(jié)省了計(jì)算量。

圖5-3

5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

例5.2請(qǐng)計(jì)算考慮縱向?qū)岬奶坠軗Q熱器內(nèi)管各點(diǎn)溫度分布微分方程：解：首先根據(jù)前面的知識(shí)，將所求的方程離散化，先假設(shè)以下各式：代入微分方程并化簡(jiǎn)得：分析式（5-3）可知，如果知道了某一時(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è)溫度可利用給定的邊界條件離散化求得：有了以上各式，上面的微分方程就可以求解了。

5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

和前面不考慮熱傳導(dǎo)的情況比較，可以發(fā)現(xiàn)溫度有細(xì)微的變化，如果導(dǎo)熱系數(shù)足夠大，則溫度的變化會(huì)更大。如導(dǎo)熱項(xiàng)的系數(shù)為0.2時(shí)，其計(jì)算公式變?yōu)椋?/p>

由于導(dǎo)熱的緣故，已經(jīng)加熱的向前流動(dòng)的流體卻要向后方向進(jìn)行熱傳導(dǎo)，從而降低了總體傳熱效率，使在相同時(shí)刻、相同位置點(diǎn)的溫度比沒有熱傳導(dǎo)時(shí)要低。

前

后

5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

3、穩(wěn)態(tài)導(dǎo)熱/擴(kuò)散方程

在化工導(dǎo)熱及擴(kuò)散過程中，沒有物流的流動(dòng)，僅靠導(dǎo)熱及擴(kuò)散進(jìn)行熱量及質(zhì)量的傳遞。如果此時(shí)系統(tǒng)達(dá)到穩(wěn)定狀態(tài)，也就是說系統(tǒng)中每一個(gè)控制單元的各項(xiàng)性質(zhì)如溫度、濃度等不再隨時(shí)間的改變而改變，系統(tǒng)中的各種性質(zhì)只與其所處的位置有關(guān)，利用化工知識(shí)，我們可以得到下面二維、三維的穩(wěn)態(tài)導(dǎo)熱或擴(kuò)散偏微分方程：二維：三維：下面我們主要介紹二維的求解方法，二維的穩(wěn)態(tài)導(dǎo)熱或擴(kuò)散偏微分方程又稱調(diào)和方程，其方程示意圖見圖5-4所示。常見有三種邊界條件：第一類邊界條件：

第二類邊界條件：

第三類邊界條件：

圖5-45.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

在化工中碰到較多的是第一類邊界條件，下面以第一類邊界條件為例，說明其求解方法：

首先利用5.2中提供的離散化公式（不考慮u中的上標(biāo)變量n），可得下面離散化公式：

取，經(jīng)化簡(jiǎn)得：對(duì)于每一個(gè)邊界內(nèi)的離散點(diǎn)均可列出這樣的五點(diǎn)格式。若中有邊界點(diǎn)，用邊界值代入。若靠邊界很近，也可以看作邊界節(jié)點(diǎn)，從靠它最近的邊界點(diǎn)上的值來取代。由于此計(jì)算格式不存在時(shí)間上的遞推問題，它只是不同空間位置上變量的求解問題，而已知條件僅僅知道邊界上的值，這樣要求邊界內(nèi)點(diǎn)的值只能通過離散化的偏微分方程來求解，幸好有多少個(gè)內(nèi)節(jié)點(diǎn)就有多少個(gè)離散化的方程，構(gòu)成了一個(gè)未知數(shù)個(gè)數(shù)與方程個(gè)數(shù)相等的稀疏方程組，既可直接求解，也可迭代求解，一般用迭代法解比較好。

5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

下面介紹3種迭代格式：

(1)同步迭代：

(2)異步迭代：

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

5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

例5.3

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

試列出其傳熱微分方程，并求出各點(diǎn)的溫度分布（間隔以），并畫出溫度分布圖。

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

1xy10(1,1)5.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

根據(jù)上圖的導(dǎo)熱表達(dá)圖可得：化簡(jiǎn)得：由已知條件可知，，則有：

此偏微分方程和前面介紹的調(diào)和方程一致，可用五點(diǎn)格式同步迭代計(jì)算，其中需要計(jì)算的內(nèi)點(diǎn)共有16點(diǎn)。也可列出16個(gè)線性方程，組成方程組利用第3章介紹的方法進(jìn)行求解，本書介紹利用五點(diǎn)格式同步迭代計(jì)算的VB程序

[+（）]－[＋]=05.15.45.35.2總目錄5.3幾種常見方程的離散化計(jì)算

計(jì)算結(jié)果：

圖示：

x

y00.20.40.60.81.003040.873.2127.2202.83000.284104.53134.89178.09234.133000.4138158.45183.73216.13255.653000.6192207.54225.46247.06272.343000.8246254.26263.5274.3286.663001.03003003003003003005.15.45.35.2總目錄VB調(diào)用5.4吸附床傳熱傳質(zhì)模型中偏微分方程求解實(shí)例

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)入或離開吸附器。

圖5-5吸附器結(jié)構(gòu)示意圖

5.15.45.35.2總目錄5.4吸附床傳熱傳質(zhì)模型中偏微分方程求解實(shí)例

2.吸附床外流體傳熱的一些基本假設(shè)

１．忽略流體在環(huán)隙寬度δ上的溫度梯度；２．忽略熱損失；3.忽略吸附器壁厚δb上的溫度梯度，用集中參數(shù)法求取吸附器壁面溫度。

３.吸附床內(nèi)傳熱傳質(zhì)的一些基本假設(shè)

１.

吸附床內(nèi)的吸附質(zhì)氣體處于氣滯狀態(tài)；２.

忽略蒸發(fā)器、冷凝器和吸附床之間的壓力差；３.

吸附床內(nèi)各計(jì)算微元內(nèi)達(dá)到吸附平衡。吸附量可利用回歸方程計(jì)算；４.

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

5.15.45.35.2總目錄5.4吸附床傳熱傳質(zhì)模型中偏微分方程求解實(shí)例

5.4.2流體傳熱模型的建立

在軸方向上取一環(huán)隙微元，見圖５－７，作能量分析如下：１.流體通過流動(dòng)流入環(huán)隙微元的能量去

其中f為流體的密度uf

為環(huán)隙的流體速度，

Sf為環(huán)隙的橫截面積，Cpf

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

處的熱量導(dǎo)入

圖５－6流體傳熱微元模型

5.15.45.35.2總目錄5.4吸附床傳熱傳質(zhì)模型中偏微分方程求解實(shí)例

4.流體熱傳導(dǎo)在x+x處的熱量導(dǎo)入

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

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

7總能量平衡方程

其中,為流體的橫截面積。

5.15.45.35.2總目錄5.4吸附床傳熱傳質(zhì)模型中偏微分方程求解實(shí)例

5.4.3吸附床內(nèi)吸附劑傳熱傳質(zhì)模型的建立

吸附床內(nèi)發(fā)生著熱量和質(zhì)量的傳遞，但質(zhì)量的傳遞是建立在熱量傳遞基礎(chǔ)上的，故只要建立熱量傳遞方程，就可以根據(jù)平衡吸附量方程求出各處的吸附量。吸附床內(nèi)的熱量傳遞主要以熱傳導(dǎo)為主，既有經(jīng)向的熱傳導(dǎo)，也有軸向的熱傳導(dǎo)，為了便于建模分析，我們選取如圖5.8的吸附床微元體，具體分析如下：

x+x

r

r+r

xx

圖5.7吸附床內(nèi)傳熱傳質(zhì)微元體

1.軸向熱量導(dǎo)入:5.15.45.35.2總目錄5.4吸附床傳熱傳質(zhì)模型中偏微分方程求解實(shí)例

2.軸向熱量導(dǎo)出

3.經(jīng)向熱量導(dǎo)入

4.徑向熱量導(dǎo)出

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

其中,為吸附床層內(nèi)的有效比熱。6.總能量平衡方程

其中

5.15.45.35.2總目錄5.4吸附床傳熱傳質(zhì)模型中偏微分方程求解實(shí)例

5.4.4吸附器壁面溫度軸向分布方程

和前面的分析方法一樣，通過微元能量平衡方程并作適當(dāng)化簡(jiǎn)可得：

其中為吸附器壁面的橫截面積。

5.15.45.35.2總目錄5.4吸附床傳熱傳質(zhì)模型中偏微分方程求解實(shí)例

5.4.5吸附器內(nèi)/外無因子化方程根據(jù)前面推導(dǎo)已經(jīng)得到的方程，并對(duì)變量作以下無因子化處理：

通過以上的無因子化處理，可得吸附器內(nèi)、外無因子化傳熱傳質(zhì)方程如下：