偏微分方程式之求解

1、第六章 偏微分方程式之求解在化工的领域中，有不少程序之动态是由以偏微分方程式 (Partial differential equation；PDE)所描述的，例如热与质量在空间中的传递等。这些用以描述实际 问题的 PDE，除非具有某些特定的方程式型态及条件，否则甚难以手算的方式 找出解析解。而在数值求解方面，最常被采用的方法为有限差分法 (finite difference)何有限元素法 (finite element)。然对于某些不熟悉数值分析及程序编写 的化工人而言，欲充分了解以偏微分方程式所描述之系统动态是相当不容易的， 更遑论进一步的设计与分析了。值得庆幸的是， MATLAB 的环境中

2、提供了一个求解 PDE 问题的工具箱，让 使用者得以利用简单的指令或图形接口工具输入欲解的 PDE，并求解。使得 PDE 之数值解在弹指之间完成， 使用者不在为数值法所苦恼， 轻松掌握偏微分方程式 系统的动态，并可进一步进行后续之设计工作。本章将以循渐进的方式，介绍 PDE 工具箱及其用法，并以数个典型的化工 范例进行 示范，期能使初学者很快熟悉 PDE工具箱， 并使用它来设计与分析 以偏微方方程式所描述的程序系统。6.1 偏微分方程式之分类偏微分方程式可根据其阶数 (order)，线性或非线性型态，以及边界条件进 行分类。6.1.1依阶数的分类偏微分方程式是以偏微分项中之最高次偏微分来定义其

3、阶数，例如：一阶偏微分方程式：u u 0 xy二阶偏微分方程式：32y2u u 0 xy三阶偏微分方程式：3u3x2uxyu u 0 xy6.1.2 依非线性程度分类偏微分方程式亦可以其线性或非线性情况，区分为线性 (linear)，似线性 (quasilinear)，以及非线性三类。 例如，以下之二阶偏微分方程式 (Constantinides and Mostoufi,1999)22a()b( ) u c( ) u2 d() 0x y x可依系数 ( )之情况，进行如下表之归类类别情况线性 似线性系数 ( )为定值，或仅为 (x,y)函数系数 ( )为依变数 (dependent vari

4、able)u或其比方程式中之偏微 分低阶之偏微分项的函数，如 ( ) (x,y,u, u x, u y)非线性系数 ( )中，具有与原方程式之偏微分同阶数之变数，如 () (x,y,u, 2u x2 , 2u y2, 2u x y)另外，对于线性二阶偏微分方程式，可进一步将其分类为椭圆型 (elliptic) ，拋 物线型 (parabolic)，以及双曲线型 (hyperbolic) 。具体上来说，此类偏微分方程 式二阶线性之一般式为xy22uuc2dy2xfu g 0系数a,b,c,d,e和 f 是定值或为 u 的函数。若 g=0，则上式为其次是偏微分方程 式。式子 ( )之分类及代表性例

5、子，请见下表方程式类别 判断式椭圆型b2 4ac 0拋物线型b2 4ac 0代表性范例Laplace方程式，Poisson方程式，x22y20f (x,y)(c u) au f热传导或扩散方程式2u2xd u(c u) au f双曲线型22b2 4ac 0 波动方程式 2 u2u2xt2 u d 2(c u) au ft2注：二维系统之运操作数 之定义为 i j xy6.1.3 起始条件和边界条件的分类为了能获得偏微分方程式之解答， 为三类。现以一维之动态热传递方程式其起始条件和边界条件可依其特性区分(拋物线型偏微分方程式 )2Tx2为例，进一步说明如何区分这些边界条件及起始条件 (Const

6、antinides and Mostoufi ,1999)。(i) 第一类： Dirichlet Condiction若依变量 (T)本身， 在某个独立变量值时， 被指定， 则此条件称为 Dirichlet Condiction，亦称为 essential边界条件。下图为一典型的 Dirichlet 条件示 意图TT1,t0T001T f (t), t0图 6.1 平板 Dirichlet Condiction 示意图由图中很清楚的显示，该平板之边界条件为T f (t) ,at x0,t0T T1 ,at x0,t0T T0 ,at x0,t0此边界条件依定义，即为Dirichlet Cond

