CN103065015B - 一种基于内力路径几何形态的承载结构低碳节材设计方法 - Google Patents

Abstract

本发明公开了一种基于内力路径几何形态的承载结构低碳节材设计方法，该方法将承载结构的材料模量（如弹性模量）作为以低碳节材为导向的设计变量，并利用数值方法求解得到的内力路径几何形态（包括疏密度与弯曲度）对单元弹性模量进行惩罚，通过循环迭代获得承载结构的最佳传力路径，并以此为依据定义结构的基本承载架构与空间要求。该方法将内力路径的相关理论应用到结构拓扑优化设计中，为重大装备承载结构的低碳节材设计提供了一种新的思路。

Description

一种基于内力路径几何形态的承载结构低碳节材设计方法

技术领域

本发明属于结构拓扑优化设计领域，涉及一种基于内力路径几何形态的承载结构低碳节材设计方法。

背景技术

重型、超重型制造装备是融合了多种学科技术的复杂集成系统，随着产品结构的复杂化与服役环境的极限化趋势越来越显著，客户对于产品在工作过程中的“全局高刚性”要求也越来越高。与此同时，在当前“低碳”发展的形势要求下，大型制造装备的“轻量化表现”将逐渐成为衡量“资源与能源耗费”的重要考核指标。在这样的发展趋势下，越来越多的设计人员将优化技术应用到承载结构的低碳节材设计上来。结构优化可分为尺寸优化、形状优化、形貌优化和拓扑优化。其中形状优化和尺寸优化经过数十年的发展，已达到较为成熟的水平。但是这两种方法是对一个已有的拓扑结构进行优化，并不能改变结构的拓扑关系，即使形状和尺寸达到最优但是拓扑结构并不一定最优，导致最终结果并不是最优的。而拓扑优化设计是在结构初始的拓扑关系未知的情况下，寻求材料在空间中的最佳分布从而得到最优的拓扑结构，初始拓扑结构设计的好坏对整个设计过程有着重要的影响。因此，近年来拓扑优化设计成为国内外结构优化领域的研究热点。但是，由于拓扑优化的描述及算法繁琐复杂，一些比较关键的技术还处在探索和发展阶段。

发明内容

本发明的目的在于提供一种基于内力路径几何形态的承载结构低碳节材设计方法。

为达到上述目的，本发明采用了以下技术方案。

1)设定一种弹性模量可变的单元并将单元的弹性模量作为设计变量，预先设定弹性模量的变化范围；

2)建立承载结构的有限元模型，为有限元模型中各单元设定弹性模量初始值，然后施加边界条件并求解；

3)经过步骤2）后，从有限元分析结果中提取单元节点应力值，通过插值得到各单元质心处应力值，根据各单元质心处应力值得出各单元质心处的内力传递矢量，模值最大的内力传递矢量代表内力路径在单元质心处的疏密程度，以此为基础求解内力路径在单元质心处的曲率，曲率值代表了内力路径在单元质心处的弯曲程度；

4)利用步骤3）中得到的各个单元质心处内力路径的几何形态指标量对单元当前的弹性模量值进行惩罚，并获得新的单元弹性模量值，几何形态指标量包括疏密程度与弯曲程度；

5)对步骤4）中得到的各单元弹性模量值进行光滑平均化处理；

6)返回步骤2）将步骤5）中得到的单元弹性模量值赋给有限元模型中相应的各个单元，并进行循环迭代计算，直至单元的平均弹性模量达到预先设定的数值，这样便可得到承载结构的最佳传力路径，也即材料在设计空间中的最佳分布形式。

弹性模量的最小值设定为一个接近0但不等于0的值，以避免有限元分析过程中刚度矩阵的奇异，弹性模量的最大值设定为承载结构选用材料的实际弹性模量值，并将最大值作为初始值赋给各个单元。

所述对单元当前的弹性模量值进行惩罚，包括以下步骤：

