偏微分方程式之求解

1、第六章 偏微分方程式之求解在化工的領(lǐng)域中,有不少程序之動(dòng)態(tài)是由以偏微分方程式(Partial differential equation ;PDE所描述的,例如熱與質(zhì)量在空間中的傳遞等。這些用以描述實(shí)際問(wèn)題的PDE ,除非具有某些特定的方程式型態(tài)及條件,否則甚難以手算的方式找出解析解。而在數(shù)值求解方面,最常被采用的方法為有限差分法(finitedifference何有限元素法(finite element。然對(duì)于某些不熟悉數(shù)值分析及程序編寫的化工人而言,欲充分了解以偏微分方程式所描述之系統(tǒng)動(dòng)態(tài)是相當(dāng)不容易的,更遑論進(jìn)一步的設(shè)計(jì)與分析了。值得慶幸的是,MATLAB 的環(huán)境中提供了一個(gè)求解PDE

2、問(wèn)題的工具箱,讓使用者得以利用簡(jiǎn)單的指令或圖形接口工具輸入欲解的PDE ,并求解。使得PDE 之?dāng)?shù)值解在彈指之間完成,使用者不在為數(shù)值法所苦惱,輕松掌握偏微分方程式系統(tǒng)的動(dòng)態(tài),并可進(jìn)一步進(jìn)行后續(xù)之設(shè)計(jì)工作。本章將以循漸進(jìn)的方式,介紹PDE 工具箱及其用法,并以數(shù)個(gè)典型的化工范例進(jìn)行 示范,期能使初學(xué)者很快熟悉PDE 工具箱,并使用它來(lái)設(shè)計(jì)與分析以偏微方方程式所描述的程序系統(tǒng)。6.1 偏微分方程式之分類偏微分方程式可根據(jù)其階數(shù)(order,線性或非線性型態(tài),以及邊界條件進(jìn)行分類。偏微分方程式是以偏微分項(xiàng)中之最高次偏微分來(lái)定義其階數(shù),例如:一階偏微分方程式:0=+y u x u 二階偏微分方程式:

3、032222=+ +y ux u y u x u 三階偏微分方程式:02233=+ y ux u y x u x u偏微分方程式亦可以其線性或非線性情況,區(qū)分為線性(linear,似線性(quasilinear,以及非線性三類。例如,以下之二階偏微分方程式(Constantinides and Mostoufi,19990(22222=+d x uc y x u b y u a可依系數(shù)(之情況,進(jìn)行如下表之歸類類別 情況線性 系數(shù)(為定值,或僅為(x,y函數(shù)似線性 系數(shù)(為依變數(shù)(dependent variableu 或其比方程式中之偏微分低階之偏微分項(xiàng)的函數(shù),如,(y u x u u y

4、x =非線性 系數(shù)(中,具有與原方程式之偏微分同階數(shù)之變數(shù),如,(22222y x u y u x u u y x =另外,對(duì)于線性二階偏微分方程式,可進(jìn)一步將其分類為橢圓型(elliptic,拋物線型(parabolic,以及雙曲線型(hyperbolic。具體上來(lái)說(shuō),此類偏微分方程式二階線性之一般式為022222=+g fu y ue x u d yu c y x u b x u a系數(shù)f e d c b a 和,是定值或?yàn)閡 的函數(shù)。若g=0,則上式為其次是偏微分方程式。式子( 之分類及代表性例子,請(qǐng)見下表方程式類別 判斷式 代表性范例橢圓型 042<-ac b Laplace 方

5、程式,02222=+yux u拋物線型 042=-ac b Poisson 方程式,(2222y x f yux u =+f u a u c (=+-熱傳導(dǎo)或擴(kuò)散方程式t uxu =22f u a u c tu d (=+-雙曲線型 042>-ac b 波動(dòng)方程式22222tux u = f u a u ctu d (22=+-注:二維系統(tǒng)之運(yùn)操作數(shù)之定義為j yi x +=為了能獲得偏微分方程式之解答,其起始條件和邊界條件可依其特性區(qū)分為三類?，F(xiàn)以一維之動(dòng)態(tài)熱傳遞方程式(拋物線型偏微分方程式22xTt T = 為例,進(jìn)一步說(shuō)明如何區(qū)分這些邊界條件及起始條件(Constantinides

6、 andMostoufi ,1999。(i 第一類:Dirichlet Condiction若依變量(T本身,在某個(gè)獨(dú)立變量值時(shí),被指定,則此條件稱為DirichletCondiction ,亦稱為essential 邊界條件。下圖為一典型的Dirichlet 條件示意圖T 0,1>=t T T 0,(>=t t f T 01x圖6.1 平板Dirichlet Condiction 示意圖由圖中很清楚的顯示,該平板之邊界條件為0,0,(>=t x at t f T 0,0,1>=t x at T T 0,0,0>=t x at T T此邊界條件依定義,即為Diri

7、chlet Condiction 。同時(shí),若再起始時(shí),各處之溫度分布可以位置之函數(shù)表示,即10,0,(=x t at t f T 此亦屬Dirichlet 型之邊界條件。(ii 第二類:Neumann conditionNeumann condition 系指依變量之變化率之邊界條件為定值,抑或獨(dú)立變量之函數(shù)之情況。例如 0,1,0=t x at xT或10,0,(=x t at t f xTNeumann 型邊界條件,亦稱為natural boundary condition 。 (iii 第三類:Robbins condition若依變量之變化率之邊界條件,為自身之函數(shù)(非獨(dú)立變量之函數(shù)時(shí)

8、,被稱為Robbins condition 。例如,0,0,(=-=t x atT T h xTkf上式之邊界條件,當(dāng)發(fā)生在固液相間之傳遞上,亦即熱流通量(heat flux正比于固液兩端之溫差,其示意圖如下:1x固體平板Tt液相,fT 0,(>-=t T T h xTkf(iv Cauchy conditionCauchy condition 系指問(wèn)題中同時(shí)存在Dirichlet 和Neumann 邊界條件。例如下圖01x絕緣(0t f T t =>1,00=>x at xTt即0,(>=t t f T “Dirichlet condition” 0,1,0>=

9、t x xT“Neumann condition”6.2 MATLAB PDE 工具箱MATLAB 提供了一個(gè)指令pdepe ,用以解以下的PDE 方程式,(,(,(xu u t x s x u u t x f x x x t u x u u t x c m m +=- 其中時(shí)間介于f t t t 0之間,而位置x 則介于b a 有限區(qū)域之間。m 值表示問(wèn)題之對(duì)稱性,其可為0,1或2,以表示平板(slab,圓柱(cylindrical或球體(spherical之對(duì)秤情況。因而,如果如果m >0,則a 必等于b ,也就是說(shuō)其具有圓柱或球體之對(duì)稱關(guān)系。同時(shí),式中,(x u u t x f 一

10、項(xiàng)為流通量(flux,而,(x u u t x s 為來(lái)源(source項(xiàng)。,(x u u t x c 為偏微分方程式之對(duì)角線系數(shù)矩陣。若某一對(duì)角線元素為0,則表示該相應(yīng)偏微分方程式為橢圓型偏微分方程式,若為正值(不為0,則為拋物線型偏微分方程式。請(qǐng)注意c 之對(duì)角線元素必不全為0。偏微分方程式之起始值可表示為(,(00x v t x u = 而邊界條件之形式為0,(,(,(=+x u u t x f t x g u t x p 其中x 為兩端點(diǎn)位置,即a 或b 。用以解含上述起始值及邊界值條件的偏微分方程式( 之指令pdepe 的用法如下:,(o p t i o n s t s p a n x

11、 m e s h b c f u n i c f u n p d e p e m p d e p e s o l = 其中m 為問(wèn)題之對(duì)稱參數(shù)xmesh 為獨(dú)立變量x 之網(wǎng)取點(diǎn)(mesh位置向量,即N x x x x xmesh 21=,其中a x =0(起點(diǎn),b x N =(終點(diǎn)。tspan 為獨(dú)立變量t (時(shí)間之向量,即M t t t tspan 10=,其中0t 為起始時(shí)間,M t 為終點(diǎn)時(shí)間。pdefun 為使用者提供之pde 函數(shù)文件。其檔案之格式如下:,(,d u d x u t x p d e f u n s f c = 亦即,使用者僅需提供偏微分方程式中之系數(shù)向量。c ,f 和

12、s 均為行(column向量,而向量c 即為矩陣c之對(duì)角線元素。i c f u n 提供解u 之起始值,其格式為(x icfun u =值得注意的是u 為行向量。b c f u n使用者提供之邊界條件檔,格式如下:(t ur xr ul xl bcfun ql pr ql pl ,=其中ul 和ur 分別表示左邊界(b xl =和右邊界(a xr =之u 的近似解。輸出變量中,pl 和ql 分別表示左邊界p 和q 之行向量,而pr 和qr 則為右邊界p 和q 之行向量。s o l 為解答輸出。sol 為多維的輸出向量,:(:,i sol 為i u 的輸出,即:,(:,i sol u i =。另

13、,元素,(,(i k j sol k j u i =表示在(j tspan t =和(k xmesh x =時(shí)i u 之答案。o p t i o n s 為解答器之相關(guān)解法參數(shù)。詳細(xì)請(qǐng)見odeset 之使用法。注:1. MATLAB PDE 解答器pdepe 之算法,主要事將原一組的橢圓型和拋物線型偏微分方程式轉(zhuǎn)化為一組常微分方程式。此轉(zhuǎn)換的過(guò)程系基于使用者所指定之mesh 點(diǎn),以二階空間離散化(spatialdiscretization技術(shù)為之(Keel and Berzins,1990,然后以ode15s 的指令求解。采用ode15s 的ode 解法,主要是因?yàn)樵陔x散化的過(guò)程中,橢圓型偏微

14、分方程式為被轉(zhuǎn)為一組代數(shù)式,而拋物線型的偏微分方程式則被轉(zhuǎn)為一組聯(lián)立的微分方程式。因而,原偏微分方程式被離散化后,遂變成一組同時(shí)伴有微分方程式與代數(shù)式之微分代數(shù)方程式系統(tǒng),故以ode15s 使可順利求解。2. x 之取點(diǎn)(mesh位置影響解的精確度甚鉅,若pdepe 解答器回應(yīng)出 "has difficulty finding consistent initial considition "之訊息時(shí),使用者可進(jìn)一步將mesh 點(diǎn)取密一點(diǎn),即增加mesh 點(diǎn)數(shù)。另外,若狀態(tài)u 再某些特定點(diǎn)上有較快速之變動(dòng)時(shí),亦需將此處之點(diǎn)取密集些,以增加精確度。值得注意的是pdepe 并不會(huì)

15、自動(dòng)做并不會(huì)自動(dòng)做xmesh 的自動(dòng)取點(diǎn),使用者必須觀察解的特性,自行作曲點(diǎn)的動(dòng)作。一般而言,所取之點(diǎn)數(shù)至少需大于3以上。3. tspan 之選取主要是基于使用者對(duì)那些特定時(shí)間之狀態(tài)有興趣而選定。而間距(step size之控制由程序自動(dòng)完成。4. 若要獲得特定位置及時(shí)間下之解,可配合以pdeval 的指令。跟pdepe 指令之格式如下:,(,xout ui xmesh m pdeval duoutdx uout =其中m 代表問(wèn)題之對(duì)稱性。m =0,代表平板;m =1,圓柱體;m =2表示球體。其意義同pdep e 中之自變量m 。x m e s h 使用者在pdepe 中所指定之輸出點(diǎn)之位

16、置向量。N x x x xmesh 10=ui 即:,(i j sol 。也就是說(shuō)其為pdepe 輸出中第i 個(gè)輸出變數(shù)ui 在各點(diǎn)位置xmesh 處,時(shí)間固定為(j tspan t j =下之解。x o u t為所欲內(nèi)插輸出點(diǎn)位置向量。此為使用者重新指定之位置向量。u o u t為基于所指定位置xout ,固定時(shí)間f t 下之相對(duì)應(yīng)輸出。d u o u t d x 為相對(duì)應(yīng)之dx du 輸出值。ref. Keel,R.D. and M. Berzins,”A Method for the Spatial Discritization ofParabolic Equations in One

17、Space Variable”,SIAM J. Sci. and Sat.Comput.,V ol.11,pp.1-32,1990.以下吾人將以數(shù)個(gè)范例,詳細(xì)說(shuō)明pdepe 之用法。范例6.1 試解以下之偏微分方程式222xu t u = 其中10x ,且滿足以下之條件限制式 (i起始值條件IC :sin(0,(x x u =(ii邊界條件BC1:0,0(=t uBC2:0,1(=+-xt u e t 注:本問(wèn)題之解析解為sin(,(x e t x u t -=解答:以下吾人將分述求解之步驟與過(guò)程。當(dāng)以下之各步驟完成后,可進(jìn)一步將其匯整為一主程序ex6_1.m ,然后求解之。步驟1 將欲解之偏

18、微分方程式改寫成如式子( 之標(biāo)準(zhǔn)式。 0002+ =x u x x x t u 此即(2,=x u u t x c(xu x u u t x f =, (0,=x u u t x s和 0=m步驟2 編寫偏微分方程式之系數(shù)向量檔function c,f,s=ex6_1pdefun(x,t,u,DuDxc=pi2;f=dudx;s=0;步驟3 編寫起始值條件function u0=ex6_1ic(xu0=sin(pi*x;步驟4 編寫邊界條件。在編寫之前,先將邊界條件改寫成標(biāo)準(zhǔn)形式,如式( ,找出相對(duì)之(p 和(q 的函數(shù),然后寫成邊界條件檔,例如,原邊界條件可寫成BC1:0,0,0(0,0(=

19、+x t x t u BC2:1,0,1(1=+-x t xu e t 即和,0,0(=ql t u pl1,=-qr e pr t 故而,邊界條件檔可編寫成function pl,ql,pr,qr=ex6_1bc(xl,ul,xr,ur,tpl=ul;ql=0;pr=pi*exp(-t;qr=1;步驟5 取點(diǎn)。例如x=linspace(0,1,20; %x 取20點(diǎn)t=linspace(0,2,5; %時(shí)間取5點(diǎn)輸出步驟6 利用pdepe 求解m=0; %依步驟1之結(jié)果sol=pdepe(m,ex6_1pdefun,ex6_1ic,ex6_1bc,x,t; 步驟7 顯示結(jié)果u=sol(:,:

20、,1;surf(x,t,utitle('pde 數(shù)值解'xlabel('位置'ylabel('時(shí)間' zlabel('u'若要顯示特定點(diǎn)上之解,可進(jìn)一步指定x 或t的位置,以便繪圖。例如,吾人欲了解時(shí)間為2(終點(diǎn)時(shí),各位置下之解,可輸入以下指令(利用pdeval指令:figure(2; %繪成圖2M=length(t; %取終點(diǎn)時(shí)間xout=linspace(0,1,100; %輸出點(diǎn)位置uout,dudx=pdeval(m,x,u(M,:,xout;plot(xout,uout; %繪圖title('時(shí)間為2時(shí),各處之解

21、'xlabel('x'ylabel('u'綜合以上各步驟,吾人可寫成一程序求解范例6.1。其參考程序如下:ex6_1.mfunction ex6_1%范例6_1%m=0;x=linspace(0,1,20; %xmesht=linspace(0,2,20; %tspan%以pde求解%sol=pdepe(m,ex6_1pdefun,ex6_1ic,ex6_1bc,x,t;u=sol(:,:,1; %取出答案%繪圖輸出%figure(1surf(x,t,utitle('pde數(shù)值解'xlabel('位置x'ylabel(&#

22、39;時(shí)間t' zlabel('數(shù)值解u'%與解析解做比較%figure(2surf(x,t,exp(-t'*sin(pi*x;title('解析解'xlabel('位置x'ylabel('時(shí)間t' zlabel('數(shù)值解u'%t=tf=2時(shí)各位置之解%figure(3M=length(t; %取終點(diǎn)時(shí)間xout=linspace(0,1,100; %輸出點(diǎn)位置uout,dudx=pdeval(m,x,u(M,:,xout;plot(xout,uout; %繪圖title('時(shí)間為2時(shí),各處

23、之解'xlabel('x'ylabel('u'%pde函數(shù)文件%function c,f,s=ex6_1pdefun(x,t,u,DuDxc=pi2;f=DuDx;s=0;%起始條件檔%function u0=ex6_1ic(xu0=sin(pi*x;%邊界條件檔%function pl,ql,pr,qr=ex6_1bc(xl,ul,xr,ur,t pl=ul;ql=0;pr=pi*exp(-t;qr=1;執(zhí)行結(jié)果: 范例6_2 試解以下之聯(lián)立pde 系統(tǒng)(024.0212121u u F xu t u -=(170.0212222u u F xu t

24、u -+= 其中(46.11exp(73.5exp(212121u u u u u u F -=- 且10x 和0t 。此聯(lián)立pde 系統(tǒng)滿足以下之條件限制式(i起始值條件10,(1=x u00,(2=x u(ii邊界值條件0,0(1=t xu 0,0(2=t u 1,1(1=t u0,1(2=t xu 解答:步驟1: 改寫PDE 方程式為標(biāo)準(zhǔn)式-+=11c -=(2121u u F u u F s 和0=m 。另外,左邊界條件(處0=x 。寫成=20u pl =01ql 同理,右邊界條件(處1=x 為=-=011u pr =10qr 步驟2: 編寫pde 函數(shù)文件function c,f,s

25、=ex6_2pdefun(x,t,u,DuDxc=1 1'f=0.024 0.170'.*DuDx;y=u(1-u(2;F=exp(5.73*y-exp(-11.47*y;s=-F F'步驟3:編寫起始條件function u0=ex6_2ic(xu0=1 0'步驟4:編寫邊界條件檔function pl,ql,pr,qr=ex6_2bc(xl,ul,xr,ur,t pl=0 ul(2'ql=1 0'pr=ur(1-1 0'qr=0 1'步驟5: 取點(diǎn)由于此問(wèn)題之端點(diǎn)均受邊界條件之限制,且時(shí)間t 小時(shí)狀態(tài)之變動(dòng)甚鉅(由多次求解后之

26、經(jīng)驗(yàn)得知,故在兩端點(diǎn)處之點(diǎn)可稍微密集些。同時(shí)對(duì)于t 小(t small處亦可取密一些。例如, x=0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.9951;t=0 0.005 0.01 0.05 0.1 0.5 1 1.5 2;以上幾個(gè)主要步驟編寫完成后,事實(shí)上就可直接完成主程序來(lái)求解。此問(wèn)題之參考程序如下:ex6_2.mfunction ex6_2%范例6_2%m=0;x=0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1;t=0 0.005 0.01 0.05 0.1 0.5 1

27、 1.5 2;%利用pde求解%sol=pdepe(m,ex6_2pdefun,ex6_2ic,ex6_2bc,x,t;u1=sol(:,:,1; %第一個(gè)狀態(tài)之?dāng)?shù)值解輸出u2=sol(:,:,2; %第二個(gè)狀態(tài)之?dāng)?shù)值解輸出%繪圖輸出%figure(1surf(x,t,u1title('u1之?dāng)?shù)值解'xlabel('x'ylabel('t'%figure(2surf(x,t,u2title('u2之?dāng)?shù)值解'xlabel('x'ylabel('t'%pde函數(shù)%function c,f,s=ex6_2

28、pdefun(x,t,u,DuDxc=1 1'f=0.024 0.170'.*DuDx;y=u(1-u(2;F=exp(5.73*y-exp(-11.47*y;s=-F F'%起始值條件%function u0=ex6_2ic(xu0=1 0'%邊界值條件%function pl,ql,pr,qr=ex6_2bc(xl,ul,xr,ur,tpl=0 ul(2'ql=1 0'pr=ur(1-1 0'qr=0 1'執(zhí)行結(jié)果 有別于pdepe指令所能處理之一維偏微分方程式,MATLAB之PDE圖形接口工具箱,可允許使用者處理某幾類的二維

29、空間之偏微分方程式問(wèn)題。利用此工具箱,使用者僅需要再命令六窗口中鍵入 透過(guò)此圖形接口其解法步驟大致上為:(1定義PDE問(wèn)題,包括二維空間范圍,邊界條件以及PDE系數(shù)等(2產(chǎn)生離散化之點(diǎn),并將原PDE方程式離散化。(3利用有限元素法(finite element method;FEM求解并顯示答案在詳細(xì)說(shuō)明此解法工具之前,吾人將先介紹此工具所能處理之PDE問(wèn)題形式。此PDE圖形接口之選單下方,存在一排主要功能的圖標(biāo)(icon按鈕 透過(guò)這些按鈕,使用者可輕松地完成PDE方程式之求解?，F(xiàn)并將這些按鈕之主要功能敘述如下:前五個(gè)按鈕為PDE系統(tǒng)之邊界范圍繪制功能,由左至右之用法為: 方形需以按住右鍵的方

30、式繪制之。 標(biāo)左鍵繪矩形,右鍵繪正方形。 橢圓,而右鑒用來(lái)繪制餅圖形。 :以中心點(diǎn)向外的方式繪制橢圓或圓。同樣地,鼠標(biāo)左及又鍵,分別用以繪制橢圓及圓形的區(qū)域。:用以繪制多邊型等不規(guī)則區(qū)域,欲關(guān)閉此功能需按鼠標(biāo)右鍵。 在這些繪制按鈕之后之按鈕功能依序如下:用以給定邊界條件。在此功能選定后,使用者可在任一圖形 邊界上按住鼠標(biāo)左鍵兩次,然后在對(duì)話框上將邊界條件給入。:用以指定PDE 問(wèn)題及相關(guān)參數(shù)。 :產(chǎn)生圖形區(qū)域內(nèi)離散化之網(wǎng)點(diǎn)。 :用以進(jìn)一步將離散化之網(wǎng)點(diǎn)再取密一點(diǎn)(refine mesh。 :在指定PDE 系統(tǒng),邊界條件及圓形區(qū)域后,按此鈕即開始 解題。:用以指定顯示結(jié)果繪制方式。 :放大縮小之

31、功能,便于圖形繪制及顯示。 pdetool 圖形接口工具所能處理之PDE 問(wèn)題的基本型式為以下之橢圓型方程式f au u c =+-(其中yj xi +為L(zhǎng)aplacian 運(yùn)操作數(shù)(operator。c ,a ,f 和u 為定義在有限平面幾下之純量或復(fù)變函數(shù)值函數(shù)(complex valued function。除了這個(gè)基本問(wèn)題外,此工具箱亦可解以下的PDE 問(wèn)題。拋物線型f au u c t u d =+-( 雙曲線型f au u c tu d =+-(22 特征值問(wèn)題du au u c =+-(其中為未知特征值以上這些問(wèn)題之邊界條件型式可為Dirichlet條件hu=,rGenerali

32、zed Neumann條件+(ugn=quc其中n 為向外之單位法向量(Outward unit normal,g,q,h和r為復(fù)數(shù)值函數(shù)。要利用pdetool接口求解之前,需先定義PDE問(wèn)題,其包含三大部份:(1利用繪圖(draw模式,定義欲解問(wèn)題的空間范圍(domain(2利用boundary模式,指定邊界條件(3利用PDE模式,指定PDE系數(shù),即輸入c,a,f和d等PDE模式中之系數(shù)在定義PDE問(wèn)題之后,可依以下兩個(gè)步驟求解(1在mesh模式下,產(chǎn)生mesh點(diǎn),以便將原問(wèn)題離散化(2在solve模式下,求解(3最后,在Plot模式下,顯示答案現(xiàn)在吾人將以sPoisson'方程式f

33、-之求解為范例,詳細(xì)說(shuō)明pdetoolu=之用法。此問(wèn)題之幾何圖形及相關(guān)邊界條件,將于求解過(guò)程中加以說(shuō)明。步驟1:在命令列窗口中鍵入pdetool以進(jìn)入GUI(graphical userinterface界面。選取Options中之Grid功能,以顯示網(wǎng)格線。 步驟2:利用Draw功能,畫出問(wèn)題之幾何圖形。請(qǐng)注意:使用者可利用內(nèi)定對(duì)象”多邊型”,”矩形”,”正方形”,”圓形”,及”橢圓型”,予以組合,例如(i先選取”矩形/正方形”對(duì)象,移動(dòng)由標(biāo)至所欲輸入左上角之點(diǎn),如坐標(biāo)(-1,0之點(diǎn),按住鼠標(biāo)左鍵,往右下角拉至坐標(biāo)為(1,-0.4處,即行成代號(hào)為R1之矩形。其余圖形C1,R1和C2可選取適

34、當(dāng)之對(duì)象制作之,以形成如下圖之圖形區(qū)域。以代數(shù)公式而言,其為R1+C1+R2+C2 值得注意的是,圓形之區(qū)域需以按住鼠標(biāo)右鍵的方式來(lái)制作(非左鍵。同時(shí),如欲進(jìn)一步修改各圖形對(duì)象之大小及位置數(shù)據(jù),可在該圖上按鼠標(biāo)左鍵兩次,然后在對(duì)象對(duì)話框上輸入數(shù)據(jù)。(ii若所欲形成之圖形區(qū)域,需將C2去除,則可在公式列中直接輸入R1+C1+R2-C2即可 步驟3:選取PDE功能項(xiàng),以輸入PDE方程式的系數(shù)及類型。因問(wèn)題為f-,故此為橢圓型的問(wèn)題,且其標(biāo)準(zhǔn)形式u=-(比較得知,c=1,a=0和f=10,所以對(duì)話框輸+auuc=f入之情況如下: 步驟4:選取Boundary功能,以輸入邊界條件。假設(shè)邊界條件為Neu

35、mann 形,且為5u。其中在弧形部份與標(biāo)準(zhǔn)式知,g=-5且q=0。n=但直線部份其邊界條件則在Dirichlet type使h=0,r=0。故而 步驟5:選取Mesh功能,產(chǎn)生網(wǎng)點(diǎn)。使用者亦可進(jìn)一步利用將網(wǎng)點(diǎn)取在密一點(diǎn)(refine mesh。 步驟6:選取solve功能,解此PDE。 之”parameters”功能,選擇適當(dāng)之結(jié)果顯示圖形及數(shù)據(jù)。例如, 2.另外,若使用者欲將結(jié)果輸出至命令列窗口中,以供后續(xù)處理,可利用solve功能項(xiàng)下之”export solution”指定變量名稱來(lái)完成。6.3 化工實(shí)例演練以外部熱交換式的管形固定層觸煤反應(yīng)裝置,進(jìn)行苯加氫反應(yīng)產(chǎn)生環(huán)己烷。此反應(yīng)系統(tǒng)之質(zhì)

36、量平衡及熱平衡方程式如下:0(122=-+ +-p r B A p e GC H r r T r r T GC k L T 01022=+ +-Gy M r r f r r f u D L f ar B A e 其中T 為溫度(,f 為反應(yīng)率,L 為軸向距離,r 為徑向距離。此系統(tǒng)之邊界條件為0=L , (0r T T = , (0r f f = 0=r ,0=rfr T w r r = , (w w e T T h rTk -=- w r r = ,0=rf此外,式中之相關(guān)數(shù)據(jù)及操作條件如下: (i反應(yīng)速率式 23H 氫4331(C C B B H H BH B H A P K P K P

37、K P P K kK r +=其中P 表示分壓(atm,而速率參數(shù)為 R RT K 3.3212100ln +-= R RT K H 9.3115500ln -= R RT K B 1.2311200ln -= R RT K C 4.198900ln -=上式中,下標(biāo)B ,H 及C 分別代表苯,氫及環(huán)己烷。R 為理想氣體常數(shù)(1.987cal/mol .K(ii操作條件及物性數(shù)據(jù)總壓 a t m P t 25.1= 反應(yīng)管管徑 cm r w 5.2=壁溫100=w T 質(zhì)量速度hr m kg G =2631苯對(duì)氫之莫耳流量比 30=m 反應(yīng)管入口的苯之莫耳分率 0323.00=y 反應(yīng)氣體之平

38、均分子量 47.4=av M觸煤層密度 31200m kg P B =流體平均比熱 m o l kg kcal C P -=74.1 反應(yīng)熱 mol kg kcal H r -=49250 整體傳熱系數(shù) C hr m kcal h 0208.65=進(jìn)料溫度 1250(=T 反應(yīng)管管長(zhǎng) cm L 45=流速hr m u 03.8=有效熱傳導(dǎo)系數(shù) C hr m kcal K e 065.0= 壁境膜傳熱系數(shù) C hr m kcal h w 02112= 有效擴(kuò)散系數(shù)hr m D e 2755.0=題意解析:因反應(yīng)速率式A r 與分壓有關(guān),而分壓又與反應(yīng)率f 有關(guān)。故需進(jìn)一步將A r 由反應(yīng)率f 表

39、示之,方能求解偏微分方程式?；谝韵轮磻?yīng)式2 3H m -f -3f f 1-fm-3f f則各分壓與總壓之關(guān)系為 f m fm P P t H 313-+-=f m fm P P t B 31-+-=fm fP P tC 31-+=將上式( ( ,連同反應(yīng)速率式( ,帶入平衡方程式( 中,配合邊界條件( ,可利用pdepe 求解之。MATLAB 程序設(shè)計(jì)將原方程式改寫成如( 之標(biāo)準(zhǔn)式-+=-01(1001Gy M r GC H r r f u D r r T GC k r r r L f L T av B A p r B A e p e 因此=11c =r f uD r T GC k f

40、e pe 和-=0(Gy M r GC H r s av B A p r B A 和1+=m (圓柱。另外,左邊界條件(0=r 處寫成=+00*1100f即=00pl , =11ql同理右邊界條件(w r r =可寫成=+-00*10(f GC T T h p w w 即-=0(w w T T h pr +=0p GC qr 根據(jù)以上的分析,吾人可編寫一MATLAB 程序求解此PDE 問(wèn)題,其參考程序如下:ex6_3_1.mfunction ex6_3_1% 范例6_3_1 觸媒反應(yīng)器內(nèi)溫度及轉(zhuǎn)化率的分布%global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u

41、 R ke hw De % 給定數(shù)據(jù)%Pt=1.25; %總壓(atmrw=0.025; %管徑(mTw=100+273; %壁溫(G=631; %質(zhì)量流率(kg/m2hrM=30;y0=0.0323;Mav=4.47;rho_B=1200;Cp=1.74;dHr=-49250;h0=65.8;T0=125+273;Lw=1;u=8.03;R=1.987;ke=0.65;hw=112;De=0.755;%m=1;% 取點(diǎn)%r=linspace(0,rw,10;L=linspace(0,Lw,10;% 利用pdepe求解%sol=pdepe(m,ex6_3_1pdefun,ex6_3_1ic,e

42、x6_3_1bc,r,L; T=sol(:,:,1; %溫度f(wàn)=sol(:,:,2; %反應(yīng)率% 繪圖輸出%figure(1surf(L,r,T'-273title('temp'xlabel('L'ylabel('r'zlabel('temp (0C'%figure(2surf(L,r,f'title('reaction rate'xlabel('L'ylabel('r'zlabel('reaction rate'% PDE 函數(shù)%function c

43、1,f1,s1=ex6_3_1pdefun(r,L,u1,DuDrglobal Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De T=u1(1;f=u1(2;%k=exp(-12100/(R*T+32.3/R;Kh=exp(15500/(R*T-31.9/R;Kb=exp(11200/(R*T-23.1/R;Kc=exp(8900/(R*T-19.4/R;%a=1+M-3*f;ph=Pt*(M-3*f/a;pb=Pt*(1-f/a;pc=Pt*f/a;%rA=k*Kh3*Kb*ph3*pb/(1+Kh*ph+Kb*pb+Kc*pc4; %c1=1

44、 1'f1=ke/(G*Cp De/u'.*DuDr;%s1=ke/(G*Cp*r*DuDr(1-rA*rho_B*dHr/(G*Cp-2*h0*(T-Tw/(rw s1=-rA*rho_B*dHr/(G*CprA*rho_B*Mav/(G*y0;% 起始值條件%function u0=ex6_3_1ic(xu0=125+273 0'% 邊界條件檔%function pl,ql,pr,qr=ex6_3_1bc(rl,ul,rr,ur,Lglobal Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw Depl=0 0'ql

45、=1 1'pr=hw*(ur(1-Tw 0'qr=G*Cp 1'結(jié)果: 范例 6_3_2 擴(kuò)散系統(tǒng)之濃度分布(C&M,p.403參考如圖之裝置。管中儲(chǔ)放靜止液體B ,高度為L(zhǎng)=10,其暴置于充滿A 氣體的環(huán)境中。假設(shè)與B 液體接觸面之濃度為3001.0m mol C A ,且此濃度不隨時(shí)間改變而改變,即在操作時(shí)間內(nèi)(10=h 天維持定值。氣體A 在液體B 中之?dāng)U散系數(shù)為s m D AB 29102-=。試決定以下兩種情況下,氣體A 溶于液體B 中之流通量(flux。(aA 與B 不管發(fā)生反應(yīng)(bA 與B 發(fā)生以下之反應(yīng)C B A + , A A kC r =-

46、其反應(yīng)速率常數(shù)17102-=s k 。 10cm z 氣體A 圖 氣體A 在液體B 中之?dāng)U散題意解析:(a因氣體A 與液體B 不發(fā)生反應(yīng),故其擴(kuò)散現(xiàn)象之質(zhì)量平衡式如下:22zC D t C A AB A = 依題意,其起始及邊界條件為I.C. 00,(=z C A , 0>zB.C. 0,0(A A C t C = , 0t0=Lz Az C , 0t(b在氣體A 與液體B 會(huì)發(fā)生一次反應(yīng)之情況下,其質(zhì)量平衡式需改寫為A A AB A kC zC D t C +=22 而起始及邊界條件同上。在獲得濃度分布后,即可以Ficks law(=-=z A AB Az z C D t N 計(jì)算流通量。MATLAB 程序設(shè)計(jì):此問(wèn)題依舊可以pdepe 迅速求解。現(xiàn)并就各狀況之處理過(guò)程簡(jiǎn)述之(a將式( 與標(biāo)準(zhǔn)式( 比