1、3 弹性力学平面问题的有限元法本章包括以下的内容：3.1 弹性力学平面问题的基本方程3.2 单元位移函数3.3 单元载荷移置3.4 单元刚度矩阵3.5 单元刚度矩阵的性质与物理意义3.6 整体分析3.7 约束条件的处理3.8 整体刚度矩阵的特点与存储方法3.9 方程组解法3.1 弹性力学平面问题的基本方程弹性力学是研究弹性体在约束和外载荷作用下应力和变形分布规律的一门学科。在弹性力学中针对微小的单元体建立基本方程，把复杂形状弹性体的受力和变形分析问题归结为偏微分方程组的边值问题。弹性力学的基本方程包括平衡方程、几何方程、物理方程。弹性力学的基本假定如下：1）完全弹性，2）连续，3）均匀，4）各
2、向同性，5）小变形。3.1.1 基本变量弹性力学中的基本变量为体力、面力、应力、位移、应变，各自的定义如下。体力体力是分布在物体体积内的力，例如重力和惯性力。面力面力是分布在物体表面上的力，例如接触压力、流体压力。应力物体受到约束和外力作用，其内部将产生内力。物体内某一点的内力就是应力。图 3.1如图 3.1 假想用通过物体内任意一点 p 的一个截面 mn 将物理分为、两部分。将部分撇开，根据力的平衡原则，部分将在截面 mn 上作用一定的内力。在 mn 截面上取包含 p 点的微小面积 ，作用于 面积上的内力为 。AQ令 无限减小而趋于 p 点时， 的极限 S 就是物体在 p 点的应力。SQA0
3、lim应力 S 在其作用截面上的法向分量称为正应力，用 表示；在作用截面上的切向分量称为剪应力，用 表示。显然，点 p 在不同截面上的应力是不同的。为分析点 p 的应力状态，即通过 p 点的各个截面上的应力的大小和方向，在 p 点取出的一个平行六面体，六面体的各楞边平行于坐标轴。图 3.2将每个上的应力分解为一个正应力和两个剪应力，分别与三个坐标轴平行。用六面体表面的应力分量来表示 p 点的应力状态。应力分量的下标约定如下：第一个下标表示应力的作用面，第二个下标表示应力的作用方向。，第一个下标 x 表示剪应力作用在垂直于 X 轴的面上，第二个下标 y 表示剪应力xy指向 Y 轴方向。正应力由于
4、作用表面与作用方向垂直，用一个下标。 表示正应力作用于垂直于 X 轴x的面上，指向 X 轴方向。应力分量的方向定义如下：如果某截面上的外法线是沿坐标轴的正方向，这个截面上的应力分量以沿坐标轴正方向为正；如果某截面上的外法线是沿坐标轴的负方向，这个截面上的应力分量以沿坐标轴负方向为正。剪应力互等： xzzyzyx,物体内任意一点的应力状态可以用六个独立的应力分量 、 、 、 、 、xyzxyz来表示。zx位移位移就是位置的移动。物体内任意一点的位移，用位移在 x，y，z 坐标轴上的投影u、v、w 表示。应变物体的形状改变可以归结为长度和角度的改变。各线段的单位长度的伸缩，称为正应变，用 表示。两
5、个垂直线段之间的直角的改变，用弧度表示，称为剪应变，用 表示。物体内任意一点的变形，可以用 六个应变分量表示。zxyxzyx 、3.1.2 平面应力和平面应变问题弹性体在满足一定条件时，其变形和应力的分布规律可以用在某一平面内的变形和应力的分布规律来代替，这类问题称为平面问题。平面问题分为平面应力问题和平面应变问题。1）平面应力问题设有很薄的等厚薄板，只在板边上受到平行于板面并且不沿厚度变化的面力，体力也平行于板面且不沿厚度变化。图 3.3设板的厚度为 t，在板面上：， ， 02tz02tzx02tzy由于平板很薄，外力不沿厚度变化，因此在整块板上有， ， zzxzy剩下平行于 XY 平面的三
6、个应力分量 未知。xyx、2）平面应变问题设有很长的柱形体，支承情况不沿长度变化，在柱面上受到平行于横截面而且不沿长度变化的面力，体力也如此分布。图 3.4以柱体的任一横截面为 XY 平面，任一纵线为 Z 轴。假定该柱体为无限长，则任一截面都可以看作对称面。由对称性， ，0zxzy0w由于没有 Z 方向的位移，Z 方向的应变 。0z未知量为平行于 XY 平面的三个应力分量 ，物体在 Z 方向处于自平衡状xyx、态。3.1.3 平衡方程弹性力学中，在物体中取出一个微小单元体建立平衡方程。平衡方程代表了力的平衡关系，建立了应力分量和体力分量之间的关系。对于平面问题，在物体内的任意一点有，（3-1）
7、0YxyXyxx3.1.4 几何方程由几何方程可以得到位移和变形之间的关系。对于平面问题，在物体内的任意一点有，（3-2）xvyuxxy刚体位移由位移 u=0，v=0 可以得到应变分量为零，反过来，应变分量为零则位移分量不为零。应变分量为零时的位移称为刚体位移。刚体位移代表了物体在平面内的移动和转动。由 0xvyux可以得到刚体位移为以下形式， xv0由 可得，,yu)(),(21xfvf将 代入 可得，2,f0yudxff)()(1积分后得到， xvfyu021)(得到位移分量， xvy0当 时，物体内任意一点都沿 x 方向移动相同的距离，可见 代0,u 0u表物体在 x 方向上的刚体平移。
8、当 时，物体内任意一点都沿 y 方向移动相同的距离，可见 代,0v v表物体在 y 方向上的刚体平移。当 时，可以假定 ，此时的物体内任意一点 P（x，y）的位0,u0移分量为， xvy,P 点位移与 y 轴的夹角为 ，tgtgP 点合成位移为， ryxyvu 222)(r 为 P 点到原点的距离，可见 代表物体绕 z 轴的刚体转动。3.1.5 物理方程弹性力学平面问题的物理方程由广义虎克定律得到。1）平面应力问题的物理方程 yxxE（3-3）xyyx)1(2平面应力问题有， 0zyxzE2）平面应变问题的物理方程 yxx 12（3-4）xyyE2xyxy)1(平面应变问题有， 0zyxz在平
9、面应力问题的物理方程中，将 E 替换为 、 替换为 ，可以得到平面211应变问题的物理方程；在平面应变问题的物理方程中，将 E 替换为 、 替换为2)(，可以得到平面应力问题的物理方程。1图 3.5求解弹性力学平面问题，可以归结为在任意形状的平面区域 内已知控制方程、在位移边界 上约束已知、在应力边界 上受力条件已知的边值问题。然后以应力分量为基uSS本未知量求解，或以位移作为基本未知量求解。如果以位移作为未知量求解，求出位移后，由几何方程可以计算出应变分量，得到物体的变形情况；再由物理方程计算出应力分量，得到物体的内力分布，就完成了对弹性力学平面问题的分析。3.2 单元位移函数根据有限元法的
10、基本思路，将弹性体离散成有限个单元体的组合，以结点的位移作为未知量。弹性体内实际的位移分布可以用单元内的位移分布函数来分块近似地表示。在单元内的位移变化可以假定一个函数来表示，这个函数称为单元位移函数、或单元位移模式。对于弹性力学平面问题，单元位移函数可以用多项式表示， .26524321 yaxayxau（3-5）bbbv多项式中包含的项数越多，就越接近实际的位移分布，越精确。具体取多项，由单元形式来确定。即以结点位移来确定位移函数中的待定系数。图 3.6如图 3.6 所示的 3 结点三角形单元，结点 I、J、M 的坐标分别为 、 、),(iyx),(j，结点位移分别为 、 、 、 、 、
11、。六个节点位移只能确定六个多),(myxiuivjjmuv项式的系数，所以 3 结点三角形单元的位移函数如下，（3-6）yaxu65421v将 3 个结点上的坐标和位移分量代入公式（3-6）就可以将六个待定系数用结点坐标和位移分量表示出来。将水平位移分量和结点坐标代入（3-6）中的第一式， mmjjj iii yaxu321写成矩阵形式，（3-7）3211ayxumjjiimji令 ，Tmjjiiyx则有 （3-8）mjiua1321T*1，A 为三角形单元的面积。2T的伴随矩阵为，（3-9）T*Tijjiijji miimmjjjj xyxy令 （3-10）mjijimjjjiii cbac
12、baT*则 （3-11）mjijijiucaA213同样，将垂直位移分量与结点坐标代入公式（3-6）中的第二式，可得，（3-12）mjijijivcbaAa21654将（3-11） 、 （3-12）代回（3-6）整理后可得， )()()( mmjjjjiii uycxbauycxauyxu 21jjjjiii vvbvcbaAv令 （下标 i，j ，m 轮换）)(yxNiiii 可得 （3-13） mjiji mji vuvuNNvu00单元内的位移记为 vuf单元的结点位移记为 mjijievu单元内的位移函数可以简写成，（3-14）eNf把N称为形态矩阵， Ni 称为形态函数。选择单元位移
13、函数应满足以下条件：1）反映单元的刚体位移与常量应变。2）相邻单元在公共边界上的位移连续，即单元之间不能重叠，也不能脱离。由（3-6）可以将单元位移表示成以下的形式， yaaxu2353521 xyv353564反映了刚体位移和常应变。单元位移函数是线性插值函数，因此单元边界上各点的位移可以由两个结点的位移完全确定。两个单元的边界共用两个结点，所以边界上的位移连续。形态函数 Ni 具有以下性质：1）在单元结点上形态函数的值为 1 或为 0。2）在单元中的任意一点上，三个形态函数之和等于 1。用 来计算三角形面积时，要注意单元结点的排列顺序，当三个结点 i，j ，m 取逆时T针顺序时， ；当三个
14、结点 i，j ，m 取顺时针顺序时， 。0T21A 0T21A例题：如图 3.7 所示等腰三角形单元，求其形态矩阵N。解： 由 jmji yxajibjmixc在公式中轮换下标可以计算得， ，00ayxajji aybmji 0jmic， ，yxiij imjyacmij 0，2yxaijji aybjim0cijm三角形积为 2aA形态函数为 axaycxbNiiii )0(1)(12yaAjjjj2ayxayxaycxbaANmm 1)(1)(212形态矩阵为 ayxayx100三角形面积的计算公式可得， 2102112ayxAmjjii 如果把三个结点按顺时针方向排列，即 i（a ，0）
15、 ，j（0， 0） ，m （0，a）212112yxmjjii 3.3 单元载荷移置有限元法的求解对象是单元的组合体，因此作用在弹性体上的外力，需要移置到相应的结点上成为结点载荷。载荷移置要满足静力等效原则。静力等效是指原载荷与结点载荷在任意虚位移上做的虚功相等。单元的虚位移可以用结点的虚位移 表示为，e*（3-15）eNf*令结点载荷为 mjieYXYR1）集中力的移置如图 3.7 所示，在单元内任意一点作用集中力 yxP图 3.8由虚功相等可得，* PNRTeeT由于虚位移是任意的，则 （3-16）RTe例题 1：在均质、等厚的三角形单元 ijm 的任意一点 p（x p，y p）上作用有集
16、中载荷。yxxmjjiimji PNYXYp),(002）体力的移置令单元所受的均匀分布体力为 yxp由虚功相等可得， tdxNRTeT*（3-17）tdxype3）分布面力的移置设在单元的边上分布有面力 ，同样可以得到结点载荷，TYXP,（3-18）sTetdsNR例题 2：设有均质、等厚的三角形单元 ijm，受到沿 y 方向的重力载荷 qy 的作用。求均布体力移置到各结点的载荷。 tdxyqNYXYmjjiimji 00,0jidxytqtdxyiii AycxbaAycAxbaNiiiiii iiii 31)(2121)( tqYyi3同理， tqYymyj 31,例题 3：在均质、等厚
17、的三角形单元 ijm 的 ij 边上作用有沿 x 方向按三角形分布的载荷，求移置后的结点载荷。sxmjjiimji tdyqNNYXY0取局部坐标 s，在 i 点 s=0，在 j 点 s=l，L 为 ij 边的长度。在 ij 边上，以局部坐标表示的插值函数为， ，Li1j 0m载荷为 sqx qtLsLtdXLi 61)32()1( 020 tsqtsLj 02303.4 单元刚度矩阵根据单元的位移函数， mjiji mji vuvuNNvu00由几何方程可以得到单元的应变表达式，（3-19） mjijji mji vuvubcbcAxvyux0021记为 ，B 矩阵称为几何矩阵。eBB矩阵可
18、以表示为分块矩阵的形式 mjiB（3-20）iiibcA021由物理方程，可以得到单元的应力表达式，（3-21）eBDD称为弹性矩阵，对于平面应力问题，210)1(2E定义 为应力矩阵。BDS将应力矩阵分块表示为，mji（3-22） iiiiii bcAEBDS21)1(22应用虚功原理可以建立单元结点位移与结点力的关系矩阵，单元刚度矩阵。虚功原理：在外力作用下处于平衡状态的弹性体，如果发生了虚位移，则所有外力在虚位移上做的虚功等于内应力在虚应变上做的虚功。单元的结点力记为 Tmjjiie VUVUF单元的虚应变为 eB*单元的外力虚功为，eTF*单元的内力虚功为，tdxyT*由虚功原理可得，
19、（3-23）tdxyFTeT*TeBB)(eDSeTeeT tdxyF*（3-24）eetdxyB定义 为单元刚度矩阵。DKT在 3 结点等厚三角形单元中B和D的分量均为常量，则单元刚度矩阵可以表示为，（3-25）tATe单元刚度矩阵表示为分块矩阵：mjmi jj iijieKK（3-26）sTrrsBD（3-27）syrsxryrsK,3.5 单元刚度矩阵的性质与物理意义（一）单元刚度矩阵的物理意义假设单元的结点位移如下： Te 001由 ，得到结点力如下：eeKF（3-28）ixmyixjyixyjiKVU,表示 i 结点在水平方向产生单位位移时，在结点 i 的水平方向上需要施加的结点力。
20、ix,表示 i 结点在水平方向产生单位位移时，在结点 i 的垂直方向上需要施加的结点力。iy,选择不同的单元结点位移，可以得到单元刚度矩阵中每个元素的物理含义：表示 s 结点在水平方向产生单位位移时，在结点 r 的水平方向上需要施加的结点力。sxrK,表示 s 结点在水平方向产生单位位移时，在结点 r 的垂直方向上需要施加的结点力。sry,表示 s 结点在垂直方向产生单位位移时，在结点 r 的水平方向上需要施加的结点力。srx,表示 s 结点在垂直方向产生单位位移时，在结点 r 的垂直方向上需要施加的结点力。syrK,因此单元刚度矩阵中每个元素都可以理解为刚度系数，即在结点产生单位位移时需要施
21、加的力。（二）单元刚度矩阵的性质1）对称性利用分块矩阵的性质证明如下： sTrrsBDKrssr )( rssTrsTrTrTsTsr KBD即 eeK2）奇异性即单元刚度矩阵的行列式为零， 。0eK将定单元产生了 x 方向的刚体移动， ，此时对应的单元Te 0101结点力为零。 0101eK可以得到，在单元刚度矩阵中 1，3，5 列中对应行的系数相加为零，由行列式的性质可知， 。eK同样如果假定单元产生了 y 方向上的刚体位移 ，可以得Te 1010到，在单元刚度矩阵中 2，4，6 列中对应行的系数相加为零。3.6 整体分析得到了单元刚度矩阵后，要将单元组成一个整体结构，根据结点载荷平衡的原
22、则进行分析，即整体分析。整体分析包括以下 4 个步骤：1）建立整体刚度矩阵，2）根据支承条件修改整体刚度矩阵，3）解方程组，求出结点的位移，4）根据结点位移，求出单元的应变和应力。在这里把结点位移作为基本未知量求解。如何得到整体刚度矩阵？基本方法是刚度集成法，即整体刚度矩阵是单元刚度矩阵的集成。图 3.9如图 3.9 所示，一个划分为 6 个结点、4 个单元的结构。得到了每个单元的单元刚度矩阵后，要集成为整体刚度矩阵。3.6.1 刚度集成法的物理意义由单元刚度矩阵的物理意义可知，单元刚度矩阵的系数是由单元结点产生单位位移时引起的单元结点力。在如图 3.9 所示的结构中，使结点 3 产生单位位移
23、时，在单元（1）中的结点 2 上引起结点力。由于结点 2、3 同时属于单元（1） 、 （3） ，在单元（2）中的结点 2 上同样也引起结点力，因此，在整体结构中当结点 3 产生位移时，结点 2 上的结点力应该是单元（1） 、（2）在结点 2 上的结点力的叠加。刚体集成法即结构中的结点力是相关单元结点力的叠加，整体刚度矩阵的系数是相关单元的单元刚度矩阵系数的集成。结点 3 在整体刚度矩阵的对应系数，应该是单元（1） 、（3） 、 （4）中对应系数的集成。3.6.2 刚度矩阵集成的规则1）将单元刚度矩阵中的每个分块放到在整体刚度矩阵中的对应位置上，得到单元的扩大刚度矩阵。单元刚度矩阵系数取决于单元
24、结点的局部编号顺序，必须知道单元结点的局部编号与该结点在整体结构中的总体编号之间的关系，才能得到单元刚度矩阵中的每个分块在整体刚度矩阵中的位置。将单元刚度矩阵中的每个分块按总体编码顺序重新排列后，可以得到单元的扩大矩阵。假定单元结点的局部编号与整体的对应关系如下：单元编号 单元结点局部编号 单元结点整体编号1 i 31 j 11 m 22 i 52 j 22 m 43 i 53 j 33 m 24 i 34 j 54 m 6单元（2）的单元扩大矩阵 的分块矩阵形式如下，只列出非零的分块：)2(K局部编号 1j11i整体编号1 2 3 4 5 61j1m2 )2(jK)2(jmK)2(ji1i3
25、4 )2(mj )2()2(i5 i im62）将全部单元的扩大矩阵相加得到整体刚度矩阵。 )4()3()2()1( KK整体刚度矩阵如下所示：整体编号1 2 3 4 5 61 )(j )1(jm)1(ji2 )1(mjK)(+ 2j+ )3(m)(miK+ 3j )2(jm)2(jiK+ 3m3 )1(ij 1i)3(j )1(i+ 3j+ )4(iK)3(ji+ 4)4(im4 )2(mj )2(m)2(miK5 )2(ij+ 3m)3(ij+ 4)2(i )2(i+ 3+ )4(j)4(jm6 )4(miKmK)4(m3.7 约束条件的处理图 3.9 所示的结构的约束和载荷情况，如图 3
26、.10 所示。结点 1、4 上有水平方向的位移约束，结点 4、6 上有垂直方向的约束，结点 3 上作用有集中力（P x，P y） 。图 3.10整体刚度矩阵K 求出后，结构上的结点力可以表示为：KF根据力的平衡，结点上的结点力与结点载荷或约束反力平衡。用 表示结点载荷和P支杆反力，则可以得到结点的平衡方程：（3-29）P这样构成的结点平衡方程组，在右边向量P中存在未知量，因此在求解平衡方程之前，要根据结点的位移约束情况修改方程(3-29)。先考虑结点 n 有水平方向位移约束，与 n 结点水平方向对应的平衡方程为：（3-30）122,112,12,1,2 nnnnn PvKuvKu根据支承情况，
27、方程（3-30）应该换成下面的方程：（3-31）0n对比公式（3-30）和（3-31） ，在式（3-29）中应该做如下修正：在K矩阵中，第 2n-1 行的对角线元素 改为 1，该行中全部非对角线元素改2,1nK为 0；在P中，第 2n-1 个元素改为 0。为了保持K矩阵的对称性，将第 2n-1 列的全部非对角元素也改为 0。同理，如果结点 n 在垂直方向有位移约束，则（3-29）中的第 2n 个方程修改为，nv在K矩阵中，第 2n 行的对角线元素改为 1，该行中全部非对角线元素改为 0；在P中，第 2n 个元素改为 0。为了保持K矩阵的对称性，将第 2n 列的全部非对角元素也改为0。（3-32
28、） nnPPvuvu21432121001000对图 3.9 所示结构的整体刚度在修改后可以得到以下的形式，（3-33）210*100*00*01 Et称对如果结点 n 处存在一个已知非零的水平方向位移 ，这时的约束条件为，*nu（3-34）*u在K矩阵中，第 2n-1 行的对角线元素 乘上一个大数 A，向量P中的对应换12,nK成 ，其余的系数保持不变。*12,nAK方程改为，（3-35）*12,2,112,12,1,2 nnnnn uAKvKuAvKu A 的取值要足够大，例如取 1010。只有这样，方程（3-35）才能与方程（3-34）等价。如果结点 n 处存在一个已知非零的垂直方向位移
29、 ，这时的约束条件为，*nv。*v也可以采用同样的方法修改整体刚度矩阵。3.8 整体刚度矩阵的特点与存储方法用有限元方法分析复杂工程问题时，结点的数目比较多，整体刚度矩阵的阶数通常也是很高的。那么，是否在进行计算时要保存整体刚度矩阵的全部元素？能否根据整体刚度矩阵的特点提高计算效率？整体刚度矩阵具有以下几个显著的特点：对称性，稀疏性，非零系数带形分布。1）对称性由单元刚度矩阵的对称性和整体刚度矩阵的集成规则，可知整体刚度矩阵必为对称矩阵。利用对称性，只保存整体矩阵上三角部分的系数即可。2）稀疏性单元刚度矩阵的多数元素为零，非零元素的个数只占较小的部分。如图 3.11 所示的结构，结点 2 只和
30、通过单元联接的 1、3、4、5 结点相关，结点 5 只和通过单元联接的2、3、4、6、8、9 结点相关。由单元刚度矩阵的物理意义和整体刚度矩阵的形成方式可知，相关结点 2、3、4、6、8、9 及结点 5 本身产生位移时，才使结点 5 产生结点力，其余结点产生位移时不在该结点处引起结点力。在用分块形式表示的整体矩阵中，与相关结点对应的分块矩阵具有非零的元素，其它位置上的分块矩阵的元素为零，如图 3.12所示。图 3.11 图 3.12 整体刚度矩阵的分块矩阵示意3）非零元素带形分布在图 3.12 中，明显可以看出，整体刚度矩阵的非零元素分布在以对角线为中心的带形区域内，这种矩阵称为带形矩阵。在包
31、括对角线元素的半个带形区域内，每行具有的元素个数叫做半带宽，用 hbd 表示。2)1( 差 值相 邻 结 点 的 编 码 的 最 大d图 3.11 所示结构的相邻结点编码的最大差值为 4，所以半带宽为 10。二维等带宽存储设整体刚度矩阵K 为一个 的矩阵，最大半带宽为 d。利用带形矩阵的特点和对n称性，只需要保存以 d 为固定带宽的上半带的元素，称为二维等带宽存储。进行存储时，把整体刚度矩阵K 每行中的上半带元素取出，保存在另一个矩阵K *的对应行中，得到一个 矩阵K *。dn把元素在K矩阵中的行、列编码记为 r、s，在矩阵K *中的行、列编码记为 r*、s *，对应关系如下：r*=rs*=s
32、-r+1图 3.13(a) 图 3.13(b)如图 3.13(a)所示的最大半带宽为 d 的整体刚度矩阵K，采用二维等带宽存储后得到如图 3.13(b)所示的矩阵K *。用新的方法存储后，K矩阵中的对角线元素保存在新矩阵中的第 1 列中，K 矩阵中的 r 行元素仍然保存在新矩阵的 r 行中，K矩阵中的 s 列元素则按照新的列编码保存在新矩阵的不同列中。采用二维等带宽存储，需要保存的元素数量与K矩阵中的总元素数量之比为 。所nd存储的元素数量取决于最大半带宽 d 的值，d 的值则由单元结点的编码方式决定。虽然在采用二维等带宽存储时，仍然会保存一些零元素，但是采用这种方法时元素寻址很方便。图 3.
33、14(a) 图 3.14(b)对于同样的有限元单元网格，按照图 3.14(a)的结点编码，最大的半带宽为 14；按照图3.14(b)的结点编码，最大的半带宽为 18；按照图 3.11 的结点编码，最大的半带宽为 10。3.9 线性方程组解法由于有限元分析需要使用较多的单元，线性方程组的阶数很高，有限元求解的效率很大程度上取决于线性方程组的解法。利用矩阵的对称、稀疏、带状分布等特点提高方程求解效率是关键。线性方程组的解法分为两大类：直接解法，迭代解法。直接解法以高斯消去法为基础，包括高斯消去法、等带宽高斯消去法、三角分解法，以及适用于大型方程组求解的分块算法和波前法等。迭代算法有高斯-赛德尔迭代
34、、超松弛迭代和共轭梯度法等。在方程组的阶数不是特别高时，通常采用直接解法。当方程组的阶数过高时，为避免舍入误差和消元时有效数损失等对计算精度的影响，可以选择迭代方法。这里给出了用 Fortran 语言编写的等带宽高斯消去法的代码，其中 NROW 为矩阵行的数目，NHBW 为最大半带宽。SUBROUTINE SOLVERB(NROW,NHBW,STIFF,DISPL)C Band elimination methodC C : SOLVE EQUATIONS KS*H=(Q) :C :C KS is stored with the half band width methodC Input:C
35、NROW - the quantity of rowsC NHBW - half band widthC STIFFNROW,NHBWC - coefficient matrix stored with the half width methodC DISPLNROW - the right hand vectorC Output:C DISPL - results of variablesC Warning:C array STIFF is changed.C*ADD:DPR*IMPLICIT DOUBLE PRECISION ( A-H,O-Z )C*END:DPR*CCINTEGER N
36、ROW,NHBWDIMENSION STIFF(NROW,NHBW),DISPL(NROW)IN = NROWKD = NHBWDO K = 1, IN-1IF(K+KD-1).GE.IN) THENIM = INELSEIM = K + KD - 1ENDIFDO I = K + 1, IML = I - K + 1CC = STIFF(K,L)/STIFF(K,1)DO J = 1,KD-L+1M = J+I-KSTIFF(I,J) = STIFF(I,J)-CC*STIFF(K,M)ENDDODISPL(I) = DISPL(I)-CC*DISPL(K)ENDDOENDDODISPL(I
37、N) = DISPL(IN)/STIFF(IN,1)DO I = IN-1,1,-1IF(KD.GT.(IN-I+1) THENJM = IN-I+1ELSEJM = KDENDIFDO J = 2,JMIH1 = J+I-1DISPL(I) = DISPL(I)-STIFF(I,J) * DISPL(IH1)ENDDODISPL(I) = DISPL(I) / STIFF(I,1)ENDDORETURNENDANSYS 提供了多种求解器供选择，如图 3.15 所示，分为直接解法和迭代解法。直接解法包括：波前法（Frontal Solver） ，稀疏法（Sparse Direct Sovler） 。迭代解法包括：雅可比共轭梯度法（Jacobi Conjugate Gradient Solver， JCG） ，不完全共轭梯度法（Incomplete Cholesky Conjugate Gradient Solver， ICCG）预处理共轭梯度法（Preconditioned Conjugate Gradient Solver ， PCG）代数多格法（Algebraic Multigrid Solver，AMG）区域分割法（Distributed Domain Solver， DDS）图 3.15 ANSYS 提供的方程组求解方法
