结构有限元程序设计相关问题

1、结构有限元程序设计相关问题有限元结构分析程序理论及模型化1 有限单元法2 弹性理论的基本公式3 假定位移场的有限元公式推导4 结构模型化及程序模块划分5 程序的单元分析模型6 约束与运动分析的模型化7 结构总刚度阵8 荷载模式与荷载工况9 数据管理与工程集成系统问题的求解模型描述 连续体模型：弹性理论是一种基本的模型，采用偏微分方程精确而简洁，但由于荷载和结构的边界条件的复杂性，使得偏微分方程的解析求解成为不可能，工程师无法得到实际应用。 离散化模型：求近似解往往应用有限自由度的结构来代替连续体。这对于计算机是十分适合，对于求解有限多未知数的题目非常有效，有限元在计算工具强有力的支持下得到了广

2、泛的应用。有限元法的求解思路从单元的角度观察，每一个单元只通过节点与结构的其他单元相连接。这些节点为单元的出口节点，出口节点，而其有关的位移为单元的出出口为移口为移。对于任一个单元，可以写出如下形式的方程组1 有限单元法11112211211222221122mmmmmmm mmmkdkdkdrkdkdkdrkdkdkdr其中 是单元出口自由度数出口自由度数，即结构与单元的相互作用是通过这 个出口为移表示，这些位移即为结构与单元连接节点位移 ， 是节点反力，与出口 为移相对应。 m(1,2,)id imid(1,2,)ir imm方程系数 称为刚度系数刚度系数，如写成矩阵形式，就是单元刚度阵，

3、记作 ijk 1112112mmmmmkkkkkkk进一步单元刚度阵与出口位移、反力关系，可表示如下 kdr有两个出口节点，每一个节点3个自由度，2个线变位，1个角变位。出口位移向量为 111222Tduvuv相应的单元刚度阵为说明 上述位移和反力都是在单元局部坐标系里描述的，相应于单元刚度阵为单元局部刚度阵。 当坐标进行变换时，单元刚度阵也需要进行相应坐标变换。 如何计算单元出口刚度阵，以及构造新的单元等等是有限元的一个研究领域。对结构中所有单元进行组装，就得到结构的正则方程，为 KDR式中 是整个结构，即有限单元拼装组成的结构总刚度总刚度阵阵。 为结构的未知位移向量，称为总位移向量总位移向

4、量。 为结构总外力向量总外力向量。 K D R基于有限元的位移法求解思路当结构的正则方程建立之后，除了总位移向量，所有的量都为已知。采用标准的线性代数方程求解器求解。计算出总位移向量之后，各个有限元出口位移就可以由总位移向量中的各个分量及结构的约束条件计算出。各个单元内的位移及内力可以根据其出口位移用插值方法，或者通过弹性力学的变分原理求出。 2 弹性理论的基本公式 平衡方程，应变-位移关系、应力-应变关系在三维弹性条件下平衡方程（静平衡、动平衡有惯性项） xyxxzxxyyyzyyzxzzzfxyzfxyzfxyz 在三维坐标系下，应变-位移关系为 000000000 xyzxyyzzxxy

5、uzvyxzxzx 其中 ， ， 为正应变； ， ， 为剪切应变； ， ， 为位移分量。 xyzxyyzzxuvw应力-应变关系为 0()E将应力-应变、应变-位移关系带入平衡方程，不计初始应变，可以得到 2221/0121/0121/012xyzeufGvxevfGvyefGvz其中 为三维拉普拉斯算子 为体积变化率。这套方程称为拉梅方程，是以位移作为未知数表示的平衡方程。 22222222,xyz uvexyz力学中两个重要要原理在所有满足内部连续条件及边界上的位移约束条件的位移中，满足平衡方程的位移使得总势能取驻值。如果驻值极小，则平衡是稳定的。如果一个内力分布与外力相平衡，则对于任一个

6、可能位移，内力与外力作功之和为零；反之，如果对于任一个可能位移，内、外力作功之和恒为零，则该力系处于平衡之中。最小总势能原理是与弹性力学位移法相对应的变分原理。弹性力学的变分原理 总势能由结构的变形能 以及外力势能 两个部分相加而成 UVUV 根据最小总势能原理，在其平衡位置附件取一个可能位移的变分，有 0UV此为虚位移原理的一个表达形式。 弹性体总势能的具体表达式 01()2TTVTTVSEdVfu dVTu ds 小结：上述弹性理论公式是对于连续体，弹性体内任一点位移都是独立未知数，因此具有无穷多的未知数。3 假定位移场的有限元公式推导 有限元将弹性体剖分成单元，将未知数选择为节点处位移。

