八节点等参单元平面有限元.doc

八节点等参单元平面有限元 分析程序 土木工程学院 2011 2 目录目录 1 概述 1 2 编程思想 2 2 1 八节点矩形单元介绍 2 2 2 有限元分析的模块组织 5 3 程序变量及函数说明 6 3 1 主要变量说明 6 3 2 主要函数说明 7 4 程序流程图 8 5 程序应用与 ANSYS 分析的比较 9 5 1 问题说明 9 5 2 ANSYS 分析结果 10 5 3 自编程序分析结果 12 5 4 结果比较分析 12 参考文献 14 附录 程序源代码 15 第 1 页 计算力学 课程大作业 1 概述 通常情况下的有限元分析过程是运用可视化分析软件 如 ANSYS SAP 等 进行前处理 和后处理 而中间的计算部分一般采用自己编制的程序来运算 具有较强数值计算和处理 能力的 Fortran 语言是传统有限元计算的首选语言 随着有限元技术的逐步成熟 它被应用 在越来越复杂的问题处理中 但在实际应用中也暴露出一些问题 有时网格离散化的区域 较大 而又限于研究精度的要求 使得划分的网格数目极其庞大 结点数可多达数万个 从 而造成计算中要运算的数据量巨大 程序运行的时间较长的弊端 这就延长了问题解决的 时间 使得求解效率降低 因为运行周期长 不利于程序的调试 特别是对于要计算多种 运行工况时的情况 同时大数据量处理对计算机的内存和 CPU 提出了更高的要求 而在 实际应用中 单靠计算机硬件水平的提高来解决问题的能力是有限的 因此 必须寻找新 的编程语言 随着有限元前后处理的不断发展和完善 以及大型工程分析软件对有限元接口的要求 有限元分析程序不应只满足解题功能 它还应满足软件工程所要求的结构化程序设计条件 能够对存储进行动态分配 以充分利用计算机资源 它还应很容易地与其它软件如CAD 的实体造型 优化设计等接口 现在可编写工程应用软件的计算机语言较多 其中C语言 是一个较为优秀的语言 很容易满足现在有限元分析程序编程的要求 C 语言最初是为操作系统 编译器以及文字处理等编程而发明的 随着不断完善 它 已应用到其它领域 包括工程应用软件的编程 近年来 C 语言已经成为计算机领域最普 及的一个编程语言 几乎世界上所有的计算机都装有 C 的编译器 从 PC 机到巨型机到超 巨型的并行机 C 与所有的硬件和操作系统联系在一起 用 C 编写的程序 可移植性极好 几乎不用作多少修改 就可在任何一台装有 ANSI C 编译器的计算机上运行 C 既是高级 语言 也是低级语言 也就是说 可用它作数值计算 也可用它对计算机存储进行操作 第 2 页 2 编程思想 本程序采用 C 语言编程 编制平面四边形八节点等参元程序 用以求解平面结构问题 程序采用二维等带宽存储整体刚度矩阵 乘大数法引入约束 等带宽高斯消去法求解位移 在有限元程序中 变量数据需赋值的可分为节点信息 单元信息 载荷信息等 对于 一个节点来说 需以下信息 节点编号 整型 节点坐标 实型 节点已知位移 实 型 节点载荷 实型 边界条件 实型 和节点温度 实型 等 同样 对于一个单 元来说 需以下信息 单元的节点联接信息 整型 单元类型信息 桁架 梁 板 壳 等 整型 单元特性信息 厚度 内力矩等 实型 材料信息 弹性模量 泊松 比等 实型 等 在 FORTRAN 程序中 以上这些变量混合在一起 很难辨认 使程序的可读性不好 如需要进行单元网络的自适应划分 节点及单元的修改将非常困难 在进行 C 语言编译过 程中 采用结构 struct 使每个节点信息存储在一个结构体数组中 提高程序的可读性 使 数据结构更趋于合理 八节点矩形单元介绍 八节点矩形单元编号如图 2 1 所示 八节点矩形单元的位移函数为 2222 12345678 uxyxxyyx yxy MERGEFORMAT 2 1 2222 910111213141516 vxyxxyyx yxy MERGEFORMAT 2 2 其形函数为 1 122334455667788 uN uN uN uN uN uN uN uN u MERGEFORMAT 2 3 1 1223 3445 56 6778 8 vN vN vN vN vN vN vN vN v 图 2 1 第 3 页 MERGEFORMAT 2 4 式 MERGEFORMAT 2 3 和式 MERGEFORMAT 2 4 中 并且采用等参 ii NN 变换 则单元的坐标变换式可取为 1 122334455667788 1122334455667788 XN xN xN xN xN xN xN xN x YN yN yN yN yN yN yN yN y MERGEFORMAT 2 5 单元应变矩阵为 x y xy u x u y uv yx MERGEFORMAT 2 6 式 MERGEFORMAT 2 6 一般简写为 B MERGEFORMAT 2 7 其中的子块矩阵为 B i i i ii N x N B y NN yx MERGEFORMAT 2 8 由于是 的函数 在中的 要按照复合函数来求导 即 i N i Bxy iii iii NNNxy xx J NNNxy yy MERGEFORMAT 2 9 从而有 第 4 页 1 ii ii NN x J NN y MERGEFORMAT 2 10 因此 单元应力矩阵为 DB MERGEFORMAT 2 11 单元刚度矩阵为 T e A KBDB hdxdy MERGEFORMAT 2 12 其中积分采用三点高斯积分 33 11 11 111 nip ijijii iji fd dfW f MERGEFORMAT 2 13 高斯积分点的总数 和或是加权系数 和是单元内的坐标 对 2 nipn i j i W i j 于三点高斯积分 高斯积分点的位置 11 0 6 5 0 9 0 22 0 0 8 0 9 0 33 0 6 5 0 9 0 单元等效节点荷载 T e S PNP ds MERGEFORMAT 2 14 结构刚度矩阵 e e KK MERGEFORMAT 2 15 结构结点荷载列阵 e e PP MERGEFORMAT 2 16 注意 对于式 MERGEFORMAT 2 15 和式 MERGEFORMAT 2 16 中的理解 e 不是简单的叠加而是集成 第 5 页 总刚平衡方程 KP MERGEFORMAT 2 17 从式 MERGEFORMAT 2 17 求出 1 KP MERGEFORMAT 2 18 将回代入式 MERGEFORMAT 2 7 和式 MERGEFORMAT 2 11 得到和 有限元分析的模块组织 一个典型的有限元分析过程主要包括以下几个步骤 1 读输入数据 定义节点及单元数组 2 由边界条件计算方程个数 赋值荷载列阵 3 读入在带状存储的总刚度矩阵中单元和载荷信息 4 定义总刚度阵数组 5 组装总刚度阵 6 解方程组 第 6 页 其流程图可见图 2 2 输入边界条件 对称条件 形成各荷载工况的节点荷载阵 总刚分解 回代求出位移及输出 计算应变 应力 形成单元刚阵 单刚向总刚投放 坐标变换 输入原始参数 计算总刚规模 形成总刚方程 向总节点荷载阵投放 形成单元荷载阵 调整几何 弹性矩阵 调整单元位移列阵 图 2 2 第 7 页 3 程序变量及函数说明 主要变量说明 主要程序变量说明 wide分析模型的宽 hight分析模型的高 wsize宽度方向网格划分尺寸 hsize高度方向网格划分尺寸 npoin节点总数 nelem单元总数 struct node 500 节点结构体变量 struct element 500 单元结构体变量 ifpre 500 节点约束信息 posgp 3 高斯积分点 weigp 3 权系数 gpcod 2 9 高斯积分点总体坐标 bmatx 3 16 单元变形矩阵 dmatx 3 3 单元弹性矩阵 xjacm 2 2 雅可比矩阵 xjaci 2 2 雅可比矩阵的逆矩阵 djacb雅可比矩阵 行列式的值 estif 136 单元刚度矩阵 shape 8 单元形函数 deriv 2 8 形函数对局部坐标的导数 cartd 2 8 形函数对整体坐标的导数 A 30000 总体刚度矩阵 eload 16 等效节点荷载 gpcod 2 9 高斯积分点的总体坐标 第 8 页 主要函数说明 主要函数说明 void meshing 网格划分 void coordinate 节点整体坐标计算 void input 信息输入 void stif 形成单元刚度矩阵 void sfr2 计算形函数对当前积分点及形函数 对局部坐标的到数值 void jacobi2 形成雅可比矩阵 void modps 形成弹性矩阵 void bmatps 形成应变矩阵 void load 计算等效节点荷载 void asem 形成整体刚度矩阵和整体荷载列阵 void solve 求解整体方程 void stre 计算单元应力 第 9 页 4 程序流程图 任何一个有限元程序处理过程 都可以划分为三个阶段 1 前处理 读入数据和生成数据 2 处理器 有限元的矩阵计算 组集和求解 3 后处理 输出节点位移 单元应变 等 为了更清楚地 描述程序 我们 给出了程序的流程图 输入原始数据 包括模型尺寸 网格尺寸 边界条件 材料类型及属性等 计算相关数值 包括网单元总 数 节点总数 节 点坐标等 void meshing 函数 void coordinate 函数 形成单元刚度矩阵 void stif 函数 集成单元刚度矩阵 及等效节点荷载列 阵 void load 函数 void asem 函数 求解整体平衡方程 void solve 函 数 计算单元应力 void stre 函数 第 10 页 5 程序应用与 ANSYS 分析的比较 为了验证程序的正确性 我们取了一个简单框架结构 分别用自编程序和 ANSYS 进 行分析和比较 问题说明 图 5 1 所示一简支梁 高 3m 长 18m 承受均布荷载 2 10 N m 10 2 10EPa 取m 作为平面应力问题 0 167 1t 由于结构和荷载对称 只对右边一半进行有限单元计算 如图 5 2 所示 而在 y 轴各 节点布置水平连杆支座 图 5 1 图 5 2 第 11 页 5 2 ANSYS 分析结果 ANSYS 分析中 采用 plane82 单元 在板单元上边施加均布荷载 在 y 轴上 2 10 N m 的各结点约束 UX 方向的自由度 在板单元右下角的结点约束 UX 自由度 结点布置 网 格划分方案如图 5 3 所示 图 5 3 图 5 4 位移图和图 5 5 各单元图分别给出了结构的位移图和应力云图 X X 第 12 页 图 5 5 各单元图 X 图 5 4 位移图 第 13 页 从图 5 4 位移图和图 5 5 各单元图可以看到 跨中的位移和正应力最大 与简支 X 梁承受均布荷载 跨中挠度和正应力最大的力学常识相符合 可以初步认为模型是正确的 表格 1 给出了的截面上的正应力和表格 2 给出了处各横截面方0 75xm X 1 5ym Y 向位移 其中 的单位为 的单位为 X aPy m 表格 1 考查点的 y m1 51 00 50 0 5 1 0 1 5 正应力 X 270 20 178 25 88 570 0088 57178 25270 21 表格 2 考查点的 x m01 53 04 56 07 59 0 方向位移 10 6 Y 0 3485 0 3380 0 3069 0 2571 0 1914 0 1133 0 5 3 自编程序分析结果 用自编程序进行分析时 采用了整个结构模型进行计算 其他条件不变 因编程水平 有限 在后处理阶段没有给出节点处的位移与应力 而只能得到高斯积分点处的位移与应 力 得到积分点处的应力后 根据数值分析相关知识采用了插值进行处理 从而得到 的截面上的正应力和 处各横截面方向位移 分别列出表格如0 75xm X 1 5ym Y 下 表格 3 考查点的 y m1 51 00 50 0 5 1 0 1 5 正应力 X 297 2 196 06 93 250 0093196 08297 23 表格 4 考查点的 x m01 53 04 56 07 59 0 方向位移 10 6 Y 0 3481 0 3376 0 3065 0 2568 0 1910 0 1125 0 5 4 结果分析比较 为了更直观的比较自编程序和 ANSYS 的计算结果 将表格 1 和表格 3 的数据绘入图 5 6 将表格 2 和表格 4 的数据绘入图 5 7 第 14 页 图 5 6 应力比较 图 5 7 位移比较 自编程序所得结果与 ANSYS 分析结果进行比较发现 应力偏大而位移偏小 分析 其中的原因 一方面是编程水平有限 自编程序有很多不完善之处 有很多因素没有考虑 温度 变形 非线性等 所得结果可信度不是很高 只能得到应力以及位移大概的分布 第 15 页 情况 另一方面是在后处理阶段 在对高斯积分点结果进行处理时 误差难以避免 还有 一方面的原因是在进行求解时保留数据精度不够 造成误差较大 并且进行数值运算时可 能会造成大数吃小数 从而引起结果的差异 第 16 页 参考文献 1 美 史密斯 Smith I m 著 王崧等译 有限单元法编程 第三版 M 北京 电子工业出版社 2003 2 王勖成 有限单元法 M 北京 清华大学出版社 2003 3 宋建辉 涂志刚 MATLAB 语言及其在有限元编程中的应用 J 湛江师范学院学报 2003 6 24 101 105 4 郑大玉 连宏玉 周跃发 有限元编程与 C 语言 J 黑龙江商学院学报 1997 3 13 23 28 5 王伟 刘德富 有限元编程中应用面向对象编程技术的探讨 J 三峡大学学报 2001 2 23 124 128 6 汪文立 吕士俊 二级 C 语言程序设计教程 M 北京 中国水利水电出版社 2006 7 赵翔龙 FORTRAN 90 学习教程 M 北京 北京大学出版社 2002 第 17 页 附录 程序源代码 include include The gemotry of the model float mesh 2 float wide hight float wsize hsize young poiss thick int i j k knelem npoin elem 500 ielem float coord 2 1 props 3 1 lnods 8 500 ifpre 1 float posgp 3 weigp 3 estif 136 elcod 2 8 int npoin nelem kevab igaus jgaus kgasp ngash djacb float gpcod 2 9 bmatx 3 16 dmatx 3 3 float shape 8 deriv 2 8 float xjacm 2 2 xjaci 2 2 elcod 2 8 float cartd 2 8 float bmatx 3 16 dmatx 3 3 float nload 1 float press 3 2 pgash 2 dgash 2 struct node int nodenu The number of the node float cor 2 The coordinate of the node float displacement 2 The displacement of the node float load 2 The load of the node float boundary 2 node 500 struct int elementnu The number of element int localnu 8 Local number int localcorx 8 int localcory 8 Local coordinate of node float qx qy t The stress and strain element 500 第 18 页 void meshing float wide float hight float mesh 2 printf Please input the meshing size n scanf f f getchar mesh 0 wide wsize mesh 1 hight hsize The global coordinate of node void coordinate int i j 1 float x y for i 1 i 2 mesh 1 1 i x 0 while x wide if i 2 0 node j cor 1 wsize 2 node j cor 2 i 1 hsize j else node j cor 1 wsize node j cor 2 i 1 hsize j main int top 500 bottom 500 The number of top and bottom element void input void stif void jacobi2 void load void asem void solve void stre 第 19 页 npoin 8 for i 1 i 8 i lnods i 1 i for i 1 i 3 i scanf f printf The EX is e nThe uo is 8f nThe bt is 8f props 1 1 props 2 1 props 3 1 getch printf Please input the gemotry of the model n scanf f f getchar printf The wide is 0 3f the hight is 0 3f wide hight getchar meshing wide hight mesh printf The wide size is f the hight size is f n mesh 0 mesh 1 input getchar nelem mesh 0 mesh 1 for i 1 i nelem i The element number element i elementnu i for i 1 i mesh 0 i npoin 5 for i 1 i mesh 1 i npoin 3 mesh 0 2 printf d npoin for i 1 i npoin i node i nodenu i for i 1 i nelem i for j 1 j 8 j element i localnu j j for i 1 i mesh 0 i bottom i i j 1 for i mesh 0 mesh 1 1 1 i nelem i top j i 第 20 页 for i 1 i j i printf the top number is d n top i printf d n element 1 localnu 8 coordinate stif jacobi2 load asem solve stre getchar data input void input int nnode 8 int notreadchar input element node number printf input element node number n for ielem 1 ielem nelem ielem for i 1 i nnode i scanf d Input restriction message printf double digit is symbol of the restriction 1 representates stable 2 representates gived displacement n i 1 while i scanf d if i 0 scanf d scanf d else break Element stiffness matrix 第 21 页 void stif int kevab nstre ievab istre lnode jstre jevab float kgasp exisp etasp dvolu float btdbm dbmat 3 void sfr Gauss point posgp 1 sqrt 0 6 posgp 2 0 posgp 3 sqrt 0 6 weight coefficient weigp 1 5 0 9 0 weigp 2 8 0 9 0 weigp 3 5 0 9 0 number of stress nstre 3 form elastic matrix for ielem 1 ielem nelem ielem printf The number of element is d n ielem for i 1 i 8 i lnode lnods i ielem for j 1 j 2 j coord j lnode node lnode cor j elcod j i coord j lnode young props 1 1 poiss props 2 1 thick props 3 1 modps young poiss The initial value of k kevab 0 for i 1 i 16 i for j 1 j i j kevab 第 22 页 estif kevab 0 0 The gauss point shape function and deriv kgasp 0 for igaus 1 igaus 3 igaus for jgaus 1 jgaus 3 jgaus kgasp kgasp 1 exisp posgp igaus etasp posgp jgaus printf The number of gauss point is d n kgasp sfr2 exisp etasp jacob2 ielem djacb kgasp elcod dvolu djacb weigp igaus weigp jgaus thick bmatps The initial value of elastic matrix D kevab 0 for ievab 1 ievab 16 ievab for istre 1 istre nstre istre dbmat istre 0 0 for jstre 1 jstre nstre jstre dbmat istre dbmat istre bmatx jstre ievab dmatx jstre istre for jevab 1 jevab ievab ievab kevab kevab 1 btdbm 0 0 for istre 1 istre nstre istre btdbm btdbm dbmat istre bmatx istre jevab estif kevab estif kevab btdbm dvolu float sfr2 float s float t 第 23 页 float s2 t2 ss tt st stt sst st2 Shape function these variables are to simplify the formula s2 s 2 0 t2 t 2 0 ss s s tt t t st s t stt s t t sst s s t st2 s t 2 0 shape function shape 1 1 0 st ss tt sst stt 4 0 shape 2 1 0 t ss sst 2 0 shape 3 1 0 st ss tt sst stt 4 0 shape 4 1 0 s tt stt 2 0 shape 5 1 0 st ss tt sst stt 4 0 shape 6 1 0 t ss sst 2 0 shape 7 1 0 st ss tt sst stt 4 0 shape 8 1 0 s tt stt 2 0 differential coefficient to local coordinate deriv 1 1 t s2 st2 tt 4 0 deriv 1 2 s st deriv 1 3 t s2 st2 tt 4 0 deriv 1 4 1 0 tt 2 0 deriv 1 5 t s2 st2 tt 4 0 deriv 1 6 s st deriv 1 7 t s2 st2 tt 4 0 deriv 1 8 1 0 tt 2 0 deriv 2 1 s t2 ss st2 4 0 deriv 2 2 1 0 ss 2 0 deriv 2 3 s t2 ss st2 4 0 deriv 2 4 t st deriv 2 5 s t2 ss st2 4 0 第 24 页 deriv 2 6 1 0 ss 2 0 deriv 2 7 s t2 ss st2 4 0 deriv 2 8 t st Jacobi matrix void jacobi2 xjacm 2 2 jacobi matrix xjaci 2 2 j 1 elcod 2 8 local coordinate float cartd 2 8 gpcod 2 9 int i j k for i 1 i 2 i gpcod i kgasp 0 0 for j 1 j 8 j gpcod i kgasp gpcod i kgasp elcod i j shape j for i 1 i 2 i for j 1 j 2 j xjacm i j 0 0 for k 1 k 8 k xjacm i j xjacm i j deriv i k elcod j k Caculate the value of Jacobi matrix djacb xjacm 1 1 xjacm 2 2 xjacm 1 2 xjacm 2 1 printf The det of jacobi matrix is f n djacb if djacb 1e 6 printf The jacobi det of d is less or equal 0 n ielem The element of j 1 xjaci 1 1 xjacm 2 2 djacb xjaci 2 2 xjacm 1 1 djacb xjaci 1 2 xjacm 1 2 djacb xjaci 2 1 xjacm 2 1 djacb The deriv of shape function to global coordinate for i 1 i 2 i 第 25 页 for k 1 k 8 k cartd i k 0 0 for j 1 j 2 j cartd i k cartd i k xjaci i j deriv j k Form elastic matris float modps float bmatx 3 16 dmatx 3 3 float constant int i j y p y props 1 1 p props 2 1 for i 1 i 3 i for j 1 j 3 j dmatx i j 0 0 If Ntype 1 it is plane stress question constant y 1 0 p p dmatx 1 1 constant dmatx 2 2 constant dmatx 1 2 constant p dmatx 2 1 constant p dmatx 3 3 constant 1 0 p 2 0 Form strain matrix void bmatps float cartd 2 8 shape 8 deriv 2 8 bmatx 3 16 dmatx 3 3 int npoin nelem ngash mgash float gpcod 2 9 The initial value of B for i 1 i 16 i for j 1 j 3 j bmatx j i 0 0 第 26 页 Form B ngash 0 for i 1 i 8 i mgash ngash 1 ngash mgash 1 bmatx 1 mgash cartd 1 i bmatx 1 ngash 0 0 bmatx 2 mgash 0 0 bmatx 2 ngash cartd 2 i bmatx 3 mgash cartd 2 i bmatx 3 ngash cartd 1 i equal load of node void load float nload 500 float elcod 2 3 eload 16 posgp 3 weigp 3 float press 3 2 pgash 2 dgash 2 noprs 3 int iedge nedge kount float sgmar tau s dvolu pxcom pycom n m t int lnode number nloca The number of load side element load The position of gauss point printf input the number of load side n scanf d posgp 1 sqrt 0 6 posgp 2 0 0 posgp 3 sqrt 0 6 The weight of gauss point weigp 1 5 0 9 0 weigp 2 8 0 9 0 weigp 3 5 0 9 0 第 27 页 The initial load of element for i 1 i nelem i nload i 0 0 for i 1 i number i The initial value of equal node load for i 1 i 16 i eload i 0 scanf d for iedge 1 iedge nedge iedge for i 1 i 3 i scanf f f if sgmar 0 input node load q 1 2 3 and t 1 2 3 if fabs sgmar 1e 8 j 3 j for i 1 i 2 i scanf f else for i 1 i 3 i press i 1 sgmar press i 2 tau The coordinate of the load node for i 1 i 3 i input load node printf input node load number n scanf d lnode noprs i for j 1 j 2 j coord j lnode node lnode cor j elcod j i coord j lnode t 1 0 for igaus 1 igaus 3 igaus s posgp igaus sfr s t 第 28 页 for j 1 j 2 j pgash j 0 0 dgash j 0 0 current point q t and the value of deriv pgash j pgash j press i j shape i dgash j dgash j elcod j i deriv 1 i dvolu weigp igaus pxcom dgash 1 pgash 2 dgash 2 pgash 1 pycom dgash 1 pgash 1 dgash 2 pgash 2 for i 1 i 8 i nloca lnods i ielem if nloca noprs 1 j i 2 break j i 2 kount 0 for k i k8 then it is 1 if k 8 n 1 m 2 sum to get the equal node load eload n eload n shape kount pxcom dvolu eload m eload m shape kount pycom dvolu 第 29 页 print The element load matrx is n for i 1 i 16 i printf f n eload i To form global stiffness matrix and load matrix void asem float v 1 a 1 float lnods 8 1 ifpre 1 nload 1 maxa 1 float ncodf 2 estif 136 eload 16 float dlarge x int ipoin kmin lnode khigh kdofn nn nnm nwk iwk in int inodi inode icodf idofn ievab icoln jnode lnodj jdofn jevab jrow ldis lnodi maxa 1 1 for ipoin 1 ipoin npoin ipoin kmin npoin 1 for ielem 1 ielem nelem ielem for i 1 i 8 i lnode lnods i ielem if lnode ipoin break for i 1 i kmin break kmin lnode khigh ipoin 1 2 j 1 for j 1 j 2 j kdofn ipoin 1 2 j 1 maxa kdofn maxa kdofn 1 khigh j 第 30 页 nn npoin 2 nnm nn 1 nwk maxa nnm 1 printf The sum is d n nwk for iwk 1 iwk nwk iwk a iwk 0 0 for in 1 in nn nn v in 0 0 dlarge 1 0e 15 for ielem 1 ielem nelem ielem printf The number of element is d n ielem for inode 1 inode 8 inode inodi lnods inode ielem icodf ifpre lnodi ncodf 1 icodf 10 ncodf 2 icodf ncodf 1 10 for idofn 1 idofn 2 idofn ievab inode 1 2 idofn icoln lnodi 1 2 idofn for jnode 1 jnode 8 jnode lnodj lnods jnode ielem for jdofn 1 jdofnicoln break if jevab ievab kevab jevab jevab 1 2 ievab else kevab ievab ievab 1 2 jevab ldis maxa icoln icoln jrow a ldis a ldis estif kevab if nload ielem 0 第 31 页 v icoln v icoln eload ievab if ncodf idofn 0 break kevab ievab ievab 1 2 ldis maxa icoln a ldis a ldis dlarge estif kevab printf Enter gived displacement n for i 1 i npoin i icodf ifpre i ncodf 1 icodf 10 ncodf 2 icodf ncodf 1 10 for j 1 j 2 j if ncodf j 2 scanf f icoln i 1 2 j v icoln a ldis x To solve the equation void solve float v 1 a 1 maxa 1 int npoin nelem ntype nmats int l k nn nnm nwk n ku kl kh kn float ic b c klt ki nd kk m1 m2 kul nn npoin 2 nnm nn 1 for n 1 n nn n kn maxa n kl kn 1 第 32 页 ku maxa n 1 1 kh ku kl if kh 0 if a kn 0 printf Stop run Coefficient matrix is not illegal Non positive pivot for equation d pivot f n n a kn break else if kh 0 k n b 0 0 for kk kl kk0 k n kh ic 0 klt ku for j 1 j kh j ic ic 1 klt klt 1 ki maxa k