实际工程结构中的设计变量多是在离散集中取值的, 由于设计变量的离散性, 使得传统的连续变量优化设计的理论及种种算法已不再适用。人们对离散变量优化问题提出了许多方法, 大致可分为三类[1]:1) 连续优化解的圆整法, 其优点是可以利用比较成熟的连续变量优化方法进行离散变量优化, 缺点是其解不是离散变量的最优解, 有时相差很远。2) 组合算法, (0, 1) 规划问题就属于这类算法。该算法大大缩减了原问题的许用集合, 减少了迭代次数, 孙焕纯、柴山等[2,3,4,5]提出了寻求离散变量结构优化设计0-1规划可行集的差商向量法, 并且发展了离散变量结构优化设计的组合算法。该算法的缺点是不能用连续变量优化的方法, 并且计算量比较大。3) 附加离散性条件的连续变量优化设计方法, 1983年, Templeman等[6]用这种方法对桁架结构离散变量优化进行了很好的尝试, 将每根杆多节单元的各段面积映射为多节单元的长度, 构造成以分节杆杆长为设计变量的连续变量线性规划问题, 求得各单元长度的最优值后再反过来确定整根杆件的截面。该法的优点是把离散变量问题映射为连续变量问题, 但是他没有解决变量连接的问题, 而且该方法只适合于桁架等内力单元;1987年, 文[7,8,9]在Templeman的基础上首次发展了该理论, 采用长度的无量纲化设计变量解决了变量合将该法推广到变内力单元的离散变量优连接问题, 并采用两节无限小单元的无穷组化中。此外, Stocki, Kolanek等[10]进行了基于可靠性离散优化技巧的研究, 尚凌云等[11]以杆件横截面积及工程中常用的两种球节点体积为设计变量, 采用基于离散变量的两级优化方法对空间网壳结构进行了截面优化设计, 张安宁等[12]提出了将连续变量视为离散变量、并进一步将混合变量转换为单一的纯整型变量, 在整型空间内对全离散或混合离散变量、约束非线性问题进行优化设计的新方法, 邓华[13]用相对差商法对空间网格结构杆件截面尺寸的离散变量优化设计问题进行了研究。本文叙述了离散变量通过力学映射的优化方法, 将文[7,8,9]的方法对膜结构进行了离散变量截面优化, 并在大型有限元分析软件MSC.Nastran上二次开发了相应的优化模块, 数值算例验证了该方法的正确性和二次开发的可行性、有效性。

因为膜单元每点的内力都是变化的, 所以借助于建立微分方程的途径, 在长度微元dx上, 视内力为常数, 在dx上构造分节单元, 每一节都对应一个离散厚度。具体做法为, 设计变量为膜的截面厚度, 映射模型将其化为厚度对应的连续表面积变量, 将表面积无量纲化以实现不同表面积单元的变量连接, 考虑到变内力单元的特点, 把膜单元剖分成无穷个无限小单元, 采用构造相同的微元再拼装成原来大小的整体单元, 其微观结构如图1所示。微元面积为dFi, 包含tkddk与tkd+1d+1k两种厚度, 分别对应于 (1-αk) dFi与αkdFi, 这里k为设计变量号, i∈k表示i单元被k号变量控制, 本来无量纲αk应取0或者1, 求解时暂时放松为0≤αk≤1。

在不考虑离散集的条件下先进行连续变量模型求解, 得到连续最优解t*k和动态尺寸下限tˉk‚ktˉk‚k表示第k号独立变量。由于t*是最优解, 考虑了位移、应力、尺寸三种约束, 而tˉktˉk只考虑了应力和尺寸约束, 因此

tˉk≤t∗k (1)tˉk≤tk*(1)

在离散集中按从小到大的顺序选离散厚度:

a) 若tˉktˉk和t*k靠得较近以至于中间没有离散的厚度, 则取刚好大于或等于t*k的离散厚度作为两节模型的下限tkddk, 即

tˉk≤t∗k≤tkd (2)tˉk≤tk*≤tdk(2)

b) 若tˉktˉk和t*k中间有离散的厚度, 就取紧邻t*k的离散厚度作为tkddk, 即

tˉk≤tkd≤t∗k (3)tˉk≤tdk≤tk*(3)

两节模型的上限tkd+1d+1k为离散集中刚大于tkddk的厚度。这样两节模型上、下限都满足了应力约束。

