CN101697176B - 多组件结构系统布局优化设计方法 - Google Patents

Abstract

本发明公开了一种多组件结构系统布局优化设计方法，用于多组件结构系统布局优化设计中几何参数的半解析灵敏度求解，通过将多组件结构系统的有限元模型划分为结构网格、组件网格和过渡网格三部分，分别评估不同几何参数扰动对于三部分网格及其刚度矩阵和质量矩阵的影响，从而获得刚度矩阵和质量矩阵的半解析灵敏度，进一步通过推导获得各种设计指标的灵敏度。对比现有布局优化设计所采用的有限差分法，本发明采用的半解析方法有效避免了几何参数扰动更新后的有限元求解过程，节省了灵敏度分析与优化设计时间，缩短了多组件结构系统布局优化设计周期，提高了设计效率。

Description

多组件结构系统布局优化设计方法

技术领域

本发明涉及一种优化设计方法，特别是多组件结构系统布局优化设计方法。 

背景技术

文献“Integrated layout design of the multi-component system，Zhu J.H.，Zhang W.H.，BeckersP.，International Journal for Numerical Methods in Engineering.2009，78(6)：631-651”公开了一种多组件结构系统布局优化设计方法，这种方法结合了结构拓扑优化技术和装填布局优化技术，将描述结构拓扑的伪密度变量和描述组件位置、方向的几何参数同时定义为设计变量，以结构的刚度为设计指标，实现了一定装填区域内组件的位置布局和支撑连接结构形式同时设计。 

由于这种拓扑优化设计方法的设计变量定义在每个单元上，因此当组件移动导致设计区域有限元网格变化时，无法进行设计变量的迭代更新。因此文献公开的密度点技术和嵌入式网格划分技术，将传统定义在单元上的拓扑设计变量转而定义在设计区域内的固定几何点上，另外采用组件网格和结构网格分别划分和嵌入的手段最大限度地减少网格重新划分的工作量。 

文献公开的结构应变能对于组件位置的有限差分灵敏度求解方法，首先利用嵌入式网格划分技术，将组件嵌入到设计区域内合适的位置和方向，进行有限元分析获得结构的整体应变能。然后对上述几何参数进行摄动，获得组件新的位置和方向，重新利用嵌入式网格划分进行有限元建模，并分析获得摄动后的结构应变能，继而利用有限差分法求解灵敏度。这种灵敏度求解方法需要在每次几何参数摄动后进行嵌入式网格划分和应变能求解，涉及到反复的几何布尔运算和部分有限元网格的重新划分，以及有限元分析。当设计问题中包含大量几何参数时，该方法需要耗费大量的计算时间，无法实现复杂多组件结构系统的快速设计。 

例如在对包含3600拓扑设计变量、六个几何设计变量(其中四个平动位置变量，两个转动位置变量)的二维组件进行优化设计时，以上建模与有限元分析的过程需要5分钟以上(计算机配置奔腾IV3.0，2G内存)。而一次优化设计迭代除了对原始模型的求解之外，还需要对每个几何设计变量进行差分，总共需要七次建模与有限元分析，因此一次迭代将耗时35分钟以上，按照优化设计平均需要迭代七十次达到收敛来计算，则需要40小时才能获得优化结果。而且一旦设计问题包含大量的几何参数，求解时间将成倍增加，这无法满足现代结构系统设计的敏捷性需求。 

发明内容

为了克服现有技术多组件结构系统布局优化设计耗时大的不足，本发明提供一种多组件 结构系统布局优化设计方法，通过对有限元节点位置和刚度矩阵的摄动，避免重复的有限元建模和分析过程，可以提高设计效率，减少多组件结构系统布局优化设计过程的耗时。 

本发明解决其技术问题所采用的技术方案：一种多组件结构系统布局优化设计方法，其特点是包括以下步骤： 