7、单元内部的位移由节点处的位移插值得到。插值函数依赖于有限个参数(节点处的位移)。与瑞利李兹法相似，未知函数在整个区域依赖于有限个参数。 通过变分原理可以推导出求解这些参数的代数方程组。有限位移法就相当于对最小总势能原理使用李兹法。利用上述思路建立结构正则方程。设单元内任一点位移由下式给出 TuuvNd其中 是单元出口位移向量， 为单元形函数矩阵。 d N利用应变-位移关系有 Bd其中 为应变出口位移映射矩阵。 B利用上述关系，任意一个单元的变形能为 01()()2eeTTTTeVVUdBEB dVddBdV eeTTTTeVsVdNf dVdNT ds单元的外力势能为对于有 个单元的结构，结构

8、的总势能为 em 1emTeeDP edTD其中 为转换矩阵。 eT对于任一个单元，出口位移与结构总位移的关系可由下式给出 结构的总势能成为 1011()2()eeeeeumTTTeeeVmTTTTeeVVTTTTuSSDTBEB dv TDDTBdvNf dvDNT dsDPDNT ds 此式即为结构有限元模型下结构总势能的表达式，该模型总势能由节点独立自由度(位移)所决定。 式中，最高二次项为 11()2emTTeeeDTk TD其中 eTVkBEB dv为单元在局部坐标系中的出口刚度阵。 整理结构总势能的表达式得到 12TTDKDDR 式中 为结构的总刚度矩阵，由 个单元刚度阵相加而成

9、Kem 1emTeeeKTk T由最小总势能原理，处于平衡位置的位移 满足下面的方程组 D120,0,0nDDD整理后的结构总势能的模型按上式求偏微商之后，便得 KDR这就是结构有限元位移法的正则方程。 4 结构模型化及程序模块划分 模型化的角度 模型化的方法 模型化的内容模型化的两方面内容工程实际对程序功能要求，计算模型能与实际结构工作状况一致。程序应具有的结构模型处理能力，这是程序的结构模型化所要解决的问题。（结构的模拟、载荷的模拟等等）对于已有的结构分析通用程序，在程序功能范围内抽象手中具体一个结构计算模型。这是具体结构模型化的问题。（动力分析模型、静力分析模型）（高低频考虑、单元划分的

10、稠密等等）工程实际对结构分析程序功能要求 通用性，尽可能处理各种复杂结构计算。 可读性，具有二次开发的接口。 提供各种模型化手段，避免方程病态保证求解精度。 模型的建立方便性。命令流、图形建模、与其它有限元模型的兼容等等。 可扩充各种新的单元。 易于扩展到不同分析领域，如优化、稳定、控制等等。 容错性、对于错误尽可能给出报告信息。 计算结果的动态可视化显示等等。多目标程序的核心问题 数据库管理数据库管理 如：数据文件就是一个简单的数据库，同一有限元程序的不同的版本其输入文件数据库的构造可能有所不同不兼容等等有限元程序都的数据管理方式 在有限元程序中，对于计算机外存（硬盘、软盘）开设一个虚拟的、

11、很大的可以任一处随机读写的“数组”BANK库。对于BANK的读写是通过调用子例程进行访问的。有效信息区有效信息区未使用的部分未使用的部分BANK:BANK:IPIPIPTRGIPTRG（1 1）* *512512IPTRGIPTRG（2 2）有限元程序都的数据管理方式 对于计算机内存采用动态管理，应用内存文件管理系统进行，并完成同BANK的数据交换。内存文件管理系统是数据管理及存储管理的一个辅助程序系统。它同有限元分析程序没有关系。 文件构造图：文件构造图： 1 12 23 34 45 56 67 7NRECNRECNRECNRECNRECNREC7 7个字个字1f1f2f2f3f3f导向器区

12、导向器区LDIRECLDIRECLENLEN文件总长文件总长记录区记录区注：文件在内存注：文件在内存IDTIDT数组中存区的基址为数组中存区的基址为IBSWIBSW。IDT(IBSW+1)IDT(IBSW+1)：文件的第：文件的第1 1个字，其中存文件在内存中占有的长度个字，其中存文件在内存中占有的长度LENLEN，以，以2 2字节为单位。字节为单位。IDT(IBSW+2)IDT(IBSW+2)：文件的第：文件的第2 2个字，记为个字，记为IPTIPT。内存中前面。内存中前面IPTIPT个字不是立即可个字不是立即可用的资源；而后面用的资源；而后面LENLENIPTIPT个字是立即可用的资源。个

