通过密度方法进行拓扑优化

2019年 1月 4日

工程师在设计飞机和空间应用中的轻量化结构部件时有很大的自由度,因此使用能够开发自由度的方法很有意义。拓扑优化是在早期设计阶段普遍使用的方法。拓扑优化方法通常需要进行正则化和特殊的插值函数才能获得有意义的设计,这对于新手和有经验的仿真用户而言都比较困难。为了简化拓扑优化问题的解决方案,COMSOL® 软件提供了一种密度拓扑功能。

拓扑优化和密度方法

顾名思义,拓扑优化是一种能够针对给定的目标函数和约束条件为工程结构找出新的更好拓扑的方法。该方法通过引入一组设计变量来描述这些新拓扑,即描述设计空间中材料是否存在。这些变量被定义在网格的每个单元内或网格的每个节点上。因此,更改这些设计变量类似于更改拓扑。这意味着结构中的孔可以出现、消失和合并,并且边界可以采用任意形状。此外,控制参数在某种程度上是自动定义的,并与离散化相关。

从 COMSOL Multiphysics® 软件 5.4 版本开始,附加的优化模块具有密度拓扑功能,可以提高拓扑优化的易用性。该功能作为密度方法使用(参考文献 3),这意味着控制参数可以通过插值函数更改材料参数。固体和流体力学的插值函数已经内置到该功能中,并应用在 COMSOL Multiphysics 案例库的所有示例模型中。

支架几何的拓扑优化对比图
支架的几何形状经过拓扑优化后,仅保留了 50% 的材料,这些材料对刚度的贡献最大。

打印的支架几何形状。
打印的支架几何形状

密度方法包括定义控制变量场 \theta_c 的范围在 0~ 1。在固体力学中,\theta_c=1 对应于构建结构的材料,而 \theta_c=0 对应于非常柔软的材料。默认情况下,空域杨氏模量为固体杨氏模量的 0.1%。在流体力学中,惯例规定 \theta_c=1 对应于流体,而 \theta_c=0 是具有反渗透系数 \alpha 的(微)渗透材料,即将阻尼项添加到纳维-斯托克斯方程(Navier-Stokes equation)中:

\rho \frac{D\mathbf{v}}
{Dt}
= -\mathbf
{\nabla}p+\eta\nabla^2\mathbf{v}-\alpha(\theta_c)\mathbf{v}

 
在流体域中,阻尼项为 0,而在固体域中,则使用较大的值。这些不同的值近似给出了不同域之间界面上的无滑移边界条件。

密度模型功能简介

密度模型 功能通过亥姆霍兹方程(Helmholtz equation)支持正则化(参考文献 1),并使用过滤器半径引入最小长度比例:

\theta_f = R_\mathrm{min}^2\mathbf{\nabla}
^2\theta_f + \theta_c.

 
式中,\theta_c 是经优化程序修正的原始控制变量,\theta_f 是过滤后的变量。网格边尺寸是过滤器半径的默认值。尽管这在优化问题正则化方面效果很好,但为了获得独立于网格的结果,设置一个固定长度(大于网格边尺寸)非常重要。

亥姆霍兹过滤器的解析方程原理图
 
MBB 梁的原始控制变量和拓扑优化过滤后的控制变量
上图:亥姆霍兹过滤器方程可被解析求解为一维 Heaviside 函数。下图:此图来自 MBB 梁的拓扑优化教程模型。图左侧为原始控制变量,右侧为过滤后的控制变量。

亥姆霍兹过滤会产生明显的灰度,该过程没有清晰的物理解释。通过施加一个平滑的阶梯函数可以降低灰度,这个函数在拓扑优化中被称为投影。投影降低了灰度,但也使优化器收敛变得更加困难。密度拓扑函数支持基于双曲正切函数的投影,并且可以通过投 \beta 来控制投影量。

\theta = \frac{\tanh(\beta(\theta_f-\theta_
{\beta}))+\tanh(\beta\theta_{beta}
)}{\tanh(\beta(1-\theta_
{\beta}))+\tanh(\beta\theta_{beta}
)}

 
式中,\theta_{\beta} 是投影点。

过滤区域与投影区域的对比图
图示左侧为已过滤区域,右侧为已投影区域。