(a)将多组件结构有限元模型划分为结构网格、过渡网格和组件网格三部分，对多组件原始结构进行有限元分析，获得过渡网格和组件网格两部分的刚度矩阵M_(T)和K_(C)l、质量矩阵M_(T)和K_(C)l、结构位移矢量u以及目标函数数值； 

(b)利用组件的节点移动完成几何参数的摄动，根据组件的平动或转动由坐标变换公式求解出节点移动值； 

(c)重新组装过渡网格的刚度矩阵、质量矩阵参数信息，在组件转动情况下利用转动变换矩阵求解更新后的组件的刚度矩阵、质量矩阵参数信息； 

(d)在组件平动情况下，仅对过渡网格的刚度矩阵、质量矩阵参数信息进行有限差分计算以近似代替其偏导数；在组件转动情况下，对过渡网格和组件网格的刚度矩阵、质量矩阵参数信息分别进行有限差分计算以近似代替其偏导数；从而获得目标函数对当前几何参数的偏导数，即灵敏度信息； 

(e)对下一个几何参数进行摄动，直到多组件结构系统布局优化设计完成。 

本发明的有益效果是：由于半解析方法，省去了摄动后重新建模和有限元分析的过程，用系统部分结构的刚度矩阵和质量矩阵的差分代替现有的有限差分法求解灵敏度数据，在配置相同的计算机上对相同结构的多组件结构系统布局进行优化设计，所需时间由现有技术的40小时减少到6分钟左右，节省了大量的计算时间，提高了分析和优化设计效率，从而缩短了组件结构系统的设计周期。 

下面结合具体实施方式对本发明作详细说明。 

具体实施方式

实施例1：结构系统整体应变能对几何参数的半解析灵敏度。 

典型的静力学环境下多组件结构系统的布局优化问题定义如下： 

$\mathrm{find}:\mathrm{\η}=({\mathrm{\η}}_1,{\mathrm{\η}}_2,...,{\mathrm{\η}}_n),s=\left[\begin{array}{ccc}{x}_1& y_1& {\mathrm{\θ}}_1\\ x_2& y_2& {\mathrm{\θ}}_2\\ & .& \\ & .& \\ & .& \\ x_m& y_m& {\mathrm{\θ}}_m\end{array}\right]---\left(1\right)$

min： $C=\frac{1}{2}{u}^T\mathrm{Ku}$

其中，η为材料布局的拓扑设计变量矢量，s为与m个组件相关的平动、转动位置几何设计变量矩阵，C为结构应变能，u为结构位移矢量，K为结构刚度。在满足材料用量、组件几何位置不干涉等约束条件的前提下，使整个结构系统获得最小的应变能即最大的刚度。优 化设计的关键在于获取目标函数对于设计变量的灵敏度数值。 

拓扑设计变量的灵敏度求解格式可以通过解析法推导获得。而对于几何设计变量的灵敏度，由于涉及到几何参数以及有限元网格的变化，无法推导出解析的灵敏度格式。因此，为解决这一问题，本发明首先提出将经过嵌入式网格划分的结构系统有限元网格分成结构网格、组件网格和过渡网格三部分。结构的整体刚度矩阵可以表示为这三部分刚度矩阵的叠加： 

$K=K_{\left(S\right)}+K_{\left(T\right)}+\underset{l=1}{\overset{n_c}{\mathrm{\Σ}}}{K}_{\left(C\right)l}---\left(2\right)$

其中，K为整体刚度矩阵，K_(S)为结构刚度矩阵，K_(T)为过渡刚度矩阵，K_(C)l为第l个组件的刚度矩阵。 

本实施例中半解析灵敏度的求解步骤具体为： 

1)有限元分析。 

对原始结构进行有限元分析，获得优化设计的目标函数C、结构位移矢量u，以及组件网格和过渡网格两部分子结构的刚度矩阵K_(C)l、K_(T)； 

2)根据几何变量摄动有限元网格。 

