一种基于应力惩罚和自适应体积的结构拓扑优化设计方法与流程

本发明属于结构优化设计相关技术领域，更具体地，涉及一种基于应力惩罚和自适应体积的结构拓扑优化设计方法。

背景技术：

在工程实际应用中，结构的强度性能要求需要被满足，因此应力约束常作为结构设计中重要的一个考虑因素。但是现有的基于应力约束下的结构拓扑优化方法仍具有一些缺点：(1)采用局部应力法，即考虑设计域内各单元上的应力，每一个设计单元均需要添加一个应力约束，虽然可以较精确地控制结构各处的应力，但是导致了庞大的约束数目，造成昂贵的计算代价，计算效率低；(2)若采用全局应力法，即利用一个应力评价函数来考虑结构整体的应力，如结构的最大应力，虽然可以获得较高的计算效率，但是无法控制局部应力，并且会造成优化的不稳定和参数依赖性；(3)现有的应力控制方法往往是通过改变整个结构来降低应力，此时结构的其他性能(如刚度)，会大幅下降。

另一方面，由于汽车、航天航空等领域的发展，轻量化的设计需求被提出。但是轻量化的实施通常是在保证稳定提升结构性能的基础上和保证结构强度的前提下进行的。因此在结构的优化过程中，结构性能、强度要求和轻量化需要同时被考虑。在优化模型中实现结构的轻量化有两种方式，但均具有各自的缺点：(1)直接将结构体积(或材料用量)最小化设置为优化目标，但是这种处理方式会极大地降低其他的结构性能(如刚度)；(2)选取一个合适且较小的结构体积(或材料用量)作为约束条件，但是因为这个约束值往往是由设计人员凭借主观经验来设定的，所以它可能会导致体积约束和应力约束无法被同时满足，可行解很难被获得。

技术实现要素：

针对现有技术的以上缺陷或改进需求，本发明提供了一种基于应力惩罚和自适应体积的结构拓扑优化设计方法，用以解决需同时增大结构刚度和减少材料用量，并满足应力约束的结构优化问题。

为实现上述目的，按照本发明，提供了一种基于应力惩罚和自适应体积的结构拓扑优化设计方法，适用于考虑应力约束、刚度最大化和体积分数最小化的结构优化，其特征在于，所述方法包括以下步骤：

(1)利用基于应力惩罚和参数化水平集方法的优化模型，求解体积约束下基于应力的柔度最小化结构优化，获得体积约束下最优结构的材料分布和应力分布，再利用区间搜索方法调整体积约束，以缩小结构最优体积分数的搜索范围，从而获得最优结构的体积分数上限值，具体包括以下子步骤：

(1.1)拓扑优化初始化：给定结构设计域、载荷和边界条件，设定许用应力、初始体积约束和初始结构，并对优化算法的参数进行初始化；

(1.2)对结构进行有限元分析，以获得结构位移场，计算并记录当前结构的柔度值和最大应力值；

(1.3)将柔度乘以应力惩罚函数作为目标函数，计算目标函数及体积约束对设计变量的灵敏度；

(1.4)通过获得的灵敏度构建优化准则，利用优化准则更新设计变量及水平集方程；

(1.5)判断判别准则a是否满足，若满足判别准则a，记录当前结构及其性能参数，并转至步骤(1.8)，否则转至下一步骤，其中，判别准则a为体积约束下基于应力的柔度最小化拓扑优化收敛条件，所记录的当前结构及其性能参数是体积分数为当前体积分数的搜索区间上限值时所获得的结构及性能参数；

(1.6)判断判别准则b是否满足，若满足判别准则b，转至下一步骤，否则返回步骤(1.2)，其中，判别准则b为自适应应力惩罚因子调整策略执行条件；

(1.7)利用自适应调整策略调整应力惩罚因子，并返回步骤(1.2)；

(1.8)判断判别准则c是否满足，若不满足判别准则c，转至下一步骤，否则进入步骤2)，其中，判别准则c为区间搜索方法执行条件；

(1.9)利用区间搜索方法调整体积约束，确定结构最优体积分数的搜索区间，并返回步骤(1.2)；

