【文章內(nèi)容簡介】 f,s]=ex20_2pdefun(x,t,u,dudx) c=[1 1]39。 f=[ ]39。.*dudx。 y=u(1)u(2)。 F=exp(*y)exp(*y)。 s=[F F]39。 步驟 3：編寫初始條件函數(shù) function u0=ex20_2ic(x) u0=[1 0]39。 步驟 4：編寫邊界條件函數(shù) function [pl,ql,pr,qr]=ex20_2bc(xl,ul,xr,ur,t) pl=[0 ul(2)]39。 ql=[1 0]39。 pr=[ur(1)1 0]39。 qr=[0 1]39。 步驟 5： 取點。 由于此問題的端點均受邊界條件的限制，且時間 t 很小時狀態(tài)的變動很大(由多次求解后的經(jīng)驗得知)，故在兩端點處的點可稍微密集些。同時對于 t 小處亦可取密一些。例如， x=[0 1]。 t=[0 1 2]。 以上幾個主要步驟編寫完成后，事實上就可直接完成主程序來求解。此問題的參考程序如下： function ex20_2 %***************************************%求解一維偏微分方程組的一個綜合函數(shù)程序 %*************************************** m=0。 x=[0 1]。 t=[0 1 2]。 %************************************* %利用 pdepe 求解 %************************************* sol=pdepe(m,@ex20_2pdefun,@ex20_2ic,@ex20_2bc,x,t)。 u1=sol(:,:,1)。 %第一個狀態(tài)之數(shù)值解輸出 u2=sol(:,:,2)。 %第二個狀態(tài)之數(shù)值解輸出 %************************************* %繪圖輸出 %************************************* figure(1) surf(x,t,u1) title(39。u1 之數(shù)值解39。) xlabel(39。x39。) ylabel(39。t39。) % figure(2) surf(x,t,u2) title(39。u2 之數(shù)值解39。) xlabel(39。x39。) ylabel(39。t39。) %*************************************** %pde 函數(shù) %*************************************** function [c,f,s]=ex20_2pdefun(x,t,u,dudx) c=[1 1]39。 f=[ ]39。.*dudx。 y=u(1)u(2)。 F=exp(*y)exp(*y)。 s=[F F]39。 %**************************************** %初始條件函數(shù) %**************************************** function u0=ex20_2ic(x) u0=[1 0]39。 %**************************************** %邊界條件函數(shù) %**************************************** function [pl,ql,pr,qr]=ex20_2bc(xl,ul,xr,ur,t) pl=[0 ul(2)]39。 ql=[1 0]39。 pr=[ur(1)1 0]39。 qr=[0 1]39。 化工應用實例 例 4 觸煤反應裝置內(nèi)溫度及轉(zhuǎn)換率的分布 以外部熱交換式的管形固定層觸煤反應裝置，進行苯加氫反應產(chǎn)生環(huán)己烷。此反應 系統(tǒng)之質(zhì)量平衡及熱平衡方程式如下： 其中T 為溫度(℃)， f 為反應率，L 為軸向距離，r 為徑向距離。此系統(tǒng)的邊界條件為 此外，式中之相關(guān)數(shù)據(jù)及操作條件如下： (i)反應速率式 其中 P 表示分壓(atm)，而速率參數(shù)為 上式中，下標 B，H 及 C 分別代表苯，氫及環(huán)己烷。R 為理想氣體常數(shù)(K)。(ii)操作條件及物性數(shù)據(jù) 題意解析： 因反應速率式 A r 與分壓有關(guān)，而分壓又與反應率 f 有關(guān)。故需進一步將 A r 由反應 率 f 表示，方能求解偏微分方程?；谝韵碌姆磻匠?則各分壓與總壓之關(guān)系為 將上式，連同反應速率式，帶入平衡方程式中，配合邊界條件，可利用 pdepe 求解。 MATLAB 程序設計 將原方程改寫成如式（35）的標準式 因此 和m = +1(圓柱)。另外，左邊界條件( r = 0處)寫成 即 同理右邊界條件()可寫成 即 根據(jù)以上的分析，可編寫 MATLAB 程序求解此 PDE 問題，其參考程序如下： function ex60_3_1 %****************************** % 觸媒反應器內(nèi)溫度及轉(zhuǎn)化率的分布 %****************************** global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De %****************************** % 給定數(shù)據(jù) %****************************** Pt=。 %總壓(atm) rw=。 %管徑(m) Tw=100+273。 %壁溫(℃) G=631。 %質(zhì)量流率(kg/m2hr) M=30。 y0=。 Mav=。 rho_B=1200。 Cp=。 dHr=49 h0=。 T0=125+273。 Lw=1。 u=。 R=。 ke=。 hw=112。 De=。 %******************** m=1。 %******************** % 取點 %******************** r=linspace(0,rw,10)。 L=linspace(0,Lw,10)。 %*********************** % 利用 pdepe 求解 %*********************** sol=pdepe(m,@ex20_3_1pdefun,@ex20_3_1ic,@ex20_3_1bc,r,L)。 T=sol(:,:,1)。 %溫度 f=sol(:,:,2)。 %反應率 %*********************** % 繪圖輸出 %*********************** figure(1) surf(L,r,T39。273) title(39。temp39。) xlabel(39。L39。) ylabel(39。r39。) zlabel(39。temp (0C)39。) % figure(2) surf(L,r,f39。) title(39。reaction rate39。) xlabel(39。L39。) ylabel(39。r39。) zlabel(39。reaction rate39。) %************************************************* % PDE 函數(shù) %************************************************* function [c1,f1,s1]=ex20_3_1pdefun(r,L,u1,DuDr) global 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)+)。 Kh=exp(15500/(R*T))。 Kb=exp(11200/(R*T))。 Kc=exp(8900/(R*T))。 % a=1+M3*f。 ph=Pt*(M3*f)/a。 pb=Pt*(1f)/a。 pc=Pt*f/a。 % rA=k*Kh^3*Kb*ph^3*pb/(1+Kh*ph+Kb*pb+Kc*pc)^4。 % c1=[1 1]39。 f1=[ke/(G*Cp) De/u]39。.*DuDr。 %s1=[ke/(G*Cp*r)*DuDr(1)rA*rho_B*dHr/(G*Cp)2*h0*(TTw)/(rw) s1=[rA*rho_B*dHr/(G*Cp)。rA*rho_B*Mav/(G*y0)]。 %********************************** %初始條件函數(shù) %********************************** function u0=ex20_3_1ic(x) u0=[125+273 0]39。 %********************************** % 邊界條件%********************************** function [pl,ql,pr,qr]=ex20_3_1bc(rl,ul,rr,ur,L) global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De pl=[0 0]39。 ql=[1 1]39。 pr=[hw*(ur(1)Tw) 0]39。 qr=[G*Cp 1]39。 例 5 擴散系統(tǒng)之濃度分布 參考如圖 3 的裝置。管中儲放靜止液體 B，高度為 L=10 ㎝，放置于充滿 A 氣體的環(huán)境中。假設與 B 液體接觸面之濃度為