第六章 有限元程序设计中的若干问题
有限元计算注意事项
1、计算原理任何有限元模拟的第一步都是利用一个有限单元的集合离散结构的实际几何形状,每个单元（element）代表这个实际结构的一个离散部分.这些单元通过公用的节点（node）来连接.节点和单元的集合成为网格（mesh）.在一个特定网格中的单元数目称为网格密度（mesh density）.在ansys计算过程中,程序以每个节点的每个自由度建立平衡方程,以节点的位移作为未知量,利用矩阵求解节点的位移.一旦节点位移求出,整个结构的应力和应变都很容易计算出来.这种计算的过程和方法,数学上称之为隐式方法.从上叙述来看,整个计算过程中就是求解n个n元一次方程组（n表示节点数量）,当计算模型复杂而且庞大时,隐式求解方法的计算量还是相当大的.与之相对应的,显式求解方法.显式求解方法是通过动态方法从一个增量步前推到下一个增量步得到的.具体显式求解方法和隐式求解方法例子如下：（1）隐式求解（2）显式求解隐式求解中,计算的精度完全控制于计算步数,在一般的计算软件中（flac、abaqus）,软件均是利用不平衡力来控制计算步数（当不平衡力<10-5时,停止计算）.不平衡力＝A+B.A表示施加在节点上的集中力；B表示：在n步数下,根据第n步计算出来的应力,求出节点的内力.Flac软件中 B,以上公式是根据虚功原理推倒而得到.具体推倒过程见《flac 原理》.2、Ansys计算注意事项：计算单位、参数、荷载、标准值、设计值,计算过程中系数的加入. （1）b eam单元对于beam单元.Ansys软件中我们常用的有两种梁单元：beam188和beam4.这两种单元均是三维的梁单元,每个节点都具有6个自由度（ux、uy、uz、mx、my、mz）,并且单元坐标系x轴是i点指向j 点.Beam188单元是基于铁木辛哥理论的梁,beam4单元是我们常用的经典结构力学梁.（铁木辛哥理论考虑了梁的剪切变形,而我们常用的经典结构力学梁只考虑了弯矩对结构的变形影响）所以说,beam188可以更精确的计算梁单元,因此我们结构计算中,一般都采用beam188单元.当然还有beam189单元,189单元属于三维二次的梁单元（beam188属于三维一次梁单元）,精度比beam188更加高.定义beam188单元,一般采用如下形式：!定义单元/prep7 !进入前处理et,1,beam188 !定义单元188号标号为1!定义材料属性mp,ex,1,2.55e7 !定义弹性模量(kn/m2)mp,nuxy,1,0.167 !定义泊松比mp,dens,1,2.5 !定义密度（KN/N*KG/M3）nummrg,all ！合并重合节点numcmp,all ！压缩编号!定义梁截面SECTYPE,1,BEAM,RECT,A1,0 ! 1表示梁编号 ; RECT表示是矩形梁（还有其他t型等等,具体见ansys帮助）; A1 表示梁的名称 ; 0表示薄壁梁单元网格划分精细程度（0～5）. SECDATA,1,3,4,12 !1表示梁b ; 3表示梁h ; 4和12定义对应宽长等分份数.SECOFFSET,CENT !cent质心 ; shrc剪切中心 ; origin原始中心 ; user用户定义;!注意：当梁单元和壳单元一起使用时,可以设置梁单元的偏心,使梁的一面和壳的一面共面.（secoffset,user,offset－y,offset－z）,如下图：!划分网格LSEL,S,,,1 !选中编号为1的线.LATT,1,1,1,,,,1 !mp,r,et,,方向点,,SECTYPE截面号.LESIZE,ALL,0.2,,,,,,,1 !0.2是单元大小,1是确认细分规则.LMESH,ALL !用beam单元离散模型,形成网格.!对于划分网格,空间的beam单元,由于需要确定b、h的方向,ansys软件利用方向点来控制b、h的方向.方向点的编号最好定义的很大,如果定义太小,会影响后面的加载.具体方向点如何控制见上面的latt命令和ansys帮助.自己试两下就知道怎么用了. AllsFINISH!加载加约束/SOLUACEL,,,9.8 !重力加速度.注意方向,数值和整体坐标相反,比如重力指向z轴负向,则为正值.SFCUM,ALL,ADD !设置单元荷载是叠加还是替代,只对加在单元和节点上的荷载有效,对于加在面、线上的荷载,都只有替代作用（对同一个面,第二次加的荷载替代第一次加的荷载）!对于beam单元,只能根据sfbeam命令增加均布荷载①等大小的均布荷载.Lsel,s,,,1ESLL,S,1sfbeam,ALL,1,PRES,-161.5 !1表示作用在beam单元的①面上(如下图,③面表示beam单元的轴向,②面表示单元侧面,①面表示beam单元顶面),-161.5表示均布荷载大小,正负号可以控制作用力的方向.②梯形均布荷载Sfbeam命令是对每个单元进行加载.如果一根梁承受10～100的梯形均布荷载,而且这根梁被分成了10个beam单元,这样施加荷载就非常困难.因此我将这种加载过程写成命令流,让软件自动进行加载.命令流如下：LSEL,S,,,1 !选中要加载的那根梁（线）ESLL,S,1 !选中属于这根梁（线）的beam单元*GET,Nelem,ELEM,,COUNT, , , , !获得当前所选单元个数,赋予参数Nelem*GET,Ne,ELEM,,NUM,MIN, , , , !获得当前所选单元最小编号,赋予参数Ne*DO,I,1,Nelem !循环加载,循环次数＝单元个数ESEL,S,,,NeNSLE,S,1*GET,Nnode,NODE,,COUNT, , , , !获得当前所选节点个数,赋予参数Nnode*GET,Nn,NODE,,NUM,MIN, , , , !获得当前所选节点最小编号,赋予参数NnNN1X=NX(NN) !将nn节点的x坐标赋予NN1X(NX表示x坐标,NY表示y坐标) NN=NDNEXT(NN) !NN=当前所选节点的下一个编号NN2X=NX(NN) !将nn节点的x坐标赋予NN2Xsfbeam,ALL,1,PRES,-1630.76/3.23*(3.23-NN1X),-1630.76/3.23*(3.23-NN2X)!以上荷载公式应根据实际情况进行调整LSEL,S,,,1ESLL,S,1NE=ELnext(NE) !NE=当前所选单元的下一个编号*ENDDO!对于此命令流,根据不同的实际情况,ABC部分需要修改,其他不需要修改.!后处理 (XY平面) （大拇指指向y，就是my）etable,ImY,smisc,2 !显示弯距etable,JmY,smisc,15PLLS,IMY,JMYETABLE,IFX,SMISC,1 !显示轴力ETABLE,JFX,SMISC,14PLLS,IFX,JFXETABLE,IFY,SMISC,5 !显示剪力ETABLE,JFY,SMISC,18PLLS,IFY,JFY!注意：beam单元的结果输出都是以单元坐标系输出的,且拉为正、压为负.前面我们已经知道,单元坐标系x轴就是i点指向j点,其他坐标可以根据整体坐标系推出.详细内容见ansys 帮助.（2）S hell单元对于shell单元应用的范围,ansys软件并没有强制规定,只是从字面上区分了薄壳和厚壳.我以前看过一本电子教案《仿真在线》,里面说一般规定壳体的主尺寸是厚度的10倍左右,都是可以用壳体来模拟的.一般高度与跨度之比（非与单元尺寸比较）<1/15,可以当作薄壳处理,>1/15 & <1/10,可以当作厚壳来处理.shell63是薄壳单元,他包含弯曲和薄膜效应,但是忽略横向剪切变形；shell43,shell143,shell181,shell91,shell93和shell99,都属于厚壳单元,不仅有弯曲、薄膜效应,他也包含了横向剪切效应.横向剪切被表示为整个厚度上的常剪切应变.这种一阶近似只适用于中等厚度壳体.线形分析时,如果不包含横向剪切应变,使用63,163单元;如果横向剪切变形重要,则遵守以下原则:均匀材料,使用43,93,143单元,复合材料使用91,99,181.我们土木工程中,一般利用shell43计算.!定义单元/prep7 !进入前处理et,1,shell43 !定义单元43号标号为1!定义材料属性mp,ex,1,2.55e7 !定义弹性模量(kn/m2)mp,nuxy,1,0.167 !定义泊松比mp,dens,1,2.5 !定义密度（KN/N*KG/M3）!定义墙体厚度!①等厚度板R,1,2 !1表示编号,2表示厚度（m）R,2,3Asel,s,,,1 !选中1号面Aatt,1,1,1 !mp,real,typeESIZE,0.2 !定义单元大小为0.2左右MSHAPE,0,2D !规定划分单元形状,0表示四边形（1表示三角形）,2d表示划分面（3d表示划分体）MSHKEY,2 !指定是自由划分还是映射划分,2表示：尽量用映射划分,不符合要求就自动使用自由划分,具体参见ansys帮助的eshkey命令. AMESH,ALL !划分面单元.！！！注意：在网格剖面方面，最好全部用四边形，而且形状尽量规则、均匀！因为将来后处理内力提取的时候，提取出来的力和单元的大小有直接的关系。
有限元程序设计三维问题的有限元方法PPT课件
0
0
N i N i
z y
N i
x
0
如何推导B矩阵？（形函数在自然坐标系下得到）
2021/5/29
第17页/共25页
六面体单元
2021/5/29
Ni Ni x Ni y Ni z
x y z
Ni Ni x Ni y Ni z
x y z
Ni Ni x Ni y Ni z
1=i
4=l P 3=k
L1 L2 L3 L4 1 VP234 VP134 VP124 VP123 V1234
z
y
x
2=j
2021/5/29
第4页/共25页
四面体常应变单元
• 形函数
1
at the home node i
Li 0 at the the remote nodes jkl
=0
x
=1
z=Z O [x=(1)(x4xB)xB, y=(1)(y4yB)yB, z=(1)z4]
2021/5/29
第12页/共25页
四面体常应变单元
Jacobian:
x x x
J
y
y
y
z
z
z
x21 x31 det[J] y21 y31
0
x31 y31
z4
x41 x21 x31 y41 y21 y31 6V 2
x N i
y
J
N i
1
N i
Ni
Ni
z
ke
BTcBdA 1 1 1 BTcB det[J]d dd 1 1 1
Ve
1 1 1
nml
I
f (,)dd
11-有限元若干问题讨论-1
边 界 条 件 的 处 理 与 支 反 力 的 计 算
有限元分析中的若干问题讨论
例：平面问题斜支座的处理（p200）
边 界 条 件 的 处 理 与 支 反 力 的 计 算
单 元 形 函 数 矩 阵 与 刚 度 矩 阵 的 性 质
单元形状函数性质1：0/1性质
3. 考虑单元发生刚体位移的情形 设单元有刚体位移 都为 ，即 ，由于是刚体位移，则单元的位移场函数及节点位移
有限元分析中的若干问题讨论
一、形状函数矩阵的性质
单 元 形 函 数 矩 阵 与 刚 度 矩 阵 的 性 质
二、刚度矩阵的性质
单 元 形 函 数 矩 阵 与 刚 度 矩 阵 的 性 质
3. 考察位移 根据以上讨论，总结出以下性质: 单元刚度矩阵性质5：奇异性质
单元刚度矩阵性质6：行(或列)的代数和为零的性质 刚度矩阵的任一行（或列）代表一个平衡力系；当节点位移全部为线位移时 (即为C0型问题)，任一行（或列）的代数和应为零。
有限元分析中的若干问题讨论
位移边界条件BC(u)在大多数情形下有两种类型
边 界 条 件 的 处 理 与 支 反 力 的 计 算
(5-46)
有限元分析中的若干问题讨论
一、处理边界条件的直接法
边 界 条 件 的 处 理 与 支 反 力 的 计 算
(5-49) (5-50) (5-51)
有限元分析中的若干问题讨论
有限元分析中的若干问题讨论
单 元 的 节 点 编 号 与 总 刚 度 阵 的 存 储 带 宽
由于刚度矩阵是对称的，可以看出，若节点的DOF数为λ，则每一
个单元在整体刚度矩阵的半带宽(semi bandwidth)为：
其中n为整个结构系统的单元数。显然，对于2D问题，有λ＝2， 对于3D问题，有λ＝3。 因此在计算机中，一般都采用二维半带宽存储刚度矩阵的系数，为 等带宽存储。
c++面向对象的有限元程序设计
《C++面向对象的有限元程序设计》一、引言在计算机科学和工程中，有限元方法是一种数值分析技术，广泛应用于工程设计和科学研究领域。
C++作为一种流行的编程语言，在有限元程序设计中也扮演了重要角色。
本文将从深度和广度两个方面对C++面向对象的有限元程序设计进行全面评估，并撰写一篇有价值的文章，以帮助读者更全面、深刻地理解这一主题。
二、C++面向对象的有限元程序设计的基本概念1. 有限元方法的基本原理有限元方法是一种数值计算方法，用于求解偏微分方程和积分方程。
通过将求解区域分割为有限个单元，建立单元之间的联系，将连续的问题转化为离散的代数问题，从而得到数值解。
在有限元程序设计中，需要考虑如何有效地表示和处理单元、节点、边界条件等信息。
2. 面向对象的程序设计思想面向对象的程序设计思想强调将现实世界中的问题抽象成对象，通过封装、继承和多态等机制构建模块化、可复用的代码结构。
在C++中，类和对象是面向对象程序设计的核心概念，有限元程序设计可以通过抽象出单元、节点、网格等对象来实现。
三、深入探讨C++面向对象的有限元程序设计1. C++语言特性在有限元程序设计中的应用在C++语言中，有丰富的特性可以用于实现面向对象的有限元程序设计。
类的封装可以用于表示单元和节点对象的属性和行为，继承可以用于构建具体单元类型的层次结构，多态可以实现对不同单元类型的统一处理。
2. 优化设计思路下的C++面向对象有限元程序设计针对大规模的有限元计算，优化的设计思路是必不可少的。
C++中提供了丰富的性能优化手段，如模板元编程、内联函数、移动语义等，可以在面向对象的有限元程序设计中发挥重要作用。
四、总结和回顾在本文中，我们对C++面向对象的有限元程序设计进行了全面评估，并撰写了一篇有价值的文章。
通过深入探讨原理、语言特性和优化设计思路，帮助读者更全面地理解了这一主题。
从我的个人观点看，C++面向对象的有限元程序设计是一个值得深入研究的领域，它不仅涉及到程序设计技术，还涉及到数值计算和工程应用等多个领域的知识。
有限元单元法程序设计
有限元单元法程序设计有限元单元法是一种用于工程结构分析和设计的计算方法，它将大型结构分解为许多小的离散单元，通过分析单元之间的相互作用来预测结构的力学行为。
有限元单元法程序设计是指针对特定工程问题，编写计算机程序来实现有限元分析的过程。
下面将介绍有限元单元法程序设计的基本流程和关键要点。
一、问题建模和网格划分有限元单元法程序设计的第一步是对工程结构进行合理的建模和网格划分。
建模的目的是将实际结构抽象为适用于有限元分析的数学模型，包括定义结构的几何特征、材料属性、边界条件等。
网格划分是将结构分解为许多小的单元，每个单元具有一定的形状和尺寸，以便于数值计算。
常用的单元形状包括三角形、四边形、四面体、六面体等，根据结构的特点选择合适的单元形状和尺寸。
二、单元刚度矩阵和载荷矩阵的求解在有限元单元法程序设计中，需要编写算法来求解每个单元的刚度矩阵和载荷矩阵。
单元刚度矩阵描述了单元内部的力学性能，包括刚度、弹性模量、泊松比等，它们通常通过数学公式或有限元理论推导得到。
载荷矩阵描述了单元受到的外部荷载，可以是均匀分布载荷、集中载荷或者边界条件引起的约束力。
通过合适的数值积分方法，可以计算得到每个单元的刚度矩阵和载荷矩阵。
三、组装全局刚度矩阵和载荷向量在有限元单元法程序设计中，需要将所有单元的刚度矩阵和载荷向量组装成整个结构的全局刚度矩阵和载荷向量。
这涉及到单元之间的连接关系以及边界条件的处理。
采用适当的组装算法，可以将各个单元的刚度矩阵和载荷向量叠加在一起，形成整个结构的刚度矩阵和载荷向量。
四、求解位移和应力有限元单元法程序设计的最后一步是求解结构的位移和应力。
通过斯蒂芬-泰勒算法或者其他迭代算法，可以得到整个结构的位移分布，然后根据位移场计算各个点的应变和应力。
这一过程涉及到对整个结构刚度矩阵的求解和对位移的后处理。
有限元单元法程序设计是一个复杂而又精密的工作，需要深入理解有限元原理、结构力学知识和数学方法。
《高等有限元方法-张年梅》第6章误差估计及收敛new
第六章误差，误差估计及收敛在科学研究和工程计算中，需要对有限元近似解提出一定的精度要求。
因此，仅仅在初始网格上进行一次有限元计算往往是不够的，必须对计算结果进行误差估计，以判断是否满足特定精度的要求；根据上述估计控制计算过程，以最小的代价获得满足要求的结果。
以下主要针对与时间无关的线性问题进行讨论。
§1 误差源一般来说，凡用FEA计算的结果都包含误差。
此处的“误差”是指FEA结果与数学模型精确解之间的不一致。
首先，假设计算软件适合于所处理的任务而不存在缺陷；指定的几何体、边界条件、载荷以及物理模型的材质属性适合手边的问题；用户在把数据输入到软件中时没有犯任何的错误；选择了适用的一般类型的单元（比如空间比平面可能更恰当）。
因此在描述建模误差、用户误差、软件缺陷之后，再分析计算误差的来源。
可能的误差源可以按如下的方法归类。
1、建模误差建模误差指的是物理系统和它的数学模型的差别。
FEA要分析的不是实际问题，而是简化了的数学模型，建模的过程中略去了实际问题中的好多细节，保留下来的用可以接受的数学公式来描述（比如弹性力学的平面理论，薄板理论，热传导方程，等等）。
主要来自于1）几何形状的简化紧固件的细节、小孔、其它的几何不规则性，以及材料属性的微小不均匀性都可能被忽略，至少在初始的分析中是这样的。
2）载荷被简化边界条件被理想化，比如认为支撑是刚性的。
经常把问题表示为平面的而不是三维的，或线性的而不是非线性的，或静态的而不是动态的，等等。
总体来说，建模误差指的是有意的合理的和经过考虑的近似，而不是错误，常常还有载荷及边界条件实际性能方面的不确定误差。
2、用户误差用户误差指的是在理解了物理问题，决定了要分析回答的问题，以及创建了合适的数学模型之后软件用户所犯的错误。
用户误差包括1）选错一般的单元类型（可能需要壳体单元的地方选择成了平板单元）；2）选择不合适的单元尺寸和形状；3）数据输入中的直接错误以至于使所描述的模型并不是想要的模型。
结构有限元程序设计相关问题
进一步单元刚度阵与出口位移、反力关系，可表示如下
k1m
kmm
kdr
对于平面内梁单元：有两个出口节点，每一个节点3个自由度，2个线变位，1个角变位。出口位移向量为
d u 1 v 1
1 u 2 v 2
T
2
相应的单元刚度阵为
说明
• 上述位移和反力都是在单元局部坐标系里描述的，相应于单元刚度阵为单元局部刚度阵。 • 当坐标进行变换时，单元刚度阵也需要进行相应坐标变换。 • 如何计算单元出口刚度阵，以及构造新的单元等等是有限元的一个研究领域。
V
UV
根据最小总势能原理，在其平衡位置附件取一个可能位移的变分，有
UV0
此为虚位移原理的一个表达形式。
弹性体总势能的具体表达式
V(12T E0T )dV
f T udV TT uds
V
S
小结：上述弹性理论公式是对于连续体，弹性体内任一点位移都是独立未知数，因此具有无穷多的未知数。
3 假定位移场的有限元公式推导
2 弹性理论的基本公式
• 平衡方程，应变-位移关系、应力-应变关系
在三维弹性条件下平衡方程（静平衡、动平衡有惯性项）
x x
xy y
xz z
fx
xy x
y y
yz z
fy
xz x
yz y
z z
fz
在三维坐标系下，应变-位移关系为
x
0
0
x
0
y
0
y
设单元内任一点位移由下式给出
u uvTN d
其中 是d单元出口位移向量， 为单元形函数矩阵。N
利用应变-位移关系有
Bd
其中 为B应变－出口位移映射矩阵。
西工大19秋《有限元及程序设计》在线作业参考答案
西工大19秋《有限元及程序设计》在线作业参考答案
西工大19春《有限元及程序设计》在线作业
试卷总分:100 得分:100
一、单选题 (共 10 道试题,共 20 分)
1.三角形单元三个节点编码应按（）编排。
A.顺时针
B.非线性
C.逆时针
D.线性
答案:C
2.下列不属于薄板小挠度弯曲理论基本假定的是（）。
A.直法线假定
B.板内无挤压假定
C.曲线法假定
D.中面位移假定
答案:C
3.属于不规则单元的有（）。
A.正四面体单元
B.正三棱体单元
C.任意曲边单元
D.任意六面体单元
答案:D
4.空间问题的基物力方程有（）个。
A.8
B.7
C.6
D.5
答案:C
5.不属于规则单元的是（）。
A.正四面体单元
B.正六面体单元
C.正三棱体单元
D.曲边单元
答案:D
6.在应力函数上任意增减一个（），对应力分量无影响。
A.非线性项
B.边界项
C.线性项
D.体力项
答案:C
7.总体刚度矩阵的形成方法有直接刚度法和（）。
A.置大数法。
《有限元程序设计》课件
有限元程序设计的前景展望
广泛应用
随着计算机技术的不断发展，有 限元程序设计将在更多领域得到 广泛应用，为工程设计和科学研 究提供有力支持。
技术创新
未来有限元程序设计将不断涌现 出新的技术和方法，推动该领域 不断发展壮大。
国际化发展
随着国际化交流的加强，有限元 程序设计将实现国际化发展，推 动国际合作和共同进步。
求解
求解整体方程组得到近似解。
有限元方法的应用领域
01
02
03
04
结构力学
用于分析各种结构的力学行为 ，如桥梁、建筑、机械零件等
。
流体动力学
用于模拟流体在各种介质中的 流动行为，如流体动力学、渗
流等。
热传导
用于分析温度场在各种介质中 的分布和变化。
电磁场
用于分析电磁场在各种介质中 的分布和变化，如电磁场、电
磁波等。
02
有限元程序设计的关键技术
网格生成技术
网格生成技术是有限元分析中 的重要步骤，它涉及到将连续 的物理空间离散化为有限个小 的单元，以便进行数值计算。
网格的生成需要满足一定的规 则和条件，以保证计算的精度
和稳定性。
常见的网格生成方法包括结构 化网格、非结构化网格和自适 应网格等。
网格生成技术需要考虑的问题 包括网格大小、形状、方向和 连接方式等。
02
详细描述
弹性地基板的有限元分析是一 个二维问题，需要考虑复杂的 边界条件和非线性方程的求解 。通过将地基板划分为若干个 四边形单元，可以建立非线性 方程组进行求解。
03
计算过程
04
首先将地基板划分为若干个四边 形单元，然后根据每个单元的物 理性质和边界条件建立非线性方 程组。最后通过迭代方法求解非 线性方程组得到每个节点的位移 和应力。
有限元单元法程序设计
有限元单元法程序设计有限元单元法程序设计是一种广泛应用于工程领域的数值计算方法，它能够模拟复杂结构的受力情况并计算出相应的应力、变形等物理量。
本文将从有限元单元法的基本原理、程序设计流程、关键步骤等方面入手，为您详细介绍有限元单元法程序设计的相关内容。
一、有限元单元法基本原理有限元单元法是一种工程结构分析的数值计算方法，它基于弹性力学原理，将结构划分为有限个小单元（有限元）进行离散化处理，通过对各个单元的力学行为进行分析来描述整个结构的受力情况。
有限元单元法的基本原理可以总结为以下几个步骤：1. 将结构离散化为有限个小单元，每个单元内的应力、变形等物理量满足弹性力学理论。
2. 建立每个单元的位移与节点力之间的关系，通常采用单元刚度矩阵来描述。
3. 根据整个结构的连接条件和边界条件，组装各个单元的刚度矩阵，形成整个结构的刚度矩阵。
4. 应用外载荷和边界条件，求解整个结构的位移场，并由此计算出应力、变形等物理量。
二、有限元单元法程序设计流程有限元单元法程序设计通常包括以下几个关键步骤，我们将逐步介绍其设计流程：1. 确定结构的几何形状和材料性质，将结构进行离散化处理，确定有限元的类型和数量。
2. 建立单元刚度矩阵的表达式，通常采用弹性力学理论和数值积分方法来进行推导和计算。
3. 将各个单元的刚度矩阵组装成整个结构的刚度矩阵，考虑节点之间的连接关系以及边界条件的处理。
4. 应用外载荷和边界条件，求解整个结构的位移场，并计算出节点处的应力、变形等物理量。
5. 对程序进行稳定性和准确性的验证，包括收敛性分析、误差估计等。
6. 编写相应的有限元单元法程序，实现结构的建模、求解和结果输出等功能。
三、有限元单元法程序设计的关键步骤在有限元单元法程序设计中，有几个关键的步骤需要特别重视：1. 单元选择和刚度矩阵的建立：选择适合结构特点的有限元类型，建立单元的刚度矩阵表达式，考虑单元的形函数、应变-位移关系等。
2. 结构刚度矩阵的组装：将各个单元的刚度矩阵通过节点的连接关系组装成整个结构的刚度矩阵，考虑节点自由度的排序和边界条件的处理。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
第六章有限元程序设计中的若干问题基本步骤：ⅰ.结构离散化,输入或生成结点信息-结点坐标单元信息-单元结点编号ⅱ.计算单元刚度矩阵,形成体刚度矩阵包括计算[]Bⅲ.形成结点载荷向量ⅳ.引入约束条件ⅴ.解线性方程组ⅵ.求出结点位移ⅶ.计算单元的应力并输出§6-1 约束条件的处理1．对称性与反对称性(1)对称结构承受对称载荷作用时(2)对称结构承受反对称载荷作用2. 约束位移的引入主元置1法主元赋大值§6-2 总刚度矩阵的存贮法1．半带宽存贮法2．一维压缩存贮法考虑到总体刚度矩阵中各行的带宽并不相等，有时由于结构的几何形状的原因，使总体刚度矩阵某些行的带宽特别大。
这种情况下如采用半带宽存贮法，就可能把许多零元素也包含了进去，这对节省计算机的存贮量是很不利的。
一维压缩存贮法是将总体刚度矩阵的夏三角形中每一行从第一个非零元素开始按行将元素排成一序列，存放于一维数组)SK中。
但(L是为了确定SK中的元素在[K]中的行列号，还需要将[K]中各行对角线的元素在伊维数组中的序号存放于另一辅助数组KD（N2）中（N2是总刚度矩阵的阶数）。
现举例说明这一存贮法：设有一系数阵⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡----1.30.00.00.06.00.00.00.04.87.10.00.00.00.00.01.50.00.07.10.01.50.00.00.00.00.00.02.100.03.10.03.52.03.12.05.4 在一维数组SK （13）中依次存放的是[]1.3007.16.04.81.52.1003.13.52.05.4--而辅助数组KD （6）中存放的是[]1398631KD （6）其实就是[K]中对角元素在一维数组SK （13）中的地址。
将一结构离散化后，对结点进行编号，就能依据单元号确定出总刚度矩阵[K]各行的带宽，由它依次累加就可得出其对角线元素一维存贮中的序号。
显然，形成了数组Kd ，就确定了[K]中被存贮的元素分布情况以及SK 和[K]中元素的对应关系，例如可求出[K]中第I 行带宽为)1()(--=I KD I KD L i也可确定出[K]中第I 行左边第一个非零元素在[K]中的列号 1)1()(+-+-=I KD I KD I M i此外，也能立即确定出单元刚度矩阵e K ][中的某子矩阵[]⎥⎦⎤⎢⎣⎡=⨯)4()3()2()1(22I SK I SK I SK I SK K eij 组集到一维数组存贮总刚度矩阵SK 中的地址)(2)12(1J I I KD I -⨯--⨯=，112+=I I1)(2)12(1--⨯--⨯=J I I KD I ，134+=I I结构总刚度矩阵的以上两种存贮方法，一般应用于直接解法中。
附录A平面问题有限元应力分析程序本程序是用FORTRAN语言编写的教学使用程序，采用常应变三角形单元，用于解决弹性力学平面应力问题.它极易由学员扩充为求解平面应力和平面应变问题的通用应力分析程序。
程序中总刚度矩阵按一维压缩存贮，线性代数方程组用三角分解直接法求解。
为使教学程序易读易懂，计算时输入计算机的载荷是结点载荷，任何其他载荷都要事先换算成等效的结点载荷.程序所处理的约束仅是X或Y轴方向上的零位移约束。
计算结果除输出结点位移外，还输出单元形心处的σx，σy，τxy，主应力σ1,σ2及主方向角θ。
程序的结构便于修改和扩充，易于连接图形库扩充为包含前、后处理程序（网格自动生成及计算结果的图形显示）的更完整的程序系统.1. 输入数据说明（以输入的先后为序，自由格式）（1）NN，NE，KU，KV，KRX，KRY，EO，PO————6个整型数，2个实型数，其中NN 结点总数（≤500）；NE 单元总数（≤700）；KU x方向位移受约束的结点数（≤50）；KV y方向位移受约束的结点数（≤50）；KRX 在x方向受结点载荷作用的结点数（≤60）；KRY在y方向受结点载荷作用的结点数(≤60)；EO 材料的弹性模量；PO 材料的泊桑比；（2） X（NN）———— NN个实型数，为结点的x坐标. （3） Y（NN）———— NN个实型数，为结点的y坐标. （4） IJM（NE，3）———— 3*NE个整型数，按行输入，为单元.按逆时针向的结点编号（5） JU (KU) ———— KU个整型数，为x方向位移受约束的结点号.（6） JV（KV）———— KV个整型数，为y方向位移受约束的结点号.（7） NRX(KRX) ———— KRX个整型数，为在x方向受结点载荷作用的结点号.（8） NRY(KRY) ———— KRY个整型数，为在y方向受结点载荷作用的结点号.（9） RX(KRX)———— KRX个实型数，为x方向结点载荷的大小.（10）RY(KRY) ———— KRY个实型数，为y方向结点载荷的大小.2. 其他标识符说明NF（=2*NN）方程的总数；LK 总刚度矩阵下三角形一维存贮的总长；D1,D2,D3 弹性矩阵中的元素；B(3),C(3) Bi,Ci(i=1～3)的工作单元；DEL 三角形单元面积的工作单元；U（NF）开始存放结点载荷，求解后存放结点位移；SK(LK) 按一维存贮的总刚度矩阵；EK(21) 单元刚度矩阵下三角形元素按一维存贮的数组；BM（3,6）用于存放2Δ〔B〕矩阵的元素；CM（3,6）用于存放2Δ﹝D﹞〔B〕矩阵的元素；JD（NF）总刚度矩阵下三角形对角线元素在一维存贮中的序号指示数组；JLL（6）单元刚度矩阵元素和总刚度矩阵元素行列对应关系的工作数组；SGM（3）存放单元形心处的应力σx，σy，τxy；SGM 1，SGM 2 存放单元形心处的主应力σ1,σ2；CET A 存放主应力的方向角.3. 各子程序功能SKDD 计算总刚度矩阵下三角形对角线元在一维存贮中的序号SHAPE（N）用于计算第N个单元的有关常数.FEK 计算单元刚度矩阵，并存放于EK（21）中.SKKE 单元刚度矩阵向总刚度矩阵传送.FIXD 用于处理约束条件.SOLVE 解线性代数方程组.STRESS 用于计算单元形心处的应力，主应力和主向，并输出.4. 出错信息的处理本程序中若干数组是半动态数组,如求解问题的规模超过所设数组下标的界限，则停机，等待处理（1）求解问题网格中的结点数大于500时，则停机.此时需将数组X（500），Y（500），U(1000),JD(1000)下标的大小作相应的扩大.（2）求解问题网格中的单元数大于700时，则停机.此时需将JE(700,3)数组中的前一个下标的大小作相应的扩大.（3）当求解问题的x方向或y方向的约束自由度数超过50时，则停机.此时需将JU（50）或JV（50）下标的大小作相应的扩充. （4）当求解问题的x方向或y方向的结点载荷数超过60时，则停机.此时需将NRX（60）,RX(60)或NRY(60),RY（60）下标的大小作相应的扩充.（5）当一维存贮的总刚度矩阵长度大于12000时，则停机.此时应将数组SK(12000)下标的大小作相应的扩大.5．主程序框图6.C P.FEM---A STRESS ANALYSIS FINITE ELEMENT CC PROGRAM FOR PLANE PROBLEM C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC MAIN PROGRAM C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCOMMON /CONT/ NN,NE,NF,LKCOMMON /MAT/ D1,D2,D3COMMON /ELEM/ B(3),C(3),DELCOMMON /DATA1/ X(500),Y(500),U(1000),JE(700,3) COMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6) DIMENSION JU(50),JV(50),NRX(60),NRY(60),RX* (60),* RY(60)OPEN(7,FILE= ′DATA1′,STATUS= ′OLD′ )OPEN(8,FILE= ′DATA2′,STATUS= ′NEW′ )READ(7,* ) NN,NE,KU,KV,KRX,KRY,EO,POIF(NN.LE.500) GOTO 10WRITE( *,901)STOP10 IF(NE.LE.700) GOTO 20WRITE( *,902)STOP20 IF(KU.LE.50) GOTO 30STOP30 IF(KV.LE.50) GOTO 40WRITE( *,904)STOP40 IF(KRX.LE.60) GOTO 50WRITE( *,905)STOP50 IF(KRY.LE.60) GOTO 55WRITE( *,906)STOP55 READ(7，*) (X(I),I=1,NN)READ(7，*) (Y(I),I=1,NN)READ(7，*) ((JE(I,J),J=1,3),I=1，NE) READ(7，*) (JU(I),I=1,KU)READ(7，*) (JV(I),I=1,KV)READ(7，*) (NRX(I),I=1,KRX)READ(7，*) (NRY(I),I=1,KRY)READ(7，*) (RX(I),I=1,KRX)READ(7，*) (RY(I),I=1,KRY)NF=2 * NNCALL SKDDWRITE(8,1100) NN,NE,KU,KV,KRX,KRY,EO,PO WRITE(8,* )′X(NN)= ′WRITE(8,1200)(X(I),I=1,NN)WRITE(8,* )′Y(NN)= ′WRITE(8,1200)(Y(I),I=1,NN)WRITE(8,* )′JE(NE,3)= ′WRITE(8,1300)((JE(I,J),J=1,3),I=1，NE) WRITE(8,* )′JU(KU)= ′WRITE(8,1300)(JU(I),I=1,KU)WRITE(8,* )′JV(KV)= ′WRITE(8,1300) (JV(I),I=1,KV)WRITE(8，*)ˊNRX(KRX)=ˊWRITE(8,1300) (NRX(I),I=1,KRX)WRITE(8,*)ˊNRY(KRY)=ˊWRITE(8,1300) (NRY(I),I=1,KRY)WRITE(8,*)ˊRX(KRX)=ˊWRITE(8,1200) (RX(I),I=1，KRX)WRITE(8,*)ˊRY(KRY)=ˊWRITE(8，1200) (RY(I),I=1，KRY)D1=EO/(1.-PO*PO)D2=D1*POD3=D1*(1.-PO)/2DO 60 I=1，LK60 SK(I)=0.DO 100 I=1，NF100 U(I)=0.DO 120 I=1，KRXK1=2*NRX(I)-1120 U(K1)=RX(I)DO 140 I=1,KRYK1=2*NRY(I)140 U(K1)=RY(I)DO 160 M=1,NEWRITE(*,*)ˊELE=ˊ,MCALL SHAPE(M)CALL FEK160 CALL SKKEDO 180 I=1,KU180 CALL FIXD(2*JU(I)-1)DO 200 I=1,KV200 CALL FIXD(2,*JV(I))CALL SOLVEWRITE(*,*)ˊSOLVE PASSEDˊWRITE(8,3000)WRITE(8,4000) (I,U(2*I-1),U(2*I),I=1，NN)WRITE(8,5000)DO 240 M=1，NECALL SHAPE(M)CALL STRESS(M)240 CONTINUESTOP901 FORMAT(5X,ˊNUMBER OF NODES EXCEEDSˊ* ˊ ALLOWABLEˊ)902 FORMAT(5X,ˊNUMBER OF ELEMENTˊ* ˊEXCEEDA ALLOWABLEˊ)903 FORMAT(5X,ˊNUMBER OF FIXED FREEDOM AT X-ˊ, * ˊDIRECTION EXCEEDS ALLOWABLEˊ)904 FORMAT(5X,ˊNUMBER OF FIXED FREEDOM AT Y-ˊ, * ˊDIRECTION EXCEEDS ALLOWABLEˊ)905 FORMAT(5X,ˊNUMBER OF NODAL FORCES AT X-ˊ, * ˊDIRECTION EXCEEDS ALLOWABLEˊ) 906 FORMAT(5X,ˊNUMBER OF NODAL FORCES AT Y-ˊ, * ˊDIRECTION EXCEEDS ALLOWABLEˊ) 1000 FORMAT (//10X,ˊPRIMARY DATAˊ//,1X,ˊCONˊ, * ˊSTANTS=ˊ)1100 FORMAT (6I5,2F12.2)1200 FORMAT (5F12.4)1300 FORMAT (15I5)3000 FORMAT (//10X,ˊNODAL DISPLACEMENTˊ/,*ˊNODEˊ,9X,ˊUˊ,13X,ˊVˊ,4X,ˊNODEˊ,9X,* ˊUˊ,13X,ˊVˊ/ˊ)4000 FORMAT (2(I5，2E14.4))5000 FORMAT (//2X,ˊELE.ˊ,2X,ˊELEMENT STESSˊ,/8X, * ˊSGMXˊ,8X,ˊSGMYˊ,9X,ˊTXYˊ,8X,ˊSGM1ˊ, * 8X,ˊSGM2ˊ,10X,ˊCETAˊ/)ENDCCSUBROUTINE SKDD CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALCULATE ORDERS OF DIAGONAL CC ELEMENTS OF MATRIX K C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCOMMON /CONT/ NN,NE,NF,LKCOMMON /DATA1/ X(500),Y(500),U(10001),JE(700，3)COMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6)JD(1)=1JD(2)=3DO 200 N=2，NNKK=0DO 100 I=1，NEIO=JE(I,1)JO=JE(I,2)KO=JE(I,3)IF(IO.NE.N.AND.JO.NE.N.AND.KO.NE.N)GOTO 100 M1=N-IOM2=N-JOM3=N-KOIF(M1.GT.KK) KK=M1IF(M2.GT.KK) KK=M2IF(M3.GT.KK) KK=M3100 CONTINUEKK=2 * KKJD(2*N-1)=JD(2*N-2)+KK+1JD(2*N)=JD(2*N-1)+KK+2200 CONTINUELK=JD(NF)WRITE(*,*)ˊLK=ˊ,LKIF(LK.LE.12000) GOTO 210WRITE(*,1000)STOP1000 FORMAT(5X,ˊNUMBER OF ELEMENTS OF ONE-DIˊ,*ˊMENSIONAL GLOBAL MATRIX EXCEEDS ALLOWˊ,* ˊABLEˊ)210 RETURNENDCCSUBROUTINE SHAPE(N) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EVALUATE ELEMENT PARAMETERS C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCOMMON /DATA1/ X(500),Y(500),U(1000),JE(700，3) COMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6) COMMON /ELEM/ B(3),C(3),DELDIMENSION XO(6),YO(6)DO 100 I=1，3K1=JE(N,I)K2=I+3XO(I)=X(K1)XO(K2)=XO(I)YO(I)=Y(K1)YO(K2)=YO(I)K2=2*I-1K3=K2+1JLL(K2)=K1*2-1100 JLL(K3)=K1*2DO 200 I=1,3K1=I+1K2=I+2B(I)=YO(K1)-YO(K2)200 C(I)=XO(K2)-XO(K1)DEL=(B(1)*C(2)-B(2)*C(1))/2.RETURNENDCCSUBROUTINE FEK CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC CALCULATE STIFFNESS MATRIX CC OF ELEMENT C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCOMMON /STIF/SK(12000),EK(21),JD(1000),JLL(6) COMMON /ELEM/B(3),C(3),DELCOMMON /MAT/D1,D2,D3DIMENSION BM(3,6),CM(3,6)DO 100 I=1,3DO 100 J=1,6100 BM(I,J)=0.DO 200 I=1,3K2=I+IK1=K2-1BM(1,K1)=B(I)BM(3,K2)=B(I)BM(2,K2)=C(I)BM(3,K1)=C(I)CM(1,K1)=D1*B(I)CM(1,K2)=D2*C(I)CM(2,K1)=D2*B(I)CM(2,K2)=D1*C(I)CM(3,K1)=D3*C(I)200 CM(3,K2)=D3*B(I)DO 400 I=1,6K1=I*(I-1)/2DO 400 II=1,IZ=0.DO 300 JJ=1,3300 Z=Z+BM(JJ,I)*CM(JJ.II)EK(K1+II)=Z/DEL/4400 CONTINUERETURNENDCCSUBROUTINE SKKE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC ASSEMBLE GLOBAL STIFFNESS MATRIX C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC COMMON /STIF/SK(12000),EK(21),JD(1000),JLL(6) DO 200 I=1,6K1=JLL(I)KK=I*(I-1)/2DO 200 J=1,IK2=JLL(J)IF(K2.LT.K1) GOTO 100K=JD(K2)-K2+K1GOTO 200100 K=JD(K1)-K1+K2200 SK(K)=SK(K)+EK(KK+J)RETURNENDCCSUBROUTINE FIXD(K) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC READ BOUNDARY CONDITIONS C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCOMMON /CONT/ NN,NE,NF,LKCOMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6) COMMON /DATA1/ X(500),Y(500),U(1000),JE(700,3) L=JD(K)IF(K.EQ.1) GOTO 200M=JD(K-1)IA=M+1IB=L-1DO 100 I=IA,IB100 SK(I)=0.200 IA=K+1DO 300 I=IA,NFIF(JD(I)-JD(I-1).LT.I-K+1) GOTO 300M=JD(I)-I+KSK(M)=0.300 CONTINUEU(K)=0RETURNENDCCSUBROUTINE SOLVE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SOLVE FOR SOLUTION OF EQUATIONS C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCOMMON /CONT/ NN,NE,NF,LKCOMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6) COMMON /DATA1/ X(500),Y(500),U(1000),JE(700,3) DO 300 I=1,NFIG=JD(I)-IIF(I.EQ.1) GOTO 10MI=JD(I-1)-IG+1GOTO 2010 MI=120 DO 200 J=MI,IIP=IG+JJG=JD(J)-JIF(J.EQ.1) GOTO 30MJ=JD(J-1)-JG+1GOTO 4030 MJ=140 IJ=MIIF(MJ.GT.MI) IJ=MJJ1=J-1DO 100 K=IJ,J1JK=JD(K)100 SK(IP)=SK(IP)-SK(IG+K)*SK(JK)*SK(JG+K) IF(I.EQ.J) GOTO 210JJ=JD(J)SK(IP)=SK(IP)/SK(JJ)U(I)=U(I)-SK(IP)*SK(JJ)*U(J)200 CONTINUE210 JI=JD(I)U(I)=U(I)/SK(JI)300 CONTINUEDO 400 I=1，NFJ=NF-I+1IG=JD(J)-JIF(J.EQ.1) GOTO 310MI=JD(J-1)-IG+1GOTO 320310 MI=1320 JI=J-1DO 400 K=MI,J1U(K)=U(K)-SK(IG+K)*U(J)400 CONTINUERETURNENDCCSUBROUTINE STRESS(N) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC EVALUATE STRESSES OF ELEMENTS CC AND OUTPUT THEM C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCOMMON /DATA1/ X(500),Y(500),U(1000),JE(700,3) COMMON /STIF/ SK(12000),EK(21),JD(1000),JLL(6)COMMON /MAT/D1,D2，D3COMMON /ELEM/B(3),C(3),DELDIMENSION SGM(3),UE(6),CM(3，6) DO 100 I=1，6K1=JLL(I)100 UE(I)=U(K1)DO 150 I=1，3K2=I+IK1=K2-1CM(1,K1)=D1*B(I)CM(1,K2)=D2*C(I)CM(2,K1)=D2*B(I)CM(2,K2)=D1*C(I)CM(3,K1)=D3*C(I)CM(3,K2)=D3*B(I)150 CONTINUEDO 200 I=1，3SGM(I)=0.DO 200 J=1，6200 SGM(I)=SGM(I)+CM(I,J)*UE(J)/DEL/2 S1=SGM(1)S2=SGM(2)S3=SGM(3)A1=SQRT(((S1-S2)/2.)**2+S3**2) A2=(S1+S2)/2SM1=A2+A1SM2=A2-A1IF(ABS(S3).GT.1.E-4) GOTO 500IF(S1.GT.S2) GOTO 400CETA=90.GOTO 600400 CETA=0.GOTO 600500 CETA=ATAN((SM1-S1)/S3)*57.29578 600 WRITE(8,1000) N,SGM,SM1,SM2,CETA RETURN1000 FORMAT(15,5F12.4,E14.6)END。