空间刚架梁有限元程序及优化设计--liqifu40页

1、 中国仿真科技论坛电子期刊 NO.9空间刚架梁有限元程序及优化设计Simwe会员 liqifu摘要： 随着计算机技术的普及和计算速度的不断提高，经历了三十多年的发展历史，理论和算法都已经日趋完善。有限元分析在工程设计和分析中得到了越来越广泛的重视，已经成为解决复杂的工程分析计算问题的有效途径。有限元的核心思想是结构的离散化，可以解决很多实际工程需要解决而理论分析又无法解决的复杂问题。空间刚架梁结构是工程中经常遇见的结构，有限元方法是工程实际中经常运用的方法，有限元方法可以很好地模拟和分析结构，达到工程需要的要求。对结构的模拟和分析，一个主要的目的就是对结构进行优化，以达到经济效益。优化理论在结

2、构分析中的作用是越来越突出。运用结构优化理论，对空间刚架梁结构通过有限元分析，在计算机中编制有限元程序，可以达到对结构进行优化的目的。关键字：空间刚架梁 有限元 程序设计 结构优化1.引言1.1 论文背景及研究意义人们进行力学分析的方法有很多种，但归结起来可分为两类，即解析法和数值法。由于实际结构物的形状和所受荷载往往比较复杂，除了少数简单的问题之外，按解析法求解的非常困难，所以数值法已成为不可替代的广泛应用的方法，并得到了不断发展。有限单元法就是伴随着电子计算机技术的进步而发展起来的一种新兴数值分析方法。它的数学逻辑严谨，物理概念清晰，易于理解和掌握，应用范围广泛，能够灵活地处理和求解各种复

3、杂的问题，特别是它采用矩阵形式表达基本公式，便于运用技术编程运算。这些优点赋予了有限元强大的生命力。随着计算机技术的普及和计算速度的不断提高，经历了三十多年的发展历史，理论和算法都已经日趋完善。有限元分析在工程设计和分析中得到了越来越广泛的重视，已经成为解决复杂的工程分析计算问题的有效途径，如今它已成为广大科技工作者的有力工具，解决了大量科学研究和工程技术问题。空间刚架梁问题是工程实际中经常遇见的问题，可以用有限元方法对空间刚架梁结构进行模拟和分析。对空间刚架梁结构进行模拟和分析的一个重要的目的就是对结构进行优化，以产生经济效益。1.2 有限元概论1.2.1有限元发展与简介有限元法是大型复杂结

4、构或多自由度体系分析的有力工具，近20年来已广泛地用于：工程结构、传热、流体运动、电磁等连续介质的力学分析市，并在气象、地球物理、医学等领域得到应用和发展。电子计算机的出现和发展，使有限元法在许多实际问题的应用变成现实，并且有广阔的前景。有限元法的基本思路是将结构物看成由有限个划分的单元组成的整体，以单元结点的位移或结点力作为基本未知量求解。有限元法基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。在有限元方法

5、中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。1.2.2 有限元方法的基本思路和解题步骤1.建立积分方程，根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。 2.区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之

6、间的关系之外，还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。 3.确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插值条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元具有规则的几何形状，在选取基函数时可遵循一定的法则。 4.单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数(即单元中各节点的参数值)的代数方程组，称为单元有限元方程。 5.总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则进行累加，形成总

7、体有限元方程。 6.边界条件的处理：一般边界条件有三种形式，分为本质边界条件(狄里克雷边界条件 )、自然边界条件(黎曼边界条件)、混合边界条件(柯西边界条件)。对于自然边界条件，一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法则对总体有限元方程进行修正满足。 7.解有限元方程：根据边界条件修正的总体有限元方程组，是含所有待定未知量的封闭方程组，采用适当的数值计算方法求解，可求得各节点的函数值。 1.3结构的优化设计结构优化设计又名结构的综合，这是因为一个结构设计只有在综合考虑力学原理，数学方法与手段，材科与工艺，经济与使用等各个方面因素的情况下，才能得到令人满意的结

8、果。因此，怎样从工程实际问题中选择评价设计方案优劣的目标函数，确定描述设计方案的设计变量，给出反映问题的性能与状态等要求的约束条件，从而形成优化设计的计算模型是一个十分重要的步骤。它不仅直接影响到一项工程设计的技术经济效果，而且与结构优比设计能否在工程中推广应用有着非常密切关系。(一)选取合理的截面型式和刚度配置，将内力与位移的分布相对于材料的配置调整得愈是合理，结构的重量就愈轻。 改变截面的尺寸，使在一定荷载作用下，各横截面的最大纤维应力相同，从而达到充分发挥材料的承载能力，减轻结构重量的目的。(二)合理地组合各种材料，结合结构的受力特性和几何形状，愈是能充分利用不同材料的优良性能，结构重量