13、字是立即可用的资源。IDT(IBSW+3)IDT(IBSW+3)：文件的第：文件的第3 3个字，记为个字，记为IVIV。表示文件尚有。表示文件尚有IVIV个字的资源可供个字的资源可供使用。使用。 IDT(IBSW+4)IDT(IBSW+4)：文件的第：文件的第4 4个字，记为个字，记为ISTATISTAT。冻结标记，当为。冻结标记，当为1000010000时表示时表示冻结。此时只可使用立即可用的资源。所有记录原位置不动。冻结。此时只可使用立即可用的资源。所有记录原位置不动。IDT(IBSW+5)IDT(IBSW+5)：文件的第：文件的第5 5个字，为本文件的记录数个字，为本文件的记录数NREC

14、NREC。IDT(IBSW+6)IDT(IBSW+6)：文件的第：文件的第6 6个字，记为个字，记为JFILBSJFILBS。文件的库基址。由此可以找。文件的库基址。由此可以找出本文件在库出本文件在库BANKBANK中导向器存区的基址。中导向器存区的基址。 IDT(IBSW+7)IDT(IBSW+7)：文件的第：文件的第7 7个字，本文件导向器的长度，记为个字，本文件导向器的长度，记为LDIRECLDIREC。说明：说明：IBSWIBSW是本文件在内存是本文件在内存IDTIDT数组中存区的基址。所谓基址是指向存区之数组中存区的基址。所谓基址是指向存区之前的地址数。前的地址数。 有限元分析数据库

15、（模型/载荷） 节点坐标 单元管理 约束条件 物理性质等等 附加质量 荷载管理5 程序的单元分析模型 单元分析与结构体系之间的接口原则 结构体系提供的功能：按单元随机调取单元表示，即单元的连接节点及单元的性质，并且能够完成单元的局部坐标到结构总体坐标系的坐标变换，提供单元各出口节点在本身局部坐标系中的坐标值。 单元向结构体系提供：单元局部坐标系中表示的单元出口刚度阵，单元内荷载引起的单元出口节点上的反力。 结构体系在求出总位移向量后，向单元提供在其局部坐标内的单元出口位移值。 单元分析在得到出口位移后，可计算单元内力分布。 结构体系将内力、变形按格式输出。结构体系方面单元表示、单元性质、局部坐

16、标内节点坐标单元内荷载总外力向量总刚度阵总位移向量单元出口位移输出单元分析方面单元出口刚度阵单元出口反力单元内力 有限元离散后，未知数一般是节点位移。由于结构构造方面的要求，往往有些单元刚度很大而其它一些单元刚度很小，相比之下这些刚度大的单元不能发生变形。因此这些刚度大的单元出口位移就不是相互独立的了（具有相关性）。结构的自由度与独立位移就成为一个重要问题。 复杂结构结构常常需要多种单元组合进行模拟，由于不同单元出口位移性态不同，如何处理不同单元之间（如梁单元和杆单元）的连接也是一个重要问题。有时不同结构部件相互连接并不是完全固接，其计算模型的选取也是需要处理的问题。 6 约束与运动分析的模型

17、化以上这些问题的核心是不同节点位移之间的关系（节点位移之间的相关性），对于这些问题的处理可采用节点位移规格数及主从关系进行。6.1 结构节点位移规格数 在采用有限元位移法时，结构的变形状态是以节点位移作为基本未知数的。 对于一个三维空间结构，计算模型可能由多种单元组装而成，节点位移允许有6个自由度（3个线变位，3个角变位），即Txyzu 由于有不同的结构形式以及不同的边界条件或约束条件，节点任一个位移都可能有以下5种规格，分别以0，1，2，3，4的代码表示规格 几何可动 独立 相关 指定为零 指定值 代码 0 1 2 3 4 在这5种位移规格中，只有独立位移才进入结构的总未知数向量。 相关位移

18、虽然也是未知数，但它是由独立位移和指定位移所决定的，不是独立未知数，不计入结构总位移向量之中。规定结构总位移向量是由独立未知数所组成。 根据总位移向量与指定位移，可以计算出一个节点的任一个位移，称总位移向量及指定位移为控制位移。 上述节点位移规格数都是在结构全局坐标系中描写的，这个坐标系称作绝对坐标系。结构节点坐标是由绝对坐标系来定位的。 位移规格数的引入一方面可以处理非齐次节点位移，另一方面，它不仅能协助把指定为零及几何可动位移排除在总位移向量之外，还可以处理指定位移和相关位移，并也将其排除在总位移向量之外。 这对扩大解题范围、减少独立未知数，避免方程病态，从而提高求解精度有很大作用。6.3

