有限元-程序设计课件

第2章有限元程序设计方法2.1程序基本框图1、输入基本数据（结构描述）：（1）控制数据：如结点总数、单元总数、约束条件总数等；（2）结点数据：如结点编号、结点坐标、约束条件等；（3）单元数据：如单元编号、单元结点序号、单元的材料特性、几何特性等；（4）载荷数据：包括集中载荷、分布载荷等。开始输入基本数据计算单元刚度矩阵形成总体刚度矩阵形成结点荷载向量引入约束条件求解方程组，输出结点位移计算单元应力，输出结果结束第2章有限元程序设计方法2.1程序基本框图开始输入基本12、单元分析（1）各单元的bi,ci(i,j,m)，面积A；（2）应变矩阵[B]，应力矩阵[S]；（3）单元刚度矩阵[k]；（4）单元等价载荷列向量[F]。开始输入基本数据计算单元刚度矩阵形成总体刚度矩阵形成结点荷载向量引入约束条件求解方程组，输出结点位移计算单元应力，输出结果结束3、系统分析（1）整体刚度矩阵[K]的组装；（2）整体载荷列阵{P}的形成；[K]的存储；约束引入；求解2、单元分析（1）各单元的bi,ci(i,j,m)，面2总刚存贮全矩阵存贮法：不利于节省计算机的存贮空间，很少采用。K[i,j]对称三角存贮法：存贮上三角或下三角元素。半带宽存贮法：存贮上三角形（或下三角形）半带宽以内的元素。一维压缩存贮法：半带宽存贮中仍包含了许多零元素。存贮每一行的第一个非零元素到主对角线元素。总刚存贮全矩阵存贮法：不利于节省计算机的存贮空间，很少采用。3等带宽形式UBWUBW行号1→IR

→N→1列号JC行号1→IR→N→1JC-(IR-1)方阵形式（1）半带宽存贮法等带宽形式UBWUBW行号1→IR→N→1列号JC行4方阵存贮和半带宽存贮地址关系存贮方式行号列号方阵存贮IRJC等带宽存贮IRJC-IR+1半带宽计算：设结构单元网格中相邻结点编号的最大差值是d，则最大半带宽为UBW：结点编号：欲使最大半带宽UBW最小，必须注意结点编号方法，使直接联系的相邻节点的最大点号差最小。方阵存贮和半带宽存贮地址关系存贮方式行号列号方阵存贮IRJC5例：计算下图半带宽。结点数N=91，总刚[K]中的元素总数为：82（91×2）×（91×2）=33124最大半带宽UBW=(7+1)×2=16，半带宽存储矩阵元素总数为182×16=2912，约方阵元素的8.8%。例：计算下图半带宽。结点数N=91，总刚[K]中的元素总数为6（2）变带宽存贮（一维压缩存贮）等带宽存贮虽然已经节省了不少内存，但认真研究半带宽内的元素，还有相当数量的零元素。在平衡方程求解过程中，有些零元素只增加运算工作量而对计算结果不产生影响。如果这些零元素不存、不算，更能节省内存和运算时间，采用变带宽存贮可以实现（也称一维数组存贮）。变带宽存贮编程技巧要求较高，程序较长。（2）变带宽存贮（一维压缩存贮）等带宽存贮7对称方阵形式的刚度矩阵[K]UBW=4顶线顶线以上零元素无须存贮，仅顶线以下元素。对称方阵形式的刚度矩阵[K]UBW=4顶线顶线以上零8124610121618MAXA22一维数组[A]存贮刚度矩阵[K]124610121618MAXA22一维数组[A]存贮刚度9

变带宽存贮：按列存贮方式。从左到右，逐列存放；对每一列，先存主对角线元素，然后由下而上顺序存放，直到顶线下第一个元素为止。为避免混淆，我们把存贮[K]的一维数组称为[A]。实现变带宽存贮的关键问题是：总刚中元素Kij在一维数组A中的地址是什么？为此，需要知道主元Kii在A中的位置和相应列高hi。主元位置：采用一个一维数组MAXA存主元在A中位置。MAXA=[1,2,4,6,10,12,16,18，22]。列高hj：第j行的左带宽。变带宽存贮：按列存贮方式。从左到右，逐列存放10从第j列的主对角线元素起到该列上方第一个非零元素为止，所含元素的个数称为第j列的列高，记为hj；如果把第j列上方第1个非零元素的行号记为mj，则第j列的列高为

hj=j-mj+1其实，hj就是第j行的左带宽，因而必有UBW=max(hj)

j=1,2,…,N利用节点位移信息数组ID（去约束后节点位移自由度编码），可容易地确定刚度矩阵[K]任何一列的列高。