通过直接移动单元节点的方式来体现几何变量的摄动即组件位置的平移或转动，保证了有限元网格拓扑形式和自由度数目保持不变，这也是对刚度矩阵进行摄动和差分的前提条件。 

对于第l个组件的平动Δx或者Δy，仅需要将属于当前组件的所有有限元节点平动Δx或者Δy，为了获得可信的的灵敏度数据，通常摄动量在单元平均尺寸的1/50到1/10之间。 

对于第l个组件相对于自身参考点的转动Δθ，可以建立组件每个节点到参考点的相对坐标，再利用坐标变换公式获得转动后的新的节点坐标，从而完成转动的摄动。 

3)获得摄动后的结构整体刚度矩阵。 

对于第l个组件的平动，整体结构中的组件刚度矩阵K_(C)l以及结构刚度矩阵K_(S)保持不变，仅有过渡刚度矩阵K_(T)由于组件移动的扯动效应而变化，因此仅需要重新组装摄动后的过渡刚度矩阵K^* _(T)。 

对于第l个组件相对于自身参考点的转动，整体结构中的结构刚度矩阵K_(S)保持不变，转动后的组件刚度矩阵可以通过坐标变换公式获得： 

$K_{\left(C\right)l}^{*}=A{\left(\mathrm{\Δ\θ}\right)}^T{K}_{\left(C\right)l}A\left(\mathrm{\Δ\θ}\right)---\left(3\right)$

A为转动矩阵，是转动角Δθ的函数。因此，仍然只需要重新组装摄动后的过渡刚度矩阵K^* _(T)。 

4)利用刚度矩阵差分获得灵敏度。 

目标函数C的灵敏度求解的推导过程如下，首先写出目标函数偏导数与刚度矩阵偏导数 的函数关系： 

$\frac{\∂C}{\∂s_{\mathrm{\ϵ}}}=-\frac{1}{2}{u}^T\·\frac{\∂K}{\∂s_{\mathrm{\ϵ}}}\·u---\left(4\right)$

s_ε是某一个几何变量。下一步写出刚度矩阵的偏导数的差分求解方法： 

$\frac{\∂K}{\∂s_{\mathrm{\ϵ}}}\≈\frac{\mathrm{\ΔK}}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}=\frac{K^{*}-K}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}---\left(5\right)$

K^*是摄动后的整体刚度矩阵，由于结构刚度矩阵K_(S)和组件刚度矩阵K_(C)l不受平动摄动的影响，公式可以简化为： 

$\frac{\∂K}{\∂s_{\mathrm{\ϵ}}}\≈\frac{K^{*}-K}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}=\frac{K_{\left(T\right)}^{*}-K_{\left(T\right)}}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}---\left(6\right)$

转动摄动状况下，仅有结构刚度矩阵K_(S)不受平动摄动的影响，公式简化为： 

$\frac{\∂K}{\∂s_{\mathrm{\ϵ}}}\≈\frac{K^{*}-K}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}=\frac{K_{\left(T\right)}^{*}+K_{\left(C\right)l}^{*}-K_{\left(T\right)}-K_{\left(C\right)l}}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}---\left(7\right)$

5)返回第二步，对下一个几何变量进行摄动并求解直到完成所有灵敏度的求解。 

由此可见，对比完全差分的重新划分网格、组装整体刚度矩阵、有限元分析、差分计算的步骤，半解析灵敏度的求解仅需要摄动部分网格、组装小部分刚度矩阵、差分计算三步即可获得几何变量的灵敏度数据，不仅节省了网格划分、组装刚度矩阵和有限元分析的大量计算时间，也提高了灵敏度求解的精度，尤其对于大型复杂结构、大量几何设计变量的情况效果更加明显。 

实施例2：结构系统一阶频率对几何参数的半解析灵敏度。 

典型的动力学环境下多组件结构系统的布局优化问题定义如下： 

