有限元和边界元方法

1、计算物理计算物理3/lesson/ComputationalPhysics3/lesson/ComputationalPhysics有限元和边界元方法有限元和边界元方法有限元和边界元方法有限元和边界元方法n物理问题的变分原理物理问题的变分原理n泊松方程的有限元方法泊松方程的有限元方法n扩散方程的有限元方法扩散方程的有限元方法n波动方程的有限元方法波动方程的有限元方法n边界积分方程边界积分方程n边界元近似边界元近似n单一边界下的边界元法单一边界下的边界元法n两种介质的边界元方法两种介质的边界元方法物理问题的变分原理物理问

2、题的变分原理(1/(1/3 3) )n有限元方法有限元方法n基于变分原理的离散化方法基于变分原理的离散化方法部分逼近地离散化部分逼近地离散化n划分整体区域为有限个基本块划分整体区域为有限个基本块( (单元单元) )n在单元上插值逼近，得到结构简单的函数集在单元上插值逼近，得到结构简单的函数集( (有限元空有限元空间，是泛函间，是泛函 J( (y) ) 的定义域的子集的定义域的子集) )n将边值问题转化为泛函的极值问题将边值问题转化为泛函的极值问题n在有限元空间中寻找泛函在有限元空间中寻找泛函 J( (y) ) 的极小值，作为近似解的极小值，作为近似解n物理中的变分物理中的变分n例：力学体系的最