19、主从关系及其应用对于如图的空间结构，假设结构中间有一块的地方的刚性远较其它杆件的刚性大，在图中用斜线表示。现在要分析这个结构的受力变形状态或者振动固有值问题。该结构为空间杆系，每个节点的未知位移的个数为6。 这个刚性块上与其它杆件相连共有10个节点，图中ABCDEFGHIJ来代表。直观的处理方法是把10个节点的位移都认为是独立未知数，为了实现刚性，则在这些节点相互间用一些其刚性远较其余杆件大的杆件人为地连接起来，组成新的体系来模拟，之后进行计算。这个处理方法虽然可行，但具有明显的缺点。显然，每个节点6个独立未知数，刚体部分10个节点就有60个独立未知数；然而从刚体运动学可知，一个刚体只有6个独

20、立位移。这样，54个独立未知数是多余的了。未知数的增加伴随而来的结构总刚度阵的带宽也随之增加，无疑对计算资源增加负担。新的解决方案 中间刚性块只有6个独立位移，应该可以把其中54个位移未知数看成这6个独立位移的相关位移，而在生成总刚度阵时将这些相关位移排除在总位移向量之外。这样总刚度阵的维数就减少了。 在计算出独立位移之后，根据位移的相关性质，可以反过来计算出这些相关位移的数值。 好处：可以大量节省总刚度阵的存区，减少总刚度阵三角化分解时间，而且避免了不必要的方程病态，改善了计算精度。 上所讲述的各节点位移主从关系是完全相关的情形。对于其它许多实际课题，并不是从节点的全部6个位移都是相关位移，

21、也就是部分相关的情况。 7 结构总刚度阵 结构分析的关键一步是组装成整个结构的正则方程 KDR 从公式上看，方程是非常简洁的，但由于复杂结构未知数很多，因此总刚度阵一般是一个稀疏矩阵，对于怎样组织稀疏矩阵以及有效地执行计算，是一个重要问题。特别是当节点自由度很高的情况下。可以将此划分为三个部分 总刚度阵的组织、未知数的编排 总刚度阵的装配 总刚度阵的求解总刚度阵的组织、未知数的编排 总刚度阵的体积很大，所以要用特殊的数据结构以及存储格式并且选择有效的算法。 算法的选择与问题有关 对于非线性迭代来说，用迭代法求解总刚度阵为好（刚度阵在每一迭代步上可能要重新生成） 对于线性问题，要分析多种复杂工况

22、的课题，则一般都要选择直接法求解总刚度阵。3.23.2 -1.0-1.04.34.3 0.80.82.02.0 4.04.00.10.18.88.8 4.14.14.04.0 2.12.1 -0.8-0.81.21.20.30.30 00 00 02.82.8 1.51.50 00 03.13.1 0.70.70 03.73.7 -0.2-0.20 00 0-1.2-1.20 0-2.1-2.14.14.17.27.2总刚度矩阵组织安排 由于总刚度阵是对称正定阵，可以选择修正的平方根方法，或LDLT方法求解。首先将总刚度阵分解为 的形式，其中 是对角阵，而 则是对角元素等于1.0且上半矩阵全为

23、零元素的下半三角阵。 TKLDL D L可采用“一维变带宽”的方案进行总刚度阵的组织。即如图324中的矩阵，在对角线左面有部分区域有不为零的元素，但在左侧的边界线之左还有一区全为零元素。由于是对称矩阵。因此只要存储左半部分的元素就可以了。一维变带宽矩阵的实际存储方式是将各行首尾相连，成为一个一维数组。该一维数组当然很长，常用的方法是用对角元地址数组IDIAG1:ND辨别每一行的分界线，如图325所示，这样可以将带外的零元素排除在外。但这个一维数组的总刚度阵在内存一般是放不下的，因此需要放在BANK库中。则计算机整型字的限制32767可能不够对角元地址之用。在这里不采用IDIAG数组，而是存放每

24、一行的半带宽（不包括对角元在内）。1 13 35 58 81010111119192323262634340 01 11 12 21 10 06 63 32 27 7对角元地址数组半带宽数组1 13 35 58 81010111119192323262634340 01 11 12 21 10 06 63 32 27 7对角元地址数组半带宽数组3.23.2 -1.0-1.04.34.3 0.80.82.02.0 4.04.00.10.18.88.8 4.14.14.04.0 2.12.1 -0.8-0.81.21.20.30.30 00 00 02.82.8 1.51.50 00 03.13.