(2)利用基于应力惩罚和参数化水平集方法的优化模型，求解体积约束下基于应力的柔度最小化结构优化，获得体积约束下最优结构的材料分布和应力分布，再利用局部搜索方法调整体积约束，以获得最优结构的体积分数，并输出最优结构。

优选地，步骤(2)包括以下子步骤：

(2.1)重新进行优化初始化：将结构及其性能参数设置为步骤(1.5)中记录的结构及其性能参数；

(2.2)利用局部搜索方法调整体积约束，然后对结构进行有限元分析，以获得结构位移场，计算并记录当前结构的柔度值和最大应力值；

(2.3)将柔度乘以应力惩罚函数作为目标函数，使其最小，计算目标函数及体积约束对设计变量的灵敏度；

(2.4)通过获得的灵敏度构建优化准则，利用优化准则法更新设计变量及水平集方程；

(2.5)判断判别准则a是否满足，若满足判别准则a，转至下一步骤，否则返回步骤(2.2)；

(2.6)判断判别准则d是否满足，若不满足判别准则d，转至下一步骤，若满足判别准则d，结束优化并输出最优拓扑结构，其中，判别准则d是局部搜索方法执行条件；

(2.7)利用局部搜索方法调整体积约束，并返回步骤(2.2)。

优选地，步骤(1)中的适用于考虑应力约束、刚度最大化和体积分数最小化的结构优化模型如以下公式(1)所示：

其中，f是优化目标，是体积约束值，j(u,φ)是柔度，其用于评价结构刚度性能，ω是结构设计域，u和v分别表示实位移场和虚位移场，u是运动学允许的位移空间，u0是dirichlet边界上的位移，φ是水平集方程，并且是由局部径向基函数的形状方程组成的矩阵，β是由设计变量βi组成的向量，σv和分别是vonmises应力和许用应力，是容忍因子，结构的最大应力σv,max被期望接近许用应力以减少材料用量。g(φ)是结构体积，h(φ)是heaviside函数，βi,min和βi,max分别是设计变量βi的上下限，βi为参数化水平集方法中第i个点上的扩展系数，a(u,v,φ)＝l(v,φ)是弹性平衡条件的弱形式，能量双线性形式a(u,v,φ)和载荷线性形式l(v,φ)分别表示如下：

a(u,v,φ)＝∫ωε^t(u)dε(v)h(φ)dω

其中，ε是应变场，ε^t(u)＝(bu)^t，ε(v)＝bv，b是形状函数，t代表矩阵的转置，d是弹性刚度，p是结构体积力，τ为边界上的牵引力，δ(φ)为dirac函数。

优选地，用于求解体积约束下基于应力的柔度最小化结构优化的基于应力惩罚和参数化水平集方法的优化模型的具体表达式如公式(2)所示：

其中，是同时优化柔度和应力的优化目标，ω是结构设计域，ε是应变场，d是弹性刚度，α是应力惩罚因子，在优化过程中由自适应应力惩罚因子调整策略调整，hobj(·)是heaviside函数，u和v分别表示实位移场和虚位移场，u是运动学允许的位移空间，u0是dirichlet边界上的位移，φ是水平集方程，并且是由局部径向基函数的形状方程组成的矩阵，β是由设计变量βi组成的向量，σv和分别是vonmises应力和许用应力，g(φ)是结构体积，h(φ)是heaviside函数，是体积约束值，βi,min和βi,max分别是设计变量βi的上下限，βi为参数化水平集方法中第i个点上的扩展系数，a(u,v,φ)＝l(v,φ)是弹性平衡条件的弱形式，能量双线性形式a(u,v,φ)和载荷线性形式l(v,φ)分别表示为：

a(u,v,φ)＝∫ωε^t(u)dε(v)h(φ)dω

其中，p是结构体积力，τ为边界上的牵引力，δ(φ)为dirac函数。

优选地，自适应应力惩罚因子调整策略如公式(3)所示：

α＝α+χα,当满足和和时(3)