从第j列的主对角线元素起到该列上方第一个非零111234①②③④1○○○○○○6783452xY例：求图示框架结构h7=?。利用ID数组得各单元的连接数组LM（假定小号为i）（1）ID数组节点号：1234按列，遇1变0，遇0加1。1234①②③④1○○○○○○6783452xY例：求图示框12连接数组：1号单元：

LM=[0,0,1,0,0,2]2号单元：LM=[0,0,2,3,4,5]

3号单元：LM=[3,4,5,6,7,8]4单元：LM=[0,0,1,6,7,8]1234①②③④1○○○○○○6783452xY1234连接数组：1234①②③④1○○○○○○6783452x13a)如果ID（i,j)=0则表明j号节点第i个自由度受有约束。b)如果ID(i,j)≠0则j号节点第i个自由度不受约束。并且，j号节点第i个位移分量在非约束节点位移列向量f中的序号就是:

ID(i,j)

a)如果b)如果14主元在一维数组[A]中的地址数组MAXA的长度是[K]的行或列数加1（N+1）。[K]的任何一个主对角元在一维数组A中的地址：第j列主对角线元素Kjj在一维数组A中的地址等于前(j-1)列的列高之和加1，即确定第j列列高的办法是：从1号单元起，对所有单元逐个进行检查。其中，与7号位移分量同在一个连接数组中的最小非零号码就是m7。显然有m7=1第7列的列高为：h7=j-mj+1=7–1+1=7主元在一维数组[A]中的地址数组MAXA的长度是15MAXA(j)=h1+…+hj-2-hj-1+1=（h1+…+hj-2+1）+hj-1=MAXA（j-1）+hj-1因为永远有MAXA(1)=1，MAXA(2)=2故计算主元地址的公式可写为

MAXA(j+1)=MAXA(j)+hj（2-1）式中，j=2，3，…，N；hj——刚度矩阵[K]第j列的列高。一维数组A的总长度(S)，即刚度矩阵K按变带宽存贮的总存贮量S=MAXA(N+1)-MAXA(1)MAXA(j)=h1+…16Ki,j在一维数组[A]中的地址

