《计算力学》课程设计说明书.doc

有限元程序ADINA之模块AD32(1)代码分析1 课程设计的目的 本课程的教学目的及要求是使学生通过本课程的学习，掌握有限元方法的基本理论，并能够学以致用，把所学的知识结合编程上机运算来数值求解固体力学的一些实际问题，对应力分析、模态分析、热分析等实际计算固体力学问题计算步骤有一定了解。2程序模块功能描述 有限元程序ADINA之模块AD32是一个比较庞大的程序模块，它的功能就是进行梁的有限元分析，但它的分析比较齐全，基本包括梁的分析的各个方面。它不仅包括梁的静力学分析，还包括梁的动力学分析；不仅考虑了梁的弹性分析，还考虑到了梁的非线性分析（包括材料非线性和几何非线性分析）。同时它还考虑了两种典型的不同的梁的有限元模型（无剪梁单元和Timoshenko梁单元）。并且在考虑梁的材料非线性分析时，考虑了基于两种不同的强化准则（等向强化和随动强化）的弹塑性梁的模型。而且在进行不同的梁模型的分析或进行不同性质的分析时，都考虑了截面形式比较常见的矩形截面梁及管形截面梁。程序模块AD32（1）只包括了一个子程序，即BMEL子程序，该子程序计算梁的单元矩阵及应力，它主要负责数据的通信，调用其他的子程序进行整体刚度矩阵及整体质量矩阵的装配，以及计算积分点位置，求解积分点上的应力应变等等。它是AD32程序模块中对其它子程序进行综合调度与处理的子程序。其运行的总体框图如图1所示。图1 子程序BMEL执行总体框图将数据写入文件，通过文件来传输数据，以实现在不同子函数中实现数据的通信。然后开始初始化工作向量，通过调用各种函数（包括线性单元刚度矩阵、非线性单元刚度矩阵，质量矩阵、及荷载列向量）转配整体刚度矩阵。由于涉及非线性分析（包括材料非线性和几何非线性）和动态响应分析，处理起来较复杂。下面分别简要说明。对于非线性分析：刚度矩阵在具体处理时初始矩阵只包括线性刚度矩阵，在求解迭代运算时，不但用非线性项进行修正，即利用NewtonRaphson方法或其改进方法。对于动态响应分析，由于结构具有加速度，必须考虑动力效应，可等效用附加惯性力的平衡体系来考虑。此时需要系统的质量矩阵，而梁结构为连续系统，为将其转化为有限自由度，并求其数值解，故需将其离散化，变为等效的多自由度振动系统来考虑。而对于梁上的集中质量块，则应特殊考虑，将其分配与整合。装配好质量矩阵与刚度矩阵及外荷载列向量，就可以利用数值解法，求得动力响应的各阶主振动（包括主振动频率及主振型）3程序的理论基础有限元程序ADINA之模块AD32几乎包括有限元分析的所有分析过程，详述其理论基础显然是不可能的，下面就自己个人理解，简述与程序模块息息相关的一些理论。3.1 有限元分析的基本步骤有限元分析中对于具体不同的分析，其步骤有一定的不同，但其最基本的分析步骤基本一致，如图2所示。（1）选择单元及单元形函数（2）将实体模型离散化，即进行网格划分（3）计算各单元的单刚及单元节点荷载列向量（4）装配整刚及整体节点荷载列向量（5）引入约束（强制约束，一般为位移约束）（6）方程组的求解（7）由节点位移计算应力，应变等后处理图2 有限元分析的基本步骤3.2经典的梁单元（无剪梁单元）和Timoshenko梁单元(1) 无剪梁单元这种单元基于经典的梁理论，即梁的变形服从平面假设：梁的横截面变形后仍保持为平面，其仍垂直于梁的中线。对于无剪梁单元，以2节点梁单元为例，可推出其单元刚度矩阵为：（1）（2） Timoshenko梁单元Timoshenko 梁理论考虑了横向剪力对梁变形的影响，梁的变形不再服从平面假设。横截面的转角与挠度不再存在导数关系。对于如图3 所示的二次抛物线等参元为例，介绍Timoshenko梁单元 图3 三节点 Timoshenko梁节点自由度取为梁的挠度和转角 单元自由度为：（2）对应的节点力（3）单元挠度和转角分别独立插值，位移模式为（4）式中（5）单元的坐标采用等参形式，即（6）单元的广义应变可表示为（7）单元刚度阵为（8）式中坐标变换的Jacobian行列式。单位长度上分布压力q，装换为等效节点力（9）3.3 结构的非线性分析结构非线性包括材料非线性，几何非线性和边界非线性问题，在AD32程序模块中涉及到了材料非线性与几何非线性问题。材料非线性问题中应力与应变的关系不再是线性的，对于本程序模块处理的材料非线性中梁的弹塑性分析。在有限元分析中，对材料非线性问题处理较为简单，不需重新列出整个问题的表达格式，只要将材料本构关系线性化，就可用线性问题的表达格式对其求解。这类问题一般通过迭代解法求解。先形成试探解，然后通过迭代，不断修正刚度矩阵，从而得到最终问题的解。较常见求解方法有NewtonRaphson方法。3.4 结构的动力响应分析对于动态响应分析，由于结构具有加速度，必须考虑动力效应，可等效用附加惯性力的平衡体系来考虑。此时需要系统的质量矩阵，而梁结构为连续系统，为将其转化为有限自由度，并求其数值解，故需将其离散化，变为等效的多自由度振动系统来考虑。而对于梁上的集中质量块，则应特殊考虑，将其分配与整合。装配好质量矩阵与刚度矩阵及外荷载列向量，就可以利用数值解法，求得动力响应的各阶主振动（包括主振动频率及主振型）4 程序设计思路及框图4.1变量说明：XYZ： 梁单元两端节点的坐标XLT： 计算所得的梁单元的长度ISECT：指示梁截面类型，其中ISECT=1 表示矩形截面梁ISECT=2 表示管形截面梁TOL： 所设置的阈值，判断计算精度NPT： 选取的积分点的个数SR 数组存储高斯消去法的系数JNSR 数组SR中条目的个数NTABLE=-1 计算节点力NTABLE=0 计算所有高斯点上的应力IMAX： 指示质量类型，其中IMAX=1时，表示集中质量；IMAX=2时，表示连续质量4.2子程序BMEL调用的外部子程序有：DIRCOS,LENGTH,TRANSF,ENDREL,SRIFNL,STIFF,ECHECK,MASS,ATKA,ADDBAN,ADDMA,COLHT等由于有些子程序不能在程序模块AD32中找到，故以下就对于一些能找到，且经常被调用的外部过程作一个简短介绍，即只介绍其调用格式及程序功能。 1.LENGTH：（1）格式说明：LENGTH (XLT,XYZ)（2）程序功能：计算在三维空间，给定的两点间的距离，即为已知两端节点的梁单元的长度。2.TRANSF：（1） 格式说明：TRANSF (XYZ,AS,BS,T,PDISP,GAMA,ITONLY)（2） 程序功能：将局部坐标系下的质量和刚度矩阵转化为整体坐标系的质量与刚度矩阵3.ENDREL：（1）格式说明：ENDREL (AS,IMOMNT,SREL,DISP,NDRL,IENDRL)（2）功能说明：（a）计算压缩位移 （b）计算并储存必需矩阵以追踪压缩位移 （c）压缩刚度矩阵 4.STIFNL：（1）格式说明：STIFNL(DISP,PROP,WA,AS,RE,ICS,ISHEAR,SR,RERIT)（2）功能说明：计算几何或材料非线性刚度矩阵和非平衡荷载列向量 5.STIFF：（1）格式说明：STIFF (AS,XLT,E,G,XI,YI,ZI,AREA)（2）功能说明：计算线弹性刚度矩阵6.MASS：（1）格式说明：MASS (AS,XLT,XI,YI,ZI,AREA,DEN)（2）功能说明：计算常质量矩阵4.3子程序BMEL各子步骤的具体实现过程由前面图1已知子程序BMEL总的执行框图，可根据功能将该子程序分解为几个更小的程序模块。1. 存储单元信息2.计算整体坐标系下的积分点3.装配线性刚度矩阵4.装配质量矩阵5.转配最终的非线性刚度矩阵并向整体坐标系转换6.单元应力计算以下对其几个重要程序模块进行详细的代码分析。1.存储单元信息单元信息：节点i,j,k坐标的坐标分别为i(x(i),y(i),z(i)，j(x(j),y(j),z(j)，k(x(k),y(k),z(k)，将其存储在XYZ数组中，以便进一步处理。其过程如图4所示。图4 节点坐标信息单元坐标系信息：由于空间各杆件(梁)单元的局部坐标系不相同，在进行结构分析时，需要建立统一的总体坐标系，程序中通过节点坐标系来指示单元坐标系。其具体执行过程如5图所示图5 节点坐标系信息2计算整体坐标系下的积分点 （对于线性单元不适用）计算积分点个数 NPT= NPT=INTX*INTY*INTZ设置阈值 TOL=1.0D-82.1计算在RST坐标系在整体坐标系中的单位向量实现原理：如图6所示的三节点所确定的RST坐标系。图6 三节点所确定的RST坐标系(1)计算节点与节点组成的杆件的单位向量，即R方向的单位向量。节点与节点的坐标差值为： （10）节点与节点的距离为：（11）故由节点与节点组成的梁单元的单位向量即方向余弦为（12）（2） 计算平面外的单位向量，即T方向的单位向量节点与节点的坐标差值为：（13） 计算向量（14）（15）令，根据向量叉乘的几何意义，乘积所得向量与向量、同时垂直。故将乘积所得的向量，进行单位化后，即得平面外的单位向量，也即是T方向的单位向量。（16）故有：（17）单位化后得到平面外的方向余弦为（18）（3）S方向的单位向量的计算由于RST坐标系为正交坐标系，故可方便地由已知的R方向与T方向的单位向量的乘积得到（乘积向量与两相量垂直，且两单位向量的向量积的仍为单位向量）。即有：（19）进而得到单位向量的各分量为：（20）流程图如图7 所示。图7 求RST方向的单位向量的流程图2.2 分别就截面几何形状为矩形截面或为管形截面计算在整体坐标系下的积分点的位置采用三节点梁单元，故需取节点与节点的中点作为梁单元的中间节点。即有： 对于矩形截面梁或管形截面梁，它们的几何参数可由以下程序段得到，如图8所示，其程序流程图如图9所示。DSECT1=PROP(3,MTYPE) DSECT2=PROP(4,MTYPE) DRINT=XLT/DBLE(FLOAT(IX) IF (ISECT.EQ.2) GO TO 213 DSINT=DSECT1/DBLE(FLOAT(IY)DTINT=DSECT2/DBLE(FLOAT(IZ) GOTO 214 213 RMEAN=0.5*(DSECT1 + DSECT2)DRAD=0.5*(DSECT1 - DSECT2) DANG=4.*DATAN(1.D0)/DBLE(FLOAT(IZ) IF (INTZ.EQ.8) DANG=DANG*2.214 KINTP=0 图8 计算矩形截面与管形截面梁的几何参数程序段图9 矩形截面梁或管形截面梁几何参数计算的流程图2.3 计算积分点的位置坐标程序分别就矩形截面梁或管形截面梁计算积分点的位置坐标，其流程图如图10 所示图10 积分点位置坐标的计算流程图接下来的程序段，控制输出积分点的位置坐标（写入文件）矩形截面梁 若INTLM0 输出积分点位置圆管形截面梁 若INTLM0 输出积分点位置将数据写入文件，通过文件来传输数据，以实现在不同子函数中实现数据的通信。然后开始初始化工作向量，通过调用各种函数（包括线性单元刚度矩阵、非线性单元刚度矩阵，质量矩阵、及荷载列向量）转配整体刚度矩阵。由于涉及非线性分析（包括材料非线性和几何非线性）和动态响应分析，处理起来较复杂。下面分别简要说明。对于非线性分析：刚度矩阵在具体处理时初始矩阵只包括线性刚度矩阵，在求解迭代运算时，不但用非线性项进行修正，即利用NewtonRaphson方法或其改进方法。对于动态响应分析，需考虑系统的质量，故需形成单元质量矩阵，然后组装为系统整体的质量矩阵。3.装配线性刚度矩阵将局部坐标系下的单元刚度矩阵，转换为整体坐标系下的单元刚度矩阵，再根据节点自由度的整体编码形成的转换矩阵，将单元刚度矩阵组装为整体刚度矩阵。这里只考虑材料的弹性部分，装配的刚度矩阵为线性。其程序流程图如图11 所示图 11 装配线性刚度矩阵程序的流程图4.装配质量矩阵分别对集中质量和连续质量进行离散化。对于连续质量的离散化又分截面形式为矩形截面或管形截面形式。变量IMAX 指示质量类型，其中IMAX=1时，表示为集中质量；IMAX=2时表示连续质量。程序执行流程图如图12所示。图12 质量矩阵装配的程序流程图 图12续 质量矩阵装配的程序流程图5.转配最终的非线性刚度矩阵并向整体坐标系转换 刚度矩阵包括线性部分与非线性部分，每次迭代时，更新非线性部分刚度矩阵，从而修正整体刚度矩阵。程序流程图略。6.应力计算应力计算包括几何上与材料上的线性应力计算和几何上与材料上的非线性应力计算。由于所求得的应力为整体坐标系下的应力，故欲求得各梁单元上的应力还需向局部坐标系下进行转换。程序流程图略。5小结: 本次课程设计在理论与实践之间搭起了一座桥梁，让我们有机会在学习计算力学理论知识时，运用FEAP（有限元分析程序）进行实际计算操作。虽然所给程序代码有些过于复杂，而且程序语言为FORTRAN77,里面有过多的非结构化的程序结构，另外程序中的数据的数据通信大量是通过读写文件，和开辟内存共用区来实现，而且变量个数非常多，注释又不明晰。这些都给本次课程设计带来了不少挑战。所以本此课程设计不仅是考察知识与能力，而且更是考察耐力，它要求你一遍又一遍的分析，一遍又一遍地去了解背景知识及程序实现原理。通过本次课程设计，我不仅查阅了很多相关书籍，加深了自己对有限元法的理解，增广了自己的见识。而且为了更好的读程序，自己又把FORTRAN77的语法又啃了一遍，这样一遍又一遍地查阅资料、分析代码，感觉自己获益匪浅。同时在翻阅资料时，网上有许多FEAP的现成软件，有些还是可视化软件，还可查阅里面所有代码，这给自己扩展视野很是有好处。虽然此次课程设计中多少遇到了不少困难，但我相信这个过程将有助于培养我们解决实际问题的能力，及科研的精神。参考文献1 王勖成. 有限单元法. 清华大学出版社, 2003.072杨庆生.现代计算固体力学，科学出版社，20073曹志远、张佑启.半解析数值方法，国防工业出版社，19924傅永华.有限元分析基础.武汉大学出版社，20035张洪武 关振群.有限元分析与CAE技术基础.2004.026谭浩强 天淑清.FORTRAN 77结构化程序设计.清华大学出版社2001.067屈均利 韩江水.工程结构的有限单元法.西北工业大学出版社2004.09