有限元精选课件.ppt

《有限元精选课件.ppt》由会员分享，可在线阅读，更多相关《有限元精选课件.ppt（77页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、关于有限元第一页，本课件共有77页求数值解方法求数值解方法n差分方法差分方法n有限元方法有限元方法nMATLABMATLAB的的pedpepedpe函数函数nMATLABMATLAB的的PDEtoolPDEtool工具箱工具箱第二页，本课件共有77页 偏微分方程分类偏微分方程分类n椭圆偏微分方程椭圆偏微分方程n双曲线偏微分方双曲线偏微分方程程n抛物线偏微分方抛物线偏微分方程程第三页，本课件共有77页4椭圆偏微分方程特例椭圆偏微分方程特例拉普拉斯方程拉普拉斯方程 拉普拉斯方程是最简单的椭圆偏微分方程，以下以拉普拉斯方程为例，讲述椭圆偏微分方程的的数值解法。拉普拉斯方程形式如下：第四页，本课件共有

2、77页5椭圆偏微分方程边界条件椭圆偏微分方程边界条件椭圆偏微分方程边界条件有以下三种提法：其中第一种提法最为普遍，下面以第一种边界条件，拉普拉斯方程为例介绍椭圆偏微分方程常用的五点差分格式和工字型差分格式的解法。第五页，本课件共有77页6五点差分格式五点差分格式五点差分格式最常用的格式，其形式如下：注意：这里的边界为矩形区域。注意：这里的边界为矩形区域。第六页，本课件共有77页7五点差分格式算法五点差分格式算法 注意注意：要保证：要保证x方向和方向和y方向上的网格步长相等才方向上的网格步长相等才能使用上面的公式。能使用上面的公式。第七页，本课件共有77页8 五点差分格式在五点差分格式在MATL