9、就愈轻。结构优化设计是随着最优化方法的发展而发展起来的。在电子计算机出现之前，由于受到计算手段的限制，人们只能利用以松弛变量、拉格朗日乘子、等式约束消元等措施的微分法与变分法，或者出自于“同步失效”的概念，把问题变成求解一组联立方程式的方法，寻求单独杆件的最轻设计。 “同步失效”是假定最优结构就是在外荷载作用下，能使所有可能的破坏模式同时实现的结构，从单个杆件的优化设计推而广之。如果把整个结构看作一个元件，而结构中的每一根杆件的破坏可以看作一种可能的破坏模式，则自1968年以后进一步发展起来的准则法，实际上就是“同步失效”设计的演变与发展。从60年代开始，由于有限元方法的出现，使得复杂结构的分

10、析有了可能，以及数学规划的渗透和被引入，大多数结构优比设计的问题可以抽象非线性规划的一个命题，从而使结构优化设计，发展非常迅速，理论研究与工程应用所涉及的范围也愈来愈广。目前，国内外采用的结构优化设计方法种类繁多。最为常用的优化方法，一般来说，可以分为两大类：(一)准则方法准则方法是从结构力学原理出发，建立一些最优准则，从而寻求用解析形式表达的结构设计的参数，或者通过直观的迭代运算决定结构各单元的截面参数。常用的准则方法有满应力准则法、位移准则法、能量准则法等。(二)数学规划方法数学规划方法是从解极值问题酌数学原理出发，运用数学规划中各种方法，求得一系列设计参数的最优解。当前已被引用的几种规划

11、方法有线性规划、非线性规划中的约束最优化方法与无约束最优化方法、动态规划、几何规划等。数学规划方法与准则方法两者互相渗透和互相借鉴，促使两条原来不同途径的汇合，出现了准则法与规划法联合应用的方法，把力学概念和优化方法密切地结合在一起，显著地提高了问题的解算效率，这是当前结构优化设计发展中的一个趋势。程序框图如图1.1所示图1.1程序框图有限元解的精度和收效性取决于求解的微分方程或采用的变分式和选取的单元。有三种误差来源，(a)由于域的近似性产生的误差，(b)由于解的近似方法产达的误控，(c)由于数值计算产生的误差，即在计算机中数值的积分和告入误差。通常确定这些误差不是一个简单的事情。但在一定的

12、条件下，对于一个约定的单元和问题则可以估算。我们需要对有限元解进行反复的比较，避免误差过大。该空间刚架梁有限元程序设计及优化设计的结果很大程度上依赖与约束条件，约束条件选择的合理，优化结果也就合理。1.4 本文主要内容及方法本文主要是对空间刚架梁结构进行有限元程序的优化设计，即用有限元的方法分析空间刚架梁（自行车刚架），求出目标函数，以达到节约材料的目的。利用虚位移原理建立空间梁单元的单元刚度矩阵，利用各结点编号，形成总刚度矩阵。再利用条件形成载荷矩阵，求出各杆端位移。根据临界载荷形成一个约束条件，选取杆件直径为设计变量，以总重量最轻为目标函数。将其编制成有限元优化程序，在Power Stat

13、ion 平台上运行。在约束条件下，通过循环语句对设计变量进行线性搜索，找到最优化直径，以求出目标函数。2.有限元方法对刚架结构的分析2.1 虚位移原理推导单元刚度矩阵2.1.1 虚位移原理有限元的核心就是建立单元刚度矩阵，有了单元刚度矩阵，加以适当组合，可以得到平衡方程组，这以后剩下的就是一些代数运算。虚位移可以是任意无限小的位移，它在结构的边界上必须满足运动学边界条件。如上图所示的物体，它受到外力 等的作用，记 ，在这些外力作用下，物体的应力为 ：（2.1）现假设物体发生了虚位移，在外力作用处各个外力相应方向的虚位移为 记 ：（2.2）由虚位移所产生的虚应变为 ：（2.3）在产生虚位移时，外

14、力已作用于物体，而且在虚位移过程中，外力保持不变，因此，外力在虚位移上所做的虚功是：（2.4）在物体的单位体积内，应力在虚应变上的虚应变能是：（2.5）整个物体应变能是：（2.6）虚位移原理表明，如果在虚位移发生之前，物体处于平衡状态，那么在虚位移发生时，外力所做虚功等于物体的虚应变能，即：（2.7）2.1.2 单元位移 对于空间应力问题有：（2.8）其中分别是结点沿方向的位移分量。用表示单元全部结点位移所构成的向量：（2.9）单元内任一点的位移可表示如下：（2.10）式中等为形函数，具有如下特性：在结点 ：（单位阵）在其他结点：（矩阵） 为坐标的函数。2.1.3 单元应变与应力当单元内任一点

