偏微分方程及其求解实例ppt课件.ppt

当B2-4AC0时，方程为椭圆型（elliptic）PDE当B2-4AC=0时，方程为抛物型（parabolic）PDE当B2-4AC0时，方程为双曲型（hyperbolic）PDE,第5章 化工中的偏微分方程及其求解,(1) 导热方程:,(2) 拉普拉斯方程: 如稳态静电场和稳态温度分布模型,(3) 波动方程: 一维弦振动模型,偏微分方程的边界条件,(1) Dirichlet边界条件(第一类边界条件)-在边界上给定u(x,y)的值(2) Neumann边界条件(第二类边界条件)-在边界上给定u(x,y)的外法向导数值(3)Robin边界条件(混合边界条件, 第三类边界条件)-在边界上给定u(x,y)与其外法向导数值的线形组合的值,即,偏微分方程数值求解方法,(1) 有限差分法(finite difference method) 将求解域内划分一定数目的网格(交叉点为网格点或节点),在所有节点上,PDE方程用有限差分近似代替，得到代数方程组。这种方法数值稳定性好，但是不能获得节点之间的解，对不规则几何形状或导热系数突变的情况不适合。,偏微分方程数值求解方法,(2) 正交配置法(method of weighted residuals) 首先选择一个多项式作为试函数，将此试函数代入微分方程，求出多项式的根作为配置点，令在各配置点试函数代入微分方程后的残差（或余量）为零，得到关于多项式系数的代数方程组，然后求解此方程组得到多项式中各项的系数，得到的多项式即为微分方程的近似解析解。,偏微分方程数值求解方法,(3) MOL法(method of lines) 将一个自变量当成连续变量，对其余的自变量用有限差分法或者正交配置法进行离散，从而把偏微分方程转变为常微分方程组，然后用龙格库塔法积分求解。常用于求解一维动态和二维稳态PDE方程。,偏微分方程数值求解方法,(4) 有限元法(finite element methods,FEM) 把求解域先划分为大量的单元（elements），其中任意大小和方向的三角形网格尤其适用于二维的情况。三角形顶点称为节点（nodes），并与相邻单元相连接。将PDE离散为代数方程组求解。优点是易于处理复杂几何区域，易于与各种边界条件组合使用，可同时提供节点和整个求解域内的解。,有限差分法求解偏微分方程,差分格式：,有限差分法求解偏微分方程,差分格式：,有限差分法求解偏微分方程,差分格式：,1、显式差分法2、隐式差分法,有限差分法求解偏微分方程一维动态PDE模型的求解,初始条件： t=0和0 x10，T=0 边界条件： x=0 cm和所有t，T=100 x=10cm和所有t，T=0 ,j,t,x,i-1,i,i+1,j+1,j-1,显式差分格式,function PDE1_FDM_im% 显式差分法求解一维热传导方程clear all;clc n=6;% 空间节点数m=8;% 时间节点数T(1,1:m)=100;T(n,1:m)=0; % 边界条件T(1:n,1)=0; % 初始条件alpha=2; % 热扩散系数,m/sa=10; % for x坐标长度,cmb=8; % for 时间长度,sh=a/(n-1); % 空间步长k=b/(m-1); % 时间步长r=alpha.*k./h.2;,for j=2:m % for time for i=2:n-1 % for x T(i,j+1)=(1-2.*r).*T(i,j)+r.*(T(i+1,j)+T(i-1,j); endendT8=T(1:n,m),j,t,x,i-1,i,i+1,j+1,j-1,Crank-Nicolson差分格式,i=1,i=2,i=3,偏微分方程的求解实例1: 恒靠近速度时两等直径液滴形成的液膜内流 体排液速率的模拟问题描述(h(r,t)t ?),靠近方向,排液方向,排液方向,偏微分方程的求解实例1: 恒靠近速度时两等直径液滴形成的液膜内流 体排液速率的模拟模型建立,初始液膜厚度: h(r,t=0)=h0c+r2/rb边界条件: r=ra, V=2e-6 m/s,h0c=2.8e-6mrb=0.0015m=0.03N/mc=0.3Pas,偏微分方程的求解实例1: 恒靠近速度时两等直径液滴形成的液膜内流 体排液速率的模拟方程离散,r,i=1,i=n:P=0,function CVFDclear all;clc; format long eh0c=2.8e-6; % 初始液膜中心厚度, mmu=0.3; % 液相流体的粘度, Pa.srb=0.0015; % 液滴半径,mtheta=0.03; % 液膜的表面张力, N/mra=0.5*rb; % 膜出口压力为零区域V=2e-6; % 膜出口ra处靠近速度m/sm=100; r=0,0,linspace(0,ra,m); n=m+2;dr=ra./(m-1); % 节点的步长,r,m h0=h0c+r.2./rb; % 初始膜厚度mh0(1)=h0(5);h0(2)=h0(4);,hold onk=10; tf=12;for i=1:k t,h=ode15s(Non,0,tf,h0,mu,theta,rb,V,dr,n) k=length(t); % 不同时刻的液膜厚度 plot(r(3:n)./ra,h(k,3:n)./h0c) h0=h(k,:);endhold off function dhdt=Non(t,h,mu,theta,rb,V,dr,n)h(n+1)=(4.*dr.2./rb+2.*h(n)-h(n-1).*(1-1./(2.*(n-3)./(1./(2.*(n-3)+1);dhdt(3)=-theta./9./mu.*(h(3).3.*(2.*h(5)-8.*h(4)+6.*h(3)./dr.4);h(1)=h(5);h(2)=h(4);,for i=4:n-1Dhr(i)=(h(i+1)-h(i-1)/(2*dr);D2hr(i)=(h(i+1)-2*h(i)+h(i-1)/(dr2);D3hr(i)=(h(i+2)-2.*h(i+1)+2.*h(i-1)- h(i-2)./(2.*dr.3);D4hr(i)=(h(i+2)-4.*h(i+1)+6.*h(i)- 4.*h(i-1)+h(i-2)./(dr.4);H1(i)=D4hr(i)+2./(i-1).*dr).*D3hr(i)- 1./(i-3).2.*dr.2).*D2hr(i)+ 1./(i-3).3.*dr.3).*Dhr(i);H2(i)=Dhr(i).*(D3hr(i)+1./(i-3).*dr).* D2hr(i)-1./(i-3).2.*dr.2).*Dhr(i);dhdt(i)=-theta./8./mu.*(h(i).3.*H1(i)./3+ h(i).2.*H2(i);enddhdt(n)=-V;dhdt=dhdt(:);,% 不同时刻液膜表面的压力分布 p(k,3)=1-rb.*(h(k,4)-h(k,3)./dr.2; hn1=(4./rb.*dr.2+2.*h(k,n)-h(k,n-1).*(1- 1./(2.*(n-3)./(1./(2.*(n-3)+1); for j=4:n-1 p(k,j)=1-rb./4.*(h(k,j+1)-h(k,j-1)./(j- 3).*2.*dr.2)+(h(k,j+1)-2.*h(k,j)+ h(k,j-1)./dr.2); p(k,n)=1-rb./4.*(hn1-h(k,n-1)./(n-3).*2.*dr.2)+ (hn1-2.*h(k,n)+h(k,n-1)./dr.2); end plot(r(3:n)./ra,p(k,3:n).*theta.*2./rb),function PDE1Dd_CrankNicolson% 使用Crank-Nicolson有限差分方法求解一维动态传热模型 c1 = 100;c2 = 0;a = 10;b = 8;alpha = 2;n = 6;m = 8; U = CrankNicolson(ic,c1,c2,a,b,alpha,n,m) % -function f = ic(x)f = 0;,偏微分方程的求解实例2:,function U = CrankNicolson(f,c1,c2,a,b,alpha,n,m)% 使用Crank-Nicolson有限差分方法求解一维动态传热模型% Initialize parameters and Uh = a/(n-1);k = b/(m-1);r = alpha*k/h2;s1 = 2+2/r;s2 = 2/r-2;U = zeros(n,m);% Boundary conditionsU(1,1:m) = c1;U(n,1:m) = c2; % Generate first row (that is, the first column in matrix U)U(2:n-1,1) = feval(f, h:h:(n-2)*h);,% Form the diagonal and off-diagonal elements of A and% the constant vector B and solve tridiagonal system AX=BVd(1,1:n) = s1*ones(1,n);Vd(1) = 1;Vd(n) = 1;Va = -ones(1,n-1);Va(n-1) = 0;Vc = -ones(1,n-1);Vc(1) = 0;Vb(1) = c1;Vb(n) = c2;for j=2:m for i=2:n-1 Vb(i) = U(i-1,j-1)+U(i+1,j-1)+s2*U(i,j-1); end X = TDMA(Va,Vd,Vc,Vb); U(1:n,j)=X;endU = U,function X = TDMA(A,D,C,B)% Tridiagonal Matrix Algorithm（三对角阵算法）% Function: to solve the Tridiagonal system CX=B, where C is a Tridiagonal matrix% Input - A is the subdiagonal of the coefficient matrix% - D is the main diagonal of the coefficient matrix% - C is the superdiagonal of the coefficient matrix % - B is the constant vector of the linear system% Output - X is the sulution vectorN = length(B);for k=2:N mult = A(k-1)/D(k-1); D(k) = D(k)-mult*C(k-1); B(k) = B(k)-mult*B(k-1);endX(N) = B(N)/D(N);for k=N-1:-1:1 X(k) = (B(k)-C(k)*X(k+1)/D(k);end,