其中，α是应力惩罚因子，χα是设定的应力惩罚因子调整值，是在第k次迭代后结构所对应的最大应力，ξ是一个极小的正数，是许用应力。

优选地，对应用于求解体积约束下基于应力的柔度最小化结构优化的基于应力惩罚和参数化水平集方法的优化模型，其目标函数及体积约束对设计变量的灵敏度可以分别如公式(4)和公式(5)所示：

其中，是同时优化柔度和应力的优化目标，ω是结构设计域，α是应力惩罚因子，在优化过程中由自适应应力惩罚因子调整策略调整，hobj(·)是heaviside函数，u和v分别表示实位移场和虚位移场，u是运动学允许的位移空间，u0是dirichlet边界上的位移，φ是水平集方程，并且是由局部径向基函数的形状方程组成的矩阵，是局部径向基函数的形状方程，x是设计域坐标，β是由设计变量βi组成的向量，σv和分别是vonmises应力和许用应力，δ(φ)为dirac函数，g(φ)是结构体积，h(φ)是heaviside函数，是体积约束值，βi是参数化水平集方法中第i个点上的扩展系数，同时是本方法中的设计变量，ke＝b^tdb，ce＝b^td^tvdb，d是弹性刚度，b是应变-位移矩阵，v为设置的3×3矩阵。

优选地，体积确定由区间搜索和局部搜索方法共同求解，其中，区间搜索方法表示如公式(6)所示：

其中，是第i次体积约束更改后的体积分数，χ是区间搜索步长，和分别是在第i-1次和第i次体积约束更改后体积分数为的最优结构所对应的最大应力，是许用应力；区间搜索方法的终止条件表示如公式(7)所示：

当区间搜索方法的终止条件被满足时，若则可获得合适的体积约束值的搜索区间若则可获得合适的体积约束值的搜索区间

优选地，局部搜索方法表示如公式(8)所示：

其中，和是分别是第i+1次和第i次体积约束更改后的体积分数，j代表当局部搜索方法被实施后体积约束第j次改变，为初始体积约束值，是在第i次体积约束更改后体积分数为的最优结构所对应的最大应力，是许用应力。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：

(1)在所述的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法中，自适应体积约束算法被采用，将以应力为约束，以刚度最大化和体积分数最小化为目标的结构优化问题分解成一个体积约束下基于应力的柔度最小化问题和一个体积确定问题，运用循环的方式进行求解，简化了优化问题的求解过程，并且避免了直接将结构体积最小化作为优化目标或人为主观选取体积约束带来的弊端；

(2)在所述的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法中，用于求解体积约束下基于应力的柔度最小化结构优化的基于应力惩罚和参数化水平集方法的优化模型被提出，其在可以控制结构局部应力的同时具有较高的计算效率，并且避免了像传统处理应力约束的方法一样为了使应力约束满足而去改变整个结构，导致结构的其他性能(如刚度)大幅下降。参数化水平集方法被用于描述和更新拓扑结构，可以保证获得的结构具有清晰光滑的边界，保证应力计算和结构描述的准确性；

(3)在所述的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法中，区间搜索和局部搜索相结合的方法被用于求解体积确定问题，提高了优化效率；

(4)所述的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法适用于连续体结构，适用范围广，简单易行；

(5)采用所述的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法，优化后获得的结构具有高刚度、轻量化的优点，并且结构的强度要求被满足，应力集中的问题被缓解。

附图说明

图1是按照本发明所构思的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法的基本流程图；

图2是用于示范性显示l型梁的载荷与边界条件示意图；

图3a和图3b分别是用于示范性显示l型梁初始结构的材料分布图和其对应的应力分布图；

图4a、图4b分别是用于示范性显示l型梁在优化过程中，迭代步数为151，对应的结构体积分数为0.5时的材料分布图和应力分布图；

图5a、图5b分别是用于示范性显示l型梁在优化过程中，迭代步数为303，对应的结构体积分数为0.4；

图6a、图6b分别是用于示范性显示l型梁在优化过程中，迭代步数为455，对应的结构体积分数为0.3；