15、的位移已知时，通过适当的微分运算可求出单元内任一点的应变，一般可表示如下：（2.11）根据广义胡克定理，单元应力可表示如下：（2.12）其中为弹性矩阵，将应变表达式代入上式得：（2.13）式中 是应力矩阵。2.1.4 结点力与单元刚度矩阵用结构力学方法求解连续介质的应力，以作用于单元结点上的等效集中力代替分布于单元边界上的应力，称为结点力。用 表示i点的结点力。结点力的个数和方向必须与结点位移保持一致，对于空间应力问题:（2.14）其中, , 分别是沿方向作用于结点的结点力。用表示单元全部结点力所组成的向量：（2.15）用虚功原理推导结点力的表达式。假想在单元中发生了虚位移，相应的结点虚位移为

16、 ，则有：（2.16）单元内发生的虚应变为：（2.17）结点力所做的虚功等于每个结点力分量与相应的结点位移分量的乘积之和：（2.18）用矩阵表示，即为：（2.19）在整个单元内，应力在虚应变上的虚应变能是：（2.20）把代入上式，得到：（2.21）根据虚功原理故：（2.22）上式对于任何虚位移都必须成立，由于虚位移可以是任意的，所以矩阵也是任意的，上式两边与它相乘的矩阵应当相等，于是得到： （2.23）将应力表达式代入得：（2.24）令：（2.25）于是：（2.26）上式建立了结点力和结点位移之间的关系。矩阵称为单元刚度矩阵，它的元素表示当单元发生一定的结点位移时，所对应的结点力。单元刚度矩阵

17、决定于该单元的形状，大小，方向和弹性常数，而与单元的位置无关，即不随单元或坐标轴的平移而改变。2.1.5 结点荷载为了能用结构力学方法求解连续介质的应力，所有分布荷载都必须代换以等效的结点荷载。用表示结点的等效结点荷载。对于空间应力问题：（2.27）其中分别是沿方向作用于点的集中荷载。用表示单元全部结点荷载所组成的向量：（2.28）2.1.6 用虚位移原理推导各种结点荷载算式。1.分布体积力设单元体积内承受的体积力为：（2.29）当单元中发生虚位移时，体积力所做的功为：（2.30）这应当等于等效结点荷载所做的功，即：（2.31）由上式可得体积力的等效结点荷载如下：（2.32）上式右边的重积分应

18、当在单元的整个体积内进行。2分布面力设单元e是靠近边界的单元，在其边界S上作用着分布的面力 当单元发生虚位移时，面力所做的功为：（2.33）它必须等于等效结点荷载所做的功，由此得到：（2.34）上式右边的面积分是在分布荷载所作用的表面S上进行的。2.2 形成空间刚架梁单元刚度矩阵1.将结构物离散成杆单元如图2.2所示：图2.2杆单元如图取结点为刚架梁单元，用节点位移表示应力和应变。单元节点列阵：（2.35）（2.36）节点力列阵：（2.37）（2.38）有几何方程：（2.39）有物理方程：（2.40）式中的D为梁单元的弹性矩阵。（2.41）由2.1虚功原理推导出梁单元的刚度矩阵为： （2.42

19、）将上式的推导代入得：（2.43）其中2.3 单元刚度矩阵的坐标变换单元刚度矩阵均是在局部坐标系中得到的。由于局部坐标是根据单元的几何形状选取的，故单元刚度矩阵的形式和推导都比较简单但不同的单元一般具有不同的局部坐标系，因而不能进行不同单元刚度矩阵问的混合运算。即得不到整个计算模型的有限元计算格式。解决这个问题肋唯一途径就是建立整体坐标与局部坐标间的坐标变换，并通过这种变换来得到单元刚度矩阵在整体坐标系下的显式。若某个单元刚度矩阵本来就是在整体坐标系下得到的，那就根本不存在坐标变换问题了。首先考虑单元在端点的三个杆端力分量，在局部坐标系中，它们是。现在推算，与之间的关系。设轴与轴的夹角分别为，