1)根据模值最大的内力传递矢量对承载结构中所有单元进行降序排列，获得单元列表，从列表的顶层单元开始，对单元体积求和直至所得体积达到总体积的预期保留比例为止，将此时对应单元的模值最大的内力传递矢量定义为截止矢量；

2)根据式2对承载结构中各个单元的弹性模量进行惩罚更新，直至单元平均弹性模量达到β·E_max时程序结束（式5），β表示材料保留率，E_max表示弹性模量的最大值，P_i表示单元主传递矢量，C_i表示单元质心处的曲率；

$E_i^{\left(k\right)}=E_i^{(k-1)}\·F(P_i^{(k-1)},C_i^{(k-1)})$ $\left(2\right)$

$=E_i^{(k-1)}\·F_1{\left(P_i^{(k-1)}\right)}^{F_2\left(C_i^{(k-1)}\right)}$

$F_1\left(P_i^{(k-1)}\right)=\left(\frac{\left|P_i^{(k-1)}\right|}{\left|P_{\mathrm{cut}}^{(k-1)}\right|}\right)---\left(3\right)$

$F_2\left(C_i^{(k-1)}\right)={\mathrm{\λ}}^{\frac{C_i^{(k-1)}}{C_i^{(k-1)}+1}}---\left(4\right)$

$\underset{i=1}{\overset{N}{\mathrm{\Σ}}}{E}_i^{\left(k\right)}/N=\mathrm{\β}\·E_{\max }---\left(5\right)$

k表示迭代步数，λ表示优化系数，1<λ<10，P_cut表示截止矢量，是惩罚更新后单元的弹性模量计算值，是惩罚更新前单元的弹性模量计算值，当 $E_i^{\left(k\right)}>E_{\max }$ 时，取 $E_i^{\left(k\right)}=E_{\max };$ 当 $E_i^{\left(k\right)}<E_{\min }$ 时，取 $E_i^{\left(k\right)}=E_{\min },$ E_min表示弹性模量的最小值，N表示单元数量。

本发明以内力路径的相关理论为基础，设定一种弹性模量可变的单元并将其弹性模量作为设计变量，并通过各单元质心处的载荷向量以及内力路径曲率对单元弹性模量进行处理，并进行循环迭代求解。最终得到用单元弹性模量来表征的材料在空间的合理的分布结果，为承载结构的低碳节材设计提供了一种新方法。

目前对于内力路径可视化的研究是一个新兴的研究领域，通过对内力路径的认识，可以使设计人员对承载结构内部的载荷分布与传递情况有更好地理解，对结构的承载性能有更加清楚的认识，从而帮助设计人员在更短的时间内提出更加合理高效的结构设计方案。准确提取承载结构中内力路径的几何形态指标量（如疏密程度与弯曲程度），并将其引入到结构拓扑优化算法当中，将是一个有较高工程应用价值的思路与方法。

附图说明

图1为本发明拓扑优化流程图；

图2为实施例中的承载结构几何模型（a）及其有限元模型（b）；

图3为承载结构水平和竖直方向的内力传递矢量场图；a）为水平方向内力传递矢量，b）为竖直方向内力传递矢量；

图4为承载结构水平和竖直方向的内力传递矢量模值云图；（a）为水平方向，（b）为竖直方向；

图5为承载结构水平和竖直方向的内力路径；a）为水平方向内力路径，b）为竖直方向内力路径；

图6为承载结构水平和竖直方向的内力路径曲率值云图；（a）为水平方向，（b）为竖直方向；

图7为拓扑优化结果；（a）β=0.3，（b）β=0.4，（c）β=0.5，（d）β=0.6，（e）β=0.7。

具体实施方式

下面结合附图和实施例对本发明作进一步说明。

