4-偏微分方程及其求解实例

当B2-4AC0时，方程为椭圆型（elliptic）PDE当B2-4AC=0时，方程为抛物型（parabolic）PDE当B2-4AC0时，方程为双曲型（hyperbolic）PDE第5章化工中的偏微分方程及其求解22222,,,,uuuuuuuABCDEFufxyuxxyyxyxy22222,,,,uuuuuuuABCDEFufxyuxxyyxyxy22uutx22220uuxy22222uutx(1)导热方程:(2)拉普拉斯方程:如稳态静电场和稳态温度分布模型(3)波动方程:一维弦振动模型偏微分方程的边界条件(1)Dirichlet边界条件(第一类边界条件)---在边界上给定u(x,y)的值(2)Neumann边界条件(第二类边界条件)---在边界上给定u(x,y)的外法向导数值(3)Robin边界条件(混合边界条件,第三类边界条件)---在边界上给定u(x,y)与其外法向导数值的线形组合的值,即,uxy,uxyn,uquxyn偏微分方程数值求解方法(1)有限差分法(finitedifferencemethod)将求解域内划分一定数目的网格(交叉点为网格点或节点),在所有节点上,PDE方程用有限差分近似代替，得到代数方程组。这种方法数值稳定性好，但是不能获得节点之间的解，对不规则几何形状或导热系数突变的情况不适合。偏微分方程数值求解方法(2)正交配置法(methodofweightedresiduals)首先选择一个多项式作为试函数，将此试函数代入微分方程，求出多项式的根作为配置点，令在各配置点试函数代入微分方程后的残差（或余量）为零，得到关于多项式系数的代数方程组，然后求解此方程组得到多项式中各项的系数，得到的多项式即为微分方程的近似解析解。偏微分方程数值求解方法(3)MOL法(methodoflines)将一个自变量当成连续变量，对其余的自变量用有限差分法或者正交配置法进行离散，从而把偏微分方程转变为常微分方程组，然后用龙格－库塔法积分求解。常用于求解一维动态和二维稳态PDE方程。偏微分方程数值求解方法(4)有限元法(finiteelementmethods,FEM)把求解域先划分为大量的单元（elements），其中任意大小和方向的三角形网格尤其适用于二维的情况。三角形顶点称为节点（nodes），并与相邻单元相连接。将PDE离散为代数方程组求解。优点是易于处理复杂几何区域，易于与各种边界条件组合使用，可同时提供节点和整个求解域内的解。有限差分法求解偏微分方程差分格式：有限差分法求解偏微分方程差分格式：有限差分法求解偏微分方程差分格式：11,,12kkijijuuutt1、显式差分法2、隐式差分法有限差分法求解偏微分方程一维动态PDE模型的求解22TTtx初始条件：t=0和0≤x≤10，T=0℃边界条件：x=0cm和所有t，T=100℃x=10cm和所有t，T=0℃jtxi-1ii+1j+1j-1,1,1,,1,22ijijijijijTTTTTtx显式差分格式functionPDE1_FDM_im%显式差分法求解一维热传导方程clearall;clcn=6;%空间节点数m=8;%时间节点数T(1,1:m)=100;T(n,1:m)=0;%边界条件T(1:n,1)=0;%初始条件alpha=2;%热扩散系数,m/sa=10;%forx坐标长度,cmb=8;%for时间长度,sh=a/(n-1);%空间步长k=b/(m-1);%时间步长r=alpha.*k./h.^2;forj=2:m%fortimefori=2:n-1%forxT(i,j+1)=(1-2.*r).*T(i,j)+r.*(T(i+1,j)+T(i-1,j));endendT8=T(1:n,m)'jtxi-1ii+1j+1j-122222,,121,21,,1,1,1,11,122222()ijijijijijijijijijTTxxTxTTTTTTx,1,1,22(/2)ijijijTTTttCrank-Nicolson差分格式1,1,11,11,,1,2222ijijijijijijrTrTrTrTrTrT1,1,11,11,,1,2222ijijijijijijTTTTTTrr2,13,14,12,3,4,2222jjjjjjTTTTTTrr1,12,13,11,2,3,2222jjjjjjTTTTTTrr1,11,1jjTTci=1i=2i=3偏微分方程的求解实例1:恒靠近速度时两等直径液滴形成的液膜内流体排液速率的模拟—问题描述(h(r,t)~t?)Figure1.Schematicsketchofacurvefilmzrbr2hrfra2hn靠近方向排液方向排液方向偏微分方程的求解实例1:恒靠近速度时两等直径液滴形成的液膜内流体排液速率的模拟—模型建立23433223423222311128311chhhhhhtrrrrrrrhhhhhrrrrrr初始液膜厚度:h(r,t=0)=h0c+r2/rb边界条件:r=ra,V=2e-6m/s4340:9chhrhtrh0c=2.8e-6mrb=0.0015mσ=0.03N/mμc=0.3Paszrbr2hrfra2hn偏微分方程的求解实例1:恒靠近速度时两等直径液滴形成的液膜内流体排液速率的模拟—方程离散r:ninhVti=1234nn-1……..5n+12111:4121211121nnbninrhhrnhn43433354321343:94649ccihhhtrhhhhhhr11211223211233421124422222464iiiiiiiiiiiiiiiiiihhhrrhhhhrrhhhhhrrhhhhhhrri=n:P=022212bhhPrrrr1524hhhhfunctionCVFDclearall;clc;formatlongeh0c=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,mh0=h0c+r.^2./rb;%初始膜厚度mh0(1)=h0(5);h0(2)=h0(4);zrbr2hrfra2hn2h0choldonk=10;tf=12;fori=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,:);endholdofffunctiondhdt=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);35433432869chhhhhtr2114/211/211/211bnnnrrhhnhnfori=4:n-1Dhr(i)=(h(i+1)-h(i-1))/(2*dr);D2hr(i)=(h(i+1)-2*h(i)+h(i-1))/(dr^2);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(:);112iiihhhrr211222iiiihhhhrr3211233222iiiiihhhhhrr4211244464iiiiiihhhhhhrr432143223322322321221111183iiiiicihhhhHrrrrrrrhhhhHrrrrrrhhHhHt:nhinVt%不同时刻液膜表面的压力分布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);forj=4:n-1p(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);endplot(r(3:n)./ra,p(k,3:n).*theta.*2./rb)11211222222114iiiiiiibhhhrrhhhhrrrhhPrrrfunctionPDE1Dd_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)%------------------------------------------------------------------functionf=ic(x)f=0;偏微分方程的求解实例2:functionU=CrankNicolson(f,c1,c2,a,b,alpha,n,m)%使用Crank-Nicolson有限差分方法求解一维动态传热模型%Init