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

空间刚架梁有限元程序及优化设计Simwe 会员 liqifu摘要：随着计算机技术的普及和计算速度的不断提高，经历了三十多年的发展历史，理论和算法都已 经日趋完善。有限元分析在工程设计和分析中得到了越来越广泛的重视，已经成为解决复杂的工程分 析计算问题的有效途径。有限元的核心思想是结构的离散化，可以解决很多实际工程需要解决而理论 分析又无法解决的复杂问题。空间刚架梁结构是工程中经常遇见的结构，有限元方法是工程实际中经 常运用的方法，有限元方法可以很好地模拟和分析结构，达到工程需要的要求。对结构的模拟和分析， 一个主要的目的就是对结构进行优化，以达到经济效益。优化理论在结构分析中的作用是越来越突出。 运用结构优化理论，对空间刚架梁结构通过有限元分析，在计算机中编制有限元程序，可以达到对结 构进行优化的目的。关键字：空间刚架梁有限元程序设计结构优化1. 引言1.1论文背景及研究意义人们进行力学分析的方法有很多种，但归结起来可分为两类，即解析法和数值法。由 于实际结构物的形状和所受荷载往往比较复杂，除了少数简单的问题之外，按解析法求解 的非常困难，所以数值法已成为不可替代的广泛应用的方法，并得到了不断发展。有限单 元法就是伴随着电子计算机技术的进步而发展起来的一种新兴数值分析方法。它的数学逻 辑严谨，物理概念清晰，易于理解和掌握，应用范围广泛，能够灵活地处理和求解各种复 杂的问题，特别是它采用矩阵形式表达基本公式，便于运用技术编程运算。这些优点赋予 了有限元强大的生命力。随着计算机技术的普及和计算速度的不断提高，经历了三十多年的发展历史，理论 和算法都已经日趋完善。有限元分析在工程设计和分析中得到了越来越广泛的重视，已经 成为解决复杂的工程分析计算问题的有效途径，如今它已成为广大科技工作者的有力工 具，解决了大量科学研究和工程技术问题。空间刚架梁问题是工程实际中经常遇见的问题，可以用有限元方法对空间刚架梁结构 进行模拟和分析。对空间刚架梁结构进行模拟和分析的一个重要的目的就是对结构进行优 化，以产生经济效益。1.2 有限元概论1.2.1 有限元发展与简介有限元法是大型复杂结构或多自由度体系分析的有力工具，近20 年来已广泛地用于： 工程结构、传热、流体运动、电磁等连续介质的力学分析市，并在气象、地球物理、医学 等领域得到应用和发展。电子计算机的出现和发展，使有限元法在许多实际问题的应用变 成现实，并且有广阔的前景。有限元法的基本思路是将结构物看成由有限个划分的单元组成的整体，以单元结点的 位移或结点力作为基本未知量求解。有限元法基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选 择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数 的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分 方程离散求解。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单 元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基 函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的 近似解构成。1.2.2 有限元方法的基本思路和解题步骤1. 建立积分方程，根据变分原理或方程余量与权函数正交化原理，建立与微分方程初 边值问题等价的积分表达式，这是有限元法的出发点。2. 区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相 互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量 比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位 置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。3. 确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插值 条件的插值函数作为单元基函数。有限元方法中的基函数是在单元中选取的，由于各单元 具有规则的几何形状，在选取基函数时可遵循一定的法则。4. 单元分析：将各个单元中的求解函数用单元基函数的线性组合表达式进行逼近；再 将近似函数代入积分方程，并对单元区域进行积分，可获得含有待定系数 （即单元中各节 点的参数值）的代数方程组，称为单元有限元方程。5. 总体合成：在得出单元有限元方程之后，将区域中所有单元有限元方程按一定法则 进行累加，形成总体有限元方程。6. 边界条件的处理：一般边界条件有三种形式，分为本质边界条件 （狄里克雷边界条 件 ）、自然边界条件（黎曼边界条件 ）、混合边界条件（柯西边界条件）。对于自然边界条件， 一般在积分表达式中可自动得到满足。对于本质边界条件和混合边界条件，需按一定法则对总体有限元方程进行修正满足。7解有限元方程：根据边界条件修正的总体有限元方程组，是含所有待定未知量的封 闭方程组，采用适当的数值计算方法求解，可求得各节点的函数值。1.3结构的优化设计结构优化设计又名结构的综合，这是因为一个结构设计只有在综合考虑力学原理，数 学方法与手段，材科与工艺，经济与使用等各个方面因素的情况下，才能得到令人满意的 结果。因此，怎样从工程实际问题中选择评价设计方案优劣的目标函数，确定描述设计方 案的设计变量，给出反映问题的性能与状态等要求的约束条件，从而形成优化设计的计算 模型是一个十分重要的步骤。它不仅直接影响到一项工程设计的技术经济效果，而且与结 构优比设计能否在工程中推广应用有着非常密切关系。（一）选取合理的截面型式和刚度配置，将内力与位移的分布相对于材料的配置调整得 愈是合理，结构的重量就愈轻。改变截面的尺寸，使在一定荷载作用下，各横截面的最 大纤维应力相同，从而达到充分发挥材料的承载能力，减轻结构重量的目的。（二）合理地组合各种材料，结合结构的受力特性和几何形状，愈是能充分利用不同材 料的优良性能，结构重量就愈轻。结构优化设计是随着最优化方法的发展而发展起来的。在电子计算机出现之前，由于 受到计算手段的限制，人们只能利用以松弛变量、拉格朗日乘子、等式约束消元等措施的 微分法与变分法，或者出自于“同步失效”的概念，把问题变成求解一组联立方程式的方 法，寻求单独杆件的最轻设计。“同步失效”是假定最优结构就是在外荷载作用下，能 使所有可能的破坏模式同时实现的结构，从单个杆件的优化设计推而广之。如果把整个结 构看作一个元件，而结构中的每一根杆件的破坏可以看作一种可能的破坏模式，则自1968 年以后进一步发展起来的准则法，实际上就是“同步失效”设计的演变与发展。从60年 代开始，由于有限元方法的出现，使得复杂结构的分析有了可能，以及数学规划的渗透和 被引入，大多数结构优比设计的问题可以抽象非线性规划的一个命题，从而使结构优化设 计，发展非常迅速，理论研究与工程应用所涉及的范围也愈来愈广。目前，国内外采用的 结构优化设计方法种类繁多。最为常用的优化方法，一般来说，可以分为两大类：（一）准则方法准则方法是从结构力学原理出发，建立一些最优准则，从而寻求用解析形式表达的结 构设计的参数，或者通过直观的迭代运算决定结构各单元的截面参数。常用的准则方法有 满应力准则法、位移准则法、能量准则法等。（二）数学规划方法数学规划方法是从解极值问题酌数学原理出发，运用数学规划中各种方法，求得一系 列设计参数的最优解。当前已被引用的几种规划方法有线性规划、非线性规划中的约束最 优化方法与无约束最优化方法、动态规划、几何规划等。数学规划方法与准则方法两者互相渗透和互相借鉴，促使两条原来不同途径的汇合， 出现了准则法与规划法联合应用的方法，把力学概念和优化方法密切地结合在一起，显著 地提高了问题的解算效率，这是当前结构优化设计发展中的一个趋势。程序框图如图1.1所示图1.1程序框图有限元解的精度和收效性取决于求解的微分方程或采用的变分式和选取的单元。有三 种误差来源，（a）由于域的近似性产生的误差，（b）由于解的近似方法产达的误控，（c）由于 数值计算产生的误差，即在计算机中数值的积分和告入误差。通常确定这些误差不是一个 简单的事情。但在一定的条件下，对于一个约定的单元和问题则可以估算。我们需要对有 限元解进行反复的比较，避免误差过大。该空间刚架梁有限元程序设计及优化设计的结果 很大程度上依赖与约束条件，约束条件选择的合理，优化结果也就合理。1.4本文主要内容及方法本文主要是对空间刚架梁结构进行有限元程序的优化设计，即用有限元的方法分析空 间刚架梁（自行车刚架），求出目标函数，以达到节约材料的目的。利用虚位移原理建立空间梁单元的单元刚度矩阵，利用各结点编号，形成总刚度矩阵。 再利用条件形成载荷矩阵，求出各杆端位移。根据临界载荷形成一个约束条件，选取杆件 直径为设计变量，以总重量最轻为目标函数。将其编制成有限元优化程序，在Power Station 平台上运行。在约束条件下，通过循环语句对设计变量进行线性搜索，找到最优化直径， 以求出目标函数。2. 有限元方法对刚架结构的分析21虚位移原理推导单元刚度矩阵2.1.1虚位移原理有限元的核心就是建立单元刚度矩阵，有了单元刚度矩阵，加以适当组合，可以得到 平衡方程组，这以后剩下的就是一些代数运算。虚位移可以是任意无限小的位移，它在结构的边界上必须满足运动学边界条件。如上 图所示的物体，它受到外力F,F 等的作用，记F = F ,F ,F ,T，在这些外力 作用下，物体的应力为：CJ=Q Q Q ,T ,T ,T T(2.1)现假设物体发生了虚位移, 记：x y z xy yz zx在外力作用处各个外力相应方向的虚位移为5*,5*,5*123(2.2)5 *二5 *, 5 *, 5 *,123由虚位移所产生的虚应变为：e* = e*,e*,e*,y * ,y * ,y * t（2.3）x y Z xy yz zx(2.4)在产生虚位移时，外力已作用于物体，而且在虚位移过程中，外力保持不变，因此， 外力在虚位移上所做的虚功是：5V = F 5 * + F 5 * + F 5 * +. = 5 *t F 123在物体的单位体积内，应力在虚应变上的虚应变能是：(2.5)q e * +q e * +q e * +t y * +t y * +t y * = e *tqx x y y z z xy xy yz yz zx zx(2.6)整个物体应变能是：5U = ffl e * t j dxdydz(2.7)虚位移原理表明，如果在虚位移发生之前，物体处于平衡状态，那么在虚位移发生时, 外力所做虚功等于物体的虚应变能，即：5*tF = M e*tqdxdydz2.1.2单元位移对于空间应力问题有: 、U5 =i(2.8)i(2.14)iiWi其中U , V , W 分别是沿x, y,z方向作用于结点i的结点力。用Fe表示单元e全i ii部结点力所组成的向量：Fe = F,F ,F .(2.15)i j m用虚功原理推导结点力的表达式。假想在单元e中发生了虚位移r*，相应的结点虚 位移为5 *e ，则有：(2.16)(2.17)(2.18)(2.19)(2.20)(2.21)(2.22)i inWp 电丫uiiiJUu.simiiiE.izom中国仿真科技论坛电子期刊 NO.9r* = N6 *e 单元内发生的虚应变为： * = B 5 * e结点力所做的虚功等于每个结点力分量与相应的结点位移分量的乘积之和:5 V = u *U + v *V + w *W +.i ii ii i用矩阵表示，即为：5v = (5 * e )T Fe在整个单元内，应力在虚应变 *上的虚应变能是：5U =*t g dxdydz把 * = B 5 *e代入上式，得到:根据虚功原理5U =5V故：(5 *e) t F e = (5 *e) t fffBr g dxdydz上式对于任何虚位移都必须成立，由于虚位移可以是任意的，所以矩阵(5 *e )T也是 任意的，上式两边与它相乘的矩阵应当相等，于是得到：Fe = JJJBtg dxdydz(2.23)将应力表达式代入得：Fe = BJBtDBdxdydz5 e(2.24)令：K e = JJJBt D B dxdydz(2.25)于是：F e = K e 5 e(2.26)上式建立了结点力和结点位移之间的关系。矩阵Ke称为单元刚度矩阵，它的元素表 示当单元e发生一定的结点位移时，所对应的结点力。单元刚度矩阵Ke决定于该单元的 形状，大小，方向和弹性常数，而与单元的位置无关，即不随单元或坐标轴的平移而改变。2.1.5结点荷载为了能用结构力学方法求解连续介质的应力，所有分布荷载都必须代换以等效的结点 荷载。用P表示结点i的等效结点荷载。i对于空间应力问题：2.1.6 用虚位移原理推导各种结点荷载算式。=iZi其中X , Y,Z分别是沿x,y,z方向作用于i点的集中荷载。 i i i用Pe表示单元e全部结点荷载所组成的向量：Pe =P,P,P Tijm1.分布体积力设单元体积内承受的体积力为q、xq yqz当单元中发生虚位移r*时，体积力q所做的功为：ffjr*tqdxdydz = （5 *e）tfffNtqdxdydz这应当等于等效结点荷载所做的功，即：（5*e）TPe = （5*e）tNtqdxdydzq由上式可得体积力q的等效结点荷载如下：Pe = fff N t qdxdydzq上式右边的重积分应当在单元e的整个体积内进行。2分布面力设单元e是靠近边界的单元，在其边界S上作用着分布的面力p = s 它必须等于等效结点荷载所做的功，由此得到 Pe =jNT pds qs 上式右边的面积分是在分布荷载所作用的表面S上进行的。2.2 形成空间刚架梁单元刚度矩阵1.将结构物离散成杆单元如图 2.2 所示：i inW NflHuiiiJUu.simiiiE.izom如图取结点i, j为刚架梁单元，用节点位移表示应力和应变。单元节点列阵：6 = u ,v , w ,0 ,0 ,0 t(2.35)iiii xiyizi6 = u ,v ,w ,0 ,0 ,0 t(2.36)jjjj xjyjzj节点力列阵：F =N,Q,Q,M,M,Mt(2.37)iiyizixiyiziF =N,Q,Q,M,M,Mt(2.38)jjyjzjxjyjzj有几何方程：du d 2 wd 2 v d0_,k = ,zk , yk ,k t =,z，y,t(2.39)x y z x dxdx 2dx 2 dx(2.40)(2.41)有物理方程：M = N,M ,M ,M t = Dk = DB6 ey z x式中的D 为梁单元的弹性矩阵。dNdx00000dNdx0B=0 -d 2 NF000d 2 N-F0 -d 2NP00d 2 Nz3dx0d 2 Nz4dx000000dNdx0000由2.1虚功原理推导出梁单元的刚度矩阵为:0000000d 2 Nd 2 N0d 2N0z 5 dxz 6 dx0dN00dxE00E0000D =00E0000G(2.42)k = ffl BTDBdV将上式的推导代入得:k 11k 21k 12k 222.43)其中EA10012 EIz1300k =11000006 EIz1200000012 EI6EIy0y13120GJ016EI4EI40y12100006 EIz120004 EIz1EA10012EIz13伙 21二00000006EIz12EA00000012EI6EIy0y1312GJ0016EI2EIA0y12100006EI120002EIz-k 22=100000012EI13 0006 EI0012EIy1306EIy_120000006 EI1206EIy120GJ10004EIy01004 EIz112.3 单元刚度矩阵的坐标变换单元刚度矩阵均是在局部坐标系中得到的。由于局部坐标是根据单元的几何形状选取的，故单元刚度矩阵的形式和推导都比较简单但不同的单元一般具有不同的局部坐标系，因而不能进行不同单元刚度矩阵问的混合运算。即得不到整个计算模型的有限元计算格 式。解决这个问题肋唯一途径就是建立整体坐标与局部坐标间的坐标变换，并通过这种变 换来得到单元刚度矩阵在整体坐标系下的显式。若某个单元刚度矩阵本来就是在整体坐标系下得到的，那就根本不存在坐标变换问题了。首先考虑单兀e在端点i的三个杆端力分量，在局部坐标系x, y, z，中，它们是X, YZ ,。现在推算X，Y，Z与X, Y, Z ,之间的关系。设x,轴与x, y,z轴的夹角分别i i iiiii i i为xx， xy， xz，则x轴在xyz坐标系中的方向余弦为I-二cos(x, x)， l_ = cos(x, y)， I- = cos(x, z)将杆端力 X ,ixzxxz在x轴上投影，可得杆端力为:ixy同理可得:X = X l一 + Yl- + Zl-11 xx 1 xy1 xz(2.44)综合以上三式,则有:Y = Xl- + Yl- + Zl-11 yx 1 yy 1 yzZ = X l- + Yl- + Zl-i(2.45)(2.46)i zx1 zyi zzXiY=LZ-i -zz -1X ,iX, Y Z ,的转换关系式。其中两坐标的转换矩阵为：i i i这就是在端点i由整体坐标系中的杆端力XiYiZi(2.47)Y，Z推算局部坐标系中杆端力iizz -1(2.48)(2.49)(2.50)(2.51)l 一xxl- 广zx参照上述方法，同样可以推导出以M表示M，以及以Q表示Q的表达式，其转换矩 阵也是t。综合以上分析，整体坐标系中的单元杆端力分量列阵Fe与局部坐标系中单元杆端力 分量列阵Fe之间的关系，可用下式表达：Fe 二 TeFe同理，可导出整体坐标系与局部坐标系杆端位移之间的转换关系:6 e 二 Te5 e其中：_ t 000_0 t 00Te =00t0000t称为单元坐标转换矩阵；它是个正交矩阵，故有：2.52)TeT 二 TeT则空间单元在整体坐标系中的单元刚度矩阵为：Ke 二 TeTKeTe(2.53)2.4 单元刚度矩阵的组装及整体分析根据全结构的平衡方程可知，总体刚度矩阵是由单元刚度矩阵集合而成的。如果一个 结构的计算模型分成n个单元，那么总体刚度矩阵就可由各个单元的刚度矩阵Ik e组装而 成，即 ：kk 1(2.54)|k 是由每个单元的刚度矩阵k e的每个系数按其脚标编号“对号入座”叠加而成的。 这种叠加要求在同一总体坐标系下进行。如果各单元的刚度矩阵是在单元局部坐标下建立 的，就必须要把它们转换到统一的结构坐标系。将总体坐标轴分别用 X,Y,Z 表示，对某单 元有e二TS e。式中亓和 e分别是局部坐标系和总体坐标系下的单元结点位移 向量；T为坐标转换阵，仅与两个坐标系的夹角有关，于是有：Ke 二TTKeT(2.55)Ke为该单元在总体坐标系下的单元刚度矩阵。2.5 结点平衡方程与求解2.5.1 平衡方程的建立用结构力学方法建立结点平衡方程。连续介质用有限元法离散以后，取出其中任意一 个结点i,从环绕i点各单元移置而来的结点荷载为P二P eii式中工 表示对环绕结点i的所以单元求和，环绕结点i的各单元施加于结点i的结点 力为：工存因此，结点i的平衡方程可表示为：ie2.56)P戈2.57)Ef 卜=pii以k 代入平衡方程，得到以结点位移表示的结点i的平衡方程，对于每个结点，都可 列出平衡方程，于是得到整个结构的平衡方程组如下：(K %=(P(2.58)式中，k为整体刚度矩阵，为全部结点位移组成的向量，仏为全部结点荷载组 成的向量。如果各点的荷载向量也是在单元局部坐标下建立的，在合成以前，也应该把它们转换 到统一的结构坐标系下，即：2.59)式中，p是总体坐标系的结点荷载向量，k为坐标转换阵。2.5.2位移边界条件有限元法对结构进行整体分析时，建立了整体刚度矩阵，也得到了结构的刚度平衡方 程。结果刚度方程的求解相当与总刚求逆的过程。但是，从数学上看，未经处理的总刚是 对称，半正定的奇异矩阵，它的行列式值为零，不能立即求逆。从物理意义看，在进行整 体分析时，结果是处于自由状态，在结点荷载的作用下，结构可以产生任意的刚体位移。 所以，在已知结点荷载的条件下，仍不能通过平衡方程唯一地解出结点位移。为了使问题 可解，必须对结构加以足够的位移约束，也就是应用位移边界条件。首先要通过施加适当 的约束，消除结构的刚度位移，再根据问题要求设定其他已知位移。2.5.3总刚度平衡方程的求解应用有限元法，最终都是归结为解总体刚度平衡方程，它实际上是以总体刚度矩阵为 系数矩阵的大型线形代数方程组。通过对结构施加位移边界条件，消除了结构的刚体位移, 从而消除总体刚度矩阵的奇异性，解这个线形代数方程组就可求出结点位移，并同时解除 出各杆的杆端内力。3. 结构优化理论自60年代以来，运筹学和数学规划一直是非常活跃的领域，它们与工程实际问题、 计算机应用相联系，促进了多门学科的相互渗透和发展。结构优化设计的产生与发展是力 学、最优比技术和计算机应用等多个领域综合的体现。工程实践证明，在工程结构的设计 方向，充分运用最优化技术，不仅促使了结构设计发生深刻的变化，而且可以在减轻结构 重量和降低工程造价等方面收到显著的效益。传统的结构设计程序，首先是凭借经验和判断做出结构的初始方案，包扬总体布置、 材料选择、结构尺寸和制造工艺等，然后进行结构分析，最后在力学分析的基础上检验其 可行或不可行，必要时则进行一、二次修改。在这样的设计程序中，结构分析只起到一种 求其安全可行的校核作用。随着生产的迅速发展，这样的设计程序已经远远不能适应现 代化建设发展的需要，这是由于工程建设对于结构设计的要求越趋复杂，使得设计工作者 无成熟的先例可作借鉴，难以作出比较，合理的经验设计，更由于计算技术的迅速发展， 在结构设计的领域里，能够综合考虑经济、工艺、材料与使用等各个方面的因素。只要充 分应用和发展运筹思想和数学规划的方法，最大限度地利用资源、材料和设备等条件就能 达到一定的目标。结构优化设计的任务正在深刻地改变工程设计者的思想，成为结构力学 中最活跃的领域之一。3.1 一般结构优化理论的要点3.1.1 约束条件为了得到一个可行的设计，所有必须遵守的限制条件称为约束条件。在工程问题中， 约束大致上可以分为两类，即界限约束与性状约束。界限约束是限制设计变量变比范围的一种约束。它是反映了设计规范或实践经验所规 定的构造要求，例如，上述两杆构架的问题，由于受到空间的限制，构架的高度被限制在 以下范围之内，即15 H 25又如圆管直径D必须大于壁厚t,即tD这些都是显式。性状约束是反映了问题的性能和状态等要求的一种约束。上述方程式就是考虑了结构 强度、稳定性等要求的性状约束条件。在一般情况下,不同于上述两杆构架,我们不能直 接写出表示这类约束条件的显函数形式。由于结构分析普通采用有限单元法,因此,在建 立应力约束条件时,为了使优化设计流程简捷、明确,一般都是把有限单元法进行应力分 析的内容从整个计算程序中分离出来,如果设计要求限制结构的极限承载能力、需安排为 一个子程序,根据需要随时进行调用。3.1.2 设计变量一个结构设汁的方案是由若干数量来描述的,其中一部分是根据一些具体要求事先确 定的,另一部分是在设计过程中视为可以变比的量,我们称其为设计变量。随着问题的性 质不同,设计变量可以是构件的截面参数,如截面尺寸,截面面积,截面惯性矩等,也可 以是结构的几何参数和选用材料的物理参数，如结点坐标。粱的跨度与间B距，材料的弹 性模量等,以及其他各种在优化设计过程中可以变动的量。设计变量可以是连续变化的, 但是，在多数工程实际问题中，它是离散跳跃的，例初型钢的截面参数必须是型材表中能 够提供的规格，螺栓与铆钉的数目必须是整数。对于离散的变量，虽然有一定的方法来处 理，多数情况还是权宜地视为连续变量，在求得计算结果之后再选取与其最接近的离散值。在选择设计变量时，为了得到尽可能好的设计，我们需要综合地考虑更多的因素，选 取更多的设计变量，使变量具有较大的变化范围，而且希望每一设计变量的取值间隔划分 得越细越好。但是，增多变量和变量的取值数，势必会加大计算工作量，增长优化分析的 过程。因此，如何确定设计变量需要仔细推敲，既要尽可能减少设计变量的数目，缩小它 的变化范围，又要从工程实际出发，把设计变量进行分组和增大每一组设计变量的取值间 隔，从而简化优化设计的分析过程，缩短计算时间，保证能没有遗漏地找到最好的设计。3.1.3目标函数目标函数也称为直函数、评价函数。它是在许多可行方案中用来选择最好方案的标准。 目标函数是设计变量的函数，通常可以表示为F (x)或F (x , x x )究竟应采取什么内容11n作为判别方案优劣的标准，是一个十分复杂而又带有综合性因素的问题。在工程结构的设 计中，一般习惯于将结构最轻作为优化的目标，这是因为结构的重量是可以确切定量的， 并且是一个重要的经济指标，因此关于结构优化设计的研究均属最轻设计的范畴。当然如 果结构的加工制作、安装和维护等费用也能够确切定量，而且经济又是工程结构的主要矛 盾时，就应该进行最经济设计。在这种情况下，目标因数可以表述为：F(x) =a 重量+ P 价格其中八0称为权系数，随着问题的性质，两者的比例可以取不同的数值，除此以外， 如果给定材料的数量与重量，或造价，要求设计的工程结构刚度、承载力最大，或固有频 率最高，则应取结构的变位、承载力或固有频率作为目标函数。3.1.4结构优化设计的计算模型结构优化设计又名结构的综合，这是因为一个结构设计只有在综合考虑力学原理，数 学方法与手段，材科与工艺，经济与使用等各个方面因素的情况下，才能得到令人满意的 结果。因此，怎样从工程实际问题中选择评价设计方案优劣的目标函数，确定描述设计方 案的设计变量，给出反映问题的性能与状态等要求的约束条件，从而形成优化设计的计算 模型是一个十分重要的步骤。它不仅直接影响到一项工程设计的技术经济效果，而且与结 构优比设计能否在工程中推广应用有着非常密切关系。(一) 选取合理的截面型式和刚度配置，将内力与位移的分布相对于材料的配置调整得 愈是合理，结构的重量就愈轻。改变截面的尺寸，使在一定荷载作用下，各横截面的最 大纤维应力相同，从而达到充分发挥材料的承载能力，减轻结构重量的目的。(二) 合理地组合各种材料，结合结构的受力特性和几何形状，愈是能充分利用不同材 料的优良性能，结构重量就愈轻。3.2空间刚架结构的优化模型(3.1)(3.2)对于空间杆件结构，为了保证压杆正常工作，不因强度不足或因压杆失稳而破坏，必 须满足下面两个条件： QQ e式中 容许应力，压杆稳定的临界应力。e根据欧拉公式可得到此临界应力为：兀2 EIG 二一A( W )2(3.3)A为圆管的截面面积，J为截面惯性矩，l为杆的长度，D为直径。A =兀(D 2 d 2)3.4)兀(D 4 一 d 4)I =643.5)于是，这一结构优化设计问题的计算模型就可表达为以下的形式，即求设计变量 D， d 的值，使得这一问题的目标函数(构造的重量)：f (x)仝 I Aiii3.6)最小。以下是刚架梁的最大应力 轴向拉压应力：临界应力和约束表达：(u u )EQ= j Zx max l3.7)xy 平面弯曲正应力：Qz maxxz 平面弯曲正应力:3D(v v )ED(0+20)E=/ i +zjzi-12l3.8)y max3D (w w ) ED (0+ 20) E=j i + yy-12l3.9)弯曲的最大正应力：w max=、；Q 2+Q 2z maxy max3.10)扭转切应力：Tmax0.5DG(0 0 )xjxil3.11)最大主应力为：Q=J(Q+Q )2 + 3T 2Lx maxw maxmax3.12)并且要满足约束条件：3.13)考虑薄壁圆筒的稳定性，要求：3.14)有：兀2 EIQ Q =eA( W )23.15)将A =兀(D2 d2)和I =代入得:64(3.16)(3.17)(3.18)比兀2 E2 + d 2)4又有d 兀l D2i ii25 i综合上面两式有：当D取最小值D 时，w有最小值w，结构最优化。ii minii min此时：40 plQD i(3.19)i min兀 V 41E4 d D（3.20）i5i min144 p 213c(3.21)w ii1025兀 Ef (x)迓 l A 迓 w(3.22)z zi=1zi =133约束条件的比较由湖北航天科技2002年第四期关于刚架梁有限元优化设计中给出的关于杆件的 约束条件：（1）杆外径的约束根据使用要求确定杆的最大值，其约束表达式为：1 -D 0（3.23）Dmax式中：D为第i根杆的外径，D 为杆的最大值。imax（2）空心杆壁厚约束根据薄壁圆筒稳定性要求，确定壁厚的最小值，其约束表达式为：D d/、一一j 0（3.24）5min式中：D为第i根杆的外径，d为第i根杆的内径，5为杆壁厚的最小值。iimin仔细分析其条件，发现它必然满足，不能独立构成约束条件。所以根据学过的材料力学中关于压杆稳定的理论知识，即：g二上旦（q为临界应力值）（3.25）e A（W）2e用计算出来的杆内最大应力值g与临界应力值G构成约束条件，即：eq q（3.26）e另一个约束条件是满足材料的强度要求，即： g g （ q 为容许应力）（3.27）有两方面的因素，没考虑强度条件。一是问题本身对该材料的属性没明确的说明；二 是很多情况下，我们首先是在满足了稳定性条件的前提下才满足强度条件。所以选择的是 将它的稳定性条件当作它的约束条件。两个约束条件下得出的比较结果如表1.6所示：（单位：mm）表1.6* 一口 单元号12345678优化直径（前）3530303015151515优化直径（后）2637312427272722通过比较，本人的结果更符合实际些。4. 有限元程序优化与结果分析4.1分析的模型如图13所示:图1.3结构模型i inWp Y LUiuiii.siniiiiE.izom中国仿真科技论坛电子期刊 NO9结点1的(u,v, w)位移为零，结点5, 6的(v, w)位移为零，结点2载荷为F二889N,zM =270.98N - m 结点 3 的结点载荷为 F = 3556N，F =-1778N，M 二 270.98N - m，xxzxM = 180.65N -m。材料的弹性模量E = 200GPa，泊松比卩=0.3。z由关系式：得 G = 76.9231GPa。4.1.1模型分析时的数据信息基本数据如表1.1、1.2、1.3、1.4、1.5所示:(各物理量采用标准单位制)表1.1单元个数结点个数虚拟结点个 数自由度数单元类型数结点载荷数非结点载荷 数86829860表1.2结点UVW000XYZxyz1000123000.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.4M 一口 单元号12345678端点结点12 134433终占结占八、八、24345656虚拟端结点12345678虚拟端点结点7891011121314表1.5结点载荷123456位移号 数值68897-2.71E+210355612-177813270.9815180.65第i根杆其外径为Di，内径为di = 5 d9 369 369 369则有：A = 兀D 2, J =兀D 4, I =兀D 4, I =兀D 4。i 25 i i 20000 i yi 40000 izi 40000 i4.2优化程序设计流程图如图 1.4所示：子程序READ|子程序MA眄子程序calm子程序MK 子程序PE 子程序M ULV12子程序MF子程序MKEp-I子程序MR| 子程序MAKE子程序TRANp|子程序MULV|z1凝殲:门L目1r读人数赭扌山Mi1图 1.4 结构优化流程图子程序solven*1i inWp 砲丫朋HmiiJUu.simiJLiE.izoin43有限元优化结果的表达与分析4.3.1验证有限元程序的正确性图1.5为一空间刚架结构。各杆材料和几何性质相同。E=2.1E5MPa, G=9.0E4MPa,A=0.004m2,l = 2.5m,l = 2.0m ,J=2.4E-5m4,I = 1.4E 5m4,1 = 3.2E 一 5m4。12yz图1.5空间刚架结构建立数据文件SFSAP.IN如表1.7所示：EXAPLE326,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,0.,-2.40.,2.,0.,2.4,2.,0.124,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结点位移|Vt 1=1单兀号DXDYDZRXRYRZ1.0000.0000.0000.0000-.0000.00002.0000-.50029E-02.0000.22337E-02.0000-.25997E-023.0000.0000.0000.0000.0000.0000结果相符，说明该程序正确。4.3.2运行优化程序得出的结果:1.最优直径结果如表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+082.659997e-0216681.212181e+081.191198e+082.189998e-021192.结点位移（程序第一次运行得到的结果）如表1.10 所示：表1.10结点位移iVr 1 1=1 单元号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-015-.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.6353E-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.3954E-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-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.