本发明首先基于有限元软件对承载结构内部任意点的内力传递矢量进行求解，据此构建承载结构的内力传递路径，进而提取内力传递路径的几何形态指标量（疏密度与弯曲度）并以此为依据对承载结构各单元的弹性模量进行惩罚，经过循环迭代得到承载结构的最佳传力路径，也即材料在设计空间中的最佳分布形式。

参见图1，方法具体包括：

1）在建立有限元模型时，使用一种人为设定的弹性模量可变的材料单元，将单元的弹性模量作为设计变量，并用各单元弹性模量表征材料在空间的分布情况。在计算启动之前，定义相关初始化参数，包括弹性模量上限E_max与下限E_min以及材料保留率β（0＜β<1）。取E_max为零件材料的实际弹性模量值，将模型中所有单元的初始弹性模量设置为E_max；取E_min值接近于0（但又要大到一定程度使得在后续有限元求解过程中不至于将其忽略）；将外载荷作用区域以及约束区域作为模型必须保留的特征，临近上述特征的单元弹性模量设置为E_max且在后续计算过程中始终保持不变；在E_max与E_min之间均匀地划分为H个等级（比如H＝300），在后续每一步迭代计算中，每一个单元的弹性模量计算值均要元整为与其最接近的那个等级所对应的弹性模量值。

2）利用有限元软件计算承载结构在设定边界条件下的应力分布状态，提取承载结构各单元质心处的应力状态（σ_xx,σ_xy,σ_xz,σ_yx,σ_yy,σ_yz,σ_zx,σ_zy,σ_zz），并计算各单元质心处的内力传递矢量（P_x，P_y，P_z），如式1所示，式（1）中i,j,k分别表示x，y，z三个坐标轴方向的单位向量：

P_ξ=σ_xξ·i+σ_yξ·j+σ_zξ·k(ξ=x,y,z)（1）

3）在单元质心处的三个内力传递矢量中选择其中模值最大的矢量P_i作为该单元的主传递矢量，将沿该方向的内力路径L_i定义为该单元的主传递路径，并计算获取单元质心处的曲率C_i。其中|P_i|是衡量载荷传递路径疏密度的指标量，而C_i则是衡量载荷传递路径弯曲度的指标量。

4）按照单元主传递矢量P_i的大小对承载结构中所有单元进行降序排列，获得单元列表{Sequence}。从该列表的顶层单元开始，对单元体积求和直至所得体积达到总体积的预期保留比例β为止，将此时对应单元的主传递矢量P_i定义为截止矢量P_cut。

5）对承载结构中各个单元的弹性模量进行惩罚更新（如式2所示），直至单元平均弹性模量达到β·E_max时程序结束。

$E_i^{\left(k\right)}=E_i^{(k-1)}\&CenterDot;F(P_i^{(k-1)},C_i^{(k-1)})$ $\left(2\right)$

$=E_i^{(k-1)}\&CenterDot;F_1{\left(P_i^{(k-1)}\right)}^{F_2\left(C_i^{(k-1)}\right)}$

$F_1\left(P_i^{(k-1)}\right)=\left(\frac{\left|P_i^{(k-1)}\right|}{\left|P_{\mathrm{cut}}^{(k-1)}\right|}\right)---\left(3\right)$

$F_2\left(C_i^{(k-1)}\right)={\mathrm{\&lambda;}}^{\frac{C_i^{(k-1)}}{C_i^{(k-1)}+1}}---\left(4\right)$

$\underset{i=1}{\overset{N}{\mathrm{\&Sigma;}}}{E}_i^{\left(k\right)}/N=\mathrm{\&beta;}\&CenterDot;E_{\max }---\left(5\right)$

式中，k表示迭代步数，λ表示优化系数，1<λ<10，N表示单元数量，是惩罚更新后单元e_i的弹性模量计算值，是惩罚更新前单元e_i的弹性模量计算值。当 $E_i^{\left(k\right)}>E_{\max }$ 时，取 $E_i^{\left(k\right)}=E_{\max };$ 当 $E_i^{\left(k\right)}<E_{\min }$ 时，取 $E_i^{\left(k\right)}=E_{\min }.$ F是弹性模量惩罚函数，主要包括两部分即疏密度惩罚因子和弯曲度惩罚因子，其中疏密度因子作为底数惩罚项（如式3所示）而弯曲度因子作为指数惩罚项（如式4所示）。