记Ki,j在一维数组A中的地址为AIJ。则由下图可知，AIJ=MAXA[J]+J–I（2-2）其中，I=mj，mj+1，…，J。图5-12j列第i行顶线下第j列MAXA(J)(j-i)个元素AIJ第j行A中地址Kj,jKi,jKmj,j[K]中地址mjKi,j在一维数组[A]中的地址记Ki,j在一维数组A中的174、引入约束条件手算时采用去行列法，而计算机编程时采用乘大数法。即：指定结点位移对应的主对角元素乘上一个大数，同时将{P}中对应元素换为结点位移指定值与扩大了的主对角线元素的乘积。4、引入约束条件手算时采用去行列法，而计算机编程时采用乘大185、线性方程组求解求解方法常用：GAUSS消元法，QR分解法等。其程序在一些专著中列出（例如见：徐士良编。FORTRAN常用算法程序集。清华大学出版社，1992）。在此不作详细介绍，其方法参阅[数值分析]有关书籍。5、线性方程组求解求解方法常用：GAUSS消元法，QR分解法196、单元应力节点位移求单元应力。首先整体节点位移变换成单元节点位移，然后再用物理方程求单元应力。例1：对角受压的正方形薄板，载荷沿厚度均匀分布，为2N/m。由于对称性，取1/4部分作为计算对象，试用有限元程序进行计算。2N/m2N/m2m2mxy6、单元应力节点位移求单元应力。首先整体节点位移变换成单元节20例2：简支梁，梁高3m,跨度18m，厚度1m，承受均布荷载10N/m2。已知按平面应力问题进行计算。18m3mxy例2：简支梁，梁高3m,跨度18m，厚度1m，承受均布荷载121网格划分网格划分22考察点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.3考察点y(m)-1.25-0.75-0.250.250.75232.2提高计算精度的方法（1）计算结果的整理计算结果包括位移和应力两个方面。在位移方面，一般无须进行整理工作。应力结果则需要整理。通常认为计算出的应力是三角形单元形心处的应力。而相邻单元之间的应力存在突变，甚至正、负符号都不相同。为了由计算结果推算出结构内某一点的接接实际的应力，必须通过某种平均计算。通常可采用两单元平均法或绕结点平均法。2.2提高计算精度的方法（1）计算结果的整理24平均法整理单元应力两单元平均法：把两个相邻单元中的常应力加以平均，用来表示公共边界中点处的应力。绕结点平均法：把环绕某一结点的各单元常应力加以平均，用以表示该结点的应力。在内结点效果较好，而在边界结点可能很差，一般改为应由内结点的应力外推计算出来。（2）网格的细分通过网格的细分，使每个单元的面积缩小，那么尽管每个单元是应变、常应力单元，仍可较好地反映结构中的应力变化，使得到的解答收敛于问题的精确解。平均法整理单元应力两单元平均法：把两个相邻单元中的常应力加以25（3）网格合理布局根据应力梯度使网格的布局合理化。即在梯度大的区域网格密些，梯度小的区域应稀些。密、稀网格之间应逐步过渡。（4）改用高阶单元受集中力的悬臂梁，采用128个三结点三角形常应变单元，以及3个八结点四边形高阶单元结果。由图可见，采用高阶元的计算精度比常应变元高得多。（3）网格合理布局26带圆孔方板的网格划分带圆孔方板的网格划分27有限元-程序设计课件282.3通用有限分析软件1、ANSYS9结构、热、流体、电磁学、声学等。2、SAP2000土木结构分析。2.3通用有限分析软件1、ANSYS929习题1、调试教材P26-30程序FEM1。2、修改FEM1,计算P31例2-2。3、以例1为对象，研究单元细分对计算结果的影响。4、用程序完成习题3和4的分析。习题1、调试教材P26-30程序FEM1。30第2章有限元程序设计方法2.1程序基本框图1、输入基本数据（结构描述）：（1）控制数据：如结点总数、单元总数、约束条件总数等；（2）结点数据：如结点编号、结点坐标、约束条件等；（3）单元数据：如单元编号、单元结点序号、单元的材料特性、几何特性等；（4）载荷数据：包括集中载荷、分布载荷等。开始输入基本数据计算单元刚度矩阵形成总体刚度矩阵形成结点荷载向量引入约束条件求解方程组，输出结点位移计算单元应力，输出结果结束第2章有限元程序设计方法2.1程序基本框图开始输入基本312、单元分析（1）各单元的bi,ci(i,j,m)，面积A；（2）应变矩阵[B]，应力矩阵[S]；（3）单元刚度矩阵[k]；（4）单元等价载荷列向量[F]。开始输入基本数据计算单元刚度矩阵形成总体刚度矩阵形成结点荷载向量引入约束条件求解方程组，输出结点位移计算单元应力，输出结果结束3、系统分析（1）整体刚度矩阵[K]的组装；（2）整体载荷列阵{P}的形成；[K]的存储；约束引入；求解2、单元分析（1）各单元的bi,ci(i,j,m)，面32总刚存贮全矩阵存贮法：不利于节省计算机的存贮空间，很少采用。K[i,j]对称三角存贮法：存贮上三角或下三角元素。半带宽存贮法：存贮上三角形（或下三角形）半带宽以内的元素。一维压缩存贮法：半带宽存贮中仍包含了许多零元素。存贮每一行的第一个非零元素到主对角线元素。总刚存贮全矩阵存贮法：不利于节省计算机的存贮空间，很少采用。33等带宽形式UBWUBW行号1→IR

→N→1列号JC行号1→IR→N→1JC-(IR-1)方阵形式（1）半带宽存贮法等带宽形式UBWUBW行号1→IR→N→1列号JC行34方阵存贮和半带宽存贮地址关系存贮方式行号列号方阵存贮IRJC等带宽存贮IRJC-IR+1半带宽计算：设结构单元网格中相邻结点编号的最大差值是d，则最大半带宽为UBW：结点编号：欲使最大半带宽UBW最小，必须注意结点编号方法，使直接联系的相邻节点的最大点号差最小。方阵存贮和半带宽存贮地址关系存贮方式行号列号方阵存贮IRJC35例：计算下图半带宽。结点数N=91，总刚[K]中的元素总数为：82（91×2）×（91×2）=33124最大半带宽UBW=(7+1)×2=16，半带宽存储矩阵元素总数为182×16=2912，约方阵元素的8.8%。例：计算下图半带宽。结点数N=91，总刚[K]中的元素总数为36（2）变带宽存贮（一维压缩存贮）等带宽存贮虽然已经节省了不少内存，但认真研究半带宽内的元素，还有相当数量的零元素。在平衡方程求解过程中，有些零元素只增加运算工作量而对计算结果不产生影响。如果这些零元素不存、不算，更能节省内存和运算时间，采用变带宽存贮可以实现（也称一维数组存贮）。变带宽存贮编程技巧要求较高，程序较长。（2）变带宽存贮（一维压缩存贮）等带宽存贮37对称方阵形式的刚度矩阵[K]UBW=4顶线顶线以上零元素无须存贮，仅顶线以下元素。对称方阵形式的刚度矩阵[K]UBW=4顶线顶线以上零38124610121618MAXA22一维数组[A]存贮刚度矩阵[K]124610121618MAXA22一维数组[A]存贮刚度39