投影使避免灰度成为可能,但是如果优化问题对其有利,灰度仍然可以出现。如果将相同的插值函数用于质量和刚度,则灰度在体积受限的最小柔度问题中是最佳的。因此,通常使用插值函数,这些函数会产生相对于其成本的几乎没有刚度的中间值(与完全固体的值相比)。我们可以将其视为支持以下两种固体力学插值方案:固体各向同性材料罚函数(SIMP)和材料属性方法的有理逼近(RAMP)插值的材料体积因子和密度模型 中间值的 接口(如下所示)。达西插值是为流体力学提供的。插值变量称为罚材料体积因子 \theta_p ,用于插值材料参数,例如,对于 SIMP 插值,p_\textsc{simp} 指数可以增加以减小中间值的刚度,从而降低灰度。

\begin{align}
\theta_p %26= \theta_\mathrm{min}(1-\theta_\mathrm{min})\theta^{p_\textsc{simp}}\\
E_p %26= E\theta_p
\end{align}

 
式中,E 是固体材料的杨氏模量,E_p 是所有优化域中使用的罚杨氏模量。

COMSOL 中的密度模型功能界面
“密度模型”功能位于组件 > 定义的拓扑优化下。网格边的长度被用作默认过滤器半径,并且效果很好,但是必须替换为固定值才能产生与网格无关的结果。

罚杨氏模量可以定义为域变量,或者(如在支架模型案例中)可以直接在材料中定义。

屏幕截图显示了杨氏模量的材料属性
用密度方法进行拓扑优化包括在空间上改变杨氏模量。这是通过在材料属性中将固体杨氏模量乘以罚材料体积因子 dtopo1.theta_p 来实现的。

总之,密度拓扑函数添加了四个变量。过滤的材料体积因子是通过因变量隐式定义的。

符号 描述 方程
\theta_c 边界控制材料体积因子 0\leq\theta_c\leq1
\theta_f 过滤的材料体积因子 \theta_f = R_\mathrm{min}^2\mathbf{\nabla}^2\theta_f + \theta_c
\theta 材料体积因子 \theta = \frac{\tanh(\beta(\theta_f-\theta_{beta}
))+\tanh(\beta\theta_

{\beta})}{\tanh(\beta(1-\theta_{beta}
))+\tanh(\beta\theta_

{\beta}
)}

\theta_p 罚材料体积因子 \theta_p = \theta_\mathrm{min}+(1-\theta_\mathrm{min}
)\theta^{p_\textsc{simp}}
 或者 \theta_p = \frac{q_\mathrm{Darcy}(1-\theta)}{q_\mathrm{Darcy}
+\theta}

禁用过滤器后,过滤后的变量将变得不确定,投影将直接使用控制的材料体积因子。如果禁用投影,则材料体积因子仍然存在,但与投影输入相同。

应用连续性以避免局部极小值

当拓扑不太复杂时,密度拓扑功能的默认值效果很好。例如,MBB 梁的拓扑优化导出和导入拓扑优化的钢钩模型。如果最佳设计更加复杂(例如本文开头的支架示例),则可能会有很多局部最小值。为了避免这些最小值,可以在 SIMP 指数和投影斜率中使用连续性。这可以通过修改密度拓扑特征中的初始值表达式并添加参数化扫描 特征来实现如下图所示。因此,求解器使用之前案例中的最佳值作为下一个优化步骤的初始值,使指定的参数变斜。也就是说,它以较小的 SIMP 指数和投影斜率开始,然后连续到更高的值。

带有研究参考的参数化扫描设置的屏幕截图
将参数化扫描与参考研究相结合应用连续性。有关详细信息,请参见“支架-拓扑优化教程模型。

拓扑优化中的目标和约束

如果针对单一载荷工况优化了几何形状(如下左图所示),则相对于该载荷工况,最终的设计将是最佳的。这似乎很明显,但是设计人员通常会对对称性和设计拓扑进行假设。除非将这些假设用公式表示为约束条件,否则将不会遵守这些假设。因此,下面右图所示的设计使用了八个荷载工况(两个荷载组乘以四个约束组)。

单一载荷工况下,支架的优化几何形状
具有八个载荷工况的支架几何形状

左:单一载荷工况下,支架的几何形状进行了优化,从而形成了两个半松散连接的不对称设计。右:具有八个载荷工况的支架几何形状。

通常,设计人员通常有几个目标需要权衡。为了对这些目标做出明智的决定,设计人员可以使用几种具有不同权重的优化方法来跟踪帕累托最优前沿(Pareto optimal front)。

支架几何形状的帕累托最优前沿绘图
可以通过在参数化扫描中更改权重来跟踪支架几何形状的帕累托最优前沿。

 