6）为了使有限元模型中弹性模量值的分布更平滑，以使其更符合实际情况。首先定义每一个节点的弹性模量值为节点所连接的材料单元的弹性模量的平均值，再将每个材料单元的弹性模量值根据它所包含的节点上的弹性模量值的平均值重新计算。为了避免计算过程中可能出现的棋盘格现象，需要对模型中各个材料单元的弹性模量进行光滑处理。对模型中的每一个节点赋予弹性模量索引值属性，即E_{index_j}(j=1,2,…,L）,该节点的弹性模量索引值E_{index_k}定义为与其连接的所有单元弹性模量计算值的平均；对模型中的每一个单元赋予弹性模量光滑值属性，即E_{smooth_i}（i=1,2,…，M）,该单元的弹性模量光滑值E_{smooth_i}定义为其上所有节点的弹性模量索引值的再次平均。

实施例

以一端固定，另一端中点处受集中力载荷的方板为例，介绍本方法的具体实施步骤（如图2）。

1)定义有限元模型中材料弹性模量的取值范围，即最大值E_max和最小值E_min；以及优化比例，即需要保留的材料所占的比例β。在该实例中，设定材料弹性模量的最大值E_max设置为5e5Pa。设置最小值E_min接近于0，但又要足够大以使有限元分析时不至于将其忽略，取E_min=1Pa。此外，承受集中载荷和被约束的部位被视为结构的特征而且必须保留下来，因此靠近这些部位的第一层单元的弹性模量定义为最大值E_max并且不允许发生变化。

2)建立该结构的有限元模型。方板尺寸设置为220mm，有限单元类型选Shell63单元。对各节点和单元进行有序的编号以方便后续的数据处理。设置各单元初始的材料属性，其中，各单元弹性模量的值均为E_max。施加边界条件并求解（如图2）。将模型左端自由度全部约束，右端中点处施加竖直向下100N的集中力。对模型进行求解。

3)计算每一个单元质心上的各个方向上的载荷向量（对于二维结构的受力问题，只存在两个内力传递矢量P_x和P_y），如图3以及图4所示，选择其中模值最大的矢量P_i作为该单元的主传递矢量，将沿该方向的内力路径L_i定义为该单元的主传递路径，并计算获取该点的曲率C_i，如图5以及图6所示。

4)按照单元主传递矢量P_i的大小对承载结构中所有单元进行降序排列，获得单元列表{Sequence}。从该列表的顶层单元开始，对单元体积求和直至所得体积达到总体积的预期保留比例β为止，将此时对应单元的主传递矢量P_i定义为截止矢量P_cut。对承载结构中各个单元的弹性模量进行惩罚更新计算。

5)为了使有限元模型中弹性模量值的分布更平滑以使其更符合实际情况。首先定义每一个节点的弹性模量值为他所连接的单元的弹性模量的平均值，再将每个单元的弹性模量值根据它所包含的节点上的弹性模量值的平均值重新计算。

6)返回第2步，将新的弹性模量的数值赋给有限元模型的各个单元，并重新利用该算法进行迭代计算。直到模型中各个单元的弹性模量的平均值达到β·E_max时程序结束，并输出此时的各个单元的弹性模量数据。这样即可得到以各个单元的弹性模量值来表示的承载结构材料最优分布结果（如图7所示）。

Claims (2)

1.一种基于内力路径几何形态的承载结构低碳节材设计方法，其特征在于，具体按下列步骤进行：

1)设定一种弹性模量可变的单元并将单元的弹性模量作为设计变量，预先设定弹性模量的变化范围；