3、小作用量原理例：力学体系的最小作用量原理n体系的特性可以用拉格朗日体系的特性可以用拉格朗日函数函数 L( (q, , q, , t) ) 描写描写n在在时刻时刻 t1 1 和和 t2 2 之间之间体系按照以下积分取最小值的方式体系按照以下积分取最小值的方式运动运动( (即，运动轨道由泛函的的极小值决定即，运动轨道由泛函的的极小值决定) )得到运动方程，由0d),(21 SttqqLStt物理问题的变分原理物理问题的变分原理(2/(2/3 3) )n例：光学的费马原理例：光学的费马原理n光从点光从点 A 到点到点 B 的传播路径是使光程的传播路径是使光程 L 取极值取极值B A dsnLn由由

4、L 0 0 得到几何光学的折射定律和反射定律得到几何光学的折射定律和反射定律n例：电磁学的麦克斯韦方程组例：电磁学的麦克斯韦方程组n电磁场的电磁场的拉格朗日拉格朗日函数函数 L 是空间积分是空间积分：体积元：标势和磁矢势：电荷和电流密度，：电场和磁场强度导率，：介质的介电常数和磁d , ,d )(21 22AujHE AjHELn电磁学的作用量是时间积分电磁学的作用量是时间积分21 dtttLS 运动方程由泛函的的极小值决定运动方程由泛函的的极小值决定( (即即 S 0 0 ) )物理问题的变分原理物理问题的变分原理(3/(3/3 3) )n例：静电场的泊松方程例：静电场的泊松方程n第一类边界

5、条件第一类边界条件),(),( ),(02222yxuyxuyxfyuxu 等价的变分问题为求解泛函的极值问题等价的变分问题为求解泛函的极值问题),(),( ,dd)()(21)(022yxuyxuyxf uyuxuuJ 泛函的求解必须在边界条件下：条件变分问题泛函的求解必须在边界条件下：条件变分问题n第二类和第三类边界条件第二类和第三类边界条件)( ,2222 unufyuxu 等价的变分问题为求解泛函的极值问题等价的变分问题为求解泛函的极值问题s u uyxf uyuxuuJd)21(dd)()(21)(222 边界条件包含泛函中：自然边界条件边界条件包含泛函中：自然边界条件泊松方程的有限

6、元方法泊松方程的有限元方法(1/(1/1111) )n静电场中二维泊松方程的静电场中二维泊松方程的有限元方法有限元方法2102222, ),(, ),(),(D, ),(21yxyxqnuyxyxuyxuyxyxfyuxu 1 1 2 2ABD 1 1 阳极阳极 2 2 1 1 阴极阴极 2 2n例：例：阴极射线管阴极射线管( (如右图如右图) ) ，在两极上，在两极上( (边边界界 1 1) )的电势的电势 u 是已知的，在左右两侧是已知的，在左右两侧( (边界边界 2 2) )的的 q 是已知的。如果管中的自由是已知的。如果管中的自由电荷密度分布电荷密度分布 ( (x, , y) ) 已知

7、，则已知，则qnuuuyxfyxu21 , ),(),(002 以上的泊松方程以上的泊松方程等价为求解以下泛函等价为求解以下泛函 J( (u) ) 的极值问题的极值问题0D212 ,ddd )(21)(uusq uyxf uuuJ泊松方程的有限元方法泊松方程的有限元方法(2/(2/1111) )n有限元方法的具体步骤有限元方法的具体步骤n划分整体区域为有限个单元和编号划分整体区域为有限个单元和编号n划分要点划分要点n三角形的顶点相连三角形的顶点相连n避免钝角避免钝角( (因引入较大误差因引入较大误差) )n每个三角形不跨越不同的介质每个三角形不跨越不同的介质n每个三角形最多只有一条边每个三角形

8、最多只有一条边在在 2 2 上上( (方便计算方便计算) )n三角形覆盖尽量多的区域三角形覆盖尽量多的区域n编号约定编号约定n三角形单元的编号：三角形单元的编号：e , , , , , ,n顶点顶点的的编号：按逆时针为编号：按逆时针为 1, 1, 2, 2, 3 3n顶点顶点的坐标：的坐标：( (x1 1, , y1 1), ), ( (x2 2, , y2 2), ), ( (x3 3, , y3 3) )n单元的单元的泛函：泛函：Je( (u) )e123en整体的整体的泛函：泛函：1)()(eeuJuJ泊松方程的有限元方法泊松方程的有限元方法(3/(3/1111) )3 (x3, y3)

9、e1 (x1, y1)2 (x2, y2)n构造线性插值函数构造线性插值函数n假设每个单元内假设每个单元内 u( (x, , y) ) 是是 x 和和 y 的的线性线性函数函数n每个每个的的三个基函数三个基函数nu( (x, , y) ) 的的插值表达式中，插值表达式中，a, , b, , c, , d 可由三角形的顶点坐可由三角形的顶点坐标确定，只标确定，只剩余剩余 u1 1, , u2 2, , u3 3 未知未知泊松方程的有限元方法泊松方程的有限元方法(4/(4/1111) )n建立单元的矩阵建立单元的矩阵n单元的泛函单元的泛函123e2eeeesq uf uuuJ2ddd )(21)(

10、2n第一项积分与单元刚度矩阵第一项积分与单元刚度矩阵 ( (z i j) )n第二项积分与单元矩阵第二项积分与单元矩阵 ( (r f j) )n第三项积分与单元矩阵第三项积分与单元矩阵 ( (r q j) )泊松方程的有限元方法泊松方程的有限元方法(5/(5/1111) )n建立顶点和结点的建立顶点和结点的( (V n) )对应关系对应关系n单元编号：有一条边在单元编号：有一条边在 2 2 上且上且 q 0 0 的单元编号为的单元编号为 1, 1, 2, 2, , , e1 1，其，其余的单元编号余的单元编号为为 e1 11,1, e1 12,2, , e0 0 2 2 1 1e1 1e1 1

11、11e0 0n顶点编号：用顶点编号：用 V( (e, , i) ) 表示，逆时针方向，表示，逆时针方向，2 2和和 3 3在在 2 2 上上n结点编号：内部和结点编号：内部和 2 2 上的结点编号为上的结点编号为 1, 1, 2, 2, , , n1 1， 1 1 上的结点编号为上的结点编号为 n1 11,1, n1 12,2, , , n0 0n建立顶点和结点的对应关系：建立顶点和结点的对应关系：V( (e, , i) ) nn集成泛函和建立方程集成泛函和建立方程n泛函的离散化泛函的离散化nK 为总体刚度矩阵，由单元刚度矩阵为总体刚度矩阵，由单元刚度矩阵 ( (z i j) ) 合成合成nR

12、f 由由单元矩阵单元矩阵 ( (r f j) ) 合成，合成，Rq 由由单元矩阵单元矩阵 ( (r q j) ) 合成合成nJ( (u) ) 被离散化为二次多元函数被离散化为二次多元函数 J( (u1 1, , u2 2, , , , un0 0) )泊松方程的有限元方法泊松方程的有限元方法(6/(6/1111) )n有限元方程有限元方程( (关于关于 um 的线性方程组的线性方程组) )n求解方程求解方程泊松方程的有限元方法泊松方程的有限元方法(7/(7/1111) )n例：如右图的边长为例：如右图的边长为 1 1 的正方形区域的正方形区域 1 1 2 2 2 2 2 2xyOyxnuuux

13、yyuxu2 , 5 . 0) 1 , 0( , 4 . 0) 1 , 1 (2222n划分整体区域为有限个单元和编号划分整体区域为有限个单元和编号n单元：单元：n顶点：顶点：1 2 31 2 3( (2 32 3在在 2 2上上) )n结点：结点：( (先先内部和内部和 2 2 ) )1 12 21 11 11 12 22 22 23 33 33 33 3泊松方程的有限元方法泊松方程的有限元方法(8/(8/1111) )n构造线性插值函数构造线性插值函数 1 1 2 2 2 2 2 2xyO1 12 21 11 11 12 22 22 23 33 33 33 3泊松方程的有限元方法泊松方程的

14、有限元方法(9/(9/1111) )n建立单元的矩阵建立单元的矩阵n建立顶点和结点的对应关系建立顶点和结点的对应关系 V ( (e, , i) ) 1 1 2 2 2 2 2 2xyO1 12 21 11 11 12 22 22 23 33 33 33 3泊松方程的有限元方法泊松方程的有限元方法(10/(10/1111) )n集成泛函和建立方程集成泛函和建立方程n求解方程求解方程泊松方程的有限元方法泊松方程的有限元方法(11/(11/1111) )n例：如下图的环形均匀带电板，内径例：如下图的环形均匀带电板，内径 6 6，外径，外径 1010，外圈外圈 1 1 的电势为常数，内圈的电势为常数，

15、内圈 2 2 的电场为常数的电场为常数12 ,100 , 4212222nuuyuxu 1 1 2 26 61010n考虑对称性，取考虑对称性，取 1/41/4 环形环形区域以简化计算区域以简化计算 2 2 1 11 12 23 3扩散方程的有限元方法扩散方程的有限元方法(1/(1/1 1) )n二维二维扩散方程扩散方程0 ),()0 ,( ),(02uyxuyxuyxfuDtun离散化离散化ijjjiiiiiDiDiDiyxyxttyxufuDtu),( , , ),()(),(d dd为基函数n关于关于 i( (t) ) 的常微分方程组的常微分方程组DiiDjiijDjiijf fDKMf

16、KtMd ,d ,dddn初始条件初始条件DiuMd ,)0(0n求解方法：二级欧拉法求解方法：二级欧拉法波动方程的有限元方法波动方程的有限元方法(1/(1/1 1) )n二维二维波动方程波动方程0 ),( ),( ),(000222uyxRtuyxuuyxfuDtuttn离散化离散化ijjjiiiiiDiDiDiyxyxttyxufuDtu),( , , ),()(),(d dd22为基函数n关于关于 i( (t) ) 的常微分方程组的常微分方程组DiiDjiijDjiijf fDKMfKtMd ,d ,ddd22n初始条件初始条件DiDiRrurtMMd ,d ,d)0(d ,)0(0n求

17、解方法：二级欧拉法求解方法：二级欧拉法边界积分方程边界积分方程(1/(1/2 2) )n边界元方法的特点边界元方法的特点n基于基于边界积分方程的近似方法边界积分方程的近似方法n结点仅分布在区域的边界结点仅分布在区域的边界n以边界积分方程为控制方程，将边界离散插值，转化以边界积分方程为控制方程，将边界离散插值，转化为代数方程组为代数方程组n未知量的个数少，求解的计算量少未知量的个数少，求解的计算量少n面积分面积分边界积分边界积分n定理：如果定理：如果 u( (x, , y) ) 和和 v( (x, , y) ) 是定义在平面域是定义在平面域 D 上的两个上的两个任意函数，并且它们在边界外法线上的

18、导数为任意函数，并且它们在边界外法线上的导数为xyOddsnpnvqnu ,supvqvuuvd)(d)( D22那么n证明：利用格林公式证明：利用格林公式( (略略) )边界积分方程边界积分方程(2/(2/2 2) )n 函数函数n定义：定义：1d)( ,0, 10, 0)(xxxxxn性质：性质：)0(d)()(fxxxfn平面域：平面域：的距离到定点为动点 P ),(ln212irrrn边界积分方程边界积分方程nv 和和 p 是已知的，上式给出是已知的，上式给出 D 内任意一点内任意一点( (等等式左边式左边) )与与 上某点上某点( (等式右边等式右边) )的的 u 和和 q 之间的之

19、间的关系关系n上式称为边界积分方程，是边界元法的基础上式称为边界积分方程，是边界元法的基础边界元近似边界元近似(1/(1/4 4) )n边界元方程边界元方程n当定点当定点 i 在在边界边界 上上n常数边界离散化常数边界离散化n边界积分方程的右边：将边界边界积分方程的右边：将边界 分成分成 N 段，以每段中点段，以每段中点的的 u 和和 q 近似该段的函数值近似该段的函数值n例：泊松方程例：泊松方程n当定点当定点 i 不在不在边界边界 上：利用前面的结果，代入最上公式上：利用前面的结果，代入最上公式边界元近似边界元近似(2/(2/4 4) )n对角对角元元 Hii 和和 Gii 的的计算计算(

20、(定点定点 i 在在 i 上上) )n非对角非对角元元 Hij 和和 Gij 的的计算计算( (定点定点 i 不不在在 j 上上) )边界元近似边界元近似(3/(3/4 4) )nBi 的的计算：将区域计算：将区域 D 划分为有限个三角形单元划分为有限个三角形单元eeirfrfBdlndlnDn例：三角区域的电势和电量例：三角区域的电势和电量0 ,2 , 02左边斜边、底边nuyuuxOy11 泊松方程的混合边界问题泊松方程的混合边界问题n常数边界常数边界离散化离散化n对角对角元元 Hii 和和 GiiA 2 2 1 1 1 1n非对角非对角元元 Hij 和和 GijA 2 2 1 1 1 1A 2 2 1 1 1 1A 2 2 1 1 1 1A 2 2 1 1 1 1nBi ：因为：因为 f 0 0，所以，所以 Bi 0 0边界元近似边界元近似(4/(4/4 4) )n构造构造方程方程 AX Rn内部内部( (三角形中心三角形中心(1/3,(1/3, 0)0) )：u 0.6230.623单一边界下的边界元法单一边界下的边界元法(1/(1/6 6) )n例：三角区域的电势和电量例：三角区域的电势和电量0 ,2 , 02左边斜边、底边nuyuun主程序主程序xOy11单一边界下的边界元法单一边界下的边界元法(2/(2/6 6) )n边界边界离