20、则轴在xyz坐标系中的方向余弦为，将杆端力，在轴上投影，可得杆端力为：（2.44）同理可得：（2.45）（2.46）综合以上三式，则有：（2.47）这就是在端点i由整体坐标系中的杆端力，推算局部坐标系中杆端力的转换关系式。其中两坐标的转换矩阵为：（2.48）参照上述方法，同样可以推导出以表示，以及以表示 的表达式，其转换矩阵也是。 综合以上分析，整体坐标系中的单元杆端力分量列阵与局部坐标系中单元杆端力分量列阵之间的关系，可用下式表达：（2.49）同理，可导出整体坐标系与局部坐标系杆端位移之间的转换关系：（2.50）其中：（2.51）称为单元坐标转换矩阵；它是个正交矩阵，故有：（2.52）则空间

21、单元在整体坐标系中的单元刚度矩阵为：（2.53）2.4 单元刚度矩阵的组装及整体分析根据全结构的平衡方程可知，总体刚度矩阵是由单元刚度矩阵集合而成的。如果一个结构的计算模型分成n个单元，那么总体刚度矩阵就可由各个单元的刚度矩阵组装而成，即 ：（2.54）是由每个单元的刚度矩阵的每个系数按其脚标编号“对号入座”叠加而成的。这种叠加要求在同一总体坐标系下进行。如果各单元的刚度矩阵是在单元局部坐标下建立的，就必须要把它们转换到统一的结构坐标系。将总体坐标轴分别用表示，对某单元有 。式中和分别是局部坐标系和总体坐标系下的单元结点位移向量； 为坐标转换阵，仅与两个坐标系的夹角有关，于是有：（2.55）为

22、该单元在总体坐标系下的单元刚度矩阵。2.5结点平衡方程与求解2.5.1平衡方程的建立用结构力学方法建立结点平衡方程。连续介质用有限元法离散以后，取出其中任意一个结点i,从环绕i点各单元移置而来的结点荷载为：（2.56）式中表示对环绕结点i的所以单元求和，环绕结点i的各单元施加于结点的结点力为：因此，结点的平衡方程可表示为：（2.57）以代入平衡方程，得到以结点位移表示的结点的平衡方程，对于每个结点，都可列出平衡方程，于是得到整个结构的平衡方程组如下：（2.58）式中，为整体刚度矩阵，为全部结点位移组成的向量，为全部结点荷载组成的向量。如果各点的荷载向量也是在单元局部坐标下建立的，在合成以前，也

23、应该把它们转换到统一的结构坐标系下，即：（2.59）式中, 是总体坐标系的结点荷载向量, 为坐标转换阵。2.5.2位移边界条件有限元法对结构进行整体分析时，建立了整体刚度矩阵，也得到了结构的刚度平衡方程。结果刚度方程的求解相当与总刚求逆的过程。但是，从数学上看，未经处理的总刚是对称，半正定的奇异矩阵，它的行列式值为零，不能立即求逆。从物理意义看，在进行整体分析时，结果是处于自由状态，在结点荷载的作用下，结构可以产生任意的刚体位移。所以，在已知结点菏载的条件下，仍不能通过平衡方程唯一地解出结点位移。为了使问题可解，必须对结构加以足够的位移约束，也就是应用位移边界条件。首先要通过施加适当的约束，消

24、除结构的刚度位移，再根据问题要求设定其他已知位移。2.5.3总刚度平衡方程的求解应用有限元法，最终都是归结为解总体刚度平衡方程，它实际上是以总体刚度矩阵为系数矩阵的大型线形代数方程组。通过对结构施加位移边界条件，消除了结构的刚体位移，从而消除总体刚度矩阵的奇异性，解这个线形代数方程组就可求出结点位移，并同时解除出各杆的杆端内力。3.结构优化理论自60年代以来，运筹学和数学规划一直是非常活跃的领域，它们与工程实际问题、计算机应用相联系，促进了多门学科的相互渗透和发展。结构优化设计的产生与发展是力学、最优比技术和计算机应用等多个领域综合的体现。工程实践证明，在工程结构的设计方向，充分运用最优化技术

25、，不仅促使了结构设计发生深刻的变化，而且可以在减轻结构重量和降低工程造价等方面收到显著的效益。 传统的结构设计程序，首先是凭借经验和判断做出结构的初始方案，包扬总体布置、材料选择、结构尺寸和制造工艺等，然后进行结构分析，最后在力学分析的基础上检验其可行或不可行，必要时则进行一、二次修改。在这样的设计程序中，结构分析只起到一种求其安全可行的校核作用。 随着生产的迅速发展，这样的设计程序已经远远不能适应现代化建设发展的需要，这是由于工程建设对于结构设计的要求越趋复杂，使得设计工作者无成熟的先例可作借鉴，难以作出比较，合理的经验设计，更由于计算技术的迅速发展，在结构设计的领域里，能够综合考虑经济、工