3、ABMATLAB中实现中实现function u=peEllip5(nx,minx,maxx,ny,miny,maxy)%x方向的节点数：nx%求解区间x的左端：minx%求解区间x的右端：maxx%y方向的节点数：ny%求解区间y的左端：miny%求解区间y的右端：maxy%求解区间上的数值解：uformat long;hx=(maxx-minx)/(nx-1);hy=(maxy-miny)/(ny-1);u0=zeros(nx,ny);for j=1:ny u0(j,1)=EllIni2Uxl(minx,miny+(j-1)*hy);u0(j,nx)=EllIni2Uxr(maxx,min

4、y+(j-1)*hy);endfor j=1:nx u0(1,j)=EllIni2Uyl(minx+(j-1)*hx,miny);u0(ny,j)=EllIni2Uyr(minx+(j-1)*hx,maxy);end%边界条件的离散第八页，本课件共有77页9 五点差分格式在五点差分格式在MATLABMATLAB中实现中实现A=-4*eye(nx-2)*(ny-2),(nx-2)*(ny-2);b=zeros(nx-2)*(ny-2),1);for i=1:(nx-2)*(ny-2);if mod(i,nx-2)=1 if i=1 A(1,2)=1;A(1,nx-1)=1;b(1)=-u0(1,

5、2)-u0(2,1);else if i=(ny-3)*(nx-2)+1 A(i,i+1)=1;A(i,i-nx+2)=1;%注意边界节点的离散方式 b(i)=-u0(ny-1,1)-u0(ny,2);else A(i,i+1)=1;A(i,i-nx+2)=1;A(i,i+nx-2)=1;b(i)=-u0(floor(i/(nx-2)+2,1);end end else if mod(i,nx-2)=0 if i=nx-2第九页，本课件共有77页10 五点差分格式在五点差分格式在MATLABMATLAB中实现中实现 A(i,i-1)=1;%注意边界节点的离散方式 A(i,i+nx-2)=1;b

6、(i)=-u0(1,nx-1)-u0(2,nx);else if i=(ny-2)*(nx-2)A(i,i-1)=1;A(i,i-nx+2)=1;b(i)=-u0(ny-1,nx)-u0(ny,nx-1);else A(i,i-1)=1;A(i,i-nx+2)=1;A(i,i+nx-2)=1;b(i)=-u0(floor(i/(nx-2)+1,nx);end end else if i1&i(ny-3)*(nx-2)&i0,a0格式不同是为了满足差分格式的稳定性，若第一个式子a0 for j=1:(n+M)u0(j)=IniU(minx+(j-M-1)*h);%向左延拓M个节点的函数值 end

7、else for j=1:(n+M)u0(j)=IniU(minx+(j-1)*h);%向左延拓M个节点的函数值 endendu1=u0;for k=1:M if a0 for i=(k+1):n+M u1(i)=-dt*a*(u0(i)-u0(i-1)/h+u0(i);end else for i=1:n+M-k u1(i)=-dt*a*(u0(i+1)-u0(i)/h+u0(i);end 一维对流方程一维对流方程迎风格式迎风格式第十九页，本课件共有77页20一维对流方程一维对流方程迎风格式算例迎风格式算例end u0=u1;endif a0 u=u1(M+1):M+n);else u=u1

8、(1:n);endformat long;第二十页，本课件共有77页 2022/12/821一维对流方程一维对流方程迎风格式算例迎风格式算例 然后在MATLAB窗口输入下列命令：u=peHypbYF(1,0.005,101,0,1,100);第二十一页，本课件共有77页22一维对流方程一维对流方程迎风格式算例结果迎风格式算例结果t=0时，时，u,x分布图分布图t=0.5时，时，u,x分布图分布图第二十二页，本课件共有77页23一维对流方程一维对流方程拉克斯拉克斯-弗里德里希斯格式弗里德里希斯格式第二十三页，本课件共有77页24一维对流方程一维对流方程拉克斯拉克斯-弗里德里希斯格式弗里德里希斯格

9、式第二十四页，本课件共有77页25拉克斯拉克斯-弗里德里希斯格式算例弗里德里希斯格式算例第二十五页，本课件共有77页26拉克斯拉克斯-弗里德里希斯格式算例弗里德里希斯格式算例第二十六页，本课件共有77页27拉克斯拉克斯-弗里德里希斯格式算例弗里德里希斯格式算例第二十七页，本课件共有77页28一维对流方程一维对流方程拉克斯拉克斯-温德洛夫格式温德洛夫格式第二十八页，本课件共有77页29一维对流方程一维对流方程拉克斯拉克斯-温德洛夫格式温德洛夫格式第二十九页，本课件共有77页30一维对流方程一维对流方程拉克斯拉克斯-温德洛夫格式算例温德洛夫格式算例第三十页，本课件共有77页31一维对流方程一维对流

10、方程拉克斯拉克斯-温德洛夫格式算例温德洛夫格式算例第三十一页，本课件共有77页32一维对流方程一维对流方程拉克斯拉克斯-温德洛夫格式算例温德洛夫格式算例第三十二页，本课件共有77页33一维对流方程一维对流方程比姆比姆-沃明格式沃明格式第三十三页，本课件共有77页34一维对流方程一维对流方程比姆比姆-沃明格式沃明格式第三十四页，本课件共有77页35一维对流方程一维对流方程比姆比姆-沃明格式算例沃明格式算例第三十五页，本课件共有77页36一维对流方程一维对流方程多步格式多步格式 多步格式也有多种，这里只简单介绍其中几种格式。包括多步格式也有多种，这里只简单介绍其中几种格式。包括Richtmyer多

11、步格式、多步格式、拉克斯拉克斯-温德洛夫多步格式、温德洛夫多步格式、MacCormack多步格式。多步格式。第三十六页，本课件共有77页37一维对流方程一维对流方程多步格式多步格式第三十七页，本课件共有77页 2022/12/838一维对流方程一维对流方程多步格式算例多步格式算例 Richtmyer多步格式算出的结果并不理想，多步格式算出的结果并不理想，不但左边有波动，而且光滑性也不好。拉克斯不但左边有波动，而且光滑性也不好。拉克斯-温温德洛夫多步格式算出的结果比较不错，虽然左边德洛夫多步格式算出的结果比较不错，虽然左边有点小波动，但是初始函数的宽度和高度都保持有点小波动，但是初始函数的宽度和

12、高度都保持的不错。的不错。MacCormack多步格式求得的结果和拉多步格式求得的结果和拉克斯克斯-温德洛夫多步格式算出的结果差不多。温德洛夫多步格式算出的结果差不多。第三十八页，本课件共有77页39双曲线偏微分方程双曲线偏微分方程二维对流方程二维对流方程第三十九页，本课件共有77页40二维对流方程二维对流方程拉克斯拉克斯-弗里德里希斯格式弗里德里希斯格式第四十页，本课件共有77页41二维对流方程二维对流方程拉克斯拉克斯-弗里德里希斯格式弗里德里希斯格式第四十一页，本课件共有77页42二维对流方程二维对流方程拉克斯拉克斯-弗里德里希斯格式算例弗里德里希斯格式算例u=peHypb2LF(1,1,

13、0.005,101,0,1,101,0,1,100);第四十二页，本课件共有77页43二维对流方程二维对流方程拉克斯拉克斯-弗里德里希斯格式算例弗里德里希斯格式算例 结果与初始值对比，可以看出，拉克斯结果与初始值对比，可以看出，拉克斯-弗里德里希斯格式算出的结果非常好。弗里德里希斯格式算出的结果非常好。第四十三页，本课件共有77页44二维对流方程二维对流方程近似分裂格式近似分裂格式近似分裂格式也是一种不错的格式，其结果也非常接近理论值。近似分裂格式也是一种不错的格式，其结果也非常接近理论值。第四十四页，本课件共有77页45抛物线偏微分方程抛物线偏微分方程扩散方程扩散方程 在实际应用中遇到的抛物

14、线偏微分方程主要是扩散方程。扩散方程有很强在实际应用中遇到的抛物线偏微分方程主要是扩散方程。扩散方程有很强的物理背景，例如不用物质之间的扩散过程、热传递过程、波传播等过程都可的物理背景，例如不用物质之间的扩散过程、热传递过程、波传播等过程都可以用扩散过程来描述。下面以扩散方程为例介绍几种差分格式。以用扩散过程来描述。下面以扩散方程为例介绍几种差分格式。第四十五页，本课件共有77页46扩散方程扩散方程显式格式显式格式第四十六页，本课件共有77页47扩散方程扩散方程显式格式显式格式第四十七页，本课件共有77页48扩散方程扩散方程显式格式算例显式格式算例u=peParabExp(1,0.005,10

15、1,0,1,100)第四十八页，本课件共有77页49扩散方程扩散方程显式格式算例显式格式算例 我们知道显示格式虽然简单，但其精度很差，而且求得的解容易出现震荡。次算例的结果如我们知道显示格式虽然简单，但其精度很差，而且求得的解容易出现震荡。次算例的结果如下：下：数值结果震荡的非常厉害，说明显式格式在这种条件下不稳定。数值结果震荡的非常厉害，说明显式格式在这种条件下不稳定。第四十九页，本课件共有77页50扩散方程扩散方程跳点格式跳点格式相同的算例数值结果同样震荡的非常厉害，说明跳点格式在这种条件下也不稳定。相同的算例数值结果同样震荡的非常厉害，说明跳点格式在这种条件下也不稳定。第五十页，本课件共

16、有77页51扩散方程扩散方程隐式格式隐式格式第五十一页，本课件共有77页52扩散方程扩散方程隐式格式隐式格式第五十二页，本课件共有77页53扩散方程扩散方程隐式格式隐式格式第五十三页，本课件共有77页54扩散方程扩散方程隐式格式算例隐式格式算例用隐式格式求解下面扩散方程的初值问题：u=peParabImp(1,0.005,101,0,1,0,1,100)第五十四页，本课件共有77页55扩散方程扩散方程隐式格式算例隐式格式算例结果是一条比较稳定的直线。结果是一条比较稳定的直线。第五十五页，本课件共有77页56扩散方程扩散方程克拉克克拉克-尼科尔森格式尼科尔森格式第五十六页，本课件共有77页57扩

17、散方程扩散方程克拉克克拉克-尼科尔森格式尼科尔森格式第五十七页，本课件共有77页58扩散方程扩散方程克拉克克拉克-尼科尔森格式尼科尔森格式第五十八页，本课件共有77页59扩散方程扩散方程克拉克克拉克-尼科尔森格式算例尼科尔森格式算例第五十九页，本课件共有77页60扩散方程扩散方程克拉克克拉克-尼科尔森格式算例尼科尔森格式算例第六十页，本课件共有77页61扩散方程扩散方程加权隐式格式加权隐式格式 通过算例我们可以知道，通过取不同参数，用加权隐式格式算得的结果差别通过算例我们可以知道，通过取不同参数，用加权隐式格式算得的结果差别不大，只是达到稳态的时间不一样而已。不大，只是达到稳态的时间不一样而已

18、。第六十一页，本课件共有77页62差分方法小结差分方法小结 以上我们介绍了差分方法在椭圆型、抛物型和双曲型偏微分方程中的应用，用差分格式求解偏微分方程的基本步骤是一样的，首先把连续的问题离散化，建立差分格式，然后根据差分格式对求解区域进行网格剖分，最后求解方程。下面将简单介绍有限元方法。另外变分法、边界元法、混合有限元法和多重网格法等也是偏微分方程数值求解方法，有兴趣的同学可以参考相关书籍。第六十二页，本课件共有77页63有限元方法有限元方法介绍介绍 有限元方法是数值求解偏微分方程边值问题的一种方法，此方法首先于20世纪50年代初由工程师提出，并用于求解简单的结构问题。有限元方法是这一种系统的

19、数值方法，并奠定其数学基础，是在60年代中期以冯康先生为代表的中国学者与西方学者独立并行完成的。有限元方法不同于差分方法，主要有以下三大特点：(1)从数学物理问题的变分原理出发，而不是从微分方程出发，因此事从问题的整体描述而不是从问题的局部描述出发。(2)对所考虑问题的区域(以二维情形为例)作三角形（或其他简单多边形）剖分，而不是仅作矩形剖分。(3)用剖分区域上的简单函数(如分片多项式)去逼近原问题的解，而不是只在剖分节点上的数值逼近。第六十三页，本课件共有77页64有限元方法有限元方法一维边值问题算例一维边值问题算例用有限元方法求解如下一维边值问题：解：一维边值问题线性有限元数值解法的MAT

20、LAB程序如下：%线性有限元方法%参数设置N=10;X=0:1/(N+1):1;b=zeros(N+1,1);A=zeros(N+1,N+1);for i=2:N+1 F1=(x)2*(N+1)*(x-X(i-1).*sin(pi*x/2);%句柄函数 R1=quad(F1,X(i-1),X(i);F2=(x)2*(N+1)*(X(i+1)-x).*sin(pi*x/2);R2=quad(F2,X(i),X(i+1);b(i-1)=R1+R2;第六十四页，本课件共有77页65endF1=(x)1.*sin(pi*x/2);b(N+1)=quad(F1,X(N+1),X(N+2);%quad:数

21、值积分%适应度矩阵a11=(N+1)+(pi2)/(12*(N+1);a12=-(N+1)+(pi2)/(24*(N+1);for i=1:N A(i,i)=2*a11;A(i,i+1)=a12;A(i+1,i)=a12;endA(N+1,N+1)=a11;%得到初始数值解%解方程Ax=bc=Ab;x=vertcat(0,c);%垂直串联矩阵y=4/(pi2)*sin(pi*X/2);y=y;Error=x-y;%绘制图像figure(1);有限元方法有限元方法一维边值问题算例一维边值问题算例第六十五页，本课件共有77页66grid on;plot(X,y,ro-,X,x,b);title(N

22、umerical solutions vs Accurate solutions);legend(Accurate solutions,Numerical solutions,0);%添加图例 如图所示为数值解与解析解的比较，可知有限元方法对这个一维边值问题是比较好的。有限元方法有限元方法一维边值问题算例一维边值问题算例第六十六页，本课件共有77页67MATLABMATLAB的的pedpepedpe函数函数pedpepedpe函数的说明函数的说明 MATLAB软件提供了pdepe函数，该函数不但可以用来求解偏微分方程，也可以用来求解偏微分方程组，函数的调用格式为：输入的参数中 pdefun是偏

23、微分方程的描述函数，方程必须具有如下形式函数pdefun由用户自己编写，函数形式为：其输出的c,f,s即为式(1)中的三个已知函数c,f,s,它们也可以是向量值函数，x,t,u与方程(1)中的参数意义相同，du表示的是u对x的一阶倒数。(1)第六十七页，本课件共有77页68MATLABMATLAB的的pedpepedpe函数函数pedpepedpe函数的说明函数的说明pdebc是偏微分方程的边界条件描述函数，函数必须具有如下形式：函数pdebc由用户自己编写，函数形式为：其中是xa,xb,ua,ub分别表示变量x,u的下边界和上边界。pdeic是偏微分方程的初始条件描述函数，函数必须具有如下形

24、式：函数pdeic由用户自己编写，函数形式如下：函数pdepe中的m即为方程(1)中的m。x,t是偏微分方程的自变量，它们一般是多维向量。输出的sol是一个三维数组，sol(i,j,k)表示的是自变量分别取x(i),t(j)时u(k)的值。由sol可以直接通过pdeval()某个点的函数值。第六十八页，本课件共有77页69MATLABMATLAB的的pedpepedpe函数函数pedpepedpe函数的实例函数的实例求解偏微分方程组(2)解：分别编写pdefun函数、pdebc函数、pdeic函数：第六十九页，本课件共有77页70MATLABMATLAB的的pedpepedpe函数函数pedp

25、epedpe函数的实例函数的实例%目标PDE函数function c,f,s=pdefun(x,t,u,du)c=1;1;f=0.024*du(1);0.17*du(2);temp=u(1)-u(2);s=-1;1.*(exp(5.73*temp)-exp(-11.46*temp);%边界条件函数function pa,qa,pb,qb=pdebc(xa,ua,xb,ub,t)%a表示下边界，b表示上边界pa=0;ua(2);qa=1;0;pb=ub(1)-1;0;qb=0;1;%初值条件函数function u0=pdeic(x)u0=1;0;编写好以上函数之后执行：第七十页，本课件共有77

26、页71MATLABMATLAB的的pedpepedpe函数函数pedpepedpe函数的实例函数的实例x=0:0.05:1;t=0:0.05:2;m=0;sol=pdepe(m,pdefun,pdeic,pdebc,x,t);subplot(121)surf(x,t,sol(:,:,1)title(The Solution of u_1)subplot(122)surf(x,t,sol(:,:,2)title(The Solution of u_2)第七十一页，本课件共有77页72MATLABMATLAB的的pedpepedpe函数函数pedpepedpe函数的实例函数的实例第七十二页，本课件

27、共有77页73MATLABMATLAB的的PDEtoolPDEtool工具箱工具箱 MATLAB专门提供了用于求解偏微分方程的工具箱PDEtool，它可用来解各种常见的二阶偏微分方程，但不能求解偏微分方程组。工具箱对三种常见的二阶偏微分方程有一定的要求。双曲型方程 抛物型方程上两式中，d,c,a,f必须是常数；椭圆型方程式中，c,a,f为给定的函数或常数。第七十三页，本课件共有77页 74 MATLABMATLAB的的PDEtoolPDEtool工具箱工具箱打开PDEtool的交互窗口，如图所示 pdetool在MATLAB命令行窗口利用命令：可解方程类型设置边界输入方程调节网格大小求解方程与

28、绘制方程解的图形第七十四页，本课件共有77页75 MATLABMATLAB的的PDEtoolPDEtool工具箱工具箱实例实例抛物型方程定解问题：第一步：单机工具栏的“PDE”按钮，在弹出的左窗口左侧选择Parabolic(抛物型)，在右侧输入相应的值，输入完成后单机“OK”按钮。c=2,a=3,f=0,d=1.第二步：设置绘制区域，在“Option”菜单中选择“Axis Limits”,打开绘制区域窗口，设置x和y的范围均为0,1。选中“Option”里面的“grid”使得绘制区域内画出网格。第七十五页，本课件共有77页76 MATLABMATLAB的的PDEtoolPDEtool工具箱工具箱实例实例第三步：设置初始条件，选择“Solve”菜单中的“parameters”,在弹出的求解参数窗口中设置初始值。第四步：设置边界条件，首先选择“Boundary”菜单中的“Boundary Mode”进入边界条件设置模式。利用工具条左边的椭圆，矩形等按钮选择要设置的区域，选择“Boundary”菜单中的“Secify Boundary Conditions.”设置该区域的边界条件，重复这个步骤直到设置好全部的边界条件。本题是矩形上的Dirichlet边界条件。第七十六页，本课件共有77页 2022/12/8感感谢谢大大家家观观看看第七十七页，本课件共有77页