《角形单元的有限元法程序设计.ppt》由会员分享，可在线阅读，更多相关《角形单元的有限元法程序设计.ppt（23页珍藏版）》请在金锄头文库上搜索。
1、第6章 有限元程序设计方法6.1 程序基本框图程序基本框图1、输入基本数据（结构描述）、输入基本数据（结构描述）：（1）控制数据：如结点总数、单元总数、约束条件总数等；（2）结点数据：如结点编号、结点坐标、约束条件等；（3）单元数据：如单元编号、单元结点序号、单元的材料特性、几何特性等；（4）载荷数据：包括集中载荷、分布载荷等。开始输入基本数据计算单元刚度矩阵形成总体刚度矩阵形成结点荷载向量引入约束条件求解方程组，输出结点位移计算单元应力，输出结果结束2、单元分析、单元分析 （1）各单元的bi,ci(i,j,m) ， 面积A；（2）应变矩阵B，应力矩阵S；（3）单元刚度矩阵k；（4）单元等价载
2、荷列向量F。开始输入基本数据计算单元刚度矩阵形成总体刚度矩阵形成结点荷载向量引入约束条件求解方程组，输出结点位移计算单元应力，输出结果结束3、系统分析、系统分析（1）整体刚度矩阵K的组装；（2）整体载荷列阵P的形成；K的存储；约束引入；求解的存储；约束引入；求解u总刚存贮总刚存贮全矩阵存贮法：不利于节省计算机的存贮空间，很少采用。Ki,j对称三角存贮法：存贮上三角或下三角元素。半带宽存贮法半带宽存贮法 ：存贮上三角形（或下三角形）半带宽以内的元素 。一一维压缩存贮法维压缩存贮法 ：半带宽存贮中仍包含了许多零元素。存贮每一行的第一个非零元素到主对角线元素。 等带宽形式等带宽形式UBWUBW行行
3、号号1 IR N1列列 号号JC行行 号号1 IR N1JC-(IR-1)方阵形式方阵形式（1 1）半带宽存贮法）半带宽存贮法方阵存贮和半带宽存贮地址关系方阵存贮和半带宽存贮地址关系存贮方式存贮方式行号行号列号列号方阵存贮方阵存贮IRJC等带宽存贮等带宽存贮IRJC-IR+1u 半带宽计算：设结构单元网格中相邻结点编号的最半带宽计算：设结构单元网格中相邻结点编号的最大差值是大差值是d，则最大半带宽为，则最大半带宽为UBW：u结点编号：欲使最大半带宽结点编号：欲使最大半带宽UBW最小，必须注最小，必须注意结点编号方法，使直接联系的相邻节点的最大点意结点编号方法，使直接联系的相邻节点的最大点号差最
4、小。号差最小。 （2） 变带宽存贮变带宽存贮（一一维压缩存贮）维压缩存贮） 等带宽存贮虽然已经节省了不少内存，但认真等带宽存贮虽然已经节省了不少内存，但认真研究半带宽内的元素，还有相当数量的零元素。在研究半带宽内的元素，还有相当数量的零元素。在平衡方程求解过程中，有些零元素只增加运算工作平衡方程求解过程中，有些零元素只增加运算工作量而对计算结果不产生影响。如果这些零元素不存、量而对计算结果不产生影响。如果这些零元素不存、不算，更能节省内存和运算时间，采用变带宽存贮不算，更能节省内存和运算时间，采用变带宽存贮可以实现（也称一维数组存贮）可以实现（也称一维数组存贮） 。变带宽存贮编程。变带宽存贮编
5、程技巧要求较高，程序较长。技巧要求较高，程序较长。对对 称称u方阵形式的刚度矩阵方阵形式的刚度矩阵KUBW=4顶顶 线线顶线以上零元素无须存贮，仅顶线以下元素。顶线以上零元素无须存贮，仅顶线以下元素。124610121618MAXA 22u一维数组一维数组A存贮刚度矩阵存贮刚度矩阵K 变变带带宽宽存存贮贮：按按列列存存贮贮方方式式。从从左左到到右右，逐逐列列存存放放；对对每每一一列列，先先存存主主对对角角线线元元素素，然然后后由由下下而而上上顺顺序序存存放放，直直到到顶顶线线下下第第一一个个元元素素为为止止。为为避避免免混混淆淆，我们把存贮我们把存贮K的一维数组称为的一维数组称为A。 实实现现
6、变变带带宽宽存存贮贮的的关关键键问问题题是是：总总刚刚中中元元素素Kij在在一一维维数数组组A中中的的地地址址是是什什么么？为为此此，需需要要知知道道主主元元Kii在在A中的位置和相应列高中的位置和相应列高hi。u主主元元位位置置：采采用用一一个个一一维维数数组组MAXA存存主主元元在在A中中位置。位置。 MAXA =1,2,4,6,10,12,16,18，22。u列高列高hj：第：第j行的左带宽。行的左带宽。 从从第第j列列的的主主对对角角线线元元素素起起到到该该列列上上方方第第一一个个非非零零元元素素为为止止，所所含含元元素素的的个个数数称称为为第第j列列的的列列高高，记记为为hj ；如如
7、果果把把第第j列列上上方方第第1个个非非零零元元素素的的行行号号记记为为mj，则则第第j列的列高为列的列高为 hj = j - mj + 1其实，其实，hj就是第就是第j行的左带宽，因而必有行的左带宽，因而必有 UBW= max(hj) j=1,2, ,N 利利用用节节点点位位移移信信息息数数组组ID （去去约约束束后后节节点点位位移移自自由由度度编编码码），可可容容易易地地确确定定刚刚度度矩矩阵阵K任任何何一一列列的的列列高。高。 4、引入约束条件、引入约束条件手算时采用去行列法去行列法 ，而计算机编程时采用乘乘大数法大数法。即：指定结点位移对应的主对角元素乘上一个大数 ，同时将P中对应元素
8、换为结点位移指定值与扩大了的主对角线元素的乘积 。5、线性方程组求解、线性方程组求解求解方法常用：GAUSS消元法，QR分解法等。其程序在此不作详细介绍，其方法参阅数值分析有关书籍。6、单元应力、单元应力节点位移求单元应力。首先整体节点位移变换成单元节点位移，然后再用物理方程求单元应力。例例1：对角受压的正方形薄板，载：对角受压的正方形薄板，载荷沿厚度均匀分布，为荷沿厚度均匀分布，为2N/m。由。由于对称性，取于对称性，取1/4部分作为计算对部分作为计算对象，试用有限元程序进行计算。象，试用有限元程序进行计算。 2N/m2N/m2m2mxy例2：简支梁，梁高3m,跨度18m，厚度1m，承受均布
9、荷载10N/m2。已知 按平面应力问题进行计算。18m3mxy网格划分网格划分考察点y(m)-1.25-0.75-0.250.250.751.25有限元结果19711436-36-114-197弹性力学结果22513444-44-134-225误差28208-8-20-28考察点y(m)-1.25-0.75-0.250.250.751.25有限元结果16.231.237.233.720.73.6弹性力学结果10.926.734.634.626.710.9误差5.34.52.6-0.9-6.0-7.36.2 提高计算精度的方法提高计算精度的方法 （1） 计算结果的整理计算结果的整理计算结果包括位
10、移和应力两个方面。在位移方面，一般无须进行整理工作。应力结果则需要整理。通常认为计算出的应力是三角形单元形心处的应力。而相邻单元之间的应力存在突变，甚至正、负符号都不相同。为了由计算结果推算出结构内某一点的接接实际的应力，必须通过某种平均计算。通常可采用两单元平均法或两单元平均法或绕结点平均法绕结点平均法。平均法整理单元应力两单元平均法两单元平均法：把两个相邻单元中的常应力加以平均，用来表示公共边界中点处的应力。绕结点平均法绕结点平均法：把环绕某一结点的各单元常应力加以平均，用以表示该结点的应力。在内结点效果较好，而在边界结点可能很差，一般改为应由内结点的应力外推计算出来。（2）网格的细分）网
11、格的细分通过网格的细分，使每个单元的面积缩小，那么尽管每个单元是应变、常应力单元，仍可较好地反映结构中的应力变化，使得到的解答收敛于问题的精确解。（3）网格合理布局）网格合理布局 根据应力梯度使网格的布局合理化。即在梯度大的区域网格密些，梯度小的区域应稀些。密、稀网格之间应逐步过渡。（4）改用高阶单元）改用高阶单元 受集中力的悬臂梁，采用128个三结点三角形常应变单元，以及3个八结点四边形高阶单元结果。由图可见，采用高阶元的计算精度比常应变元高得多。带圆孔方板的网格划分6.3 通用有限分析软件通用有限分析软件1、ANSYS 结构、热、流体、电磁学、声学等。2、SAP2000土木结构分析。习习 题题1、调试教材程序。2、修改FEM1,计算算例。3、以算例为对象，研究单元细分对计算结果的影响。