26、艺、材料与使用等各个方面的因素。只要充分应用和发展运筹思想和数学规划的方法，最大限度地利用资源、材料和设备等条件就能达到一定的目标。结构优化设计的任务正在深刻地改变工程设计者的思想，成为结构力学中最活跃的领域之一。3.1一般结构优化理论的要点3.1.1 约束条件为了得到一个可行的设计，所有必须遵守的限制条件称为约束条件。在工程问题中，约束大致上可以分为两类，即界限约束与性状约束。界限约束是限制设计变量变比范围的一种约束。它是反映了设计规范或实践经验所规定的构造要求，例如，上述两杆构架的问题，由于受到空间的限制，构架的高度被限制在以下范围之内，即又如圆管直径D必须大于壁厚t，即 这些都是显式。性

27、状约束是反映了问题的性能和状态等要求的一种约束。上述方程式就是考虑了结构强度、稳定性等要求的性状约束条件。在一般情况下，不同于上述两杆构架，我们不能直接写出表示这类约束条件的显函数形式。由于结构分析普通采用有限单元法，因此，在建立应力约束条件时，为了使优化设计流程简捷、明确，一般都是把有限单元法进行应力分析的内容从整个计算程序中分离出来，如果设计要求限制结构的极限承载能力、需安排为一个子程序，根据需要随时进行调用。3.1.2设计变量一个结构设汁的方案是由若干数量来描述的，其中一部分是根据一些具体要求事先确定的，另一部分是在设计过程中视为可以变比的量，我们称其为设计变量。随着问题的性质不同，设计

28、变量可以是构件的截面参数，如截面尺寸，截面面积，截面惯性矩等，也可以是结构的几何参数和选用材料的物理参数，如结点坐标。粱的跨度与间B距，材料的弹性模量等，以及其他各种在优化设计过程中可以变动的量。设计变量可以是连续变化的，但是，在多数工程实际问题中，它是离散跳跃的，例初型钢的截面参数必须是型材表中能够提供的规格，螺栓与铆钉的数目必须是整数。对于离散的变量，虽然有一定的方法来处理，多数情况还是权宜地视为连续变量，在求得计算结果之后再选取与其最接近的离散值。在选择设计变量时，为了得到尽可能好的设计，我们需要综合地考虑更多的因素，选取更多的设计变量，使变量具有较大的变化范围，而且希望每一设计变量的取

29、值间隔划分得越细越好。但是，增多变量和变量的取值数，势必会加大计算工作量，增长优化分析的过程。因此，如何确定设计变量需要仔细推敲，既要尽可能减少设计变量的数目，缩小它的变化范围，又要从工程实际出发，把设计变量进行分组和增大每一组设计变量的取值间隔，从而简化优化设计的分析过程，缩短计算时间，保证能没有遗漏地找到最好的设计。3.1.3 目标函数目标函数也称为直函数、评价函数。它是在许多可行方案中用来选择最好方案的标准。目标函数是设计变量的函数，通常可以表示为或究竟应采取什么内容作为判别方案优劣的标准，是一个十分复杂而又带有综合性因素的问题。在工程结构的设计中，一般习惯于将结构最轻作为优化的目标，这

30、是因为结构的重量是可以确切定量的，并且是一个重要的经济指标，因此关于结构优化设计的研究均属最轻设计的范畴。当然如果结构的加工制作、安装和维护等费用也能够确切定量，而且经济又是工程结构的主要矛盾时，就应该进行最经济设计。在这种情况下，目标因数可以表述为：其中、称为权系数，随着问题的性质，两者的比例可以取不同的数值，除此以外，如果给定材料的数量与重量，或造价，要求设计的工程结构刚度、承载力最大，或固有频率最高，则应取结构的变位、承载力或固有频率作为目标函数。3.1.4 结构优化设计的计算模型结构优化设计又名结构的综合，这是因为一个结构设计只有在综合考虑力学原理，数学方法与手段，材科与工艺，经济与使

31、用等各个方面因素的情况下，才能得到令人满意的结果。因此，怎样从工程实际问题中选择评价设计方案优劣的目标函数，确定描述设计方案的设计变量，给出反映问题的性能与状态等要求的约束条件，从而形成优化设计的计算模型是一个十分重要的步骤。它不仅直接影响到一项工程设计的技术经济效果，而且与结构优比设计能否在工程中推广应用有着非常密切关系。（一)选取合理的截面型式和刚度配置，将内力与位移的分布相对于材料的配置调整得愈是合理，结构的重量就愈轻。 改变截面的尺寸，使在一定荷载作用下，各横截面的最大纤维应力相同，从而达到充分发挥材料的承载能力，减轻结构重量的目的。(二)合理地组合各种材料，结合结构的受力特性和几何形