find：η＝(η₁，η₂，...，η_n)， $s=\left[\begin{array}{ccc}{x}_1& y_1& {\mathrm{\θ}}_1\\ x_2& y_2& {\mathrm{\θ}}_2\\ & .& \\ & .& \\ & .& \\ x_m& y_m& {\mathrm{\θ}}_m\end{array}\right]$

max：ω₁ ²                 (8) 

s.t.：g(η，s)≤0 

其中，η为材料布局的拓扑设计变量矢量，s为与m个组件相关的平动、转动位置几何设计变量矩阵，ω₁为结构的一阶圆频率，通过结构的运动方程的特征值分析获得： 

(K-ω²M)u＝0              (9) 

其中K和M为结构的整体刚度矩阵和质量矩阵，ω为结构固有频率，u为结构位移矢量。 仍然将结构系统有限元网格分成结构网格、组件网格和过渡网格三部分。结构的整体刚度矩阵和整体质量矩阵都可以表示为这三部分的叠加： 

$K=K_{\left(S\right)}+K_{\left(T\right)}+\underset{l=1}{\overset{n_c}{\mathrm{\Σ}}}{K}_{\left(C\right)l}---\left(10\right)$

$M=M_{\left(S\right)}+M_{\left(T\right)}+\underset{l=1}{\overset{n_c}{\mathrm{\Σ}}}{M}_{\left(C\right)l}---\left(11\right)$

本实施例中半解析灵敏度的求解步骤具体为： 

1)有限元分析。 

对原始结构进行有限元分析，获得优化设计的目标函数ω₁ ²、系统位移矢量u，以及组件网格和过渡网格两部分子结构的刚度矩阵K_(T)、K_(C)l和质量矩阵M_(T)、M_(C)l； 

2)根据几何变量摄动有限元网格。 

按照实施例1中介绍的方法完成对组件平动和转动的摄动。 

3)获得摄动后的结构整体刚度矩阵和质量矩阵 

对于第l个组件的平动，整体结构中的组件刚度矩阵K_(C)l、质量矩阵M_(C)l以及结构刚度矩阵K_(S)、质量矩阵M_(S)保持不变，仅有过渡刚度矩阵K_(T)、质量矩阵M_(T)由于组件移动的扯动效应而变化，因此仅需要重新组装摄动后的过渡刚度矩阵K^* _(T)和质量矩阵M^* _(T)。 

对于第l个组件相对于自身参考点的转动，整体结构中的结构刚度矩阵K_(S)、质量矩阵M_(S)保持不变，转动后的组件刚度矩阵和质量矩阵可以通过坐标变换公式获得： 

$K_{\left(C\right)l}^{*}=A{\left(\mathrm{\Δ\θ}\right)}^T{K}_{\left(C\right)l}A\left(\mathrm{\Δ\θ}\right)---\left(12\right)$

$M_{\left(C\right)l}^{*}=A{\left(\mathrm{\Δ\θ}\right)}^T{M}_{\left(C\right)l}A\left(\mathrm{\Δ\θ}\right)---\left(13\right)$

A为转动矩阵，是转动角Δθ的函数。因此，仍然只需要重新组装摄动后的过渡刚度矩阵K^* _(T)和质量矩阵M^* _(T)。 

4)利用刚度矩阵和质量矩阵差分获得灵敏度。 

目标函数ω₁ ²的灵敏度求解的推导过程如下，首先写出目标函数偏导数与刚度矩阵偏导数的函数关系： 

$\frac{\∂{\mathrm{\ω}}_1^2}{\∂s_{\mathrm{\ϵ}}}=\frac{u^T\frac{\∂K}{\∂s_{\mathrm{\ϵ}}}u-{\mathrm{\ω}}^2{u}^T\frac{\∂M}{\∂s_{\mathrm{\ϵ}}}u}{u^T\mathrm{Mu}}---\left(14\right)$

s_ε是某一个几何变量。写出刚度矩阵和质量矩阵的偏导数的差分求解方法： 

