实际工程问题中, 结构优化设计的约束条件中通常要求既满足强度限制又满足刚度限制。Rizzi基于准则法进行了多约束结构优化的研究[1]。Schmit和Miura通过设计变量连接, 将设计变量分组, 减少独立设计变量的个数, 并在每次进行结构分析之后, 删除无效约束等近似方法获得了高效的优化方法[2]。Fleury将对偶理论用到结构重量优化问题上, 利用变量可分离的对偶规划进行求解, 取得了与最优准则法相近的计算结果[3]。Fleury和Schmit利用虚载荷方法将某些临界应力约束作为有效应力约束, 其它应力约束转化为上下限约束, 提出了混合最优性准则方法[4]。

大连理工大学课题组开发出“多单元、多工况、多约束的结构优化设计系统——DDDU”[5,6]。钱令希、钟万勰等人将序列二次规划 (SQP) 运用到工程结构优化设计中, 为解决多工况、多约束问题提供了有效途径[7,8,9]。用无量纲设计变量, 通过变量连接, 实现了对多种单元组合结构的位移约束问题的结构优化理论推导[10];文献[5]对应力约束进行零阶近似, 即将其化为动态尺寸下限处理。隋允康、阳志光等[11]利用两点积累信息和两点有理逼近对对偶优化方法进行了改进, 提出了原、倒变量展开的对偶优化方法。叶元烈、秦东晨[12]等对薄壁覆盖件结构优化设计做了研究, 选择了序列二次规划法进行优化, 并编制了OPTSHEET程序。

本文将尺寸、应力约束转化为动态尺寸约束, 然后利用虚载荷法[5,6]建立了设计变量与位移约束的显式关系, 构造了问题的数学模型。在建立这个显式关系的过程中, 摒弃了过去的应力与应变的Mohr积分形式, 而是用单元节点实力向量与节点虚位移向量的乘积和形式代替。用对偶数学规划精确映射原问题, 用泰勒展式建立对偶问题的二阶近似, 最后用二次规划求解器求解。程序中用到了射线步技术、粗选有效约束技术、主被动变量循环技术[10]。

同时满足尺寸、应力和位移约束的膜结构重量最轻问题可以表示为:

其中Ak (i=1, …, K) 是无量纲设计变量, 与膜单元的初始厚度t0i和设计厚度ti的关系为ti=Akt0i (i∈k) , w0k是第k号设计变量控制的单元的初始重量, ur是第r号位移约束对应的位移, uˉuˉk是第r号位移约束的许用值, A˜kA˜k为动态尺寸下限 (由应力约束转化来的动态尺寸约束和尺寸下限的最大值) , AˉˉˉAˉk为尺寸上限。

在使用有限元法的情况下考虑单位载荷法, 对于当前设计变量A0k (k=1, …, K) 来说, 第k号变量区域对第r号位移约束点的位移的贡献可写成

τrk=Σi∈kA0kFrTiuriτrk=Σi∈kAk0FirΤuir (2)

其中Fri为在第r号位移约束对应的实工况下第i号单元的节点力总向量, uri为在第r号位移约束对应的虚工况下第i号单元的节点虚位移总向量, i∈k表示第i号单元由第k号变量控制。

由此得位移ur对设计变量Ak (k=1, …, K) 的近似函数为[10]

ur(A)=Σk=1Kτrk/Akur(A)=Σk=1Κτrk/Ak (3)

经过上面的变量变换及公式推导, 膜结构优化问题的数学模型可写成[10]

用对偶规划精确映射上述问题, 其形式如下[10]

其中ϕ(λ)=minA˜k≤Ak≤Aˉˉˉk[Σk=1KwkAk+Σr=1Rλr(Σk=1Kτrk/Ak−uˉr)] (6)ϕ(λ)=minA˜k≤Ak≤Aˉk[Σk=1ΚwkAk+Σr=1Rλr(Σk=1Κτrk/Ak-uˉr)](6)

根据库恩—塔克条件, 得