32、状，愈是能充分利用不同材料的优良性能，结构重量就愈轻。3.2空间刚架结构的优化模型对于空间杆件结构，为了保证压杆正常工作，不因强度不足或因压杆失稳而破坏，必须满足下面两个条件：（3.1）（3.2）式中容许应力，压杆稳定的临界应力。根据欧拉公式可得到此临界应力为：（3.3）A为圆管的截面面积，J为截面惯性矩，为杆的长度，D为直径。（3.4）（3.5） 于是，这一结构优化设计问题的计算模型就可表达为以下的形式，即求设计变量D，d的值，使得这一问题的目标函数（构造的重量）：（3.6）最小。以下是刚架梁的最大应力，临界应力和约束表达：轴向拉压应力：（3.7）平面弯曲正应力：（3.8）平面弯曲正应力:（

33、3.9）弯曲的最大正应力：（3.10）扭转切应力： （3.11）最大主应力为：（3.12）并且要满足约束条件：（3.13）考虑薄壁圆筒的稳定性，要求：（3.14）有：（3.15）将 和代入得： （3.16）又有，则有：（3.17）（3.18）综合上面两式有：当取最小值时，有最小值 ，结构最优化。此时：（3.19）（3.20）（3.21）（3.22）3.3 约束条件的比较由湖北航天科技 2002年第四期关于刚架梁有限元优化设计中给出的关于杆件的约束条件：（1）杆外径的约束根据使用要求确定杆的最大值，其约束表达式为：（3.23）式中：为第根杆的外径，为杆的最大值。（2）空心杆壁厚约束根据薄壁圆筒稳

34、定性要求，确定壁厚的最小值，其约束表达式为：（3.24）式中：为第根杆的外径，为第根杆的内径，为杆壁厚的最小值。仔细分析其条件，发现它必然满足，不能独立构成约束条件。所以根据学过的材料力学中关于压杆稳定的理论知识，即：（为临界应力值）（3.25）用计算出来的杆内最大应力值与临界应力值构成约束条件，即：（3.26）另一个约束条件是满足材料的强度要求，即： （为容许应力）（3.27）有两方面的因素，没考虑强度条件。一是问题本身对该材料的属性没明确的说明；二是很多情况下，我们首先是在满足了稳定性条件的前提下才满足强度条件。所以选择的是将它的稳定性条件当作它的约束条件。两个约束条件下得出的比较结果如表

35、1.6所示：（单位：mm）表1.6单元号12345678优化直径（前）3530303015151515优化直径（后）2637312427272722通过比较，本人的结果更符合实际些。4. 有限元程序优化与结果分析4.1分析的模型如图1.3所示：图1.3结构模型结点1的位移为零，结点5，6的位移为零，结点2载荷为 ，结点3的结点载荷为， ，。材料的弹性模量，泊松比 。由关系式：得 。4.1.1 模型分析时的数据信息基本数据如表1.1、1.2、1.3、1.4、1.5所示：（各物理量采用标准单位制）表1.1单元个数结点个数虚拟结点个数自由度数单元类型数结点载荷数非结点载荷数86829860表1.2结

36、点uvwXYZ1000123000.4324567892.500.5331011121314150.410041617181920210.58500.53522002324250.840.130626002728290.84-0.130表1.3虚拟结点7891011121314坐标X0.410.410.5850.250.410.410.5850.585坐标Y00000000坐标Z000.530.53000.530.53表1.4单元号12345678端点结点12134433终点结点24345656虚拟端结点12345678虚拟端点结点7891011121314表1.5结点载荷123456位移号6

37、710121315数值889-2.71E+23556-1778270.98180.65第i根杆，其外径为，内径为。则有：， 。4.2优化程序设计流程图如图1.4所示：图1.4结构优化流程图4.3有限元优化结果的表达与分析4.3.1验证有限元程序的正确性图1.5为一空间刚架结构。各杆材料和几何性质相同。E=2.1E5MPa，G=9.0E4MPa,A=0.004,J=2.4E-5,。图1.5空间刚架结构建立数据文件SFSAP.IN如表1.7所示：EXAPLE 3,2,6,2,1,0,20,0,0,0,0,0,0.,0.,0.1,2,3,4,5,6,2.4,0.,0.0,0,0,0,0,0,2.4,

38、0.,-2.40.,2.,0.,2.4,2.,0.1,2,4,1,2,3,5,12.1E8,9E7,0.05,2.6E-5,1.2E-5,3,3.0E-51,3,1,-10.,1.2,2,1,1,-15.,-15表1.7？运行程序后，得出结果如表1.8所示：表1.8结点位移单元号DXDYDZRXRYRZ 1.0000.0000.0000.0000-.0000.0000 2.0000-.50029E-02.0000.22337E-02.0000-.25997E-02 3.0000.0000.0000.0000.0000.0000结果相符，说明该程序正确。4.3.2运行优化程序得出的结果：1.最优