7、iction 。同时，若再起始时，各处之温度分布可以位置之函数表示，即T f (t) ,at t0,0x1此亦属 Dirichlet 型之边界条件。(ii) 第二类： Neumann conditionNeumann condition 系指依变量之变化率之边界条件为定值，抑或独立变量之函数之情况。例如0 , atx 1,t0f (t) , at t 0, 0 x 1xNeumann 型边界条件，亦称为 natural boundary condition。(iii) 第三类： Robbins condition 若依变量之变化率之边界条件，为自身之函数 (非独立变量之函数 ) 时，被称为 R

8、obbins condition。例如，k T h(T Tf ) , at x 0,t 0 x上式之边界条件，当发生在固液相间之传递上，亦即热流通量 (heat flux) 正比于固液两端之温差，其示意图如下：t液相,Tf固體平板TTk h(TTf ),xt001x(iv) Cauchy conditionCauchy condition 系指问题中同时存在 Dirichlet 和 Neumann边界条 件。例如下图即T f(t) , t 0Dirichlet condition0, x 1,t 0Neumann conditiontTx00f(t)0 , at x01絕緣6.2 MATLAB

9、 PDE 工具箱6.2.1 MATLAB PDE 解答器MATLAB 提供了一个指令 pdepe，用以解以下的 PDE 方程式c(x,t,u, u) u x m (xmf (x,t,u, u) s(x,t,u, u)x t x x x其中时间介于 t0 t tf 之间，而位置 x则介于 a b 有限区域之间。 m值 表示问题之对称性，其可为 0，1 或 2，以表示平板 (slab)，圆柱 (cylindrical) 或 球体(spherical)之对秤情况。因而，如果如果 m0，则 a 必等于 b，也就是说其 具有圆柱或球体之对称关系。同时，式中 f(x,t,u, u x) 一项为流通量 (f

10、lux) ， 而s(x,t,u, u x)为来源(source)项。c(x,t,u, u x) 为偏微分方程式之对角线系 数矩阵。若某一对角线元素为 0，则表示该相应偏微分方程式为椭圆型偏微分 方程式，若为正值 (不为 0)，则为拋物线型偏微分方程式。请注意 c 之对角线 元素必不全为 0。偏微分方程式之起始值可表示为 u(x,t0) v0 (x)而边界条件之形式为p(x,t,u) g(x,t) f(x,t,u, u x) 0其中 x 为两端点位置，即 a 或 b。用以解含上述起始值及边界值条件的偏微分方程式 ( )之指令 pdepe 的 用法如下：sol pdepe(m, pdepe,icf

11、un,bcfun,xmesh,tspan,options)其中m 为问题之对称参数xmesh 为独立变量 x 之网取点 (mesh)位置向量，即xmesh x0 x1 x2 xN ，其中 x0 a(起点)，xN b(终点)。tspan 为独立变量 t(时间)之向量，即 tspan t0 t1tM ，其中t0为起始时间， tM 为终点时间。pdefun 为使用者提供之 pde 函数文件。其档案之格式如下： c, f,s pdefun(x,t,u,dudx) 亦即，使用者仅需提供偏微分方程式中之系数向量。c， f 和s均为行 (column)向量，而向量 c即为矩阵 c之对角线元素。icfun 提

12、供解 u之起始值，其格式为 u icfun (x) 值得注意的 是 u 为行向量。bcfun 使用者提供之边界条件档，格式如下：pl,ql,pr,ql bcfun xl,ul,xr,ur,t 其中 ul和ur 分别表 示左边界 (xl b) 和右边界 (xr a) 之 u 的近似解。输 出变量中， pl 和 ql 分别表示左边界 p 和 q 之行向量， 而 pr 和 qr 则为右边界 p 和 q 之行向量。sol 为解答输出。 sol 为多维的输出向量， sol(:,:i)为 ui的输出，即ui sol(:, :, i ) 。另，元素 ui (j,k) sol(j,k,i)表示在 t tspa

13、n( j)和 x xmesh(k)时 ui之答案。options 为解答器之相关解法参数。 详细请见 odeset之使用法。 注：1. MATLAB PDE 解答器 pdepe之算法，主要事将原一组的椭圆型和 拋物线型偏微分方程式转化为一组常微分方程式。 此转换的过程系 基于使用者所指定之 mesh 点，以二阶空间离散化 (spatial discretization)技术为之 (Keel and Berzins,1990)，然后以 ode15s的指 令求解。采用 ode15s的 ode 解法，主要是因为在离散化的过程中， 椭圆型偏微分方程式为被转为一组代数式， 而拋物线型的偏微分方 程式则被

14、转为一组联立的微分方程式。 因而， 原偏微分方程式被离 散化后，遂变成一组同时伴有微分方程式与代数式之微分代数方程 式系统，故以 ode15s使可顺利求解。2. x 之取点(mesh)位置影响解的精确度甚鉅，若 pdepe解答器回应出 has difficulty finding consistent initial considition 之讯息时，使 用者可进一步将 mesh点取密一点，即增加 mesh点数。另外，若 状态 u 再某些特定点上有较快速之变动时， 亦需将此处之点取密集 些，以增加精确度。 值得注意的是 pdepe并不会自动做并不会自动 做 xmesh的自动取点， 使用者必须观

15、察解的特性， 自行作曲点的动 作。一般而言，所取之点数至少需大于 3 以上。3. tspan之选取主要是基于使用者对那些特定时间之状态有兴趣而选 定。而间距 (step size)之控制由程序自动完成。4. 若要获得特定位置及时间下之解，可配合以 pdeval 的指令。跟 pdepe指令之格式如下：uout , duoutdx pdeval (m, xmesh,ui , xout )其中m 代表问题之对称性。 m=0 ，代表平板； m =1，圆柱 体；m=2 表示球体。其意义同 pdepe中之自变量 m。 xmesh 使用者在 pdepe中所指定之输出点之位置向量。 xmesh x0 x1xN

16、ui 即sol(j,:,i)。也就是说其为 pdepe输出中第 i 个输出 变数 ui 在各点位置 xmesh处，时间固定为t j tspan( j) 下之解。xout 为所欲内插输出点位置向量。 此为使用者重新指定 之位置向量。uout 为基于所指定位置 xout ，固定时间 t f 下之相对应输出。duoutdx 为相对应之 du dx 输出值。ref. Keel,R.D. and M. Berzins,” A Method for the Spatial Discritization ofParabolic Equations in One Space Variable” ,SIAM J

17、. Sci. and Sat.Comput.,Vol.11,pp.1-32,1990.以下吾人将以数个范例，详细说明 pdepe 之用法。范例 6.1 试解以下之偏微分方程式22 u utx2其中 0 x 1 ，且满足以下之条件限制式(i) 起始值条件IC： u(x,0) sin( x)(ii)边界条件BC1： u(0,t) 0BC2： e t u(1,t) 0x 注：本问题之解析解为 u(x,t) e t sin( x)解答： 以下吾人将分述求解之步骤与过程。 当以下之各步骤完成后， 可进一 步将其汇整为一主程序 ex6_1.m，然后求解之。步骤 1 将欲解之偏微分方程式改写成如式子 ( )

18、之标准式。2 u 0 0 uxxt x x此即c x,t,u, u xf x,t,u, u xxs x,t,u, u x 0和 m 0步骤 2 编写偏微分方程式之系数向量档function c,f,s=ex6_1pdefun(x,t,u,DuDx)c=pi2;f=dudx;s=0;步骤 3 编写起始值条件function u0=ex6_1ic(x)u0=sin(pi*x);步骤 4 编写边界条件。 在编写之前， 先将边界条件改写成标准形式， 如式( )，找出相对之 p( )和q( )的函数，然后写成边界条件 档，例如，原边界条件可写成BC1： u(0,t) 0 (0,t) 0 , x 0xBC

19、2： e t 1 u(1,t) 0 , x 1x即pl u(0,t) , ql 0 , 和tpr e , qr 1故而，边界条件档可编写成function pl,ql,pr,qr=ex6_1bc(xl,ul,xr,ur,t)pl=ul;ql=0;pr=pi*exp(-t);qr=1;步骤 5 取点。例如x=linspace(0,1,20); %x取 20 点t=linspace(0,2,5); %时间取 5 点输出步骤 6 利用 pdepe求解m=0; % 依步骤 1 之结果sol=pdepe(m,ex6_1pdefun,ex6_1ic,ex6_1bc,x,t);步骤 7 显示结果u=sol(

20、:,:,1);surf(x,t,u)title(pde 数值解 )xlabel(位置 )ylabel(时间 )zlabel(u)若要显示特定点上之解， 可进一步指定 x 或t 的位置，以便绘图。例 如，吾人欲了解时间为 2(终点)时，各位置下之解， 可输入以下指令 (利 用 pdeval 指令 ) ：figure(2); % 绘成图 2M=length(t); % 取终点时间xout=linspace(0,1,100); % 输出点位置 uout,dudx=pdeval(m,x,u(M,:),xout);plot(xout,uout); %绘图title( 时间为 2时, 各处之解 )xlab

21、el(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(1)surf(x,t,u)title(pde 数值解 )xlabel( 位置 x)ylabel( 时间 t )zlabel( 数值解 u) %与解析解做比较

22、%figure(2) surf(x,t,exp(-t)*sin(pi*x);title( 解析解 ) xlabel( 位置 x) ylabel( 时间 t ) zlabel( 数值解 u)%t=tf=2 时各位置之解%figure(3)M=length(t); % 取终点时间 xout=linspace(0,1,100); % 输出点位置 uout,dudx=pdeval(m,x,u(M,:),xout);plot(xout,uout); % 绘图title( 时间为 2 时 , 各处之解 ) xlabel(x)ylabel(u)%pde 函数文件%function c,f,s=ex6_1pd

23、efun(x,t,u,DuDx)c=pi2;f=DuDx;s=0;%起始条件档%function u0=ex6_1ic(x)u0=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;范例 6_2 试解以下之联立 pde 系统2u1 0.024 u21 F(u1 u2 )tx2u2 t 其中0.170u2F (u1 u2 )F (u1 u2 ) exp( 5.73(u1 u2) exp( 11.46(u1 u2 )且 0 x 1和 t 0 。此联立 pde

24、系统满足以下之条件限制式(i) 起始值条件u1 (x,0) 1u2 (x,0) 0(ii)边界值条件u1x(0,t)u2 (0,t) 0 u1 (1,t) 1u22 (1, t )解答：步骤 1： 改写 PDE 方程式为标准式1 *u10.024 xF (u1 u2 )1 ? t u2x 0.170 u2F (u1 u2)x因此1c10.0240.170u1xu2xF (u1 u2 ) sF (u1 u2 )和 m 0 。另外，左边界条件 ( x 0處 )。写成0u210?0.024u10.170 u2pl0u2ql同理，右边界条件10(x 1處 ) 为u1001?u10.024 1x0.17

25、0 u2x步骤 2：prqr编写u1001 pde函数文件function c,f,s=ex6_2pdefun(x,t,u,DuDx)c=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(x)u0=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： 取点由于此问题之端点均受边界条件之限制，

26、且时间 t 小时状 态之变动甚鉅 (由多次求解后之经验得知 )，故在两端点处之点 可稍微密集些。同时对于 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; 以上几个主要步骤编写完成后， 事实上就可直接完成主程序来 求解。此问题之参考程序如下：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.

27、9951;t=0 0.005 0.01 0.05 0.1 0.5 1 1.5 2;%利用 pde 求解% sol=pdepe(m,ex6_2pdefun,ex6_2ic,ex6_2bc,x,t);u1=sol(:,:,1); %第一个状态之数值解输出u2=sol(:,:,2); %第二个状态之数值解输出%绘图输出%figure(1)surf(x,t,u1)title(u1 之数值解 )xlabel(x)ylabel(t)%figure(2)surf(x,t,u2)title(u2 之数值解 )xlabel(x)ylabel(t)%pde 函数%function c,f,s=ex6_2pdefu

28、n(x,t,u,DuDx)c=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(x) u0=1 0;%边界值条件%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;执行结果6.2.2 PDE 图形接口工具箱有别于 pdepe指令所能处理之一维偏微分方程式， MATLAB 之 PDE 图形 接口工具箱， 可允许使用者处理某几类的二维空间之偏

29、微分方程式问题。 利用 此工具箱，使用者仅需要再命令六窗口中键入透过此图形接口其解法步骤大致上为：(1) 定义 PDE 问题，包括二维空间范围，边界条件以及 PDE 系数等(2) 产生离散化之点，并将原 PDE 方程式离散化。(3) 利用有限元素法 (finite element method;FEM)求解并显示答案 在详细说明此解法工具之前， 吾人将先介绍此工具所能处理之 PDE 问题形式 此 PDE 图形接口之选单下方，存在一排主要功能的图标 (icon) 按钮透过这些按钮，使用者可轻松地完成 PDE 方程式之求解。现并将这些按钮之 主要功能叙述如下：前五个按钮为 PDE 系统之边界范围绘

30、制功能，由左至右之用法为：以对角绘制矩形或正方形。按住鼠标左键可绘制矩形，而正 方形需以按住右键的方式绘制之。：从中心点至某一角边的方式绘制矩形或正方形。同样地，鼠 标左键绘矩形，右键绘正方形。：由周围界线的方式绘制椭圆或圆形区域。 鼠标左键用以绘制 椭圆，而右鉴用来绘制饼图形。：以中心点向外的方式绘制椭圆或圆。 同样地，鼠标左及又键，分别用以绘制椭圆及圆形的区域：用以绘制多边型等不规则区域， 欲关闭此功能需按鼠标右键。 在这些绘制按钮之后之按钮功能依序如下：用以给定边界条件。在此功能选定后，使用者可在任一图形 边界上按住鼠标左键两次，然后在对话框上将边界条件给 入。：用以指定 PDE 问题及

31、相关参数 ：产生图形区域内离散化之网点 ：用以进一步将离散化之网点再取密一点 (refine mesh)。：在指定 PDE 系统，边界条件及圆形区域后，按此钮即开始 解题。：用以指定显示结果绘制方式。：放大缩小之功能，便于图形绘制及显示 PDE 问题形式pdetool图形接口工具所能处理之 PDE 问题的基本型式为以下之椭圆 型方程式(c u) au f其中 xiyj 为 Laplacian 运操作数 (operator)。c，a，f 和 u 为定义在有限平面几下之纯量或复变函数值函数 (complex valued function)。除 了这个基本问题外，此工具箱亦可解以下的

32、 PDE 问题。拋物线型dut 双曲线型2d t22u 特征值问题 (c(c u)(c u)u) auauaudu其中 为未知特征值 以上这些问题之边界条件型式可为Dirichlet 条件hu r ,Generalized Neumann条件n (c u) qu g其中 n 为向外之单位法向量 (Outward unit normal)，g，q，h 和 r 为复数值函数。 如何利用 pdetool 接口解 PDE 问题要利用 pdetool接口求解之前，需先定义 PDE 问题，其包含三大部份： (1) 利用绘图 (draw)模式，定义欲解问题的空间范围 (domain)(2) 利

33、用 boundary 模式，指定边界条件(3) 利用 PDE 模式，指定 PDE 系数，即输入 c，a，f 和 d 等 PDE 模 式中之系数在定义 PDE 问题之后，可依以下两个步骤求解(1) 在 mesh 模式下，产生 mesh点，以便将原问题离散化(2) 在 solve 模式下，求解(3) 最后，在 Plot 模式下，显示答案现在吾人将以 Poissons 方程式 u f 之求解为范例，详细说明 pdetool 之用法。此问题之几何图形及相关边界条件，将于求解过程中加以说明。步骤 1：在命令列窗口中键入 pdetool以进入 GUI(graphical user interface)界面

34、。选取 Options 中之 Grid 功能，以显示网格线。步骤 2：利用 Draw 功能，画出问题之几何图形。请注意：使用者可 利用内定对象 ”多边型”，”矩形”，”正方形 ”，”圆形”，及 椭圆型 ”，予以组合，例如(i)先选取 ”矩形/正方形 ”对象，移动由标至所欲输入左上角之点，如坐标 (-1,0)之点，按住鼠标左键，往右下角拉至坐标为 (1,-0.4) 处，即行成代号为 R1 之矩形。其余图形 C1，R1 和 C2 可选取 适当之对象制作之，以形成如下图之图形区域。以代数公式而 言，其为 R1+C1+R2+C2值得注意的是，圆形之区域需以按住鼠标右键的方式来制作 (非左 键)。同时，

35、如欲进一步修改各图形对象之大小及位置数据，可在 该图上按鼠标左键两次，然后在对象对话框上输入数据。(ii)若所欲形成之图形区域，需将 C2 去除，则可在公式列中直接输 入 R1+C1+R2-C2 即可步骤 3：选取 PDE 功能项，以输入 PDE方程式的系数及类型。因问题 为 u f ，故此为椭圆型的问题，且其标准形式(c u) au 入之情况如下：步骤 5：选取 Mesh 功能，产生网点。使用者亦可进一步利用将网f 比较得知， c=1，a=0 和 f=10 ，所以对话框输步骤 4：选取 Boundary 功能，以输入边界条件。假设边界条件为 Neumann 形，且为 u n 5 。其中在弧形

36、部份与标准式知， g=-5 且 q=0。 但直线部份其边界条件则在 Dirichlet type 使 h=0，r=0。故而点取在密一点 (refine mesh)。步骤 6：选取 solve 功能，解此 PDE请注意： 1.MATLAB 会以图形的方式展示结果，使用者亦可点选 plot 下 之 ”parameter功s ”能，选择适当之结果显示图形及数据。例如，2.另外，若使用者欲将结果输出至命令列窗口中，以供后续处理， 可利用 solve 功能项下之 ”export solution指定”变量名称来完成6.3 化工实例演练范例 6.3.1 触煤反应装置内温度及转换率的分布 (吕， p.235

37、)以外部热交换式的管形固定层触煤反应装置， 进行苯加氢反应产生环己烷。此反应系统之质量平衡及热平衡方程式如下：Tke2T1TrA B (Hr ) 0LGCp2 rrrGCpfDe2 f 1frA BM ar0Lur 2 rrGy0其中 T 为温度 ()， f 为反应率， L 为轴向距离， r 为径向距离。此系统之 边界条件为L 0 , T T0(r) , f f 0 (r)r0r rwke Tr rhw (T Tw )此外，式中之相关数据及操作条件如下：(i)反应速率式rA3H 2氫環己 烷kKH3 K B PH3 PB(1 KH PH KBPB KCPC )4其中 P 表示分压 (atm)，

38、而速率参数为ln K 12100 RT 32.3 Rln K H 15500 RT 31.9 Rln K B 11200 RT 23.1 Rln K C 8900 RT 19.4 R上式中，下标 B，H 及 C 分别代表苯，氢及环己烷。 R 为理想气体常数(1.987cal/molK)(ii)操作条件及物性数据总压Pt 1.25atm反应管管径rw 2.5cm壁温100质量速度 苯对氢之莫耳流量比G 631kg m2 hr m 30反应管入口的苯之莫耳分率y0 0.0323反应气体之平均分子量触煤层密度流体平均比热 反应热整体传热系数进料温度 反应管管长 流速有效热传导系数壁境膜传热系数有效扩

39、散系数题意解析：因反应速率式 rA 与分压有关 一步将 rA 由反应率 f 表示之 应式1-fM av 4.47PB 1200 kg m3CP 1. 74 kcal kg molH r 49250 kcal kg mol h0 65.8 kcal m2 hr 0 C T(0) 125 L 45cmu 8.03m hrKe 0.65kcal m hr 0Chw 112 kcal m2 hr 0CDe 0.755 m 2 hrf 有关。故需进基于以下之反而分压又与反应率方能求解偏微分方程式。3H 2-3f1-f m-3f则各分压与总压之关系为PHP m 3f t 1 m 3fPBP m f t 1

40、 m 3fPC PtC t 1 m 3f将上式 ( )( )，连同反应速率式 ( )，带入平衡方程式 ( )中，配 合边)之标准式kerGCpDe r urA B ( Hr )GCprA BM avGy0界条件 ( )，可利用 pdepe求解之。MATLAB 程序设计将原方程式改写成如 (T0L1fL因此GCprDefurrA B (Hr)GCprA BMavkeTs和mGy01(圆柱)。另外，左边界条件(r 0 处 ) 写成11?plql同理右边界条件 (rrw )可写成hw (T Tw)0GCp *1Tw)hw(Tpr w 0GCpqr 0 pMATLAB 程序求解此 PDE 问题，其参考

41、根据以上的分析， 吾人可编写一 程序如下：ex6_3_1.mfunction ex6_3_1 % 范例 6_3_1 触媒反应器内温度及转化率的分布%global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw De % 给定数据%Pt=1.25; % 总压 (atm)rw=0.025; % 管径 (m)Tw=100+273; % 壁温 ( )G=631; % 质量流率 (kg/m2hr)M=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

42、.03;R=1.987;ke=0.65;hw=112;De=0.755;%m=1;% 取点%r=linspace(0,rw,10);L=linspace(0,Lw,10);% 利用 pdepe 求解%sol=pdepe(m,ex6_3_1pdefun,ex6_3_1ic,ex6_3_1bc,r,L);T=sol(:,:,1); %温度f=sol(:,:,2); %反应率% 绘图输出%figure(1)surf(L,r,T-273)title(temp)xlabel(L)ylabel(r)zlabel(temp (0C)%figure(2)surf(L,r,f)title(reaction ra

43、te)xlabel(L)ylabel(r)zlabel(reaction rate)% PDE 函数%function c1,f1,s1=ex6_3_1pdefun(r,L,u1,DuDr)global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw DeT=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

44、;pb=Pt*(1-f)/a;pc=Pt*f/a;rA=k*Kh3*Kb*ph3*pb/(1+Kh*ph+Kb*pb+Kc*pc)4; %c1=1 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*Cp)rA*rho_B*Mav/(G*y0);% 起始值条件%function u0=ex6_3_1ic(x)u0=125+273 0;% 边界条件档%function pl,ql,pr,qr=ex6_3_1bc(rl,ul,rr,ur

45、,L)global Pt rw Tw G M y0 Mav rho_B Cp dHr h0 u R ke hw Depl=0 0;ql=1 1;pr=hw*(ur(1)-Tw) 0;qr=G*Cp 1;结果：范例 6_3_2 扩散系统之浓度分布 (C&M,p.403) 参考如图之装置。管中储放静止液体 B，高度为 L=10 ，其暴置于充满A 气体的环境中。假设与 B液体接触面之浓度为 CA0 0.01mol m3 ，且此浓度s 。试决定以下两种情况下，气体 A 溶于液不随时间改变而改变，即在操作时间内 (h 10天)维持定值。气体 A 在液体 B 中之扩散系数为 DAB 2 10 9 m2 体

46、 B 中之流通量 (flux) 。(a)A 与 B 不管发生反应kCA1。(b)A 与 B 发生以下之反应 A B C ,rA其反应速率常数 k 2 10 7 s氣體Az 10液體Bcm图 气体 A 在液体 B 中之扩散题意解析：故其扩散现象之质量平衡式如下：(a)因气体 A 与液体 B 不发生反应，C A2CAt AB z2依题意，其起始及边界条件为0CA0I.C. C A(z,0)B.C. CA (0,t)CAzA z L(b)在气体 A 与液体 B 会发生一次反应之情况下，其质量平衡式需改写DAB zC2A kC AzCA t 而起始及边界条件同上。 在获得浓度分布后，即可以 Fick s lawCAN Az(t) DAB Az0z计算流通量。MATLAB 程序设计：此问题依旧可以 pdepe迅速求解。现并就各状况之处理过程简述之 (a)将式( )与标准式 ( )比较，可得 C 1，f DAB CA z，s 0， 和 m 0 。另外，经与式 ( )比较后得知，左边界及右