偏微分方程式之求解.doc

第六章 偏微分方程式之求解在化工的领域中，有不少程序之动态是由以偏微分方程式(Partial differential equation；PDE)所描述的，例如热与质量在空间中的传递等。这些用以描述实际问题的PDE，除非具有某些特定的方程式型态及条件，否则甚难以手算的方式找出解析解。而在数值求解方面，最常被采用的方法为有限差分法(finite difference)何有限元素法(finite element)。然对于某些不熟悉数值分析及程序编写的化工人而言，欲充分了解以偏微分方程式所描述之系统动态是相当不容易的，更遑论进一步的设计与分析了。 值得庆幸的是，MATLAB的环境中提供了一个求解PDE问题的工具箱，让使用者得以利用简单的指令或图形接口工具输入欲解的PDE，并求解。使得PDE之数值解在弹指之间完成，使用者不在为数值法所苦恼，轻松掌握偏微分方程式系统的动态，并可进一步进行后续之设计工作。 本章将以循渐进的方式，介绍PDE工具箱及其用法，并以数个典型的化工范例进行 示范，期能使初学者很快熟悉PDE工具箱，并使用它来设计与分析以偏微方方程式所描述的程序系统。6.1 偏微分方程式之分类偏微分方程式可根据其阶数(order)，线性或非线性型态，以及边界条件进行分类。6.1.1依阶数的分类 偏微分方程式是以偏微分项中之最高次偏微分来定义其阶数，例如： 一阶偏微分方程式： 二阶偏微分方程式： 三阶偏微分方程式： 6.1.2 依非线性程度分类 偏微分方程式亦可以其线性或非线性情况，区分为线性(linear)，似线性(quasilinear)，以及非线性三类。例如，以下之二阶偏微分方程式(Constantinides and Mostoufi,1999) 可依系数之情况，进行如下表之归类类别 情况线性 系数为定值，或仅为(x,y)函数似线性 系数为依变数(dependent variable)u或其比方程式中之偏微分低阶之偏微分项的函数，如非线性 系数中，具有与原方程式之偏微分同阶数之变数，如 另外，对于线性二阶偏微分方程式，可进一步将其分类为椭圆型(elliptic)，拋物线型(parabolic)，以及双曲线型(hyperbolic)。具体上来说，此类偏微分方程式二阶线性之一般式为 系数是定值或为u的函数。若g=0，则上式为其次是偏微分方程式。式子( )之分类及代表性例子，请见下表方程式类别 判断式 代表性范例椭圆型 Laplace方程式，拋物线型 Poisson方程式， 热传导或扩散方程式 双曲线型 波动方程式 注：二维系统之运操作数之定义为6.1.3 起始条件和边界条件的分类 为了能获得偏微分方程式之解答，其起始条件和边界条件可依其特性区分为三类。现以一维之动态热传递方程式(拋物线型偏微分方程式) 为例，进一步说明如何区分这些边界条件及起始条件(Constantinides and Mostoufi ,1999)。 (i) 第一类：Dirichlet Condiction 若依变量(T)本身，在某个独立变量值时，被指定，则此条件称为Dirichlet Condiction，亦称为essential边界条件。下图为一典型的Dirichlet条件示意图图6.1 平板Dirichlet Condiction示意图 由图中很清楚的显示，该平板之边界条件为 此边界条件依定义，即为Dirichlet Condiction。同时，若再起始时，各处之温度分布可以位置之函数表示，即 此亦属Dirichlet型之边界条件。 (ii) 第二类：Neumann condition Neumann condition系指依变量之变化率之边界条件为定值，抑或独立变量之函数之情况。例如 或 Neumann型边界条件，亦称为natural boundary condition。 (iii) 第三类：Robbins condition若依变量之变化率之边界条件，为自身之函数(非独立变量之函数)时，被称为Robbins condition。例如， 上式之边界条件，当发生在固液相间之传递上，亦即热流通量(heat flux)正比于固液两端之温差，其示意图如下： (iv) Cauchy condition Cauchy condition系指问题中同时存在Dirichlet和Neumann边界条件。例如下图 即 “Dirichlet condition” “Neumann condition” 6.2 MATLAB PDE工具箱6.2.1 MATLAB PDE解答器 MATLAB提供了一个指令pdepe，用以解以下的PDE方程式 其中时间介于之间，而位置则介于有限区域之间。m值表示问题之对称性，其可为0，1或2，以表示平板(slab)，圆柱(cylindrical)或球体(spherical)之对秤情况。因而，如果如果m0，则a必等于b，也就是说其具有圆柱或球体之对称关系。同时，式中一项为流通量(flux)，而为来源(source)项。为偏微分方程式之对角线系数矩阵。若某一对角线元素为0，则表示该相应偏微分方程式为椭圆型偏微分方程式，若为正值(不为0)，则为拋物线型偏微分方程式。请注意c之对角线元素必不全为0。偏微分方程式之起始值可表示为 而边界条件之形式为 其中x为两端点位置，即a或b。 用以解含上述起始值及边界值条件的偏微分方程式( )之指令pdepe的用法如下： 其中 为问题之对称参数 为独立变量x之网取点(mesh)位置向量，即，其中(起点)，(终点)。 为独立变量t(时间)之向量，即，其中为起始时间，为终点时间。 为使用者提供之pde函数文件。其档案之格式如下： 亦即，使用者仅需提供偏微分方程式中之系数向量。，和均为行(column)向量，而向量即为矩阵之对角线元素。 提供解u之起始值，其格式为值得注意的是u为行向量。 使用者提供之边界条件档，格式如下：其中ul和ur分别表示左边界和右边界之u 的近似解。输出变量中，pl和ql分别表示左边界p和q之行向量，而pr和qr则为右边界p和q之行向量。 为解答输出。为多维的输出向量，为的输出，即。另，元素表示在和时之答案。 为解答器之相关解法参数。详细请见odeset之使用法。 注：1. MATLAB PDE解答器pdepe之算法，主要事将原一组的椭圆型和拋物线型偏微分方程式转化为一组常微分方程式。此转换的过程系基于使用者所指定之mesh点，以二阶空间离散化(spatial discretization)技术为之(Keel and Berzins,1990)，然后以ode15s的指令求解。采用ode15s的ode解法，主要是因为在离散化的过程中，椭圆型偏微分方程式为被转为一组代数式，而拋物线型的偏微分方程式则被转为一组联立的微分方程式。因而，原偏微分方程式被离散化后，遂变成一组同时伴有微分方程式与代数式之微分代数方程式系统，故以ode15s使可顺利求解。 2. x之取点(mesh)位置影响解的精确度甚鉅，若pdepe解答器回应出has difficulty finding consistent initial considition之讯息时，使用者可进一步将mesh点取密一点，即增加mesh点数。另外，若状态u再某些特定点上有较快速之变动时，亦需将此处之点取密集些，以增加精确度。值得注意的是pdepe并不会自动做并不会自动做xmesh的自动取点，使用者必须观察解的特性，自行作曲点的动作。一般而言，所取之点数至少需大于3以上。3. tspan之选取主要是基于使用者对那些特定时间之状态有兴趣而选定。而间距(step size)之控制由程序自动完成。4. 若要获得特定位置及时间下之解，可配合以pdeval的指令。跟pdepe指令之格式如下： 其中 代表问题之对称性。=0，代表平板；=1，圆柱体；=2表示球体。其意义同pdepe中之自变量。 使用者在pdepe中所指定之输出点之位置向量。 即。也就是说其为pdepe输出中第i个输出变数在各点位置处，时间固定为下之解。 为所欲内插输出点位置向量。此为使用者重新指定之位置向量。 为基于所指定位置，固定时间下之相对应输出。 为相对应之输出值。 ref. Keel,R.D. and M. Berzins,”A Method for the Spatial Discritization of Parabolic Equations in One Space Variable”,SIAM J. Sci. and Sat. Comput.,Vol.11,pp.1-32,1990. 以下吾人将以数个范例，详细说明pdepe之用法。范例6.1 试解以下之偏微分方程式 其中，且满足以下之条件限制式(i)起始值条件IC：(ii)边界条件 BC1： BC2：注：本问题之解析解为 解答： 以下吾人将分述求解之步骤与过程。当以下之各步骤完成后，可进一步将其汇整为一主程序ex6_1.m，然后求解之。 步骤1 将欲解之偏微分方程式改写成如式子( )之标准式。 此即 和 步骤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 编写边界条件。在编写之前，先将边界条件改写成标准形式，如式( )，找出相对之和的函数，然后写成边界条件档，例如，原边界条件可写成 BC1： BC2： 即 故而，边界条件档可编写成 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(:,:,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时,各处之解)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(1)surf(x,t,u)title(pde数值解)xlabel(位置x)ylabel(时间t )zlabel(数值解u)%与解析解做比较%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_1pdefun(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系统 其中 且和。此联立pde系统满足以下之条件限制式(i)起始值条件 (ii)边界值条件 解答： 步骤1： 改写PDE方程式为标准式因此 和。另外，左边界条件()。写成 即 同理，右边界条件()为 即 步骤2： 编写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： 取点 由于此问题之端点均受边界条件之限制，且时间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.995 1;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.995 1;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_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;%起始值条件%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图形接口工具箱，可允许使用者处理某几类的二维空间之偏微分方程式问题。利用此工具箱，使用者仅需要再命令六窗口中键入 透过此图形接口其解法步骤大致上为：(1) 定义PDE问题，包括二维空间范围，边界条件以及PDE系数等(2) 产生离散化之点，并将原PDE方程式离散化。(3) 利用有限元素法(finite element method;FEM)求解并显示答案 在详细说明此解法工具之前，吾人将先介绍此工具所能处理之PDE问题形式。 此PDE图形接口之选单下方，存在一排主要功能的图标(icon)按钮 透过这些按钮，使用者可轻松地完成PDE方程式之求解。现并将这些按钮之主要功能叙述如下： 前五个按钮为PDE系统之边界范围绘制功能，由左至右之用法为：以对角绘制矩形或正方形。按住鼠标左键可绘制矩形，而正方形需以按住右键的方式绘制之。：从中心点至某一角边的方式绘制矩形或正方形。同样地，鼠标左键绘矩形，右键绘正方形。：由周围界线的方式绘制椭圆或圆形区域。鼠标左键用以绘制椭圆，而右鉴用来绘制饼图形。：以中心点向外的方式绘制椭圆或圆。同样地，鼠标左及又键，分别用以绘制椭圆及圆形的区域。：用以绘制多边型等不规则区域，欲关闭此功能需按鼠标右键。 在这些绘制按钮之后之按钮功能依序如下：用以给定边界条件。在此功能选定后，使用者可在任一图形边界上按住鼠标左键两次，然后在对话框上将边界条件给入。：用以指定PDE问题及相关参数。：产生图形区域内离散化之网点。：用以进一步将离散化之网点再取密一点(refine mesh)。：在指定PDE系统，边界条件及圆形区域后，按此钮即开始解题。：用以指定显示结果绘制方式。：放大缩小之功能，便于图形绘制及显示。 6.2.2.1 PDE问题形式pdetool图形接口工具所能处理之PDE问题的基本型式为以下之椭圆型方程式 其中为Laplacian运操作数(operator)。c，a，f和u 为定义在有限平面几下之纯量或复变函数值函数(complex valued function)。除了这个基本问题外，此工具箱亦可解以下的PDE问题。 拋物线型 双曲线型 特征值问题 其中为未知特征值 以上这些问题之边界条件型式可为 Dirichlet条件 Generalized Neumann条件 其中为向外之单位法向量(Outward unit normal)，g，q，h和r为复数值函数。 6.2.2.2 如何利用pdetool接口解PDE问题 要利用pdetool接口求解之前，需先定义PDE问题，其包含三大部份：(1) 利用绘图(draw)模式，定义欲解问题的空间范围(domain)(2) 利用boundary模式，指定边界条件(3) 利用PDE模式，指定PDE系数，即输入c，a，f和d等PDE模式中之系数在定义PDE问题之后，可依以下两个步骤求解(1) 在mesh模式下，产生mesh点，以便将原问题离散化(2) 在solve模式下，求解(3) 最后，在Plot模式下，显示答案 现在吾人将以方程式之求解为范例，详细说明pdetool之用法。此问题之几何图形及相关边界条件，将于求解过程中加以说明。步骤1：在命令列窗口中键入pdetool以进入GUI(graphical user interface)界面。选取Options中之Grid功能，以显示网格线。步骤2：利用Draw功能，画出问题之几何图形。请注意：使用者可利用内定对象”多边型”，”矩形”，”正方形”，”圆形”，及”椭圆型”，予以组合，例如(i)先选取”矩形/正方形”对象，移动由标至所欲输入左上角之点，如坐标(-1,0)之点，按住鼠标左键，往右下角拉至坐标为(1,-0.4)处，即行成代号为R1之矩形。其余图形C1，R1和C2可选取适当之对象制作之，以形成如下图之图形区域。以代数公式而言，其为R1+C1+R2+C2值得注意的是，圆形之区域需以按住鼠标右键的方式来制作(非左键)。同时，如欲进一步修改各图形对象之大小及位置数据，可在该图上按鼠标左键两次，然后在对象对话框上输入数据。 (ii)若所欲形成之图形区域，需将C2去除，则可在公式列中直接输入R1+C1+R2-C2即可步骤3：选取PDE功能项，以输入PDE方程式的系数及类型。因问题为，故此为椭圆型的问题，且其标准形式比较得知，c=1，a=0和f=10，所以对话框输入之情况如下：步骤4：选取Boundary功能，以输入边界条件。假设边界条件为Neumann形，且为。其中在弧形部份与标准式知，g=-5且q=0。但直线部份其边界条件则在Dirichlet type使h=0，r=0。故而步骤5：选取Mesh功能，产生网点。使用者亦可进一步利用将网点取在密一点(refine mesh)。步骤6：选取solve功能，解此PDE。 请注意：1.MATLAB会以图形的方式展示结果，使用者亦可点选plot下之”parameters”功能，选择适当之结果显示图形及数据。例如， 2.另外，若使用者欲将结果输出至命令列窗口中，以供后续处理，可利用solve功能项下之”export solution”指定变量名称来完成。6.3 化工实例演练 范例6.3.1 触煤反应装置内温度及转换率的分布(吕，p.235) 以外部热交换式的管形固定层触煤反应装置，进行苯加氢反应产生环己烷。此反应系统之质量平衡及热平衡方程式如下： 其中为温度()，为反应率，为轴向距离，为径向距离。此系统之边界条件为 , , , , , 此外，式中之相关数据及操作条件如下： (i)反应速率式 其中表示分压(atm)，而速率参数为 上式中，下标B，H及C分别代表苯，氢及环己烷。R为理想气体常数(1.987cal/molK) (ii)操作条件及物性数据 总压 反应管管径壁温 质量速度 苯对氢之莫耳流量比 反应管入口的苯之莫耳分率 反应气体之平均分子量 触煤层密度 流体平均比热 反应热 整体传热系数 进料温度 反应管管长 流速 有效热传导系数 壁境膜传热系数 有效扩散系数 题意解析：因反应速率式与分压有关，而分压又与反应率有关。故需进一步将由反应率表示之，方能求解偏微分方程式。基于以下之反应式 则各分压与总压之关系为 将上式( )( )，连同反应速率式( )，带入平衡方程式( )中，配合边界条件( )，可利用pdepe求解之。 MATLAB程序设计 将原方程式改写成如( )之标准式 因此 和 和(圆柱)。另外，左边界条件(处)写成 即 ， 同理右边界条件()可写成 即 根据以上的分析，吾人可编写一MATLAB程序求解此PDE问题，其参考程序如下： 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.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 rate)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;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,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液体接触面之浓度为，且此浓度不随时间改变而改变，即在操作时间内(天)维持定值。气体A在液体B中之扩散系数为。试决定以下两种情况下，气体A溶于液体B中之流通量(flux)。(a)A与B不管发生反应 (b)A与B发生以下之反应 , 其反应速率常数。图 气体A在液体B中之扩散 题意解析：(a)因气体A与液体B不发生反应，故其扩散现象之质量平衡式如下： 依题意，其起始及边界条件为 I.C. , B.C. , , (b)在气体A与液体B会发生一次反应之情况下，其质量平衡式需改写为 而起始及边界条件同上。在获得浓度分布后，即可以Ficks law 计算流通量。 MATLAB程序设计： 此问题依旧可以pdepe迅速求解。现并就各状况之处理过程简述之(a)将式( )与标准式( )比较，可得，和。另外，经与式( )比较后得知，左边界及右边界条件之系数分别为 左边界() 右边界() (b)将式( )与标准式( )比较，可得，和。而边界条件之系数同(a)之结果。 利用以上之处理结果，吾人可编写一MA