39、直径结果如表1.9所示：表1.9单元号最大内应力临界应力最优杆直径循环次数13.205854e+093.205854e+092.599997e-0216022.189872e+082.180232e+083.680006e-0226831.358561e+081.341714e+083.069996e-0220749.128415e+079.053365e+072.369998e-0213751.065981e+081.045165e+082.749997e-0217569.953480e+079.924442e+072.679997e-0216871.784100e+081.760202e+

40、082.659997e-0216681.212181e+081.191198e+082.189998e-021192.结点位移（程序第一次运行得到的结果）如表1.10所示：表1.10结点位移单元号DXDYDZRXRYRZ1.0000.0000.0000-1.067.5794E-01-.8292E-012.6553E-02.1169-.1633E-02-1.364.5753E-01-.25313.3130E-02-.7529E-01.2840E-02.3510-.3685E-02.35274.6572E-02-.4475E-01.1682E-02-.1424-.5025E-02-.6229E-0

41、15-.1963E-01.0000.0000.1988.1028.7958E-016.2588E-01.0000.0000.2018-.7985E-01.7818E-013.结点位移（程序未得到稳定结果时运行得到的某一中间结果）如表1.11所示：表1.11结点位移号DXDYDZRXRYRZ1.0000.0000.0000-.1789E-01.8662E-03-.6686E-032.1068E-03.2043E-02-.2593E-04-.2424E-01.9291E-03-.3802E-023.6493E-04-.8260E-03.4339E-04.8370E-02.2142E-03.6353

42、E-024.1092E-03-.4816E-03.2408E-04-.2843E-02-.1870E-03-.1876E-025-.1850E-03.0000.0000.4037E-02.2027E-02-.3322E-036.3129E-03.0000.0000.3626E-02-.1395E-02-.2290E-034.结点位移（程序时运行得到的稳定结果）如表1.12所示：表1.12结点位移号DXDYDZRXRYRZ1.0000.0000.0000-.9154E-02.4858E-03-.2616E-032.6138E-04.1129E-02-.1457E-04-.1405E-01.395

43、4E-03-.1865E-023.4175E-04.9928E-04.2453E-04.7991E-02.1815E-03.3522E-024.6278E-04.6910E-04.1285E-04-.4236E-02-.1820E-03-.1335E-025.7147E-04.0000.0000.3547E-02.1582E-02-.1803E-026.1003E-04.0000.0000.3042E-02-.9628E-03-.1443E-024.4应用ANSYS验证优化结果4.4.1 ANSYS分析前数据如表1.13所示表1.13单元号外径杆长截面积转动惯量10.02601.030776E

44、-17.645374E-41.324370E-820.03685.590000E-11.531608E-35.315048E-830.03075.941380E-11.065930E-32.574361E-840.02375.578315E-16.352559E-49.143434E-950.02756.027735E-18.552979E-41.657474E-860.02686.027735E-18.123096E-41.495048E-870.02664.492215E-18.002308E-41.450916E-880.02194.492215E-15.424257E-46.66641

45、5E-9材料的弹性模量，泊松比 。选择的是PIPE16号单元，建立模型，施加载荷。4.4.2ANSYS分析结果节点线位移如表1.14所示：表1.14结点线位移节点号UXUYUZUSUM10.0.0.0.20.10279E-03-0.14845E-02-0.23956E-040.14882E-0230.10911E-030.47398E-030.17306E-040.48669E-0340.86119E-040.94586E-030.34168E-040.95039E-0350.10411E-030.0.0.36411E-036-0.19513E-040.0.0.19513E-03节点角位移如表

46、1.15所示：表1.15结点角位移节点号UXUYUZUSUM10.12438E-010.86074E-03-0.25358E-030.12723E-0220.17061E-010.87655E-030.10053E-020.17113E-0330.24652E-03-0.11076E-030.18377E-020.18575E-0240.56647E-02-0.87566E-04-0.32962E-030.56749E-0250.40397E-020.68900E-03-0.29219E-020.30291E-0260.45972E-02-0.41028E-03-0.29416E-020.30

