偏微分方程及其求解实例ppt课件.ppt
22222, , ,uuuuuuuABCDEFufx y uxx yyxyxy 22222, , ,uuuuuuuABCDEFufx y uxx yyxyxy 22uutx22220uuxy22222uutx(1) 导热方程导热方程:(2) 拉普拉斯方程拉普拉斯方程: 如稳态静电场和稳态温度分布模型如稳态静电场和稳态温度分布模型(3) 波动方程波动方程: 一维弦振动模型一维弦振动模型(1) Dirichlet边界条件边界条件(第一类边界条件第一类边界条件)-在边界上给定在边界上给定u(x,y)的值的值(2) Neumann边界条件边界条件(第二类边界条件第二类边界条件)-在边界上给定在边界上给定u(x,y)的外法向导数值的外法向导数值(3)Robin边界条件边界条件(混合边界条件混合边界条件, 第三类边界条件第三类边界条件)-在边界上给定在边界上给定u(x,y)与其外法向导数值的线形组合的值与其外法向导数值的线形组合的值,即即,ux y,ux yn,uqux yn(1) 有限差分法有限差分法(finite difference method) 将求解域内划分一定数目的网格将求解域内划分一定数目的网格( (交叉点为网格点或节点交叉点为网格点或节点),),在所有节点上在所有节点上, ,PDE方程用有限差分近似代替，得到代数方程用有限差分近似代替，得到代数方程组。这种方法数值稳定性好，但是不能获得节点之间方程组。这种方法数值稳定性好，但是不能获得节点之间的解，对不规则几何形状或导热系数突变的情况不适合。的解，对不规则几何形状或导热系数突变的情况不适合。(2) 正交配置法正交配置法(method of weighted residuals) 首先选择一个多项式作为试函数，将此试函数代入微分方首先选择一个多项式作为试函数，将此试函数代入微分方程，求出多项式的根作为配置点，令在各配置点试函数代程，求出多项式的根作为配置点，令在各配置点试函数代入微分方程后的残差（或余量）为零，得到关于多项式系入微分方程后的残差（或余量）为零，得到关于多项式系数的代数方程组，然后求解此方程组得到多项式中各项的数的代数方程组，然后求解此方程组得到多项式中各项的系数，得到的多项式即为微分方程的近似解析解。系数，得到的多项式即为微分方程的近似解析解。(3) MOL法法(method of lines) 将一个自变量当成连续变量，对其余的自变量用有限差分将一个自变量当成连续变量，对其余的自变量用有限差分法或者正交配置法进行离散，从而把偏微分方程转变为常法或者正交配置法进行离散，从而把偏微分方程转变为常微分方程组，然后用龙格库塔法积分求解。常用于求解微分方程组，然后用龙格库塔法积分求解。常用于求解一维动态和二维稳态一维动态和二维稳态PDE方程。方程。(4) 有限元法有限元法(finite element methods,FEM) 把求解域先划分为大量的单元（把求解域先划分为大量的单元（elements），其中任意大），其中任意大小和方向的三角形网格尤其适用于二维的情况。三角形顶小和方向的三角形网格尤其适用于二维的情况。三角形顶点称为节点（点称为节点（nodes），并与相邻单元相连接。将），并与相邻单元相连接。将PDE离离散为代数方程组求解。优点是易于处理复杂几何区域，易散为代数方程组求解。优点是易于处理复杂几何区域，易于与各种边界条件组合使用，可同时提供节点和整个求解于与各种边界条件组合使用，可同时提供节点和整个求解域内的解。域内的解。差分格式：差分格式：差分格式：差分格式：差分格式：差分格式：11,12kki ji juuutt1、显式差分法2、隐式差分法一维动态PDE模型的求解22TTtx初始条件： t=0和0 x10，T=0 边界条件： x=0 cm和所有t，T=100 x=10cm和所有t，T=0 jtxi-1ii+1j+1j-1,1,1,1,22i ji jiji jijTTTTTtx显式差分格式显式差分格式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)jtxi-1ii+1j+1j-122222,121,21,1,1,1,11,122222()i ji ji jiji jijiji jijTTxxTxTTTTTTx,1,1,22(/ 2)i ji ji jTTTttCrank-Nicolson差分格式差分格式1,1,11,11,1,2222iji jijiji jijrTr TrTrTr TrT1,1,11,11,1,2222iji jijiji jijTTTTTTrr2,13,14,12,3,4,2222jjjjjjTTTTTTrr1,12,13,11,2,3,2222jjjjjjTTTTTTrr1,11,1jjTTci=1i=2i=3偏微分方程的求解偏微分方程的求解实例实例1 1: : 恒靠近速度时两等直径液滴形成的液膜内流恒靠近速度时两等直径液滴形成的液膜内流 体排液速率的模拟体排液速率的模拟问题描述问题描述(h(r,t)t ?)Figure1. Schematic sketch of a curve film z rb r 2h rf ra 2hn 靠近方向靠近方向排液方向排液方向排液方向排液方向偏微分方程的求解偏微分方程的求解实例实例1 1: : 恒靠近速度时两等直径液滴形成的液膜内流恒靠近速度时两等直径液滴形成的液膜内流 体排液速率的模拟体排液速率的模拟模型建立模型建立23433223423222311128311chhhhhhtrrrrrrrhhhhhrrrrrr 初始液膜厚度初始液膜厚度: h(r,t=0)=h0c+r2/rb边界条件边界条件: r=ra, V=2e-6 m/s4340:9chhrhtr h0c=2.8e-6mrb=0.0015m=0.03N/mc=0.3Pas z rb r 2h rf ra 2hn 偏微分方程的求解偏微分方程的求解实例实例1 1: : 恒靠近速度时两等直径液滴形成的液膜内流恒靠近速度时两等直径液滴形成的液膜内流 体排液速率的模拟体排液速率的模拟方程离散方程离散r:ninhVt i=1 2 3 4nn-1.5n+12111:4121211121nnbninrhhrnhn43433354321343:94649ccihhhtrhhhhhhr 11211223211233421124422222464iiiiiiiiiiiiiiiiiihhhrrhhhhrrhhhhhrrhhhhhhrri=n:P=022212bhhPrrrr1524hhhhfunction 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); z rb r 2h rf ra 2hn 2h0c 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);35433432869chhhhhtr 2114/21 1/ 211/ 211bnnnrrhhnhn 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(:);112iiihhhrr211222iiiihhhhrr3211233222iiiiihhhhhrr4211244464iiiiiihhhhhhrr432143223322322321221111183iiiiicihhhhHrrrrrrrhhhhHrrrrrrhh Hh Ht :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); 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) 11211222222114iiiiiiibhhhrrhhhhrrrhhPrrr 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 = Ufunction 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