变带宽存贮：按列存贮方式。从左到右，逐列存放；对每一列，先存主对角线元素，然后由下而上顺序存放，直到顶线下第一个元素为止。为避免混淆，我们把存贮[K]的一维数组称为[A]。实现变带宽存贮的关键问题是：总刚中元素Kij在一维数组A中的地址是什么？为此，需要知道主元Kii在A中的位置和相应列高hi。主元位置：采用一个一维数组MAXA存主元在A中位置。MAXA=[1,2,4,6,10,12,16,18，22]。列高hj：第j行的左带宽。变带宽存贮：按列存贮方式。从左到右，逐列存放40从第j列的主对角线元素起到该列上方第一个非零元素为止，所含元素的个数称为第j列的列高，记为hj；如果把第j列上方第1个非零元素的行号记为mj，则第j列的列高为

hj=j-mj+1其实，hj就是第j行的左带宽，因而必有UBW=max(hj)

j=1,2,…,N利用节点位移信息数组ID（去约束后节点位移自由度编码），可容易地确定刚度矩阵[K]任何一列的列高。

从第j列的主对角线元素起到该列上方第一个非零411234①②③④1○○○○○○6783452xY例：求图示框架结构h7=?。利用ID数组得各单元的连接数组LM（假定小号为i）（1）ID数组节点号：1234按列，遇1变0，遇0加1。1234①②③④1○○○○○○6783452xY例：求图示框42连接数组：1号单元：

LM=[0,0,1,0,0,2]2号单元：LM=[0,0,2,3,4,5]

3号单元：LM=[3,4,5,6,7,8]4单元：LM=[0,0,1,6,7,8]1234①②③④1○○○○○○6783452xY1234连接数组：1234①②③④1○○○○○○6783452x43a)如果ID（i,j)=0则表明j号节点第i个自由度受有约束。b)如果ID(i,j)≠0则j号节点第i个自由度不受约束。并且，j号节点第i个位移分量在非约束节点位移列向量f中的序号就是:

ID(i,j)

a)如果b)如果44主元在一维数组[A]中的地址数组MAXA的长度是[K]的行或列数加1（N+1）。[K]的任何一个主对角元在一维数组A中的地址：第j列主对角线元素Kjj在一维数组A中的地址等于前(j-1)列的列高之和加1，即确定第j列列高的办法是：从1号单元起，对所有单元逐个进行检查。其中，与7号位移分量同在一个连接数组中的最小非零号码就是m7。显然有m7=1第7列的列高为：h7=j-mj+1=7–1+1=7主元在一维数组[A]中的地址数组MAXA的长度是45MAXA(j)=h1+…+hj-2-hj-1+1=（h1+…+hj-2+1）+hj-1=MAXA（j-1）+hj-1因为永远有MAXA(1)=1，MAXA(2)=2故计算主元地址的公式可写为

MAXA(j+1)=MAXA(j)+hj（2-1）式中，j=2，3，…，N；hj——刚度矩阵[K]第j列的列高。一维数组A的总长度(S)，即刚度矩阵K按变带宽存贮的总存贮量S=MAXA(N+1)-MAXA(1)MAXA(j)=h1+…46Ki,j在一维数组[A]中的地址

记Ki,j在一维数组A中的地址为AIJ。则由下图可知，AIJ=MAXA[J]+J–I（2-2）其中，I=mj，mj+1，…，J。图5-12j列第i行顶线下第j列MAXA(J)(j-i)个元素AIJ第j行A中地址Kj,jKi,jKmj,j[K]中地址mjKi,j在一维数组[A]中的地址记Ki,j在一维数组A中的474、引入约束条件手算时采用去行列法，而计算机编程时采用乘大数法。即：指定结点位移对应的主对角元素乘上一个大数，同时将{P}中对应元素换为结点位移指定值与扩大了的主对角线元素的乘积。4、引入约束条件手算时采用去行列法，而计算机编程时采用乘大485、线性方程组求解求解方法常用：GAUSS消元法，QR分解法等。其程序在一些专著中列出（例如见：徐士良编。FORTRAN常用算法程序集。清华大学出版社，1992）。在此不作详细介绍，其方法参阅[数值分析]有关书籍。5、线性方程组求解求解方法常用：GAUSS消元法，QR分解法496、单元应力节点位移求单元应力。首先整体节点位移变换成单元节点位移，然后再用物理方程求单元应力。例1：对角受压的正方形薄板，载荷沿厚度均匀分布，为2N/m。由于对称性，取1/4部分作为计算对象，试用有限元程序进行计算。2N/m2N/m2m2mxy6、单元应力节点位移求单元应力。首先整体节点位移变换成单元节50例2：简支梁，梁高3m,跨度18m，厚度1m，承受均布荷载10N/m2。已知按平面应力问题进行计算。18m3mxy例2：简支梁，梁高3m,跨度18m，厚度1m，承受均布荷载151网格划分网格划分52考察点y(m)-1.25-0.75-0