47、055E-02通过比较分析，有限元解和ANSYS解基本相符合，说明该优化程序是正确的。线位移图如图1.5所示：图1.5线位移图由上图可以看出来，在2，3，4号节点处的位移变化比其它节点要大。得到的角位移图如图1.6所示：图1.6角位移图由图可知，在1，2节点的角位移变化最显著。得到的结构变形图如图1.7所示：图1.7结构变形图5.总结与展望5.1 设计变量的选择与优化设计变量是杆的内直径和外直径。这里有两个设计变量。一般在做优化问题，不论是手算还是编程，都是试图通过一系列的约束关系，将多个设计变量转化为一个设计变量。在这里，针对薄壁杆件，需要，又根据薄壁圆筒稳定性要求，最优条件下的与的关系是。

48、在求每根杆件的最大应力时，必须知道每根杆的杆端位移，而要知道每根杆的杆端位移，必须在输入文件中知道每根杆的材料常数，即：截面积，转动惯量,根据材料力学的知识，就是需要知道每根杆的直径。然而这个又是我们所要求的，所以我给每根杆预先给一个直径和一个步长，根据所给的约束条件进行线性搜索，只到找到最优解。5.2 编程时遇见的问题与解决（1）程序的嵌套。我将自己编的结构优化程序嵌套到原始程序，并且对原程序中的个别子程序做了修改，例如对子程序，将变量直径嵌入到了子程序中，发现编译和连接出现了问题，主要原因是变量的作用域模糊；对与数字输入，尤其是对在分母处的数字，没有用实数而用的是整数，导致溢出。在算杆长时

49、，用定义杆长，算出的小于1，根据fortron 默认为整形，结果导致运算时程序的崩溃。（2）输入文件SFSAP.IN的问题。要运行程序，必须先将已知的数据导到程序里去，所以要单独建立一个IN文件。导入数据时，发觉格式不匹配，程序无法解读导入的数据。主要原因是对fortron的输入与输出格式不了解。为了适应嵌套的程序，还需要对输入程序进行修改，得到自己所要的结果。（3）不能运行的问题。数据程序已经加载，可是程序不能运行。在这里，主要的问题不是程序怎么错了，而是程序在哪出错。也就是如何去调试程序。因为断点功能不能使用，通过找指导老师和其它同学，用语句去逐句得寻找，只到找到问题所在。（4）可以运行，

50、但是不能循环的问题。 程序可以运行了，却不能循环。主要原因是在语句里用了阻止了程序的循环。（5）可以循环，但是不能循环效率低下。 虽然可以循环了，但是每次循环都是只对一根杆件循环，其它的没有变化。分析了一下，这样不仅循环效率低下，产生的误差也很大。原因是每次循环回来都会对结构的内力重新计算，再重新判断约束条件。而最优化的情况是结构中的每根杆同时达到最优解。只对一根杆件循环，其它的不循环，必然导致很大的误差。（6）同时循环了，但结果出乎意料。通过改进程序，让每根杆件同时循环。发觉在直径已经不变的情况下，杆的内力却在不断的减小。调试了好久不知其解。最后选择了书上的一个小程序进行类似的分析，才知道究

51、竟。原因是循环回来时未对结构的刚度矩阵和载荷矩阵从新初始化，导致了结构的刚度矩阵逐渐增大。5.3结构优化的展望结构优化设计是随着最优化方法的发展而发展起来的。有限元方法的出现，使得复杂结构的分析有了可能，以及数学规划的渗透和被引入，大多数结构优比设计的问题可以抽象非线性规划的一个命题，从而使结构优化设计，发展非常迅速，理论研究与工程应用所涉及的范围也愈来愈广。目前，国内外采用的结构优化设计方法种类繁多。最为常用的优化方法，一般来说，可以分为两大类：一是准则方法二是数学规划方法。数学规划方法与准则方法两者互相渗透和互相借鉴，促使两条原来不同途径的汇合，出现了准则法与规划法联合应用的方法，把力学概

52、念和优化方法密切地结合在一起，显著地提高了问题的解算效率，这是当前结构优化设计发展中的一个趋势。附录一：空间刚架梁结构有限元优化程序C $debug PROGRAM SFSAPC ANALYSIS PROGRAM FOR SPACE FRAMEREAL K(100,100),KE(12,12),AKE(12,12),X(50),Y(50),Z(50),AL(50), &EAI(6,20),PJ(20),PF(2,20),R(12,12),P(50),D_D(20),E,G,D_E, !D_D为杆直径 DMAX,FE(12),D(50),ADE(12),DE(12),RT(12,12),F(6),AFE(12),AF(12)INTEGER JE(3,50),JN(6,50),JPJ(20),JPF(3,20),M(12),JEAI(40)OPEN (6,FILE=SFSAP.IN)OPEN (8,FILE=SFSAP.OUT)CALL READ(NJ,NJJ,N,NE,NM,NPJ,NPF,JN