式 (7) 为得到函数ϕ (λ) 的条件。

由于难以求得ϕ (λ) 的显式, 所以只能通过序列近似映射法求解[10]。为此需要求出∂ϕ/∂λr与∂2ϕ/ (∂λr∂λp) , 由 (7) 式及对偶理论有

∂ϕ/∂λr=Σk=1Kτrk/A∗k−uˉr∂ϕ/∂λr=Σk=1Κτrk/Ak*-uˉr (8)

∂2ϕ/(∂λr∂λp)=−Σk∈Iaτrkτpk/[2wk(A∗k)3] Ia={k|A˜k<Ak<Aˉˉˉk,k=1,⋯,K}∂2ϕ/(∂λr∂λp)=-Σk∈Ιaτrkτpk/[2wk(Ak*)3]Ιa={k|A˜k<Ak<Aˉk,k=1,⋯,Κ} (9)

其中Ia为主动变量集。

用∂ϕ/∂λr与∂2ϕ (∂λr∂λp) 构造 (5) 式的二阶近似式[10]

其中Trp=Σk∈Iaτrkτpk/[2wk(A∗k)3],br=uˉr−Σk=1Kτrk/A∗k−Σk∈Iaτrk/(2A∗k)Τrp=Σk∈Ιaτrkτpk/[2wk(Ak*)3],br=uˉr-Σk=1Κτrk/Ak*-Σk∈Ιaτrk/(2Ak*)

这样就可以用二次规划求解器求解得λ, 再根据 (7) 式求得A*k。

将上述解法付诸程序实现时, 考虑到求解效率, 宜采用射线步技术、粗选有效约束技术, 考虑到求解精度及主、被动变量循环的问题, 需要迭代求解。

由膜单元组成的结构其刚度阵是设计变量的线性函数, 所以在设计空间原点引出的射线上的内力保持不变, 利用这一特性可以对结构进行性态调整。根据结构响应 (包括应力和位移) 及其许用值求出射线步系数γ, 然后将所有设计变量皆乘以γ得到新的设计点, 在新的设计点内力保持不变, 无需结构重分析就可以把设计点从可行域内部拉到边界上, 使重量下降;也可以把设计点从可行域外拉回到边界上, 使优化迭代过程更加平稳。

采用粗选有效约束技术, 有利于提高优化算法的效率。优化模型中的约束越多, 模型的求解难度就越大, 优化的效率就越低。因此, 从模型的求解考虑, 需要根据约束的性质, 将可能无效的约束去掉, 只保留可能有效的约束, 这种技术就是粗选有效约束技术。

根据定义, 满足A˜k≤Ak≤AˉˉˉkA˜k≤Ak≤Aˉk时设计变量为主动变量, 其它的为被动变量。由二次规划求解器得到λ后, 根据式 (7) 求得的Ak有可能小于下限或大于上限, 从而成为被动变量。对于被动变量有关系式∂Ak/∂λρ=0, 在开始计算二次规划系数矩阵时, 把所有的设计变量都当成了主动变量, 这样计算结果必定存在误差。因此, 需要通过迭代求解来逼近真实解, 直到主动和被动变量集稳定为止。在这个循环中并没有改变膜单元的厚度, 不涉及到内力重新分配问题, 不需进行结构重分析。

平面膜结构的基本尺寸为120×20mm, 划分为24×4的四边形单元, F=40N作用于右边界上, 为减少应力集中, 将载荷分散在右边界的5个节点上, 左边界固支, 如图1所示。材料:弹性模量200GPa, 泊松比0.3, 许用应力120MPa, 密度7800kg/m3。取两种不同的设计变量分组方案进行计算。