图7a、图7b分别是用于示范性显示l型梁在优化过程中，迭代步数为908，对应的结构体积分数为0.17426，其为优化后获得的最优结构。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。

参照各附图，一种基于应力惩罚和自适应体积的结构拓扑优化设计方法为求解公式(1)描述的需同时增大结构刚度和减少材料用量，并满足应力约束的结构优化问题，利用自适应体积约束算法将原优化问题分解成公式(2)描述的体积约束下基于应力的柔度最小化问题和公式(3)描述的体积确定问题。优化模型表示如下

其中，f是优化目标，是体积约束值，j(u,φ)是柔度，其用于评价结构刚度性能，ω是结构设计域，u和v分别表示实位移场和虚位移场，u是运动学允许的位移空间，u0是dirichlet边界上的位移，φ是水平集方程，并且是由局部径向基函数的形状方程组成的矩阵，β是由设计变量βi组成的向量，σv和分别是vonmises应力和许用应力，是容忍因子，结构的最大应力σv,max被期望接近许用应力以减少材料用量，g(φ)是结构体积，h(φ)是heaviside函数，βi,min和βi,max分别是设计变量βi的上下限，βi为参数化水平集方法中第i个点上的扩展系数，a(u,v,φ)＝l(v,φ)是弹性平衡条件的弱形式，能量双线性形式a(u,v,φ)和载荷线性形式l(v,φ)分别表示如下：

a(u,v,φ)＝∫ωε^t(u)dε(v)h(φ)dω

其中，ε是应变场，ε^t(u)＝(bu)^t，ε(v)＝bv，b是形状函数，t代表矩阵的转置，d是弹性刚度，p是结构体积力，τ为边界上的牵引力，δ(φ)为dirac函数。

其中，是同时优化柔度和应力的优化目标，ω是结构设计域，ε是应变场，d是弹性刚度，α是应力惩罚因子，在优化过程中由自适应应力惩罚因子调整策略调整，hobj(·)是heaviside函数，u和v分别表示实位移场和虚位移场，u是运动学允许的位移空间，u0是dirichlet边界上的位移，φ是水平集方程，并且是由局部径向基函数的形状方程组成的矩阵，β是由设计变量βi组成的向量，σv和分别是vonmises应力和许用应力，g(φ)是结构体积，h(φ)是heaviside函数，是体积约束值，βi,min和βi,max分别是设计变量βi(参数化水平集方法中第i个点上的扩展系数)的上下限，a(u,v,φ)＝l(v,φ)是弹性平衡条件的弱形式，能量双线性形式a(u,v,φ)和载荷线性形式l(v,φ)分别表示为：

a(u,v,φ)＝∫ωε^t(u)dε(v)h(φ)dω(4)

其中，ε是应变场，d是弹性刚度，p是结构体积力，τ为边界上的牵引力，δ(φ)为dirac函数。

请参阅图1，通过所述的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法，上述问题的求解过程包括以下两个阶段：

阶段一，利用基于应力惩罚和参数化水平集方法的优化模型求解体积约束下基于应力的柔度最小化(即刚度最大化)问题，利用区间搜索方法调整体积约束，以缩小结构最优体积分数的搜索范围，获得最优结构的体积分数上限值，阶段一包括以下子步骤：

(1.1)拓扑优化问题初始化，给定结构设计域、载荷和边界条件，设定许用应力、初始体积约束和初始结构，并对优化算法的参数进行初始化。

(1.2)对结构进行有限元分析，以获得结构位移场，计算并记录当前结构的柔度值和最大应力值。

(1.3)将柔度乘以应力惩罚函数作为目标函数，使其最小，计算目标函数及体积约束对设计变量的灵敏度。用于求解体积约束下基于应力的柔度最小化结构优化的基于应力惩罚和参数化水平集方法的优化模型表达式为：