$\frac{\∂K}{\∂s_{\mathrm{\ϵ}}}\≈\frac{\mathrm{\ΔK}}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}=\frac{K^{*}-K}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}---\left(15\right)$

$\frac{\∂M}{\∂s_{\mathrm{\ϵ}}}\≈\frac{\mathrm{\ΔM}}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}=\frac{M^{*}-M}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}---\left(16\right)$

与实施例1相同，平动摄动的情况下，公式可以简化为： 

$\frac{\∂K}{\∂s_{\mathrm{\ϵ}}}\≈\frac{K^{*}-K}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}=\frac{K_{\left(T\right)}^{*}-K_{\left(T\right)}}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}---\left(17\right)$

$\frac{\∂M}{\∂s_{\mathrm{\ϵ}}}\≈\frac{M^{*}-M}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}=\frac{M_{\left(T\right)}^{*}-M_{\left(T\right)}}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}---\left(18\right)$

转动摄动状况下，公式简化为： 

$\frac{\∂K}{\∂s_{\mathrm{\ϵ}}}\≈\frac{K^{*}-K}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}=\frac{K_{\left(T\right)}^{*}+K_{\left(C\right)l}^{*}-K_{\left(T\right)}-K_{\left(C\right)l}}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}---\left(19\right)$

$\frac{\∂M}{{\∂s}_{\mathrm{\ϵ}}}\≈\frac{M^{*}-M}{{\mathrm{\Δs}}_{\mathrm{\ϵ}}}=\frac{M_{\left(T\right)}^{*}+M_{\left(C\right)l}^{*}-M_{\left(T\right)}-M_{\left(C\right)l}}{\mathrm{\Δ}{s}_{\mathrm{\ϵ}}}---\left(20\right)$

5)返回第二步，对下一个几何变量进行摄动并求解直到完成所有灵敏度的求解。 

同样，由于避免了对结构有限元网格的重新划分、整体刚度矩阵和质量矩阵的组装以及有限元分析过程，半解析灵敏度的分析策略极大地节省了一阶频率最大化问题的求解计算时间，提高了灵敏度分析精度。 

例如在对包含3600拓扑设计变量、六个几何设计变量(其中四个平动位置变量，两个转动位置变量)的二维组件进行优化设计时(计算机配置奔腾IV3.0，2G内存)，每次对平动位置变量的半解析灵敏度求解的计算时间只需要5秒左右，对转动位置变量的求解由于涉及到较复杂的节点位置换算，需要计算时间10秒左右。因此，包含原始模型建模与求解的一次优化迭代总共需要计算时间6分钟左右，是现有有限差分法求解时间的1/5～1/6。而且在增加几何参数数目的情况下，一次迭代时间的延长并不明显。 

Claims (1)

1.一种多组件结构系统布局优化设计方法，其特征在于包括以下步骤：

(a)将多组件结构有限元模型划分为结构网格、过渡网格和组件网格三部分，对多组件原始结构进行有限元分析，获得过渡网格和组件网格两部分的刚度矩阵K_(T)和K_(C)l、质量矩阵M_(T)和M_(C)l、结构位移矢量u以及目标函数数值；

(b)利用组件的节点移动完成几何参数的摄动，根据组件的平动或转动由坐标变换公式求解出节点移动值；

(c)重新组装过渡网格的刚度矩阵、质量矩阵参数信息，在组件转动情况下利用转动变换矩阵求解更新后的组件的刚度矩阵、质量矩阵参数信息；

(d)在组件平动情况下，仅对过渡网格的刚度矩阵、质量矩阵参数信息进行有限差分计算以近似代替其偏导数；在组件转动情况下，对过渡网格和组件网格的刚度矩阵、质量矩阵参数信息分别进行有限差分计算以近似代替其偏导数；从而获得目标函数对当前几何参数的偏导数，即灵敏度信息；

(e)对下一个几何参数进行摄动，直到多组件结构系统布局优化设计完成。