有限元方法

有限元方法2有限元法是求解偏微分方程问题的一种重要数值方法，它的基础分两个方面：一是变分原理，二是剖分插值．从第一方面看，有限元法是Ritz-Galerkin方法的一种变形．它提供了一种选取“局部基函数”的新技巧，从而克服了Ritz-Galerkin方法选取基函数的固有困难．从第二方面看，它是差分方法的一种变形．差分法是点近似，它只考虑在有限个离散点上函数值，而不考虑在点的邻域函数值如何变化；有限元方法考虑的是分段（块）的近似．因此有限元方法是这两类方法相结合，取长补短而进一步发展了的结果．在几何和物理条件比较复杂的问题中，有限元方法比差分方法有更广泛的适应性．3§7.两点边值问题的有限元方法本节以两点边值问题为例，并从Ritz法和Galerkin法两种观点出发来叙述有限元法的基本思想及解题过程.7.1基于Ritz法的有限元方程考虑两点边值问题其中，(7.2)0)(,0)((7.1),,)(buaubxafqudxdupdxdLu1,,0,,,0,,pxCabpqCabqfCab41.写出Ritz形式的变分问题与边值问题(7.1)、(7.2)等价的变分问题是：求，使其中，（7.3）式(7.3)是应用有限元法求解边值问题(7.1)、(7.2)的出发点.1*EuH1*minEuHJuJu1,,2Juauufu,,.bbaadudvauvpquvdxfufudxdxdx，52.区域剖分剖分原则与差分法相同，即将求解区域剖分成若干个互相连接，且不重叠的子区域，这些子区域称为单元．单元的几何形状可以人为选取，一般是规则的，但形状与大小可以不同．对于一维情形最为简单.将求解区间分成若干个子区间，其节点为每个单元的长度为.[,]ab01inaxxxxb1[,]iiiexx1iiihxx单元在区间中分布的疏密程度或单元尺寸的大小，可根据问题的物理性质来决定，一般来说，在物理量变化剧烈的地方，单元尺寸要相对小一些，排列要密一些.6设为的有限维子空间，它的元素为.要构造，只需构造单元基函数.构造单元基函数所遵循的原则是：其中，是单元节点序号为的节点.hV1EH)(xuhhVi(1)每个单元中的基函数的个数和单元中的节点数相同，每个节点对应一个基函数，本例中，单元有两个节点，因此基函数有两个.ie(2)基函数应具有性质,,0,,1)(kjkjxjkkjkxk73.确定单元基函数有限元法与Ritz-Galerkin方法的主要区别之一，就在于有限元方法中的基函数是在单元中选取的．由于各个单元具有规则的几何形状，而且可以不必考虑边界条件的影响，因此在单元中选取基函数可遵循一定的法则．8函数取为中的基原则，可将为线性函数，则按上述若取hiVx)()4.7(.,0,,)(,1,,2,1.,0,,,,)(1111111在别处在别处nnnnniiiiiiiiixxxhxxxnixxxhxxxxxhxxx线性组合，即的可以表示为基函数中任一函数显然，)(xuVihh)()()(2211xuxuxuunnh9),,,2,1()(,,,21niuxuuuuuiihhn在节点上的值，即是其中，表示为上，在单元)(xuehi)5.7(].,[,)()()(11111iiiiiiiiiiiihxxxhxxuhxxuxuxuxu.叠加而成由各个单元的近似函数生，全区域的近似函数合产数由单元基函数线性组可见，单元中的近似函10的集合：函数是满足下列条件的所有由以上可以看出，hhuV,0)()3();,,2,1(1)2(];,[,],[)1(2aunieubaLuubauhihhhh的多项式上是次数不超过在上连续，且在.1称为试探函数数空间，维子空间，称为试探函的一个是故hhEhVunHV11(7.6)令4.形成有限元方程，得代入问题，将式的极小上解泛函，在替代法一样，以与)3.7()4.7()3.7(Ritz1hEhVHV1()(,)(,)2hhhhJuauufu,111(,)(,)2nnijijjjijjauuuf()0hjJuu便得到确定的线性代数方程组12,,,nuuu1(,)(,),1,2,,nijijiaufjn称式（7.5）为有限元方程.12(7.8)(7.7)值得注意的是，在实际计算中，并不是按照上述步骤形成有限元方程的，而是先进行单元分析，即在单元上建立有限元特征式，然后再进行总体合成，即将各单元的有限元特征式进行累加，合成为有限元方程．具体过程如下：第一步：单元分析.注意到221()(2)2bhhhhaJupuqufudx1122111()2iiiinnxxhhhxxiipuqudxfudx作变换1iixxh13(7.10)并引入记号其中，.于是或写成(7.9)01()1,()NN可写成上，则在单元hiiiuxxe],[1101101()()()(,)ihiiiuuxNuNuNNu()ihuNu()011(,),(,)iTiiNNuuNu()11()()ihiiiuxuuhMu其中，.从而有(1/,1/)iihhM(7.11)这里(7.12)()()1,1,()()()0,1()iiiiiiiiiTTiiiiiiiaaKhpqdaaMMNN，称为单元刚度矩阵，其中(7.13)(7.16)(7.14)式中(7.15)将式（7.11）、（7.14）代入式（7.7），便有对式（7.7）右端第二项积分，有11()()()10()()(),iixiTiTihiiixfudxhfxhdNuuF101)(101)(1T)()(1)()()1)((,,dhxfhFdhxfhFFFiiiiiiiiiiiiiiiF.)(单元荷载向量为称iF111.2nnTTiiiiihiiJuuKuuF这样，我们就得到了单元有限元特征式的一般表示形式：()()()iiiKuF于是有第二步：总体合成.总体合成就是将单元上的有限元特征式进行累加，合成为总体有限元方程.这一过程实际上是将单元有限元特征式中的系数矩阵（称为单元刚度矩阵）逐个累加，合成为总体系数矩阵（称为总刚度矩阵）；同时将右端单元荷载向量逐个累加，合成为总荷载向量，从而得到关于的线性代数方程组．为此，记()()1(,)iTiiiuuuBu从而式(7.16)右端第一个和式为11111[()]222nnTiiiiiiTTTiiuKuuBKBuuKu，(7.17)其中（未标明的元素均为0）这就是总刚度矩阵．对式(7.16)右端第二个和式，有其中这就是总荷载向量.从总刚度矩阵和总荷载向量的形成过程可以看出，的计算，实际上是把中四个元素在适当的位置上“对号入座”地叠加，的计算也是如此．我们引入，只是为了叙述方便，实际上，在编制程序时并不需要．显然，方程组(7.18)的系数矩阵是对称正定的三对角矩阵，因此可采用追赶法求出在节点上的近似值.(7.18)其这样，就可将式(7.16)写成因此，有限元方程为K)(iKb)(iBKunuuu,,,2119§7.两点边值问题的有限元方法本节以两点边值问题为例，并从Ritz法和Galerkin法两种观点出发来叙述有限元法的基本思想及解题过程.7.1基于Ritz法的有限元方程7.2基于Galerkin法的有限元方程从Galerkin法出发形成有限元方程的过程与前面完全一样，针对边值问题(7.1)、(7.2)所得到的结果也是一致的．但是从Galerkin法出发形成的有限元方程更具一般性，它不仅适用于对称正定的算子方程，而且也适用于非对称正定的算子方程，所以我们今后主要是依据这一观点建立有限元方程．20与边值问题(7.1)、(7.2)等价的Galerkin变分问题是：求，使得(7.19)其中仍用分段线性函数构成的试探函数空间替代，将代入(7.19)，则得到所满足的线性代数方程组(7.20)这和方程组(7.6)是完全一样的.1EuH1,,0,EauvfvvH,,.bbaadudvauvpquvdxfufudxdxdx，hV1EH11,nnhiihiiiiuuxvvx12,,,nuuu1,,,1,2,,.nijijiaufjn与容易看出，方程组(7.20)的系数矩阵就是总刚度矩阵．在总刚度矩阵形成的过程中，注意到(7.21)而从而有即故有这就是有限元方程(7.18).由上述看出，按Galerkin法推导有限元方程更加直接方便．尤其重要的是．按这一观点推导的有限元方程，不仅适用于定常的微分方程定解问题，而且也适用于不定常的微分方程定解问题，因此具有广泛的适应性．例7.1用有限元方法解边值问题将区间[0,1]等分成4个单元.解利用上述分析结果，我们只需构造出单元刚度矩阵和单元荷载向量，然后合成为总刚度矩阵和总荷载向量．注意到(7.13)和(7.15)，并将形成单元上的中点值则不难得到其中，，单元的中点为于是有(),(),()pxqxfx1[,]iixx,,iiipqf1[,]iixx如果把单元刚度矩阵和单元荷载向量“扩大”，便得到和为类似地，可写出和.()iK()iF()iK()iF(3)K(4)K然后进行叠加，便得到总刚度矩阵和总荷载向量：依边界条件即在中划去首末两行和首末两列，在中划去首末两行，便得到如下线性代数方程组：解之，得§8.二维椭圆边值问题的有限元方法用有限元方法求解二维椭圆边值问题的过程与两点边值问题的有限元方法大体相同，只是在具体处理时比一维情形更复杂些．考虑如下椭圆型方程的第一边值问题：(8.1)(8.2)其中，是平面的一个有界域，其边值是分段光滑的简单闭曲线.以下我们从Galerkin法出发，叙述有限元求解问题(8.1)、(8.2)的全过程．与边值问题(8.1)、(8.2)等价的Galerkin变分问题是：求，使得(8.3)其中8.1区域剖分正如前章所言，对高维区域的剖分与对一维区域的剖分有很大不同.对一维区域无论作哪一种剖分，其单元仍然是一个区间，对不同的剖分只是区间长度不同而已.对高维区域而言，不同的剖分其单元的形状各异，如对二维区域，剖分后的子区域可以是三角形、矩形或四边形.限于篇幅，本书只讨论剖分后所得的子区域是三角形的情况，这种剖分称为三角形剖分.将区域划分成有限个三角形单元，剖分方法见前章，那里曾假定剖分的单元应是锐角三角形．现在我们去掉这一限制，只假定不同的单元是无重叠的内部，且单元的顶点不是其它单元边的内点．当然还要尽量避免出现大钝角的三角形．在物理量变化剧烈的地方，单元要划分得细密一些，变化缓和的地方，划分得稀一些．划分好单元之后，要对单元和节点进行编号．设是区域中的单元总数，将全区域中的单元统一编号，单元号记为．全区域中的节点也要按一定的顺序统一编号，记全区域中共有个节点，节点号记为节点编号的一般原则是尽可能使同一单元内的节点号比较接近．以后可以看到，单元内节点序号的差值决定了总体系数矩阵的带宽．8.2确定单元基函数与一维情形一样，为了构造试探函数空间我们只需在每个单元上构造插值基函数．这里，我们仅考虑三角形单元上的线性插值函数．为了便于后面积分的计算，我们先将直角坐标转换为面积坐标.1.面积坐标及有关公式(1)面积坐标的定义设是以为顶点的任意三角形单元，面积为，我们规定的次序按逆