2)建立承载结构的有限元模型，为有限元模型中各单元设定弹性模量初始值，然后施加边界条件并求解；

3)经过步骤2)后，从有限元分析结果中提取单元节点应力值，通过插值得到各单元质心处应力值，根据各单元质心处应力值得出各单元质心处的内力传递矢量，模值最大的内力传递矢量代表内力路径在单元质心处的疏密程度，以此为基础求解内力路径在单元质心处的曲率，曲率值代表了内力路径在单元质心处的弯曲程度；

4)利用步骤3)中得到的各个单元质心处内力路径的几何形态指标量对单元当前的弹性模量值进行惩罚，并获得新的单元弹性模量值，几何形态指标量包括疏密程度与弯曲程度；

5)对步骤4)中得到的各单元弹性模量值进行光滑平均化处理；

6)返回步骤2)将步骤5)中得到的单元弹性模量值赋给有限元模型中相应的各个单元，并进行循环迭代计算，直至单元的平均弹性模量达到预先设定的数值，这样便可得到承载结构的最佳传力路径，也即材料在设计空间中的最佳分布形式；

步骤3)中利用有限元软件计算承载结构在设定边界条件下的应力分布状态，提取承载结构各单元质心处的应力状态(σ_xx,σ_xy,σ_xz,σ_yx,σ_yy,σ_yz,σ_zx,σ_zy,σ_zz)，并计算各单元质心处的内力传递矢量(P_x，P_y，P_z)，如式1所示，式(1)中i,j,k分别表示x，y，z三个坐标轴方向的单位向量，ξ＝x,y,z：

P_ξ＝σ_xξ·i+σ_yξ·j+σ_zξ·k(1)

所述对单元当前的弹性模量值进行惩罚，包括以下步骤：

a)根据模值最大的内力传递矢量对承载结构中所有单元进行降序排列，获得单元列表，从列表的顶层单元开始，对单元体积求和直至所得体积达到总体积的预期保留比例为止，将此时对应单元的模值最大的内力传递矢量定义为截止矢量；

b)根据式2对承载结构中各个单元的弹性模量进行惩罚更新，直至单元平均弹性模量达到β·E_max时程序结束，β表示材料保留率，E_max表示弹性模量的最大值，P_i表示单元主传递矢量，C_i表示单元质心处的曲率；

$\begin{array}{c}{E}_i^{\left(k\right)}=E_i^{(k-1)}\&CenterDot;F(P_i^{(k-1)},C_i^{(k-1)})\\ =E_i^{(k-1)}\&CenterDot;F_1{\left(P_i^{(k-1)}\right)}^{F_2\left({C_i}^{(k-1)}\right)}\end{array}---\left(2\right)$

$F_1\left(P_i^{(k-1)}\right)=\left(\frac{\left|P_i^{(k+1)}\right|}{\left|P_{cut}^{(k-1)}\right|}\right)---\left(3\right)$

$F_2\left(C_i^{(k-1)}\right)={\mathrm{\&lambda;}}^{\frac{{c_i}^{(k-1)}}{{c_i}^{(k-1)}+1}}---\left(4\right)$

k表示迭代步数，λ表示优化系数，1＜λ＜10，P_cut表示截止矢量，E_i ^(k)是惩罚更新后单元的弹性模量计算值，E_i ^(k-1)是惩罚更新前单元的弹性模量计算值，当E_i ^(k)>E_max时，取E_i ^(k)＝E_max；当E_i ^(k)<E_min时，取E_i ^(k)＝E_min，E_min表示弹性模量的最小值。

2.如权利要求1所述的基于内力路径几何形态的承载结构低碳节材设计方法，其特征在于，弹性模量的最小值设定为一个接近0但不等于0的值，以避免有限元分析过程中刚度矩阵的奇异，弹性模量的最大值设定为承载结构选用材料的实际弹性模量值，并将最大值作为初始值赋给各个单元。