数学模型的建立就是要怎样处理约束和目标之间的关系, 是在满足约束的不可违背下的最优。但是往往已知的约束与目标之间是隐式关系, 所以建模的过程就是确定目标与约束的显式过程, 既要保证显式与隐式之间尽可能的逼近, 还要保证显式模型的可求解性。

位移约束在多工况下, 考虑动态尺寸下限形成。设实、虚载荷工况分别对应的单元应力向量与应变向量为σi (x, y, t) 和εviiv (x, y, t) 。由莫尔积分, i单元的微元对广义位移Ur的贡献为

duir=∫tkdσTiεviiv (1-αk) dFidt+∫tkd+1σTiεviivαkdFidt (4)

则i单元对位移Ur的贡献为

uir=∫Fiduir=(Δid+lr−Δidr)αk+Δidr (5)uir=∫Fiduir=(Δid+lr-Δidr)αk+Δidr(5)

其中Δid+1r=∫FidF∫tkd+1σTiεviivd (i∈k) ,

Δidr=∫FidF∫tkdσTiεviivdt (i∈k)

在使用有限元法情况下考虑载荷法, 根据虚功原理, 上式的积分形式可转换成单元节点力向量和虚位移向量的乘积和形式。这样可得到

Δidr=∑i∈kFlTiuji (6)Δidr=∑i∈kFilΤuij(6)

其中, Δidr表示第k个变量区域在第l号实工况下对第j号位移约束点的位移的贡献。Fliil为在第l号实工况下第i个单元的节点力总向量, ujiij为在第j号虚工况下第i个单元的节点虚位移总向量。Fliil和ujiij都可以从MSC.Nastran程序中提取, 所以Δidr得到显式化。在进行连续变量优化分析时已经对各个工况下各种约束进行了初选, 只保留了有效约束, 式中r为有效约束数, 这样就大大提高了计算速度。于是得到

Ur=∑k=1n∑i∈kuir=∑k=1n∑i∈k[(Δid+1r−Δidr)αk+Δidr] (7)将式(7)代入位移约束Ur≤δr‚得∑k=1nBrkαk≤δˉr (8)Ur=∑k=1n∑i∈kuir=∑k=1n∑i∈k[(Δid+1r-Δidr)αk+Δidr](7)将式(7)代入位移约束Ur≤δr‚得∑k=1nBrkαk≤δˉr(8)

其中, Brk=∑i∈k(Δid+1r−Δidr)‚δˉr=δr−∑k=1n∑i∈kΔidr‚δrBrk=∑i∈k(Δid+1r-Δidr)‚δˉr=δr-∑k=1n∑i∈kΔidr‚δr值由用户输入。

数学模型以膜结构的厚度为变量, 以应力、位移为约束, 以总体结构的重量为目标, 其简略表达式为

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪W=∑k=1nCkαk→mins.t.∑k=1nBrkαk≤δˉr (r=1‚⋯‚J∗R) 0≤αk≤1 αk={10 ∀K (k=1‚⋯‚n) (9){W=∑k=1nCkαk→mins.t.∑k=1nBrkαk≤δˉr(r=1‚⋯‚J*R)0≤αk≤1αk={10∀Κ(k=1‚⋯‚n)(9)

其中, Ck=∑i∈kρiFi(tkd+1−tkd)‚ρi‚Fi‚tiCk=∑i∈kρiFi(td+1k-tdk)‚ρi‚Fi‚ti分别为第i号膜单元的容重、面积、厚度, n是设计变量总数, J, K分别为实、虚工况总数。式 (9) 最后一个条件为单节性条件, 求解时暂不考虑它, 则式 (9) 变为线性规划问题

⎧⎩⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪W=∑k=1nCkαk→mins.t.∑k=1nBrkαk≤δˉr 0≤αk≤1 (10){W=∑k=1nCkαk→mins.t.∑k=1nBrkαk≤δˉr0≤αk≤1(10)

可通过线性规划求解器对其求解。

由线性规划求解器得到αk, 分为三种情况:1) αk=0, 则tk=tkddk, 取两节模型下限;2) αk=1, 则tk=tkd+1d+1k, 取两节模型上限;3) 对于不满足单节性条件的少数单元需要进行圆整, 即0<αk<1时, 需要判断主、被动变量:a) 对于第r号有效位移约束ur有∂ur/∂αk=Brk, 若αk在所有有效约束中对应的Brk都是非负数, 则该两节解为“下被动变量”, 简单地圆整为αk=0;若αk在所有有效约束中对应的Brk都是非正数, 则该两节解为“上被动变量”, 简单地圆整为αk=1, 这样增加了整体的刚度。b) 不属于“上、下被动变量”的两节解为主动变量, 从偏于安全考虑, 若αk≤0.3, 则取αk=0, 否则取αk=1。

