一种考虑应力约束的等几何拓扑优化设计方法与应用

1.本发明属于结构拓扑优化设计相关技术领域，更具体地，涉及一种考虑应力约束的等几何拓扑优化设计方法与应用。


背景技术：

2.拓扑优化技术作为一种高效的结构优化设计工具，其在航空航天、海洋船舶、车辆工程等众多领域都具有重要应用，同时也是当前学术研究的热点。但目前大部分关于拓扑优化技术的研究工作都聚焦于结构的刚度特性，因为结构刚度设计问题具有易于数值实现、稳定收敛等优点。在实际工程中，除需要考虑结构的刚度外，结构的强度特性更为重要，因为结构强度是保证结构正常发挥效能的前提条件。表征结构强度的指标是结构应力状态，因此需要研究基于应力的拓扑优化设计方法。
3.基于应力的拓扑优化设计分为两类，即应力最小化设计和应力约束设计，其中应力约束设计更符合工程设计的常规模式也更具挑战性。同基于刚度的拓扑优化设计相比，应力相关的设计问题存在以下三大难点：
4.1.奇异现象：奇异现象指优化问题的全局最优解位于可行域内非平滑过渡的退化子空间中，采用常规的优化求解算法无法搜寻到该类退化子空间，导致难以获得优化问题的全局最优解。
5.2.应力的局部特性：结构应力状态是与结构内任意一点空间坐标相关的场函数，且应力状态分布会随结构材料分布变化而变化。因此控制结构的最大应力水平需要在结构内任意一点上设置应力约束，这将造成优化问题中包含大量的约束条件，导致优化问题计算求解困难。
6.3.应力的高度非线性行为：结构内任意一点应力状态与该点周围的材料分布及局部几何特征高度相关，且两者间关系具有非线性特征。因此，结构应力水平对结构拓扑变化非常敏感，应力的高度非线性行为会导致优化过程难以稳定收敛。


技术实现要素：

