第五章 偏微分方程的有限元法

1、计算物理学Harbin Institute of Technology Yangkun 1/104第五章第五章 偏微分方程的有限元法偏微分方程的有限元法5.1 泛函与变分原理泛函与变分原理5.2 基于变分原理的有限元法基于变分原理的有限元法5.3 matlab有限元法工具箱有限元法工具箱计算物理学Harbin Institute of Technology Yangkun 2/104第五章第五章 偏微分方程的有限元法偏微分方程的有限元法有限元法（有限元法（FEA，Finite Element Analysis，FEM） 有限元法的基本思想是用有限元法的基本思想是用较简单的问题较简单的问题代替代

2、替复杂问复杂问题题，然后再对简单问题进行求解的数值计算方法。，然后再对简单问题进行求解的数值计算方法。 有限元法将求解域看成是由许多被称为有限元法将求解域看成是由许多被称为有限元有限元的小的互的小的互连子域组成，对每一单元假定一个较简单的近似解，然后推连子域组成，对每一单元假定一个较简单的近似解，然后推导求解这个域总的满足条件，从而得到问题的解。这个解不导求解这个域总的满足条件，从而得到问题的解。这个解不是准确解，而是近似解。有限元不仅计算精度高，而且能适是准确解，而是近似解。有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的数值计算方法。应各种复杂形状，因而成为行之有效的数值计

3、算方法。 有限元法于上世纪有限元法于上世纪50年代首先在力学领域年代首先在力学领域-飞机结飞机结构的静、动态特性分析中得到应用，随后很快广泛的应用构的静、动态特性分析中得到应用，随后很快广泛的应用于求解热传导、电磁场、流体力学等连续性问题。于求解热传导、电磁场、流体力学等连续性问题。有限元法主要用于求解有限元法主要用于求解拉普拉斯方程拉普拉斯方程和和泊松方程泊松方程所描述的所描述的各类物理场中。各类物理场中。计算物理学Harbin Institute of Technology Yangkun 3/104第五章第五章 偏微分方程的有限元法偏微分方程的有限元法有限元法有限元法-变分原理变分原理