其中，是同时优化柔度和应力的优化目标，ω是结构设计域，ε是应变场，d是弹性刚度，α是应力惩罚因子，在优化过程中由自适应应力惩罚因子调整策略调整，hobj(·)是heaviside函数，u和v分别表示实位移场和虚位移场，u是运动学允许的位移空间，u0是dirichlet边界上的位移，φ是水平集方程，并且是由局部径向基函数的形状方程组成的矩阵，β是由设计变量βi组成的向量，σv和分别是vonmises应力和许用应力，g(φ)是结构体积，h(φ)是heaviside函数，是体积约束值，βi,min和βi,max分别是设计变量βi(参数化水平集方法中第i个点上的扩展系数)的上下限，a(u,v,φ)＝l(v,φ)是弹性平衡条件的弱形式，

对于上述基于应力惩罚和参数化水平集方法的优化模型，其目标函数及体积约束对设计变量的灵敏度可以分别表示为：

其中，ke＝b^tdb，ce＝b^td^tvdb，d是弹性刚度，b是应变-位移矩阵，v为设置的3×3矩阵并且其优选定义为：

(1.4)通过获得的灵敏度构建优化准则，利用优化准则法更新设计变量及水平集方程。

(1.5)判断判别准则a(体积约束下基于应力的柔度最小化拓扑优化收敛条件)是否满足，若满足判别准则a，记录当前结构及其性能参数，并转至步骤(1.8)，否则转至下一步骤。判别准则a为:

其中，ξ和ζ都是正数，k是迭代次数，是设定的每次体积约束改变后的最大迭代次数。

(1.6)判断判别准则b(自适应应力惩罚因子调整策略执行条件)是否满足，若满足判别准则b，转至下一步骤，否则转至步骤(1.2)。判别准则b为：

其中，是在第k次迭代后结构所对应的最大应力，ξ是一个极小的正数。

(1.7)利用自适应调整策略调整应力惩罚因子，并转至步骤(1.2)。通过公式(1.2)调整应力惩罚因子：

α＝α+χα(12)

其中，χα是设定的应力惩罚因子调整值。

(18)判断判别准则c(区间搜索方法执行条件)是否满足，若不满足判别准则c，转至下一步骤，否则转至步骤(2.1)。判别准则c为：

其中，是在第i次体积约束更改后体积分数为的最优结构所对应的最大应力。

(1.9)利用区间搜索方法调整体积约束，确定结构最优体积分数的搜索区间，并转至步骤(1.2)。区间搜索方法被定义为：

其中，χ是区间搜索步长。

阶段二，利用基于应力惩罚和参数化水平集方法的优化模型求解体积约束下基于应力的柔度最小化(即刚度最大化)问题，利用局部搜索方法调整体积约束，以获得最优的结构体积分数，并输出最优结构。

进一步地，阶段二包括以下子步骤：

(2.1)重新进行拓扑优化问题初始化，将结构及其性能参数设置为步骤(1.5)中记录的体积分数为当前体积分数的搜索区间上限值时所获得的结构及其性能参数，并且利用局部搜索方法调整体积约束。局部搜索方法被定义为：

其中，其中，和是分别是第i+1次和第i次体积约束更改后的体积分数，j代表当局部搜索方法被实施后体积约束第j次改变，为初始体积约束值，是在第i次体积约束更改后体积分数为的最优结构所对应的最大应力，是许用应力。

(2.2)对结构进行有限元分析，以获得结构位移场，计算并记录当前结构的柔度值和最大应力值。

(2.3)将柔度乘以应力惩罚函数作为目标函数，使其最小，计算目标函数及体积约束对设计变量的灵敏度。用于求解体积约束下基于应力的柔度最小化结构优化的基于应力惩罚和参数化水平集方法的优化模型表达式为：

其中，a(u,v,φ)＝l(v,φ)是弹性平衡条件的弱形式，能量双线性形式a(u,v,φ)和载荷线性形式l(v,φ)分别表示为：

a(u,v,φ)＝∫ωε^t(u)dε(v)h(φ)dω

其中，p是结构体积力，τ为边界上的牵引力，δ(φ)为dirac函数。

对于上述的基于应力惩罚和参数化水平集方法的优化模型，其目标函数及体积约束对设计变量的灵敏度可以分别表示为：