支架-拓扑优化的动画。(从 COMSOL 案例库中下载 glTF™ 文件的 GLB 文件格式,自己旋转几何。)

导出和导入拓扑优化结果

我们无需重新剖分网格就可以对应力集中和屈曲进行分析来考查拓扑优化设计的结果。但是,如果要完全确定空域相不起作用,可以通过导出和导入结果设计来去除它,如下图所示。有关此过程的详细信息,请参阅上一篇博客文章

MBB梁拓扑优化设计的轮廓和插值曲线
将 MBB 梁拓扑优化设计的轮廓(左)作为插值曲线导出和导入(右)。

后续操作

了解有关解决优化问题的 COMSOL 内置工具和功能的更多信息,请单击下面的按钮,查看优化模块产品页面。

更多资源

参考文献

  1. B.S. Lazarov and O. Sigmund, “Filters in topology optimization based on Helmholtz‐type differential equations,” International Journal for Numerical Methods in Engineering, vol. 86, no. 6, pp. 765–781, 2011.
  2. F. Wang, B.S. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations in topology optimization,” Structural and Multidisciplinary Optimization, vol. 43, pp. 767–784, 2011.
  3. M.P. Bendsøe, “Optimal shape design as a material distribution problem,” Structural Optimization, vol. 1, pp. 193–202, 1989.

glTF 和 glTF 徽标是 Khronos Group Inc. 的商标。

博客分类


评论 (11)

正在加载...
泽昊 李
泽昊 李
2023-03-06

请问如何设置随迭代次数变化的投影斜率

Kristian Ejlebjærg Jensen
Kristian Ejlebjærg Jensen
2023-03-06 COMSOL 员工

The bracket_topology_optimization_stl model demonstrates this.

泽昊 李
泽昊 李
2023-03-06

参数化扫描的投影斜率并不是随着迭代次数连续变化的,是我理解错您的意思了吗?

Kristian Ejlebjærg Jensen
Kristian Ejlebjærg Jensen
2023-03-06 COMSOL 员工

No, four values are used as shown in the UI picture of “Applying Continuation to Avoid Local Minima”, but you can use more values, if you prefer.

泽昊 李
泽昊 李
2023-03-07

谢谢!

明 小
明 小
2023-03-08

请问网格边尺寸是用户定义栏默认的h吗?能否知道它的具体参数?否的话,如何确定一个固定长度(大于网格边尺寸)

屹磊 金
屹磊 金
2023-03-20 COMSOL 员工

密度模型中设置的过滤半径“h”是一个全局定义下的参数,在后面的网格剖分中也使用到了该参数。具体的数值您可以该案例中进行查看:支架 – 拓扑优化:https://cn.comsol.com/model/bracket-8212-topology-optimization-69891

孟杰 刘
孟杰 刘
2023-04-02

您好,我在对液流电池流场板进行拓扑优化设计时,需要求解二次电流分布、稀物质传递、自由与多孔介质流动等物理场,用到了三个研究步骤,但发现基于梯度的优化求解器(如MMA等)仅能优化单研究步骤,请问有什么方法可以实现多研究步骤的优化呢?

Kristian Ejlebjærg Jensen
Kristian Ejlebjærg Jensen
2023-04-03 COMSOL 员工

你好,

COMSOL 6.1 仅支持对单个研究进行基于梯度的优化,但如果物理场接口都依赖于稳态求解器,则它们可以通过单个研究步骤求解。 可以使用分离求解器来控制物理场求解的顺序(这与单向耦合特别相关)。

此致,
克里斯蒂安·詹森
技术产品经理,优化

Tom Steven
Tom Steven
2023-11-09

你好,请问我在进行电磁类超表面(吸波器)拓扑优化时,优化器常常迭代几次就结束了,得到的结果并不理想,材料的的密度并没有趋近于1或0,而是处于中间值,不知道为什么会这样,请问有相关的拓扑优化案例推荐吗?

Kristian Ejlebjærg Jensen
Kristian Ejlebjærg Jensen
2023-11-09 COMSOL 员工

Hi Tom Steven

COMSOL does not have examples of topology optimization for electromagnetic wave propagation, but the provided functionality is similar to what is used in academic literature. Problems with grayscale can be due to bad interpolation and/or the absence of a volume constraint. I suggest studying the academic literature on the subject.

Best regards,

Kristian E. Jensen

Technical Product Manager, Optimization

浏览 COMSOL 博客