4、基于变分原理的有限元法是逼近论、偏微分方程、变基于变分原理的有限元法是逼近论、偏微分方程、变分与泛函分析的巧妙结合。分与泛函分析的巧妙结合。 基于变分原理的有限元法以基于变分原理的有限元法以变分原理变分原理为基础，把所为基础，把所要求解的要求解的微分方程微分方程定解问题，首先转化为相应的定解问题，首先转化为相应的变分问变分问题题，即，即泛函求极值泛函求极值问题；它将求解域看成是由许多称为问题；它将求解域看成是由许多称为有限元的小的互连子域组成，然后利用有限元的小的互连子域组成，然后利用剖分插值剖分插值，对每一，对每一单元假定一个合适的单元假定一个合适的(较简单的）近似解，把离散化的较简单的）近

5、似解，把离散化的变分变分问题问题转化为普通多元函数的转化为普通多元函数的极值问题极值问题，然后推导求解这，然后推导求解这个域总的满足条件个域总的满足条件(边界条件），即最终归结为一组多元的边界条件），即最终归结为一组多元的代数方程组代数方程组，求解代数方程组，就得到待求边值问题的，求解代数方程组，就得到待求边值问题的数值解。数值解。计算物理学Harbin Institute of Technology Yangkun 4/104第五章第五章 偏微分方程的有限元法偏微分方程的有限元法有限元法有限元法-加权余数法加权余数法 自从自从1969年以来，某些学者在流体力学中应用加权余年以来，某些学者在流

6、体力学中应用加权余数法中的迦辽金法或最小二乘法等同样获得了有限元方程，数法中的迦辽金法或最小二乘法等同样获得了有限元方程，因而有限元法可应用于以任何微分方程所描述的各类物理因而有限元法可应用于以任何微分方程所描述的各类物理场中，而不再要求这类物理场和泛函的极值问题有所联系。场中，而不再要求这类物理场和泛函的极值问题有所联系。 加权余数法的核心思想是：近似解与解析解相比会存加权余数法的核心思想是：近似解与解析解相比会存在误差在误差R，但是可以通过一个准则使，但是可以通过一个准则使R尽量小，求解这个等尽量小，求解这个等式，就可以得到待定常数的值，也就得到了近似解。式，就可以得到待定常数的值，也就得

7、到了近似解。计算物理学Harbin Institute of Technology Yangkun 5/104第五章第五章 偏微分方程的有限元法偏微分方程的有限元法 有限元法特点有限元法特点有限元法的物理意义直观明确，理论完整可靠。有限元法的物理意义直观明确，理论完整可靠。 因为因为变变分原理分原理描述了支配物理现象的物理学中的描述了支配物理现象的物理学中的最小作用原理最小作用原理（如力学中的最小势能原理）。（如力学中的最小势能原理）。 优异的解题能力。有限元法对优异的解题能力。有限元法对边界几何形状复杂边界几何形状复杂以及以及媒媒质物理性质变异质物理性质变异等复杂物理问题求解上，有突出优点：

8、等复杂物理问题求解上，有突出优点： 不受几何形状和媒质分布的复杂程度限制。不受几何形状和媒质分布的复杂程度限制。 不必单独处理第二、三类边界条件。不必单独处理第二、三类边界条件。 离散点配置比较随意，通过控制有限单元剖分密度和离散点配置比较随意，通过控制有限单元剖分密度和单元插值函数的选取，可以充分保证所需的数值计算精度。单元插值函数的选取，可以充分保证所需的数值计算精度。计算物理学Harbin Institute of Technology Yangkun 6/1045.1 泛函与变分原理泛函与变分原理 数学上，通常自变量与因变量间的关系称为函数，数学上，通常自变量与因变量间的关系称为函数，

9、而泛函则是函数集合的函数，也就是而泛函则是函数集合的函数，也就是函数的函数函数的函数，即，即自变量为函数，而不是变量。自变量为函数，而不是变量。5.1.1 泛函的定义泛函的定义 泛函通常是指一种定义域为函数，而值域为实数的泛函通常是指一种定义域为函数，而值域为实数的“函数函数”。 设设C是函数的集合，是函数的集合，B是实数集合。如果对是实数集合。如果对C中的任中的任一元素一元素y(x)，在，在B中都有一个元素中都有一个元素J与之对应，则称与之对应，则称J为为y(x)的泛函，记为的泛函，记为Jy(x)。计算物理学Harbin Institute of Technology Yangkun 7/1

10、045.1 泛函与变分原理泛函与变分原理例例5.1.1 质点在重力作用下，沿一条光滑的从质点在重力作用下，沿一条光滑的从A点到点到B点的曲线运动，如图所示。求下落时间最短的曲线。点的曲线运动，如图所示。求下落时间最短的曲线。曲线上任一小段线元长度为：曲线上任一小段线元长度为：ABxyOx0 x122222)1 (dxdxdydydxdsdxyds)1 (2捷线问题捷线问题计算物理学Harbin Institute of Technology Yangkun 8/1045.1 泛函与变分原理泛函与变分原理线元处的质点速度为线元处的质点速度为ABxyOx0 x1gyv2dxgyyvdsdT212d

11、s线元下落时间为线元下落时间为从从A点到点到B点的下落时间为点的下落时间为)(21102xyJdxgyyTxxmin)(xyJ计算物理学Harbin Institute of Technology Yangkun 9/1045.1 泛函与变分原理泛函与变分原理5.1.2 函数的变分函数的变分 设设y(x)是泛函是泛函J定义域内任一函数，如果定义域内任一函数，如果y(x)变化为变化为新函数新函数Y(x) ，且，且Y(x)属于泛函属于泛函J的定义域，则的定义域，则Y(x)与与y(x)之差为函数之差为函数y(x)的变分。的变分。)()(xyxYy变分变分y是是x的函数，它不同于函数的的函数，它不同于

12、函数的增量增量y。性质：函数求导与求变分可以交换次序性质：函数求导与求变分可以交换次序 )()()()()(yxyxYxyxYy yy xxy x 计算物理学Harbin Institute of Technology Yangkun 10/1045.1 泛函与变分原理泛函与变分原理5.1.3 泛函的变分泛函的变分定义定义最简泛函最简泛函10),()(xxdxyyxFxyJF(x,y,y)称为泛函的称为泛函的“核函数核函数”泛函的变分泛函的变分10),(),()()(xxdxyyxFyyyyxFyJyyJJ最简泛函最简泛函: 核函数只包含自变量核函数只包含自变量 x、未知函数、未知函数y(x)

13、以及导数以及导数y(x)计算物理学Harbin Institute of Technology Yangkun 11/1045.1 泛函与变分原理泛函与变分原理利用二元函数的泰勒展开利用二元函数的泰勒展开2222222( ,)1( , ,)( , ,)( , ,)1!1( , ,)( , ,)( , ,) 22!F x yy yyF x y yF x y yF x y yyyyyF x y yF x y yF x y yyy yyyy yy 计算物理学Harbin Institute of Technology Yangkun 12/1045.1 泛函与变分原理泛函与变分原理其中其中10222

14、122! yyxxyyyyy yFyFyJdxFyFy yFyJJ 1010222122xyyxxyyyyy yxJFyFy dxJFyFy yFydx 分别称为泛函的分别称为泛函的一阶变分一阶变分和和二阶变分二阶变分。22 yyyFFFFyy()( )JJ yyJ y 计算物理学Harbin Institute of Technology Yangkun 13/1045.1 泛函与变分原理泛函与变分原理泛函取极值的必要条件：泛函取极值的必要条件：一阶变分为零一阶变分为零0J性质：对于最简泛函，变分运算可以与积分、微性质：对于最简泛函，变分运算可以与积分、微分运算交换次序分运算交换次序1100

15、( , ,)( , ,)xxxxJF x y y dxF x y y dxdyyddxdx计算物理学Harbin Institute of Technology Yangkun 14/1045.1 泛函与变分原理泛函与变分原理5.1.4 泛函的极值问题泛函的极值问题 泛函的一阶变分泛函的一阶变分10()xxFFJyy dxyy 利用利用dFdFFyydxydxyydFFydydyyxdxy1 泛函的极值问题的间接泛函的极值问题的间接解法解法 转化为微分方程：欧拉方程转化为微分方程：欧拉方程0J()yy计算物理学Harbin Institute of Technology Yangkun 15/

16、1045.1 泛函与变分原理泛函与变分原理11100011000 xxxxxxxxxxFdFFJy dxydxyydxyyFdFFydxyydxyy 对于驻定问题，对于驻定问题，两边界固定两边界固定0FdFydxy010 x xx xy 这就是最简泛函的这就是最简泛函的欧拉方程欧拉方程，等价于泛函取极值的必要条件。，等价于泛函取极值的必要条件。把把变分问题转化微分方程的定解问题（边值问题）来求变分问题转化微分方程的定解问题（边值问题）来求解解。10()xxFFJyy dxyyFdFdFyyydxydxyy计算物理学Harbin Institute of Technology Yangkun 1

17、6/1045.1 泛函与变分原理泛函与变分原理 对于例对于例5.1.1求下落时间最短的轨迹求下落时间最短的轨迹0FdFydxy利用最简泛函的利用最简泛函的欧拉方程欧拉方程。)(21102xyJdxgyyTxxmin)(xyJ212yFgy计算物理学Harbin Institute of Technology Yangkun 17/1045.1 泛函与变分原理泛函与变分原理代入欧拉方程代入欧拉方程23212 2yFyg y221Fyygyy232210212 2yFdFdyydxydxgyyg y23221012ydydxyyy212yFgy计算物理学Harbin Institute of Te

18、chnology Yangkun 18/1045.1 泛函与变分原理泛函与变分原理变换得到变换得到222101dyydxyyy进一步化简得到进一步化简得到2101ddxyy积分积分211yyc计算物理学Harbin Institute of Technology Yangkun 19/1045.1 泛函与变分原理泛函与变分原理做变量替换做变量替换cossintyt 得得22211211sinsin1ytcycty而而21112sin cos2sin(1cos2 )cossincttdtdydxctdtct dttyt计算物理学Harbin Institute of Technology Yan

19、gkun 20/1045.1 泛函与变分原理泛函与变分原理对上式积分得到对上式积分得到121(sin2 )2xc ttc这样就得到了下落时间最短曲线的参数方程这样就得到了下落时间最短曲线的参数方程12211(sin 2 )2sinxc ttcyct式中常数式中常数c c1 1和和c c2 2由始末两点位置确定由始末两点位置确定练习：画出经过练习：画出经过(0,0)和和(1,1)的下落时间最短曲线。的下落时间最短曲线。连接两个点上凹的唯一一段旋轮线连接两个点上凹的唯一一段旋轮线343sin1cosxcttcyct计算物理学Harbin Institute of Technology Yangku

20、n 21/1045.1 泛函与变分原理泛函与变分原理2 泛函的极值问题的直接解法泛函的极值问题的直接解法 基本做法：基本做法：瑞利瑞利-里兹里兹(Rayleigh-Ritz)法法(1) 选定一组具有相对完备性的基函数，构造一个线选定一组具有相对完备性的基函数，构造一个线性组合的近似函数性组合的近似函数(2) 将含有将含有n个待定系数的构造函数作为近似的极值个待定系数的构造函数作为近似的极值函数，代入泛函函数，代入泛函 12,nJ y xI a aa1 niiiiiyaa：基函数 ：待定系数10),()(xxdxyyxFxyJ计算物理学Harbin Institute of Technology

21、 Yangkun 22/1045.1 泛函与变分原理泛函与变分原理0 i=1,2,3niIa(3) 为了求泛函的极值，按照多元函数取极值的必要条件为了求泛函的极值，按照多元函数取极值的必要条件(4) 求解以上方程组，求出求解以上方程组，求出 就可以得就可以得到极值函数的近似解到极值函数的近似解 12,na aa(5) 再将含有再将含有n+1个待定系数的函数个待定系数的函数作为近似极值函数，重复作为近似极值函数，重复(2)(4)，就可以得到极值函数，就可以得到极值函数新的新的近似解近似解 。如果连续两次所得到的结果接近，就认为。如果连续两次所得到的结果接近，就认为最后得到的函数就是极值函数的近似

22、解最后得到的函数就是极值函数的近似解 。11niiiya111nniiiiiiaa计算物理学Harbin Institute of Technology Yangkun 23/1045.1 泛函与变分原理泛函与变分原理例例5.1.2 求下列泛函的极值函数。求下列泛函的极值函数。解：为了满足边界条件，取基函数为解：为了满足边界条件，取基函数为 1220(4)(0)(1)0J yyyxy dxyy(1)iixx近似函数为近似函数为11niiiya xx计算物理学Harbin Institute of Technology Yangkun 24/1045.1 泛函与变分原理泛函与变分原理当当n=1时

23、时代入泛函代入泛函11ya xx 12221111012141J yaa xa xxa xx dxI a取极值取极值1222110121221410Iaxaxxxxdxa 1220(4)J yyyxy dx计算物理学Harbin Institute of Technology Yangkun 25/1045.1 泛函与变分原理泛函与变分原理计算得到计算得到近似函数近似函数159a 5(1)9yxx 同理同理n=2时时7172 (1)36941yxxx 利用欧拉方程，得到的精确解利用欧拉方程，得到的精确解2sin2sin1xyx计算物理学Harbin Institute of Technolog

24、y Yangkun 26/1045.1 泛函与变分原理泛函与变分原理00.10.20.30.40.50.60.70.80.91-0.16-0.14-0.12-0.1-0.08-0.06-0.04-0.020 xyn=1n=2精 确 解计算物理学Harbin Institute of Technology Yangkun 27/1045.1 泛函与变分原理泛函与变分原理 泛函的极值问题可以通过变分运算产生一个微分方泛函的极值问题可以通过变分运算产生一个微分方程和相应的边界条件，即欧拉方程，其解对应于最简泛程和相应的边界条件，即欧拉方程，其解对应于最简泛函的极值函数。也就是函的极值函数。也就是泛函

25、的极值问题泛函的极值问题可以等价为可以等价为在在一定边界条件下求解微分方程问题。一定边界条件下求解微分方程问题。 变分原理变分原理 通过求解一个相应的泛函的极小函数而得到通过求解一个相应的泛函的极小函数而得到偏微分方程边值问题的解。偏微分方程边值问题的解。 有限元法有限元法正是正是里兹法里兹法与与有限差分法有限差分法相结合的成果，它相结合的成果，它取长补短地在理论上以变分为基础，在具体方法构造上又取长补短地在理论上以变分为基础，在具体方法构造上又利用了有限差分法网格离散化处理的思想。利用了有限差分法网格离散化处理的思想。计算物理学Harbin Institute of Technology Y

26、angkun 28/1045.1 泛函与变分原理泛函与变分原理 20世纪世纪60年代初首次提出结构力学计算有年代初首次提出结构力学计算有限元概念的克拉夫（限元概念的克拉夫（Clough）教授形象地将其）教授形象地将其描绘为：描绘为：“有限元法有限元法=Rayleigh Ritz法法分片函分片函数数”。 有限元法是有限元法是Rayleigh Ritz法的一种局部化法的一种局部化情况。不同于求解满足整个定义域边界条件的情况。不同于求解满足整个定义域边界条件的允许函数的允许函数的Rayleigh Ritz法（往往是困难的），法（往往是困难的），有限元法将函数定义在简单几何形状（如二维有限元法将函数定

27、义在简单几何形状（如二维问题中的三角形或任意四边形）的单元域上问题中的三角形或任意四边形）的单元域上（分片函数），且不考虑整个定义域的复杂边（分片函数），且不考虑整个定义域的复杂边界条件，这是有限元法优于其它近似方法的原界条件，这是有限元法优于其它近似方法的原因之一。因之一。计算物理学Harbin Institute of Technology Yangkun 29/1045.2 基于变分原理的有限元法基于变分原理的有限元法对于具有不同物理性质和数学模型的问题，有限对于具有不同物理性质和数学模型的问题，有限元法的基本做法是相同的，只是具体公式推导和元法的基本做法是相同的，只是具体公式推导和运算

28、求解不同。运算求解不同。有限元法基本做法有限元法基本做法首先把待求的偏微分方程边值问题转化为等价的变分问首先把待求的偏微分方程边值问题转化为等价的变分问题。题。然后通过有限单元剖分的离散处理，构造一个分片解析然后通过有限单元剖分的离散处理，构造一个分片解析的有限元子空间。的有限元子空间。通过构造近似函数，把变分问题近似地转化为有限元子通过构造近似函数，把变分问题近似地转化为有限元子空间中的多元函数极值问题，由此直接利用空间中的多元函数极值问题，由此直接利用Rayleigh Ritz法探求变分问题的近似解（极值函数解），以此作法探求变分问题的近似解（极值函数解），以此作为所求边值问题的近似解。为

29、所求边值问题的近似解。 计算物理学Harbin Institute of Technology Yangkun 30/1045.2 基于变分原理的有限元法基于变分原理的有限元法有限元法具体求解步骤有限元法具体求解步骤 建立积分方程建立积分方程根据变分原理或方程余量与权函数正交化原理，建立与微根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法分方程初边值问题等价的积分表达式，这是有限元法的出发点。的出发点。区域单元剖分区域单元剖分根据求解区域的形状及实际问题的物理特点，将区域剖分根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠

30、的单元。区域单元划分是采用为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，同时还需要列出自然边外，还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。界和本质边界的节点序号和相应的边界值。计算物理学Harbin Institute of Technology Yangkun 31/1045.2 基于变分原理的有限元法基于变分原理的有限元法确定单元基

31、函数确定单元基函数根据单元中节点数目及对近似解精度的要求，选择满足一根据单元中节点数目及对近似解精度的要求，选择满足一定插值条件的插值函数作为单元基函数。有限元方法中定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元的基函数是在单元中选取的，由于各单元 具有规则的具有规则的几何形状，在选取基函数时可遵循一定的法则。几何形状，在选取基函数时可遵循一定的法则。 单元分析单元分析将各个单元中的求解函数用单元基函数的线性组合表达式将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将进行逼近；再将 近似函数代入积分方程，并对单元区近似函数代入积分方程，并对单

32、元区域进行积分，可获得含有待定系数域进行积分，可获得含有待定系数(即单元中各节点即单元中各节点 的的参数值参数值)的代数方程组，称为单元有限元方程。的代数方程组，称为单元有限元方程。总体合成总体合成在得出单元有限元方程之后，将区域中所有单元有限元方在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进行累加，形成总体有限元方程。程按一定法则进行累加，形成总体有限元方程。计算物理学Harbin Institute of Technology Yangkun 32/1045.2 基于变分原理的有限元法基于变分原理的有限元法 边界条件的处理边界条件的处理一般边界条件有三种形式，对于第二类边

33、界条件，一般在一般边界条件有三种形式，对于第二类边界条件，一般在积分表达式中可自动得到满足。对于第一类边界条件和积分表达式中可自动得到满足。对于第一类边界条件和第三类边界条件，需按一定法则对总体有限元方程进行第三类边界条件，需按一定法则对总体有限元方程进行修正满足。修正满足。解有限元方程解有限元方程根据边界条件修正的总体有限元方程组，是含所有待定未根据边界条件修正的总体有限元方程组，是含所有待定未知量的封闭方程组，采用适当的数值计算方法求解，知量的封闭方程组，采用适当的数值计算方法求解，可求得各节点的函数值。可求得各节点的函数值。计算物理学Harbin Institute of Technol

34、ogy Yangkun 33/1045.2 基于变分原理的有限元法基于变分原理的有限元法 有限元分析可分成三个阶段：有限元分析可分成三个阶段： 前置处理、计算求解和后置处理。前置处理、计算求解和后置处理。 前置处理是建立有限元模型，完成单元网格划分；后前置处理是建立有限元模型，完成单元网格划分；后置处理则是采集处理分析结果，使用户能简便提取信息，置处理则是采集处理分析结果，使用户能简便提取信息，了解计算结果。了解计算结果。 计算物理学Harbin Institute of Technology Yangkun 34/1045.2 基于变分原理的有限元法基于变分原理的有限元法1. 求解区域离散求

35、解区域离散 离散单元基本要求：离散单元基本要求：各单元只能在顶点处相交。各单元只能在顶点处相交。不同单元在边界处相连，既不能相互分离又不能相互重不同单元在边界处相连，既不能相互分离又不能相互重叠。叠。 各单元节点编号循序应一致，一律按各单元节点编号循序应一致，一律按逆时针方向逆时针方向，从最，从最小节点号开始。同一单元节点编号相差不能太悬殊，对小节点号开始。同一单元节点编号相差不能太悬殊，对多区域的编号，按区域连续编号。多区域的编号，按区域连续编号。 把求解区域分割成有限个单元体的集合。单元体形把求解区域分割成有限个单元体的集合。单元体形状原则上是任意的，一般取有规则形体。状原则上是任意的，一

36、般取有规则形体。有限元法计算步骤有限元法计算步骤计算物理学Harbin Institute of Technology Yangkun 35/1045.2 基于变分原理的有限元法基于变分原理的有限元法 三角单元是经常使用的单元剖分方法，剖分时应注三角单元是经常使用的单元剖分方法，剖分时应注意几下几点：意几下几点： 三角形不能重叠。三角形不能重叠。不能把一个三角形的顶点取为相邻三角形的边上。不能把一个三角形的顶点取为相邻三角形的边上。剖分的三角形应该避免钝角。剖分的三角形应该避免钝角。三角形不可过于狭长，最长边一般不大于最短边的三角形不可过于狭长，最长边一般不大于最短边的3倍。倍。三角形三边之比

37、尽量接近三角形三边之比尽量接近1。不能把一个三角形跨越不同的介质。不能把一个三角形跨越不同的介质。每个三角形最多只有一个边在边界上。每个三角形最多只有一个边在边界上。三角形单元面积越小，计算精度越高三角形单元面积越小，计算精度越高计算物理学Harbin Institute of Technology Yangkun 36/1045.2 基于变分原理的有限元法基于变分原理的有限元法把求解区域划分把求解区域划分m个三角形有限单元，共有个三角形有限单元，共有n个节点个节点在有限单元在有限单元e(j,k,l)上进行分片线性插值，插值函数为上进行分片线性插值，插值函数为 e123ux,yaa xa y2

38、. 选择近似函数选择近似函数计算物理学Harbin Institute of Technology Yangkun 37/1045.2 基于变分原理的有限元法基于变分原理的有限元法在单元节点上在单元节点上jj12j3jjkk12k3kkll12l3llu x ,yaa xa yu u x ,yaa xa yuu x ,yaa xa yu求解以上方程组可以得到求解以上方程组可以得到jjjkkklll1jjkklljjkklluxyuxy1uxyaa ua ua u1xy21xy1xyjkllkkljjlljkljax yx yax yx yax yx yjjkkll1xy11xy21xy 3.

39、求解单元形函数求解单元形函数e123ux,yaa xa y123a a a计算物理学Harbin Institute of Technology Yangkun 38/1045.2 基于变分原理的有限元法基于变分原理的有限元法同理可以求出同理可以求出jjkkll2jjkklljjkkll1uy1uy11uyab ub ub u1xy21xy1xyjklkljljkbyybyybyy计算物理学Harbin Institute of Technology Yangkun 39/1045.2 基于变分原理的有限元法基于变分原理的有限元法jjkkll3jjkklljjkkll1xu1xu11xuac

40、uc uc u1xy21xy1xyjlkkjllkjcxxcxxcxxjkkj1b cb c2 计算物理学Harbin Institute of Technology Yangkun 40/1045.2 基于变分原理的有限元法基于变分原理的有限元法则插值函数可以写为则插值函数可以写为e123jjjjkkkkllllessj k lux,yaa xa yab xc y uab xc y u1 2ab xc y u u Nx,yessss1Nx,yab xc y sj,k,l2单元形函数（基函数）单元形函数（基函数）1jjkkll2jjkkll3jjkkll1aa ua ua u21ab ub u

41、b u21ac uc uc u2计算物理学Harbin Institute of Technology Yangkun 41/1045.2 基于变分原理的有限元法基于变分原理的有限元法三角元三角元e插值函数可以改写为矩阵形式插值函数可以改写为矩阵形式 jeeeejklkeeluux,yNNNuu NU jeeejklekel uuu NN NNU计算物理学Harbin Institute of Technology Yangkun 42/1045.2 基于变分原理的有限元法基于变分原理的有限元法下面以泊松方程为例讨论有限元解法下面以泊松方程为例讨论有限元解法22220uu0 x,yDxyu (

42、x,y)ux,y 22mee 1DuuJ udxdyJuxy所对应的泛函为所对应的泛函为4. 建立单元特征式建立单元特征式难点难点：寻找与微分方程对应的泛函：寻找与微分方程对应的泛函计算物理学Harbin Institute of Technology Yangkun 43/1045.2 基于变分原理的有限元法基于变分原理的有限元法在第在第e个三角元的泛函个三角元的泛函22eeeDuuJx,ydxdyxyue123ux,yaa xa y由于由于2222e22jjkkllDDu1dxdya dxdyab ub ub ux4 1jjkkll2jjkkll3jjkkll1aa ua ua u21ab

43、 ub ub u21ac uc uc u2e2u=ax计算物理学Harbin Institute of Technology Yangkun 44/1045.2 基于变分原理的有限元法基于变分原理的有限元法 jjjkjlkjkk22ejjkkllD2jjklxljljllklkljjjklkjklklljTjklklu1dxdyb ub ub ux4b1 uuu b4bbu1 b bb bb b1b bb uuu bbbb u4buuuuuuubb b4b bb bb beeeUKU改写为矩阵形式改写为矩阵形式计算物理学Harbin Institute of Technology Yangku

44、n 45/1045.2 基于变分原理的有限元法基于变分原理的有限元法其中其中jjjkjlxkjkkklljljllb bb bb b1b bb bb b4b bb bb beK2TeyDudxdyyeeeUKU同理同理jjjkjlykjkkklljljllc cc cc c1c cc cc c4c cc cc ceK jkeluuuU计算物理学Harbin Institute of Technology Yangkun 46/1045.2 基于变分原理的有限元法基于变分原理的有限元法三角元三角元e的泛函的泛函 TxTTyeJu eeeeeeeeeKUUUUKUUKeeejjjkjleeekjk

45、kkleeeljlkllKKKKKKKKKeK其中其中eerssrrsrs1KKb bc c4jjjkjlxkjkkklljljllb bb bb b1b bb bb b4b bb bb beKjjjkjlykjkkklljljllc cc cc c1c cc cc c4c cc cc ceKr sj,k,l、计算物理学Harbin Institute of Technology Yangkun 47/1045.2 基于变分原理的有限元法基于变分原理的有限元法改写改写Ke到所有到所有n个节点，即把扩充部分添零，以方便总体个节点，即把扩充部分添零，以方便总体矩阵的处理矩阵的处理 TeJu eUK

46、U其中其中eeejjjkjleeekjkkkleeeljlkll000KKK0KKK0KKK000eK123nuuUuu mee 1J uJu计算物理学Harbin Institute of Technology Yangkun 48/1045.2 基于变分原理的有限元法基于变分原理的有限元法求解区域上的总体泛函求解区域上的总体泛函 mmTTee 1e 1J uJueUKUUKU其中其中meijije 1 KKi, j1,2, n 变分问题被离散化的多元二次函数的极值问题变分问题被离散化的多元二次函数的极值问题 T12nnnijiji 1 j 1J uJ u ,u ,u K u uminUKU

47、5. 建立系统有限元方程建立系统有限元方程计算物理学Harbin Institute of Technology Yangkun 49/1045.2 基于变分原理的有限元法基于变分原理的有限元法根据多元函数极值理论根据多元函数极值理论iJ0 i1,2,nu得到第得到第i点有限元方程点有限元方程nijjj 1K u0即即 0KU求解上述有限元方程（线性代数方程组），就可以得到求解上述有限元方程（线性代数方程组），就可以得到节点上的函数值。节点上的函数值。 计算物理学Harbin Institute of Technology Yangkun 50/1045.2 基于变分原理的有限元法基于变分原理

48、的有限元法 获得有限元方程之后，就可以选择各种方法求解相应获得有限元方程之后，就可以选择各种方法求解相应的代数方程组，常用方法有高斯消去法、列元素消去法、的代数方程组，常用方法有高斯消去法、列元素消去法、迭代法等等。迭代法等等。 在变分问题中第二类、第三类边界条件已经自然包在变分问题中第二类、第三类边界条件已经自然包含在泛函达到极值的要求中，不必单独处理，称为含在泛函达到极值的要求中，不必单独处理，称为自然自然满足的边界条件满足的边界条件，只需考虑，只需考虑第一类强加边界条件第一类强加边界条件，强加边界条件的处理方法因代数方程组的解法而异。强加边界条件的处理方法因代数方程组的解法而异。6. 有

49、限元方程求解与边界条件处理有限元方程求解与边界条件处理计算物理学Harbin Institute of Technology Yangkun 51/1045.2 基于变分原理的有限元法基于变分原理的有限元法迭代法求解迭代法求解：凡是遇到边界节点所对应的方程均不迭：凡是遇到边界节点所对应的方程均不迭代，节点值始终保持给定值，不必单独处理边界。代，节点值始终保持给定值，不必单独处理边界。直接法求解直接法求解：节点节点m为边界，函数值为边界，函数值um=u0，处理方法，处理方法为，把对角元素的特征元素设置为为，把对角元素的特征元素设置为1，即，即kmm=1,然后把然后把m行与行与m列的其它元素全部设

50、置为列的其它元素全部设置为0，方程的，方程的等式右边等式右边改改为给定的函数值为给定的函数值u0，其它元素则要减去该节点处理前对，其它元素则要减去该节点处理前对应的应的m列的特征系数列的特征系数kim与与u0的乘积。的乘积。1111n11mn1nnnKKKu0KKKu0KKKu0mmmmmnnm111n110m0n1nnn0KKu001u00KKumnmK uuK u计算物理学Harbin Institute of Technology Yangkun 52/1045.2 基于变分原理的有限元法基于变分原理的有限元法 例例 5.2.1 一个边长为一个边长为1的二维正方形静电场域，电位函数的二维

51、正方形静电场域，电位函数为为(x,y)，边界条件如图所示，试用有限元法确定二维静，边界条件如图所示，试用有限元法确定二维静电场域的电位分布。电场域的电位分布。解：该二维静电场域的电位函数解：该二维静电场域的电位函数(x,y),可以用下列第一类边界条可以用下列第一类边界条件的偏微分方程描述：件的偏微分方程描述：2222010 y10 x 10 xSave As直接生成直接生成M代码。代码。计算物理学Harbin Institute of Technology Yangkun 69/1045.3 matlab有限元法工具箱有限元法工具箱PDE Toolbox求解偏微分方程类型求解偏微分方程类型 1

52、 椭圆型方程（椭圆型方程（Elliptic ）( )*div c grad ua uf2 抛物线型方程（抛物线型方程（Parabolic）*( )*d udiv cgrad ua uf3 双曲型方程（双曲型方程（Hyperbolic ）*( )*d udiv cgrad ua uf计算物理学Harbin Institute of Technology Yangkun 70/1045.3 matlab有限元法工具箱有限元法工具箱 4 特征值方程（特征值方程（Eigenmodes）( )*div cgrad ua ud u上述微分方程中上述微分方程中121212( , )( ),nnnuu x x

53、x tgrad uuxxxdivxxxc、a、d、f在在椭圆型方程椭圆型方程中可以为中可以为函数函数，但在其它方，但在其它方程中必须为程中必须为常数常数。计算物理学Harbin Institute of Technology Yangkun 71/1045.3 matlab有限元法工具箱有限元法工具箱PDE Toolbox 边界条件边界条件 1 狄里赫利条件（狄里赫利条件（Didchlet）hur2 诺依曼条件（诺依曼条件（Neumann）*( )*n c grad uq ugn为边界上的单位外法线矢量，为边界上的单位外法线矢量，h、r、q、g可以为函数可以为函数计算物理学Harbin Ins

54、titute of Technology Yangkun 72/1045.3 matlab有限元法工具箱有限元法工具箱PDE Toolbox 启动启动 1 启动启动2 界面界面计算物理学Harbin Institute of Technology Yangkun 73/1045.3 matlab有限元法工具箱有限元法工具箱PDE Toolbox 菜单菜单 Options打开或关闭栅格打开或关闭栅格调整栅格大小调整栅格大小打开或关闭捕捉栅格功能打开或关闭捕捉栅格功能绘图轴的坐标范围绘图轴的坐标范围打开或关闭绘图方轴打开或关闭绘图方轴关闭帮助信息关闭帮助信息图形缩放图形缩放选择应用模式选择应用模式

55、重新显示图形重新显示图形计算物理学Harbin Institute of Technology Yangkun 74/1045.3 matlab有限元法工具箱有限元法工具箱PDE Toolbox 菜单菜单 Draw进入绘图模式进入绘图模式对角点绘矩形对角点绘矩形固定中心绘矩形固定中心绘矩形矩形对角点绘椭圆矩形对角点绘椭圆固定中心绘椭圆固定中心绘椭圆绘多边形绘多边形旋转已选图形旋转已选图形将几何描述矩阵输出到主工作空间将几何描述矩阵输出到主工作空间计算物理学Harbin Institute of Technology Yangkun 75/1045.3 matlab有限元法工具箱有限元法工具箱P

56、DE Toolbox 菜单菜单 Boundary进入边界模式进入边界模式对已选边界输入条件对已选边界输入条件显示边界区域标识开关显示边界区域标识开关显示子区域标识开关显示子区域标识开关删除已选的子域边界删除已选的子域边界删除所有删除所有 的子域边界的子域边界将分解几何矩阵、边界条件矩阵输将分解几何矩阵、边界条件矩阵输出到主工作空间出到主工作空间计算物理学Harbin Institute of Technology Yangkun 76/1045.3 matlab有限元法工具箱有限元法工具箱PDE Toolbox 菜单菜单 PDE进入偏微分方程模式进入偏微分方程模式显示子区域标识开关显示子区域标

57、识开关调整调整PDE参数和类型参数和类型将将PDE参数输出到主工作空间参数输出到主工作空间计算物理学Harbin Institute of Technology Yangkun 77/1045.3 matlab有限元法工具箱有限元法工具箱PDE Toolbox 菜单菜单 Mesh输入网格模式输入网格模式初始化三角形网格初始化三角形网格加密当前三角形网格加密当前三角形网格优化网格优化网格退回上一步退回上一步用数字化的颜色显示网格质量，大于用数字化的颜色显示网格质量，大于0.6可接受可接受显示网格节点标识显示网格节点标识显示三角形网格标识显示三角形网格标识修改网格生成参数修改网格生成参数输出网格矩

58、阵到主工作空间输出网格矩阵到主工作空间计算物理学Harbin Institute of Technology Yangkun 78/1045.3 matlab有限元法工具箱有限元法工具箱PDE Toolbox 菜单菜单 Solve对已经定义的偏微分方程求解对已经定义的偏微分方程求解调整解调整解PDE的参数的参数输出解到主工作空间输出解到主工作空间Plot显示图形解显示图形解绘图参数设置绘图参数设置输出动画输出动画计算物理学Harbin Institute of Technology Yangkun 79/1045.3 matlab有限元法工具箱有限元法工具箱PDE Toolbox 求解步骤求解

59、步骤 求解区域设置求解区域设置应用模式设置应用模式设置输入边界条件输入边界条件微分方程参数设定微分方程参数设定网格剖分网格剖分初值和误差设置初值和误差设置解方程解方程图形解显示参数设置图形解显示参数设置File-Save As直接生成直接生成M代码代码计算物理学Harbin Institute of Technology Yangkun 80/1045.3 matlab有限元法工具箱有限元法工具箱例例5.3.1 如图带有矩形孔（如图带有矩形孔（0.1*0.8）的金属板（）的金属板（1*1.6），），金属板左侧保持在金属板左侧保持在100，右侧热量可以向环境定常流动，右侧热量可以向环境定常流动，上下侧及内孔保持绝热，初始温度为上下侧及内孔保持绝热，初始温度为0。求。求t=0.1、0.3、0.5、1.5s时金属板温度分布时金属板温度分布 解：此问题可以表示为如下定解问题解：此问题可以表示为如下定解问题222200100 1 0 0tuuutxyuununu 左边界右边界其它边界计算物理学Harbin Institute of Technology Yangkun 81/1045.3 matlab有限元法工具箱有限元法工具箱 求解区域设置求解区域设置