有限元方法6

1 §6 有限元方法 出发点:变分原理： ①虚功原理 ②最小位能原理 ( , ) ( ) u V D u v F v v V    1( ) min ( ) ( ) ( , ) ( ) 2 J u J u J v D v v F v u V      为了导出近似计算方法，核心是构造 的近似空 间。最自然的想法：多项式子空间。对于一维 问题效果较好，但是对于二维以上问题处理不 利： ①边界条件难以满足，甚至近似满足都比较困 难； ②稳定性差； ③“刚性”程度过高。 由此我们需要有限元法。 V 2 本章摘要 1．一维问题的有限元方法——线性元 2．若干理论探讨 3．二维问题 1．一维问题的有限元方法——线性元 a.研究对象 b.单元剖分及试探函数空间的构造 c.有限元方程的形成 3 a.研究对象： 等价形式 , , ( ) 0, ( ) ( ) ( ) ( ) d dup qu f a x b dx dx duu a p b b u b g b dx               ( , ) ( ) u V D u v F v v V    （更弱的 只需 ）  ' '( , ) ( ) ( ) ( ) ( ) b a b a D u v puv quv dx u b v b F v fvdx gv b         1[ , ], ( ) 0V v C a b v a    2' 2: , ( ) 0baV v v v dx v a       4 b.单元剖分及试探函数空间的构造 区间[a,b]分成N个小单元： ，长度 ， 称为该剖 分的直径。 基于此，试探函数空间 的特征如下： 0 1 Na x x x b     1[ , ]i i ie x x 1i i ih x x   max iih h hU (1) 是分段线性的; (2) 是 上的连续函数。 这就是Finite element Method. h hv U  hv hv [ , ]a b  , ( ) 0h h h hV v U v a   ( , ) ( ), h h h h h h h u V D u v F v v V     5 先来具体分析一下 的结构 ① 以局部观点考察， 在某一单元 是 怎样的。 hU h hv U ie 1 1 1 1 1 ( ) ( ) ( ) ( ) ( ) i i i h h i h ie i i i i h i i h i i x x x xv x v x v x x x x x v x N v x M            ② 以整体看，设 代表第i个顶点取值1， 其它顶点取值为0的分段线性连续函 数，则 i 0 ( ) ( ) ( ) N h h i i i i v x v x x   6 c.有限元方程的形成 有关符号： ， ， 即， 为了研究方便 设  0 1, , ,h NU span     h hv U  0 ( ) ( ) ( ) N h h i i i i v x v x x     , ( ) 0h h h hV v U v a   0( ) 0hv x     1( ), , ( ) Th h Nv x v x v 找 使 ， 即  iu 1 N h i iu u  v ( , ) ( )h h hD u v F v  , ( )i i j j j jD u v v F     7 设 ，则 以上方法称为直接形成代数问题法，在实 际工程计算中，采用所谓子结构法完成这 部分工作。 ( )ij N NA a  ( , ) ( , )Au v F v  Au F 以上方法称为直接形成代数问题法，在实际工程 计算中，采用所谓子结构法完成这部分工作。 原理：从局部着手： 因此整体二次型是由限制在每个单元上的局部二 次型进行迭加而形成的，按这种思想获得线性代 数问题即所谓子结构法。 ' ' 1 ( , ) ( ) ( ) ( ) i N h h h h h h h h i e D u v pu v qu v dx u b v b     8 (1)单元“刚度”阵和“荷载”向量的计算 (2) 总刚度矩阵和总荷载向量的形成 (3) 约束条件的处理 (1)单元“刚度”阵和“荷载”向量的计算 如何用二次型表示 ， ， ＝ ＝… ＝ ' '( ) i h h h h e pu v qu v dx ( ) ( )h hu b v b 1h i i i iu u N u M  1 1h i i i iv v N v M    ' '( ) i h h h h e pu v qu v dx   1 1 1 1 1( )( ) ( )( ) i i x i i i i i i i i i i i i i i i i x pu N uM v N vM qu N uM v N vM dx                ([ ] , ) i i i e e e K u v 9 其中 再看右端， ，   1( ) ( )i h i e h i v x v v x      1 10 0( ) ( ) , 0 1 N N h h N N N N u v u b v b u v u v                        1 ( ) ( ) i N h h h i e F v fv dx v b    其中 ，     1 1( ) ( , ) i i ii i x h i i i i e ex e fv dx f v N v M dx F v          1 1 ,i i i i i Tx x i ie x x F fN dx fM dx       0( ) , 1 Nh e v b v         10 (2) 总刚度矩阵和总荷载向量的形成 小刚度 小二次型 自然为整体二次型 大刚度阵 再做迭加即得整体大刚度阵。 同理可形成总荷载向量。   (3) 约束条件的处理 设以上形成的总刚度阵为 ，总荷载 为 ，则 ， ， 但是有很多问题存在边界约束条件，例如 本节介绍的例子： [ ]K  F     ( , ) ( , )h hD u v K u v    0 1, , TNv v v v     ( ) ( , )hF v F v 0( ) 0 0hv a v   11 其中 只需划去相应的行和列即可。同理对于 也应划去相应的行。   1 1 2 2 1 0 0 ( , ) ( , ) ( , )h h N N u v u v D u v K K u v                                               1 * K K         F 2．若干理论探讨 （1） 可解性 ， 是否可逆？ 回忆： Ax b A 0( ) 0p x p  ( ) 0q x  12 定理： 为对称正定阵。 证明：对称性： 而 ， 故 正定性： 取 ，则由构造知 A ( , )ij i ja D   ( , ) ( , )D u v D v u ij jia a  \ 0Nx R  ( , ) 0Ax x  1 N h i i i u x   ( , ) ( , )h hD u v Ax x 而 等号成立： 要使 ，必须使： 这说明 在每一个单元上都为常数，而 又 为上连续可微，故 ， 但   ' ' 1 ' 2 ' 0 ( , ) ( ) ( ) ( ) ( ) 0 i N h h h h h h h h i e b b h ha a D u v pu v qu v dx u b v b p u dz p u dz            ( , ) 0Ax x   2' '0 0ib h h ea u dz u   hu hu [ , ]a b hu C ( ) 0hu a  0 0C x    13 （2） 收敛性分析 在空间 中引入范数（可证明） ， 引理1： 存在正常数 ，使得： V   12 2'1 bav v dx  v V 1 2,M M 2 1 1 2 1 1 ( , ), ( , ) , , M v D v v v V D u v M u v u v V     第一步.给出 中 和 之关系V 0 1        2 2 2 ' 0 2' 2' 22 1 10 1 ( ) ( ) 1 ( ) ( ) ( ) 1 2 b b x a a a b x x a a a b b a a v v x dx v s ds dx ds v s ds dx x a dx v s ds b a v v c v v V                    14 第二步. 取  ' ' 22 22 2 0 0 1, ( , ) ( )b ba av V D v v pv qv dx v b p v dx p v            ' ' ' 2 3 ' ' 2 3 0 00 0 2 2 3 1 1 1 1 1 2 1 1 , , ( , ) ( ) ( ) ' ( ) ( ) ( ) ( ) ( ) b h a b b a a u v V D u v pu v quv dx u b v b c u v dx c uv dx u b v b c u v c u v u b v b c c c u v b a u v M u v                        1 0M p 引理2： 设 为原来问题之解， 为有限元 解，则 ， 、 分别满足： u hu inf h h v V u u c u v      u hu ( , ) ( ) u V D u v F v v V    ( , ) ( ) h h h h h h h u V D u v F v v V    15 证明： 由 ，故 ，hV V h hv V  ( , ) 0 ( , ) 0 ( , ) ( , ) h h h h h h h h h D u u v v V D u v v u v D u v v D v u v             特取： 知 又左边 ， 右边 h hv v u  ( , ) ( , )h h hD v u v u D u v v u      2 1 1h M v u  2 1 1h M u v v u    2 1 1 1 h Mv u u v M     16 1 1 1 1 2 1 1 1 (1 ) h h h u u u v v u u v v u Mu v M C u v               的选取：显然 的信息由 在单元顶点 上的值完全确定， 自然地： 因此 为插值算子 关键点： 如何估计 v vv 0 1, , , Nx x x  ( ) 0,1,2, ,i iv x u x i N   hu v u  1h u u 17 引理3： 设 ，则 ， 这里 。 这里 为网格剖分地最大直径， 为与 无关的常数。 2[ , ]u V C a b  41 2hu u C h u  1 22'' 2 1 N b a i u u dx      h 4C h 证明：                      ' ' 1 1 1 1 2 2 1 1 2 1' 1 1 2' ' 1 1 2 '' 1 , i i i i i i i i i Nb x h h ha x i N x i i x i i i N x i i i ix i N x x x i u u u u dx u u dx u x u x u x dx x x u x u dx x x u s ds dx                                  18  '' 1 '' 1 1 '' 1 '' 1 '' 22 1 2 1 2 22 1 1 2 2 1 22 2 2 4 1 1 1 2 2 i i i i i i i i i i i i N x x x x i N x x ix x i N x i i i ix i N x x i b a ds u ds dx u dx x dx u dx x x u dx h h u dx C h                                     定理： 有限元解 满足估 计 。 证明：由引理2、引理3即知。 举例： hu 31 2hu u M h u  '' ' ( ) 2 (0,1) (0) 0, (1) 0 u x x u u      19 第一步，变分问题： 求 使 对一切 成立。u V ( , ) ( )D u v F v v V   '1 2 20| , (0) 0V v v v dx v     1 1' ' 0 0 ( , ) , ( ) 2D u v u v dx F v vdx   第二步，有限元空间的构造 ， 有限元问题： 4, iN x ih   : , [0,1], (0) 0ih eV v v v C v    ( , ) h hh h h h h u V D u v F v v V    20 第三步，刚度矩阵及荷载向量的形成 a．局部刚度矩阵与局部荷载向量 其中     1 ' '' ' 1 1 i i i x i i i i i i i ix e u v dx u N u M v N v N dx       ' /1 1 1, , ,i ii i i i x x x xN M N M h h h h        上式    1 1 12 1 i i x i i i ix u u v v dx h        1 1 1 1 11 i i i i i i i iu v u v u v u vh           1 11 1 1ie K h      12 i i i e vdx hv hv    1 1ie h F h h            21 b．整体刚度矩阵和整体荷载的形成   1 1 0 1 1 1 0 0 ie K h          1 1 0 0 ie F h            全部加起来：   5 5 1 1 1 2 1 1 1 2 1 1 2 1 1 1 K h                 1 2 2 2 1 ie F h           22 c．约束条件的处理：将约束条件的分量对 应的行列去掉 2 1 1 2 11 1 2 1 1 1 A h           2 2 2 1 h b         第四步，求线性代数方程组 求出的是一个函数： 。 Au b 1 2 3 4 u u u u u         1 1 2 2 3 3 4 4hu u u u u       23 3． 二维问题 假设 ，且 不恒为，唯一性由极值 原理易证。 对于 ，解不一定存在， 要有解则须考虑 ？＝ 。 u f u u g n      ( ) 0x  ( )x u f u g n      ( ) 1u dxdy    fdx   左边 ， 即须满足 。 从而 是该问题有解的充分必要条 件。 1 1 u uu dxdy ds ds gds n n                 fdxdy gds      0fdx gds      24 等价变分问题：求 使u V ( , ) ( )D u v F v v V  22 2: ( , ) ( ) v vV v v dxdy x y D u v u vdxdy uvds F v fvdxdy gvds                                         a.单元划分及试探函数空间的形成 为研究方便，设 为平面多角形区域，首 先对区域 进行三角剖分，满足： 1. 任一三角形单元只能与相邻三角形单元 共边或共顶点； 2. 每一个三角形单元不可太尖或太钝，大 小亦可不相差太大； 3. 多角形区域的角点均应为单元顶点； 4. 三角形单元充满整体多角形区域。   25 对于奇异情形另当别论。 对单元编号：设只有个单元， 记为 ( )。 对顶点编号：设只有个顶点， 记为 。 三角形单元的最大直径为 。 ke 1,2, ,k NE  ( , ), 0,1, ,i i iP x y i NP   h 定义：有限元空间 为基于以上三角单元 剖分的分片线性连续函数全体所构成的 有限维空间。（ ） 下面是的具体结构： hU dim hU NP 26 1. 从整体看，设 为当顶点 取值为1， 其他顶点取值为0的分片连续函数 ，时 ， ：山形函数，非零域很小。 i iP 0,1, ,i NP    1 NP i i i v v P     i 2. 局部看，考察在任一单元上的结构，设 （逆时间方向）， ， 满足条件： i j me PP P ev ax by c   ( ) 1 ( ) 1 ( ) 1 i i i i i i i j j j j j j j m m m m m m m ax by c v P v x y a v ax by c v P v x y b v ax by c v P v x y c v                                   27 记： 为 面积， 形成单元基函数 1 1 1 2 1 i i j j m m x y e x y x y   1 1 1 1 11 , 1 , 2 2 2 1 1 i i i i i i i j j j j j j j m m m m m m m v y x v x y v a v y b x v c x y v e e e v y x v x y v       , , 为线性基函数。      1 1 12 2 2i i j j m m i i j j m m i i j j m me i i j j m m v av av a v x bv bv b v y cv cv c v e e e Nv Nv N v                 iN jN mN 28       1 2 1 2 1 2 i i i i j j j j m m m m N a x b y c e N a x b y c e N a x b y c e              性质： ，其中 重新记： 。 iN 1i j m i i j j m m i i i j m m N N N x N x N x N x y N y N y N y           1 ( ) 0l k lk l k N v l k      , , ,k k k kN a x b y c k i j m   29 b．有限元方程的形成 如何得到 ？( , ) ( )h h hD u v F v Ax b ( , ) ( ) ( )h h h h h hD u v u v dxdy u s v s ds        ( ) ( ) ( ) n h h h e F v fv dxdy g s v s ds       A.局部刚度矩阵和荷载向量的形成 B.形成整体刚度阵和整体荷载向量 C.约束条件处理 D.计算过程的若干处理 E.若干理论探讨 30 A．局部刚度矩阵和荷载向量的形成     n n n h h i i j j m m i i j j m m e e i i j j m m i i j j m m i i j j m m i i j j m me T i i j m i j m j i j m i j m m u v dxdy N u N u N u N v N v N v dxdy a u a u a u a v a v a v dxdy bu b u b u b v b v b v ua a a a a a u b b b b b b u                                                            ,n nh i j n m T Ti i i j m i j m j j n i j m i j m m m ne ee v v e v u va a a a a a u v e b b b b b b u v K u v e                                                其中：    , , n n n T i j m i j m e i j m i j m i i j je e m m a a a a a a K b b b b b b u v u u v v u v                                           , n n n n i i j j m m e e e e fvdxdy f N v N v N v dxdy F v     31 还须处理由单元边界导致的局部刚度阵和 荷载向量，只需计算位于求解区域边界 的单元上的量。 为对应单元 的边， 且有： n ne n n uvds uvds       引入局部坐标的边上点至 上的距离 为参 数考虑 则： ，     * n n i i j j m m i i j j m muvds N u N u N u N v N v N v ds          iP t 0 ,i jt P t l P    , , 0ji m i j i j i j NN Nl t t PP l PP l PP    32 由 、 、 函数的定义： 如果的边界无贡献，后一项取消。 iN jN mN               0 0 * , , n n h n n n n n n n n n l i j i j e ee l i i ee e e e e e e l t t l t tu u v v dt l l l l K u v l t tgvds g v v dt l l F v K K K F F F                                                      B.形成整体刚度阵和整体荷载向量  （略） 33 C.约束条件处理： 注：对于有约束问题，一般先对约束边界 点编号，这样处理时只需划掉开始的有 关行和列。 , 0 0 h h u f v v u       D．计算过程的若干处理   , 0 k x y T f T      ( , ) ( , ) 0 T Tk x y k x y f x x y y T                       0V v   34 单元刚度阵： ( , ) ( ) D u v k u vdxdy F v fvdxdy          ne k u vdxdy  求解积分项 ①数值积分 ②对 为多项式情形，可通过引入重心 坐标来解决 记 ( , )k x y 3 2 1 3 1 2 1 2 3 1 2 3 1 2 3 1 2 3 , ,PPP PPP PPP PP P PP P PP P          35 显然 作这种变换 把一个任意三角形化 成一个等腰直角三角形（标准参考元） 因此    1 2 3 1 21 , ,x y            1 2, ,x y     1 2 3 1 2 1 2 3 1 2( , ) ( , , ) 2 ( , , ) ne e e k x y dxdy k Jd d e k d d                ：Jacobi 行列式 ，面积之比：J 21 2 e e             31 2 1 31 2 2 3 31 2 1 2 1 2 1 2 1 1 1 1 2 1 2 20 0 1 11 1 1 10 0 1 2 3 1 2 3 1 1 1 1 ! ! ! 2 ! e d d d d d t t dt                                                     36 E.若干理论探讨 问题： 是否非奇异，要看 ①    , ,h hD u v Ax y Ax b   A  ,h hD u v 0 u f u      : 0 , ( , )V v v D u v u vdxdy        是正定的 ②    ( , ) , , 0h hD u v Ax y Ax x   ( , ) 0h hD u v   22 ( , ) 0 : D u v u vdxdyu f u V v v v dxdyn                         37 奇异    ( , ) , , , 0 ( , ) 0D u v Ax y Ax x D u u      2, 0 0 0D u u u dxdy u         A 2．矩阵带宽 看 的支集的交集情况可初 步判断 是否取零值。  ,ij i ja D   ,i j  ija