25、1 0.70.70 03.73.7 -0.2-0.20 00 0-1.2-1.20 0-2.1-2.14.14.17.27.2总刚度矩阵总未知数编排对于半带宽有重大影响 寻求未知数优化编排的方法有很多，可以采用比较优的总未知数编排，如ZXY-优先序： Z坐标值小的节点编在先在Z相同时，X坐标值小的节点先编号 如Z、X坐标值都相同，Y小的先编号总刚度阵的装配 单元分析只是在局部坐标中提供单元的出口刚度阵。将单元出口刚度阵进行转换，累加到结构总刚度阵中去是结构分析的关键之一。 单元局部坐标中的出口位移是与结构绝对坐标中的节点位移不一致的，当然用接口节点在绝对坐标中的6个位移，总可以表示出口节点的出

26、口位移，这是一个转换到关心位移的问题。 但节点在绝对坐标中的位移（关心位移）并不是结构总位移向量中的独立位移，还有转换到控制位移，即独立位移及指定位移的步骤。 在上面的基础上单元刚阵才能累加进总刚度阵去，以及将指定位移的影响累加到总外力向量中去。 对于任一个单元，出口位移与结构总位移的关系可由下式给出 edTD其中Te为转换矩阵。 在前述章节中推导有限元位移的一般公式中，只是设单元的出口位移可以由控制位移来表示，在原则上说很简明；但在数值上看是远远不够的。矩阵Te是如何构造的，它是怎么算出来的，都是必须解决的。 矩阵Te是由两个转换合成的 单元出口位移到关心位移的转换 关心位移到控制位移的转换

27、 当总位移向量编排好以后，还有计算机生成总刚度阵的半带宽数组问题，以及由此而安排总刚度阵的存区问题。这些工作也只有在将所有单元的控制转换阵Te及其构造算出后才能执行。总刚度阵的求解 当总刚度阵按一维变带宽的格式组织好后，执行三角化运算则是一个纯粹的数值分析问题。但从力学的角度应当提出一定的要求，即总刚度阵病态条件的查验。 对于多种荷载作用下的结构强度计算，振动或稳定计算，求解联立方程组要多次执行。这对于直接法是比较有利的，一次三角化之后，把结果保存下来可以多次求解。8 荷载模式与荷载工况 一个工程结构所受的载荷一般总是多种多样的，来自各种不同因素。对于许多工程领域，结构本体的计算理论可以说是比

28、较统一的，但是载荷则有很大差异。 对于结构（如建筑结构）通常要考虑的荷载有 自重 工作荷载 不同方向的风荷载 地震荷载 有限元模型介绍人字架人字架臂架臂架基座基座模型的主要结构钢索钢索工况与载荷介绍浮吊承受的主要荷载载荷类型工况类型载荷载荷工况工况DLPcSWLGgGsFwxFwzFswxFswzPHhPHbPHxPHyPHz臂架仰角臂架仰角工况类别工况类别11111171无风2111.11.11.171无风3111.11.11.1171有风4111.11.11.11171有风51111180无风6111.11.11.180无风7111.11.11.1180有风8111.11.11.11180

29、有风9111.11.11.111171有风10111.11.11.11171有风111111110放置12111.11.11.120无风13111.11.11.150无风计算工况列表DL自重载荷 Pc侧倾载荷 SWL安全工作载荷 Gg主钩重量 Gs钢丝绳重量 Fw工作风载荷 Fsw非工作风载荷 PHh回转惯性力 PHb变幅惯性力 PHx、PHy、PHz船体运动引起的惯性力 这些荷载来源不一，但是在结构实际工作中却一般是同时作用的。对于结构的线性计算理论来说，完全可以对于各个单种因素的荷载进行求解；而结构的实际工作状况则可由这些单种因素的结果线性组合叠加而得。 可将这类单种因素的荷载成为荷载模式，而将结构实际所受到的复合外力成为荷载工况。荷载工况是由荷载模式线性组合而成的。 对于用户来说，按荷载模式来组织数据是比较方便的。对于每个荷载模式的数据则只需要说明一遍，当各种荷载工况中有这个荷载模式参加组合时，只需要引用一下就行了，也就是说输入一个线性组合的系数就可以了 。 设结构受到有LDMDL个荷载模式（LoaD MoDuLes）的作用；优这些荷载模式组合成LDCS个荷载工况（LoaDing CaSes）。在数学上，线性组合可以用一个组合矩阵Q表示