7.针对现有技术的以上缺陷或改进需求，本发明提供了一种考虑应力约束的等几何拓扑优化设计方法与应用，其采用增广拉格朗日函数将原优化模型转变为无约束优化模型，实现了结构各局部区域应力状态的高效率、高精度控制，有效解决了包含应力约束的结构拓扑优化设计问题。
8.为实现上述目的，按照本发明的一个方面，提供了一种考虑应力约束的等几何拓扑优化设计方法，该方法主要包括以下步骤：
9.(1)采用改进的应力消失约束函数调换待优化结构的等几何拓扑优化设计模型中的应力约束条件，得到改进拓扑优化设计模型；改进的应力消失约束函数为：
10.11.其中，m
e,
(ρ)为评估点a处的杨氏模量；表示评估点a处的冯米塞斯应力，σ
lim
为应力约束值，nv指总的应力评估点数目；
12.(2)采用增广拉格朗日函数法构造所述改进拓扑优化设计模型的近似模型，进而推导并计算所述近似模型的灵敏度信息；
13.(3)基于得到的灵敏度信息采用无约束mma算法实现所述近似模型的设计变量的迭代更新，每经过n
l
次设计变量更新后，基于kkt条件更新一次拉格朗日乘子，直至所述近似模型的优化设计结果收敛，进而采用当前的近似模型进行结构的等几何拓扑优化设计；其中，n
l
》1。
14.进一步地，步骤(1)之前还包括以材料用量最小为目标、结构评估点处冯米塞斯应力为约束来构建待优化结构的等几何拓扑优化设计模型的步骤。
15.进一步地，优化设计变量ρ与结构密度场mv(ρ)及杨氏模量场me(ρ)之间的对应关系为：
16.ρ＝[ρ
1，1
...ρ
i，j
…
ρ
n，m
]
[0017][0018]
me(ρ)＝e
min
+(e
0-e
min
)[mv(ρ)]
p
[0019]
其中，i，j表示沿参数化坐标ξ和η方向上的控制点编号，ρ
i，j
表示i，j号控制点对应的设计变量，表示对应于i，j号控制点的nurbs基函数，r和q分别指沿ξ和η方向上的基函数阶次；指经过阈值投影处理后的设计变量，n和m分别表示沿ξ和η方向上的控制点总数；e0指实体材料杨氏模量，e
min
是设置的最小杨氏模量；p为大于1的惩罚系数，设置为3.5。
[0020]
进一步地，等几何拓扑优化设计模型的表达式为：
[0021][0022]
式中，k为结构刚度矩阵，u为控制点位移矢量，f为外载荷；表示评估点a处的冯米塞斯应力，σ
lim
为应力约束值，nc指总的应力评估点数目。
[0023]
进一步地，改进拓扑优化设计模型的表达式为：
[0024][0025]
进一步地，近似模型的表达式为：
[0026][0027]
其中，s为松弛变量，λa和μ均为拉格朗日乘子，ha(ρ)为引入松弛变量后的应力约束；ρ为优化设计变量；mv(ρ)为结构密度场；sa为松弛变量。
[0028]
进一步地，
[0029]
ha(ρ)＝ga(ρ)+sa[0030]
拉格朗日乘子λa和μ的更新准则分别为：
[0031]
μ
k+l
＝αμk，α＞1
[0032][0033]
其中，k代表拉格朗日乘子迭代次数，α为设定的更新系数。
[0034]
本发明提供了一种等几何拓扑优化设计设备，所述设备包括存储器及处理器，所述存储器储存有计算机程序，所述处理器执行所述计算机程序时执行如上所述的考虑应力约束的等几何拓扑优化设计方法。
[0035]
本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有机器可执行指令，所述机器可执行指令在被处理器调用和执行时，所述机器可执行指令促使所述处理器实现如上所述的考虑应力约束的等几何拓扑优化设计方法。
[0036]
总体而言，通过本发明所构思的以上技术方案与现有技术相比，本发明提供的考虑应力约束的等几何拓扑优化设计方法与应用主要具有以下有益效果：
[0037]
1.本发明通过采用改进的应力消失约束函数表达式替换应力约束条件能够有效提升拓扑优化设计结果的收敛速度，使得结构的最大冯米塞斯应力在优化迭代初期快速降至应力约束值附近。
[0038]
2.本发明采用的增广拉格朗日乘数法既解决了优化模型中存在大量应力约束条件所造成的计算困难，同时实现了局部应力约束控制，这种控制策略可精准约束结构最大应力水平。
[0039]
3.本发明所提供的方法能够在保证结构最大冯米塞斯应力满足应力约束条件的同时实现结构轻量化或者其他优化目标的设计。
附图说明
[0040]
图1是本发明提供的考虑应力约束的等几何拓扑优化设计方法的流程示意图；
[0041]
图2中的(a)、(b)分别是应力优化领域典型算例“l”形梁优化算例的设计域及边界条件设置及“l”形梁结构应力分布云图；
[0042]
图3中的(a)、(b)分别是包含体积约束的柔度最小化设计得到的拓扑结构及该拓扑结构对应的应力分布云图；
[0043]
图4中的(a)、(b)分别是考虑应力约束的体积最小化设计得到的拓扑结构及其对应的应力分布云图；
[0044]
图5是结构最大冯米塞斯应力和体积分数的收敛曲线。
具体实施方式
[0045]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0046]
请参阅图1，本发明提供了一种考虑应力约束的等几何拓扑优化设计方法，所述方法主要包括以下步骤：
[0047]
步骤一，以材料用量最小为目标、结构评估点处冯米塞斯应力为约束来构建待优化结构的等几何拓扑优化设计模型。
[0048]
具体地，基于待优化结构的等几何分析模型建立优化设计变量ρ与结构密度场mv(ρ)及杨氏模量场me(ρ)之间的对应关系，进而以材料用量最小化为目标，结构评估点处vonmises应力为约束来构建等几何拓扑优化设计模型。
[0049]
本实施方式中，在建立等几何拓扑优化设计模型前先基于等几何分析模型引入设计变量ρ与结构密度场mv(ρ)及杨氏模量场me(ρ)间的对应关系：
[0050]
ρ＝[ρ
1，1
…
ρ
i，j
…
ρ
n，m
]
ꢀꢀꢀ
(1)
[0051][0052]
me(ρ)＝e
min
+(e
0-e
min
)[my(ρ)]
p
ꢀꢀꢀ
(3)
[0053]
其中，i,j表示沿参数化坐标ξ和η方向上的控制点编号，ρ
i,j
表示i,j号控制点对应的设计变量，表示对应于i,j号控制点的nurbs基函数，r和q分别指沿ξ和η方向上的基函数阶次；指经过阈值投影处理后的设计变量，n和m分别表示沿ξ和η方向上的控制点总数；e0指实体材料杨氏模量，e
min
是设置的最小杨氏模量，目的是避免刚度矩阵发生奇异，在本本实施方式中设置为0.001mpa；p为大于1的惩罚系数，设置为3.5。
[0054]
由经平滑处理的设计变量阈值投影后得到，投影表达式为：
[0055][0056]
其中，β为投影曲率控制参数，β值越大则投影函数曲线曲率越大，η为投影阈值，大于η的变量经投影后趋近于1，反之趋近于0，在本实施方式中η设置为0.5；由设计变量ρ
i,j
经平滑处理后得到，其表达式为：
[0057][0058]
其中，是c4紧支撑径向基函数。
[0059]
基于式(1)-(3)即可建立等几何拓扑优化设计模型，对应的数学表达式为：
[0060][0061]
其中，k为结构刚度矩阵，u为控制点位移矢量，f为外载荷；表示评估点a处的冯米塞斯应力，σ
lim
为应力约束值，nc指总的应力评估点数目。
[0062]
具体地，评估点处的冯米塞斯应力可表示为：
[0063][0064]
其中，σa为评估点处的应力矢量，v指应力系数矩阵，对于平面应力问题σa和v的表达式分别如式(7)和(8)所示：
[0065][0066][0067]
其中，d0表示实体材料的弹性矩阵，ba表示评估点处的应变位移矩阵，表示包含评估点a的跨度单元对应的控制点位移索引矩阵。
[0068]
步骤二，采用改进的应力消失约束函数调换所述等几何拓扑优化设计模型中的应力约束条件，得到改进拓扑优化设计模型；改进的应力消失约束函数为：
[0069][0070]
其中，m
e,a
(ρ)为评估点a处的杨氏模量；表示评估点a处的冯米塞斯应力，σ
lim
为应力约束值，nc指总的应力评估点数目。
[0071]
本实施方式中，为提升优化设计结果的收敛速度，提高设计效率，本实施方式对式(6)所述优化模型中的应力约束条件进行了改进，改进后的应力约束表达式为：
[0072][0073]
改进后的应力约束乘上了评估点处的杨氏模量值，使之变为消失约束，规避了奇异现象。同时将原来的一次约束函数变为二次约束函数，该策略可在保证优化过程稳定性的同时增加对结构中高应力区域的惩罚力度，使结构应力水平快速逼近应力约束值，提高了优化设计的收敛速度。改进后的拓扑优化设计模型的表达式为：
[0074][0075]
步骤三，采用增广拉格朗日函数法构造所述改进拓扑优化设计模型的近似模型。
[0076]
通过增广拉格朗日函数法构造的上述改进拓扑优化设计模型的近似模型的表达式为：
[0077][0078]
其中，s为松弛变量，λa和μ均为拉格朗日乘子，ha(ρ)为引入松弛变量后的应力约束，其具体表达式为：
[0079]
ha(ρ)＝ga(ρ)+sa[0080]
拉格朗日乘子λa和μ的更新准则分别为：
[0081]
μ
k+1
＝αμk，α＞1
[0082][0083]
其中，k代表拉格朗日乘子迭代次数，α为设定的更新系数。
[0084]
本实施方式中，式(11)所示的优化设计模型中包含了nc个应力评估点，在本实施方式中应力评估点设置为每个跨度单元的高斯积分点，因此nc就等于总的高斯积分点数目，这就使得上述优化设计模型中包含了过多应力约束条件，易造成求解困难。为此，本实施方式采用增广拉格朗日函数法构造上述优化设计模型的近似模型，将其中的应力约束条
件整合至目标函数中。具体地，先引入松弛变量sa：
[0085]
sa≥0，a＝1，2，...，ncꢀꢀꢀ
(12)
[0086]
令
[0087]
ha(ρ)＝ga(ρ)+sa＝0
ꢀꢀꢀ
(13)
[0088]
引入拉格朗日乘子λa和μ构造惩罚函数p(ρ)
[0089][0090]
基于惩罚函数p(ρ)和目标函数构造拉格朗日函数l(ρ)得到：
[0091]
l(ρ)＝∫
ω
mv(ρ)dω+p(ρ)
ꢀꢀꢀ
(15)
[0092]
建立近似问题模型的表达式为：
[0093][0094]
在求解上述近似问题前，可先通过求解式(17)所示的惩罚项子问题消除松弛变量：
[0095][0096]
上述惩罚项子问题的解为：
[0097][0098]
将式(18)代入近似问题模型的表达式(16)即可消除松弛变量sa。
[0099]
步骤四，推导并计算所述近似模型的灵敏度信息。
[0100]
由式(16)可知，需求解的灵敏度信息为拉格朗日函数l(ρ)关于设计变量的一阶导，即：
[0101][0102]
其中，由式(2)(4)(5)可得：
[0103][0104]
由式(15)得：
[0105][0106]
其中结构总体积可由高斯积分法表述为：
[0107][0108]
式(22)中，和分别为相对空间到物理空间两个逆映射的雅可比矩阵，ωa为高斯积分权重系数，故：
[0109][0110]
由式(14)可知
[0111][0112]
由式(18)可知，当ha(ρ)取时，上式为0，当ha(ρ)取ga(ρ)时有：
[0113][0114]
由式(7)(8)可得：
[0115][0116]
由伴随法可得：
[0117][0118]
其中整体刚度矩阵k同样可由高斯积分法表示为：
[0119][0120]
故由上式可得：
[0121][0122]
步骤五，基于得到的灵敏度信息采用无约束mma算法实现所述近似模型的设计变量的迭代更新，每经过n
l
次设计变量更新后，基于kkt条件更新一次拉格朗日乘子，直至所述近似模型的优化设计结果收敛，进而采用当前的近似模型进行结构的等几何拓扑优化设计。
[0123]
本实施方式中，式(16)所示近似模型已将全部应力约束条件整合至目标函数中，故本实施方式基于近似模型的目标函数的灵敏度信息，采用无约束mma算法实现设计变量的迭代更新。该更新算法适合解决应力约束这一类具有较强非线性的优化设计问题，同时可进一步提升优化问题的求解效率，在一定程度上保证优化过程的稳定性。
[0124]
通过求解近似优化模型所得优化设计结果是原优化问题模型式(11)的近似解，为获得原优化问题的精确解，需在优化迭代过程中基于kkt条件不断更新拉格朗日乘子使近似优化模型逐渐逼近原优化问题模型直至得到稳定收敛的优化设计结果。基于kkt条件推导得出的拉格朗日乘子λa和μ的更新准则为：
[0125][0126]
其中，k代表拉格朗日乘子的迭代次数，为保证优化过程的稳定性，拉格朗日乘子的更新频率应小于设计变量的更新频率，通常设置为每经过n
l
(n
l
》1)次设计变量更新后更新一次拉格朗日乘子，在本实施方式中n
l
＝5。α为设定的更新系数，且α值越大，设计结果收敛速度越快，但算法稳定性越差，因此α值需通过经验设置，在本实施方式中α值设置为1.1。
[0127]
为了证明本发明的有效性，以下以一个具体算例对本发明进行验证。
[0128]
如图2中的(a)所示为应力优化领域典型算例“l”形梁优化算例的设计域及边界条件设置，“l”形梁上边界完全约束，在其右边界上部长度为0.06l的线段上作用有合力f＝1n的均布载荷，其中l＝1。图2中的(b)为“l”形梁结构应力分布云图，在其内侧拐角处存在明显的应力集中现象，最大冯米塞斯应力高达187.3mpa。在该算例中，所选用的nurbs基函数阶次r和q均为2次，“l”形梁被划分为9800个跨度单元，共计10296个控制点，材料杨氏模量为1mpa，泊松比为0.25。图3中的(a)为包含体积约束的柔度最小化设计得到的拓扑结构，其中体积约束设置为30％，图3中的(b)为该拓扑结构对应的应力分布云图。由图可知，柔度优化设计不能减轻结构的应力水平，在“l”形梁内侧拐角处仍然存在应力集中现象，且最大冯米塞斯应力达到了293.2mpa。
[0129]
接下来采用本实施方式构建的考虑应力约束的等几何拓扑优化设计方法对“l”形梁进行优化设计，设置冯米塞斯应力约束值为70mpa，其他优化控制参数均与前述内容保持一致。图4中的(a)展示了考虑应力约束的体积最小化设计得到的拓扑结构，图4中的(b)为其对应的应力分布云图。由图4中的(a)可以看出，本实施方式的方法所得拓扑结构完全去除了“l”形梁内侧拐角处的材料，使之变为一个圆弧形倒角，规避了应力集中现象，结构最大冯米塞斯应力值为70.0mpa，精确控制在应力约束值上，同时结构体积比为30.4％，相较于初始“l”形梁减重约70％，实现了结构轻量化设计。图5给出了结构最大冯米塞斯应力和体积分数的收敛曲线，分别为结构最大冯米塞斯应力与应力约束值的比值及结构体积分数，可以看出在整个优化过程中收敛平稳，结构最大应力水平在优化过程初期快速降至应力约束值附近，且后期迭代中无明显突变现象。
[0130]
本发明还提供了一种等几何拓扑优化设计设备，所述设备包括存储器及处理器，所述存储器储存有计算机程序，所述处理器执行所述计算机程序时执行如上所述的考虑应力约束的等几何拓扑优化设计方法。
[0131]
本发明还提供了一种计算机可读存储介质，所述计算机可读存储介质存储有机器
可执行指令，所述机器可执行指令在被处理器调用和执行时，所述机器可执行指令促使所述处理器实现如上所述的考虑应力约束的等几何拓扑优化设计方法。
[0132]
综上，本发明构建的考虑应力约束的等几何拓扑优化设计方法通过采用改进的应力消失约束函数使结构应力水平可快速降低，提高了优化设计效率。同时本发明采用增广拉格朗日函数法消除了原优化模型中的大量应力约束，避免了求解困难，通过不断更新拉格朗日乘子最终获得了原优化模型的数值解。此外，通过采用稳定性较好的无约束mma算法及合理设置优化过程控制参数实现了优化设计的稳定收敛。
[0133]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。