例1 根据实际锅炉尺寸、形状, 建立如图2所示的模型, 直径D=1600mm, 壁厚t=30mm。已知锅炉材料密度为7800kg/m3, 内压P=3.4MPa, 材料弹性模量200GPa, 泊松比0.3, 将其自重当做均布载荷处理, 按简支外伸梁计算, 取许用应力250MPa。尺寸约束上限为0.3mm, 下限为0.003mm;位移约束中间一点, 约束值1.5mm, 方向向下, 优化精度为0.01。有限元模型与设计变量分布如图3, 模型共划分为480个单元, 圆柱部分400个, 两端各40个, 各设计变量的离散集如表1所示。

在软件MSC.Nastran上已二次开发了连续变量[10]和离散变量[11]的优化模块。先对模型进行连续变量优化设计, 得到连续变量最优解和动态尺寸下限, 根据它们在离散集中确定两节模型上、下限, 进而进行优化分析。优化结果如表2所示。本文以后用统一的符号表示表格中的元素:t*k表示连续最优解, tˉktˉk表示动态尺寸下限, tkd+1d+1k表示两节模型的上限, tkddk表示两节模型的下限, 离散解用t**kk**来表示。

表2给出了离散变量的优化结果, 由表2可知, 因为在2号设计变量处加有约束, 所以此区域比较厚;4号设计变量需要圆整, 因为其对应的Brk为负, 所以离散解取两节单元上限, 这样增加了结构的刚度;从安全考虑, 工程中只能选取邻近连续解较大的尺寸, 即取圆整解, 所以离散优化后结构重量重于连续优化后的重量, 但是这样有可能造成不必要的浪费。因为有时不需要在连续解向上取值, 只需向下圆整就可以了, 如本题3、5、6号设计变量, 所以本方法给工程中在离散集中取值提供了理论依据并节省了材料。由于厚度增加使其轴向位移减少, 即离散解满足约束情况好于连续解, 可见圆整解是较为浪费的, 如表3所示。

例2 基本尺寸为120×20 (mm2) 的平面, 材料弹性模量200Gpa, 泊松比为0.3, 许用应力120MPa, 密度7800kg/m3, 厚度为5mm。划分为24×4 个矩形单元, 集中载荷F=200N作用于右边界上, 为减少应力集中, 将载荷平均分散在右边界的5个节点上, 左边界全部采用固支约束, 如图4所示。结构划分为12个设计变量, 变量区域与x的关系如表4所示。尺寸约束的上限取10mm, 下限取0.025mm。位移约束为:结构的右上角点约束值小于0.2mm, 方向向下;结构的右下角点约束值小于0.025mm, 方向水平向右。优化精度为0.01。各个设计变量的离散集如表4所示。优化结果如表5所示。

结果分析:优化结果如表5所示, 2号和7号设计变量需要圆整, 它们对应的Brk都全为负值, 所以基于增加整体刚度考虑, 这两个设计变量都取两节模型上限;从安全考虑, 工程中只能选取圆整解, 所以离散优化后结构重量重于连续优化后的重量, 为了减少不必要的浪费, 有时不需要在连续解向上取值, 只需向下圆整就可以了, 如表中的4、6、8、10、11号设计变量。由于离散优化后整体厚度增加, 所以离散解满足约束情况好于连续解, 如表6所示。

把离散优化后的结构画成立体图, 如图5所示, 设计变量1、3、5、7、9、11比较厚, 而且厚度从受约束区域向外逐渐减少, 这是符合理论解的, 中间设计变量因为受力较小, 受到位移约束的影响比较大, 所以厚度分布没有规律。整体厚度呈“工”字分布, 和工字梁有相似之处, 这在现实生活中有重要的百富策略白菜网。

本文表明力学模型映射离散变量为连续变量, 采用无量纲化技术和构造分节单元的求解途径是有效的, 为工程实际在离散集中取值提供了优化设计手段。同时利用此方法在大型有限元分析软件MSC.Nastran进行了二次开发, 通过大量算例验证了理论的正确性和二次开发的可行性。