其中，是局部径向基函数的形状方程，ke＝b^tdb，ce＝b^td^tvdb，b是应变-位移矩阵，v是被定义为：

(2.4)通过获得的灵敏度构建优化准则，利用优化准则法更新设计变量及水平集方程。

(2.5)判断判别准则a(体积约束下基于应力的柔度最小化拓扑优化收敛条件)是否满足，若满足判别准则a，转至下一步骤，否则转至步骤(2.2)。判别准则a为：

其中，ξ和ζ都是正数，k是迭代次数，是设定的每次体积约束改变后的最大迭代次数。

(2.6)判断判别准则d(局部搜索方法执行条件)是否满足，若不满足判别准则d，转至下一步骤，若满足判别准则d，结束优化并输出最优拓扑结构。判别准则d为

(2.7)利用局部搜索方法调整体积约束，并转至步骤(2.2)。局部搜索方法被定义为：

其中，j代表当局部搜索方法被实施后体积约束第j次改变，为初始体积约束值。

请参阅图2、图3a、图3b、图4a、图4b，以下以l型梁的设计来进一步说明本发明的方法，图2展示了l型梁的设计域，其中l1＝100mm，l2＝100mm，l3＝40mm，l4＝60mm。在优化过程中，结构设计域被划分成100×100的正方形网格，材料弹性模量为200gpa，泊松比为0.3的材料被用于l型梁的设计中，l型梁的顶部被固定，一个集中力f＝200kn被施加到l型梁右端的顶点上，初始体积约束值许用应力为130mpa，初始的应力惩罚因子α0＝5；优化目标为结构柔度最小化(即结构柔度最大化)和结构体积最小化(即轻量化)，结构应力被约束。

图3a、图3b展示了l型梁初始结构的材料分布图和其对应的应力分布图，从应力分布图中可以看出，在l型梁的内角处出现较高的应力，即应力集中。

基于自适应体积约束算法，优化问题被分解成一个体积约束下基于应力的柔度最小化问题和一个体积确定问题，运用循环的方式进行求解；基于应力惩罚和参数化水平集方法的优化模型被构建，用于求解体积约束下基于应力的柔度最小化问题；区间搜索和局部搜索方法被共同用于求解体积确定问题；图4a和图4b、图5a和图5b、图6a和图6b分别展示了在结构优化过程中第151、303、455次迭代后得到的l型梁的材料分布图和应力分布图。当优化结束时，l型梁的最优构型的材料分布图和应力分布图如图7a和图7b所示；图7a和图7b显示的最优结构的柔度为41096.27，通过柔度最小化问题被优化，其体积分数为0.17426，相比初始体积约束值0.5，下降了超过60％，实现结构的轻量化，其最大应力值129.71mpa，应力约束被满足，并且l型梁内角处的应力集中问题被缓解。

本发明提供的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法中，自适应体积约束算法被采用，将以应力为约束，以刚度最大化和体积分数最小化为目标的结构优化问题分解成一个体积约束下基于应力的柔度最小化问题和一个体积确定问题，运用循环的方式进行求解，简化了优化问题的求解过程，并且避免了直接将结构体积最小化直接作为优化目标或人为主观选取体积约束带来的弊端；在所述的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法中，用于求解体积约束下基于应力的柔度最小化结构优化的基于应力惩罚和参数化水平集方法的优化模型被提出，其在可以控制结构局部应力的同时具有较高的计算效率，并且避免了像传统处理应力约束的方法一样为了使应力约束满足而去改变整个结构，导致结构的其他性能(如刚度)大幅下降。参数化水平集方法被用于描述和更新拓扑结构，可以保证获得的结构具有清晰光滑的边界，保证应力计算和结构描述的准确性；在所述的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法中，区间搜索和局部搜索方法被用于共同求解体积确定问题，提高了优化效率；所述的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法适用于连续体结构，适用范围广，简单易行；采用所述的一种基于应力惩罚和自适应体积的结构拓扑优化设计方法，优化后获得的结构具有高刚度、轻量化的优点，并且结构的强度要求被满足，应力集中的问题被缓解。

本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。