方案一:从左至右每20mm作为一个设计变量, 编号分别为t1～t6, 设计精度为0.01。位移约束为最大位移 (结构右端的挠度) 不超过0.2mm, 文本程序、Nastran软件的计算结果及理论解的比较如表1所示。其中理论解是对位移与厚度关系函数利用拉格朗日乘子法求多元函数的极值得到的。很显然本程序无论从结构的总重量还是约束的满足上都比Nastran的稍好, 且优化效率要高很多。程序解与理论解相比稍大, 这是因为理论算法中忽略了剪力对变位的影响。而剪切内力在所有截面都是一样的, 所以剪力对变位的贡献在各截面是常数, 而弯矩对变位的贡献是离固支端越近越大。其次还与网格精度及优化精度也有关系。在第6号变量区域, 因为其弯矩最小, 所以忽略剪力对该区域造成了较大的误差。

方案二:设计变量与区域的对应关系如表2所示, 设计精度为0.025。

单工况、单位移约束的情况:载荷工况和位移约束同方案一。

多工况、多位移约束的情况:工况1同方案一, 工况2为左边固支, 右边界施加240N的水平拉力;位移约束为右端的挠度不超过0.15mm、水平位移不超过0.007mm。

对上述两种情况分别用本文程序与Nastran程序进行计算, 其结果比较如表3所示。第一种情况下工况和约束虽与方案一完全相同, 但由于采用了不同的变量区域划分, 其结构重量只有0.00862kg, 与方案一的结构重量0.012821kg相比结构总重量减轻非常明显。两种情况下本文程序与Nastran软件的结果基本一样, 但本文程序的迭代次数比Nastran软件少很多。

双曲型屋顶, 如图2所示, 四边的双曲线为x2/16- (y2+z2) /9=1, 其结构在y、z面投影为一正方形, 两者取值范围 (-4.31, 4.31) 。材料:弹性模量200GPa, 泊松比0.3, 密度为7800kg/m3, 许用应力100MPa。结构四边铰支, 屋顶加10MPa的压力。整个结构按图3分成三个设计变量区域。初值为0.05m, 优化精度0.025。位移约束为结构中心点的最大位移不超过0.005m。

本文程序与Nastran软件计算结果的比较如表4, 本文程序得到的结构总重量比Nastran软件的大3.3%, 但迭代次数比Nastran软件少。该题在满应力条件下的总重量为38582kg, 比位移约束下的计算结果小。

直径为10m的球面如图4所示, 初始厚度为0.1m, 受内压10MPa, 忽略重力, 用来模仿天然气储罐。材料:弹性模量200GPa, 泊松比0.3, 密度7800kg/m3, 许用应力250MPa。两种载荷工况:工况1为储满气后静止放在支架上, 支架均匀分布在球赤道上的8点固支约束;工况2为上顶点向下作用一400N的力, 位移边界同工况1。位移约束点为球顶点的最大位移为0.015m, 优化精度为0.01。球罐网格及变量分配如图2所示。

用本文程序经过7次迭代收敛, 得到了较好的优化结果, 而Nastran经过26次迭代后得到的结果不能满足应力约束要求, 其结构尺寸全都成了其尺寸下限。说明Nastran优化模块对球结构的应力约束问题不能有效解决。该题在满应力条件下的重量为248363kg, 显然比位移约束下的259509kg小。

本文对膜结构在尺寸、应力、位移约束下的优化问题进行了理论推导。首先将结构的设计参数 (厚度) 转化为无量纲设计变量, 将尺寸、应力约束转化为动态尺寸下限, 然后利用虚载荷法将位移约束转化为含设计变量的约束, 构造对等的数学问题。用对偶规划精确映射, 再用泰勒展式建立对偶规划的二阶近似。最后用二次规划求解器求解二阶近似问题。在求解的过程中, 为了进一步提高优化解法的效率又采用了射线步调整结构性态, 用准有效约束初选筛选约束;为保证程序收敛稳定, 使用了主被动变量循环。最终用MSC/Patran的PCL语言开发了与其连接的优化设计模块。通过对膜结构的单变位、多变位的结构优化问题进行计算, 并与Nastran优化模块进行比较, 考核了程序对膜结构优化设计的可靠性、高效性和稳定性, 同时通过这些结果证明了该理论算法的优越性。