CN111079279A - 一种多构型晶格结构多尺度拓扑优化设计方法 - Google Patents

Abstract

本发明属于结构拓扑优化设计领域，并具体公开了一种多构型晶格结构多尺度拓扑优化设计方法。包括：将设计区域离散成多个晶格单元，以各晶格单元的总材料占比ρ和该晶格单元中各构型材料百分比ratio为变量，构建设计变量，并以此构建设计区域的设计变量组；构建各晶格单元的总材料占比ρ、各晶格单元总材料中各构型材料的百分比以及设计变量组的刚度矩阵

的一一对应关系模型；然后建立多构型晶格结构多尺度拓扑优化模型，并对宏观尺度及微观尺度中的设计变量组进行灵敏度分析，并迭代更新宏观尺度及微观尺度中的设计变量组。本发明实现了晶格结构多尺度拓扑优化设计的可制造性，相比于单一构型的晶格结构，多构型晶格结构在力学性能上更为优异。

Description

一种多构型晶格结构多尺度拓扑优化设计方法

技术领域

本发明属于结构拓扑优化设计领域，更具体地，涉及一种多构型晶格结构多尺度拓扑优化设计方法。

背景技术

拓扑优化是通过数学建模与优化算法，在给定的设计域及一定的边界条件下优化材料分布，以实现结构性能的提升，这种方法能够通过算法自动生成结构拓扑构型(如孔洞的数量，位置及连通性)，摆脱了经验式设计的弊端，能够充分发挥材料与结构性能。面向晶格结构设计与增材制造更紧密的结合，Cheng等基于均匀化方法提出一种面向增材制造的变密度的多孔结构设计方法。该方法根据密度分布以及增材制造工艺重构多孔结构，并利用3D打印制造出来并进行力学模拟与测试。Tang等针对晶格结构提出了一种带有3D打印制造约束的设计及优化方法。其中可制造单元的概念被引入来连接设计及制造，并基于实验与人工神经网络建立了制造约束模型，最后根据制造约束利用BESO方法优化晶格结构分支的尺寸来实现优化设计与增材制造的结合。Altair公司旗下的Opti-struct在2015年推出了结合拓扑优化设计的晶格结构填充设计功能模块。该功能模块基于传统密度变量法得到的拓扑优化设计构型，对非关键承载区域进行晶格结构进一步填充设计。相比传统宏观结构拓扑构型优化设计，该宏微多级晶格结构设计在实现构件进一步减重的同时，其综合承载性能也得到一定的提升。

但上述技术中所采用的设计模式中存在以下的问题：1)填充式的、经验式的晶格结构设计宏观构型与晶格单元构型设计分离，无法充分发挥晶格结构性能；2)基于均匀化方法的晶格结构设计方法其晶格单元之间连接性不足，并且宏微观尺度间无显式的比例系数从而导致其设计结果制造性差；3)晶格结构设计方法中晶格结构单元较为单一，缺乏有效的单元构型控制和多构型的晶格结构设计方法。

因此，本领域亟待提出一种多构型晶格结构多尺度拓扑优化设计方法，从而使得设计过程中最大化考虑晶格结构性能，提高设计结果的可制造性。

发明内容

针对现有技术的以上缺陷或改进需求，本发明提供了一种多构型晶格结构多尺度拓扑优化设计方法，其中结合多构型晶格结构自身的特征及其多尺度拓扑优化设计特点，通过密度变量ρ和比例变量ratio对晶格单元的参数化，能够实现晶格结构构型在不同材料用量下的连续变化，并基于设计变量组建立了密度变量，比例变量和晶格单元刚度矩阵三者之间的关系，实现自变量组(ρ,ratio)与晶格单元构型及其刚度矩阵的一一对应，并以此建立多构型晶格结构多尺度拓扑优化模型，相比于单一构型的晶格结构，多构型晶格结构在力学性能上更为优异。因而本方法适于用来进行增材制造晶格结构拓扑优化设计的理论研究和实际制造。

为实现上述目的，本发明提出了一种多构型晶格结构多尺度拓扑优化设计方法，包括以下步骤：

S1将设计区域离散成多个晶格单元，以各晶格单元的总材料占比ρ和该晶格单元中各构型材料百分比ratio为变量，构建设计变量x＝(ρ,ratio)，并以此构建设计区域的设计变量组；

S2构建各晶格单元的总材料占比ρ、各晶格单元总材料中各构型材料的百分比以及设计变量组的刚度矩阵

的一一对应关系模型；

S3根据步骤S2中得到的一一对应关系模型，建立多构型晶格结构多尺度拓扑优化模型，然后对宏观尺度及微观尺度中的设计变量组进行灵敏度分析，并迭代更新宏观尺度及微观尺度中的设计变量组，从而确定多构型晶格结构多尺度拓扑优化设计的最优结果。

作为进一步优选的，步骤S1中，将构成晶格单元的结构划分为叉形结构和十字结构，所述叉形结构和十字结构的总材料占比为ρ；所述叉形结构的材料百分比ratio为：

其中，γ叉形结构材料占晶格单元的总材料占比ρ的材料占比。

作为进一步优选的，步骤S2具体包括以下步骤：根据所述设计变量组计算该变量组的刚度矩阵

对该刚度矩阵

进行本征特征分解，以提取不同密度下

的特征基底，采用子结构降阶近似模型对该特征基底的本征正交分解系数进行参数拟合，进而构建各晶格单元的总材料占比ρ、各晶格单元总材料中各构型材料的百分比以及设计变量组的刚度矩阵

的一一对应关系模型。

作为进一步优选的，步骤S2具体包括以下步骤：

S21采用静态压缩模型对刚度矩阵

进行静态压缩，该静态压缩模型为：

其中，K_bb为晶格单元边界节点关于边界节点的刚度矩阵，K_bi为晶格单元内部节点关于边界节点的刚度矩阵，K_ii为晶格单元内部节点关于内部节点的刚度矩阵，K_ib为晶格单元边界节点关于内部节点的刚度矩阵。将边界节点自由度与内部节点自由度分离表示，其有限元方程即可化为：

U_b、F_b分别表示边界节点处的位移向量和力向量；U_i、F_i分别表示内部节点处的位移向量和力向量。

S22对静态压缩后的刚度矩阵

进行本征特征分解，以提取不同密度下刚度矩阵

的特征基底，将静态压缩后的刚度矩阵

转化为：

其中，α₁,……,α_m为刚度矩阵

的特征基底的本征正交分解系数，[φ₁],……,[φ_m]为不同密度下刚度矩阵

的特征基底；

S23采用子结构降阶近似模型对刚度矩阵

的特征基底的本征正交分解系数进行二维参数拟合，进而构建各晶格单元的总材料占比ρ、各晶格单元总材料中各构型材料的百分比以及设计变量组的刚度矩阵

的一一对应关系模型，所述二维参数拟合的计算模型为：

其中，w_i为权重系数，q为多项式基向量，m为保留特征基底向量的个数，a＝[a₁,a₂,...]^T为插值误差最小的系数向量，α(ρ_i)为第i个晶格单元与密度相关的单元系数。

作为进一步优选的，所述一一对应关系模型为：

其中，

为本征正交分解拟合系数，m为保留特征基底向量的个数，

作为进一步优选的，步骤S3中，所述多构型晶格结构多尺度拓扑优化模型为：

其中，K表示整体刚度矩阵，U表示整体位移矩阵，F表示整体力向量，V(ρ)表示材料体积，V₀表示设计域体积，

表示预设的体积分数，ρ_min表示最小相对密度值，以避免刚度矩阵的奇异性，U_i为晶格单元的位移矩阵，

为U_i的转置。

作为进一步优选的，步骤S3中，还可采用全局收敛性GCMMA算法迭代更新宏观尺度及微观尺度中的设计变量组，从而确定多构型晶格结构多尺度拓扑优化设计的最优结果。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，主要具备以下的技术优点：

1.本发明通过密度变量ρ和比例变量ratio对晶格单元的参数化，能够实现晶格结构构型在不同材料用量下的连续变化，并基于设计变量组建立了密度变量，比例变量和晶格单元刚度矩阵三者之间的关系，实现自变量组(ρ,ratio)与晶格单元构型及其刚度矩阵的一一对应，并以此建立多构型晶格结构多尺度拓扑优化模型。

2.本发明通过密度变量ρ和比例变量ratio对晶格单元的参数化，能够实现晶格结构构型在不同材料用量下的连续变化，同时也实现了自变量与晶格构型的耦合，每一组密度变量和比例变量(ρ,ratio)能够唯一确定一个晶格单元的构型，为后续的多尺度设计奠定了基础。

3.本发明在构建双自变量且具有耦合关系的自变量更新迭代方法时，常用的OC算法在处理多种自变量问题时需要分别推导两种自变量的优化准则，过程较为复杂，难以实现准确高效的求解，因而采用具有全局收敛性的MMA算法即GCMMA算法进行双自变量的更新迭代。

附图说明

图1是本发明实施例涉及的一种多构型晶格结构多尺度拓扑优化设计方法的流程图；

图2是本发明实施例中晶格单元构型参数控制定义示意图；

图3是本发明实施例中密度变量ρ和比例变量ratio联合控制的多构型晶格单元示意图；

图4中的(a)～(f)是本发明实施例涉及的多构型晶格单元本征正交分解系数插值结果示意图。

具体实施方式

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

如图1所示，本发明一种多构型晶格结构多尺度拓扑优化设计方法包括以下步骤：定义一种多构型晶格单元，通过密度变量ρ和比例变量ratio对晶格单元的参数化，能够实现晶格结构构型在不同材料用量下的连续变化；利用上述晶格单元形成采样的设计变量组，并计算该变量组的刚度矩阵

对计算好的刚度矩阵进行本征特征分解提取不同密度下

的特征基底，根据子结构降阶近似模型，对本征正交分解系数进行二维参数拟合；建立了密度变量，比例变量和晶格单元刚度矩阵三者之间的关系，实现自变量组(ρ,ratio)与晶格单元构型及其刚度矩阵的一一对应。根据变密度法与SIMP模型，建立多构型晶格结构多尺度拓扑优化模型。将设计区域离散成N_x×N_y个子结构，每个子结构对应两个设计变量分别为密度变量ρ与比例变量ratio，设计变量通过上述本征正交分解与二维的参数拟合与子结构超单元刚度矩阵一一对应，建立多构型晶格结构多尺度拓扑优化模型；根据材料体积V对于比例变量ratio的灵敏度公式

并使用GCMMA算法求解问题模型。

按照本发明的多构型多尺度拓扑优化设计方法实现了晶格结构多尺度拓扑优化设计的可制造性，相比于单一构型的晶格结构，多构型晶格结构在力学性能上更为优异。因此，本方法适于用来进行增材制造晶格结构拓扑优化设计的理论研究和实际制造。

具体而言，即首先，将设计区域离散成多个晶格单元，以各晶格单元的总材料占比ρ和该晶格单元中各构型材料百分比ratio为变量，构建设计变量x＝(ρ,ratio)，并以此构建设计区域的设计变量组。

如图2和图3所示，将晶格单元分解成叉形结构与十字形结构，晶格单元密度为ρ，晶格单元中，叉形结构占晶格单元总材料的百分比为：

以此方式，通过密度变量ρ和比例变量ratio将晶格单元在在两个维度的连续性变化。对参数化处理的晶格单元子结构进行有限元网格划分，然后对子结构密度变量ρ和比例变量ratio分别进行均匀采样，然后将其按照以下方式形成采样的设计变量组：

然后利用公式

对子结构刚度矩阵进行静态压缩，并计算降阶的刚度矩阵

其中，K_bb为晶格单元边界节点关于边界节点的刚度矩阵，K_bi为晶格单元内部节点关于边界节点的刚度矩阵，K_ii为晶格单元内部节点关于内部节点的刚度矩阵，K_ib为晶格单元边界节点关于内部节点的刚度矩阵。

如图4所示，其中，图4中的(a)～(f)的纵坐标均为α，横坐标为ratio。本发明中，针对密度变量ρ采样个数为21，密度变量ratio采样个数为8，则通过计算能够得到168个不同密度子结构超单元的刚度矩阵

根据上文所建立的方法，将

向量化并组成采样刚度矩阵集合

对

进行本征正交分解即可提取不同密度下

的特征基底，将其转化为以下形式：

其中m＝168，由于引入了新的比例变量控制晶格单元构型，因此

和α₁(ρ,ratio),...,(ρ,ratio)均与密度变量ρ和比例变量ratio相关。

根据子结构降阶近似模型，对本征正交分解系数进行参数拟合，采用基于移动最小二乘的拟合方法，通过移动最小二乘的方法拟合系数α_k与自变量ρ的函数关系：

其中q＝[q₁,q₂,...]^T为多项式基向量，a＝[a₁,a₂,...]^T为插值误差最小的系数向量。对于二维的曲面拟合，多项式基向量以自变量ρ和r的形式表示：

线性基：q＝[1,ρ,r]^T

二次基：q＝[1,ρ,r,ρ²,ρr,r²]^T (5)

其系数向量a＝[a₁,a₂,...]^T是使得插值误差最小的值：

其中，w_i为权重系数，q为多项式基向量，m为保留特征基底向量的个数，a＝[a₁,a₂,...]^T为插值误差最小的系数向量，α(ρ_i)为第i个晶格单元与密度相关的单元系数。

其中参数α₁,α₅,α₁₂,α₂₀,α₃₀,α₄₀的对设计变量ρ和ratio的二维拟合曲面如图3所示。

通过系数的二维拟合，建立了密度变量，比例变量和晶格单元刚度矩阵三者之间的关系，实现了自变量组(ρ,ratio)与晶格单元构型及其刚度矩阵的一一对应。

根据变密度法与SIMP模型，建立多构型晶格结构多尺度拓扑优化模型。将设计区域离散成N_x×N_y个子结构，每个子结构对应两个设计变量分别为密度变量ρ与比例变量ratio，其中密度变量ρ表示每个子结构的密度值，比例变量ratio和密度变量ρ共同决定子结构构型：

其中ρ∈[ρ_min,1]，ratio∈[0,1]。

另外，设计变量通过上文中的本征正交分解与二维的参数拟合与子结构超单元刚度矩阵关联起来：

其中：[φ_k]不同密度下刚度矩阵

的特征基底，

为本征正交分解拟合系数，m为保留特征基底向量的个数。

根据以上推导，根据变密度法与SIMP模型，建立多构型晶格结构多尺度拓扑优化模型，其中目标函数为最小化柔度(即最大化刚度)，约束为体积分数：

其中，K表示整体刚度矩阵，U表示整体位移矩阵，F表示整体力向量，V(ρ)表示材料体积，V₀表示设计域体积，

表示预设的体积分数，ρ_min表示最小相对密度值，以避免刚度矩阵的奇异性，U_i为晶格单元的位移矩阵，

为U_i的转置。

为了求解已建立的多构型晶格结构多尺度拓扑优化数学模型，最小化结构柔顺度，需要分别针对密度变量ρ和比例变量ratio进行敏度分析。

灵敏度可以表示为：

其中：x表示设计变量。分别对密度变量ρ和比例变量ratio进行灵敏度分析。

对于密度变量ρ，式10可转化为：

根据刚度矩阵的表达式3，对密度变量ρ求导：

最终得到目标函数对于密度变量的灵敏度：

其中

表示拟合系数，[φ_k]表示矩阵形式的本征正交分解基底向量。由于对于系数α的拟合是采用三次样条插值的方式，所以

的值可以直接通过前文所得的参数拟合函数直接求导得到。

根据材料体积表达式可以得出材料体积V对于密度变量ρ的灵敏度：

其中：v_i表示单元体积。

另外，对于比例变量ratio，式4.10可转化为：

根据刚度矩阵的表达式4.3，对比例变量ratio求导：

最终得到目标函数对于比例变量的灵敏度：

其中：

表示拟合系数，[φ_k]表示矩阵形式的本征正交分解基底向量。同样，由于对于系数α的拟合是采用三次样条插值的方式，

的值也可以直接通过前文所得的参数拟合函数直接求导得到。

根据材料体积

表达式可以得出材料体积V对于比例变量ratio的灵敏度：

在构建双自变量且具有耦合关系的自变量更新迭代方法时，常用的OC算法在处理多种自变量问题时需要分别推导两种自变量的优化准则，过程较为复杂，难以实现准确高效的求解。本发明采用具有全局收敛性的MMA算法即GCMMA算法进行双自变量的更新迭代。

GCMMA算法中，一般非线性优化问题可表示为：

其中：x＝(x₁,...,x_n)^T为自变量，y＝(y₁,...,y_n)^T，z∈R为附加变量，

为实数，且

f₀,...,f_m为连续可微实函数。a₀,a_i,c_i,d_i均为实数，且有a₀＞0,a_i,c_i,d_i≥0，另外c_i+d_i＞0，a_ic_i＞a₀。

对于拓扑优化问题，对以上一般非线性优化问题做相应处理：a_i＝0，d_i＝0，c_i＝0即可得到结构拓扑优化标准形式：

即结构拓扑优化问题可以看作为一般非线性优化问题的特殊形式。

GCMMA算法解是否存在的判定依据在于由目标函数与约束条件组成Lagrange函数，根据KKT(Karush-Kuhn-Tucker)条件得到优化问题解的存在性判定。因此建立式19问题的Lagrange函数：

其中，λ＝(λ₁，...，λ_m)^T，ξ＝(ξ₁，...，ξ_n)^T，η＝(η₁，...，η_m)^T，μ＝(μ₁，...，μ_n)^T，ζ为非负的拉格朗日乘子。

式21的KK条件为：

如果存在一般非线性优化问题19的最优解，则存在满足上述KKT条件的Lagrange乘子λ,ξ,η,μ,ζ以及松弛因子s。

GCMMA算法加入外迭代与内迭代以保证全局收敛。其自变量迭代过程为：首先选定初始值(x⁽¹⁾,y⁽¹⁾,z⁽¹⁾)，然后根据内外迭代规则依次进行设计变量更新。假设第k次外迭代时自变量为(x^(k),y^(k),z^(k))，在本次外迭代中建立一个近似的子问题并求解，在这个子问题中利用特定凸函数

来替代原函数f_i(x)，假设此子问题的的最优解为

如果

则下一步迭代自变量

本次外迭代结束，无需进行内迭代。否则由外迭代进入内迭代，在内迭代中建立新的子问题，利用更为保守的

进行原函数的替代，假设子问题最优解为

如果

则下一步迭代自变量

此外迭代结束(仅需一步内迭代)，否则，按照重复内迭代直至满足

υ表示内迭代步数。

在内外迭代中，移动近似子问题的产生如下：

其中,

近似方程

定义如下：

其中：

其中，默认下渐近线

与上渐近线

的更新规则如下：

当外迭代步数k＝1或2时：

在之后的外迭代中，即k≥3时

其中：

式(24)中参数

严格为正，其更新规则如下：

在每一次外迭代开始时，内迭代为0即υ＝0，其初始值定义为：

在进行外迭代时，

在进行内迭代时，

其中：

在求解本章建立的多构型晶格结构多尺度拓扑优化模型时，定义设计变量x＝(ρ,ratio)^T即可进行双自变量的优化求解。

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

Claims (7)

1.一种多构型晶格结构多尺度拓扑优化设计方法，其特征在于，包括以下步骤：

S1将设计区域离散成多个晶格单元，以各晶格单元的总材料占比ρ和该晶格单元中各构型材料百分比ratio为变量，构建设计变量x＝(ρ,ratio)，并以此构建设计区域的设计变量组；

S2构建各晶格单元的总材料占比ρ、各晶格单元总材料中各构型材料的百分比以及设计变量组的刚度矩阵

的一一对应关系模型；

S3根据步骤S2中得到的一一对应关系模型，建立多构型晶格结构多尺度拓扑优化模型，然后对宏观尺度及微观尺度中的设计变量组进行灵敏度分析，并迭代更新宏观尺度及微观尺度中的设计变量组，从而确定多构型晶格结构多尺度拓扑优化设计的最优结果。

2.根据权利要求1所述的一种多构型晶格结构多尺度拓扑优化设计方法，其特征在于，步骤S1中，将构成晶格单元的结构划分为叉形结构和十字结构，所述叉形结构和十字结构的总材料占比为ρ；所述叉形结构的材料百分比ratio为：

其中，γ叉形结构材料占晶格单元的总材料占比ρ的材料占比。

3.根据权利要求1所述的一种多构型晶格结构多尺度拓扑优化设计方法，其特征在于，步骤S2具体包括以下步骤：根据所述设计变量组计算该变量组的刚度矩阵

对该刚度矩阵

进行本征特征分解，以提取不同密度下

的特征基底，采用子结构降阶近似模型对该特征基底的本征正交分解系数进行参数拟合，进而构建各晶格单元的总材料占比ρ、各晶格单元总材料中各构型材料的百分比以及设计变量组的刚度矩阵

的一一对应关系模型。

4.根据权利要求1所述的一种多构型晶格结构多尺度拓扑优化设计方法，其特征在于，步骤S2具体包括以下步骤：

S21采用静态压缩模型对刚度矩阵

进行静态压缩，该静态压缩模型为：

其中，K_bb为晶格单元边界节点关于边界节点的刚度矩阵，K_bi为晶格单元内部节点关于边界节点的刚度矩阵，K_ii为晶格单元内部节点关于内部节点的刚度矩阵，K_ib为晶格单元边界节点关于内部节点的刚度矩阵；

S22对静态压缩后的刚度矩阵

进行本征特征分解，以提取不同密度下刚度矩阵

的特征基底，将静态压缩后的刚度矩阵

转化为：

其中，α₁,……,α_m为刚度矩阵

的特征基底的本征正交分解系数，[φ₁],……,[φ_m]为不同密度下刚度矩阵

的特征基底；

S23采用子结构降阶近似模型对刚度矩阵

的特征基底的本征正交分解系数进行二维参数拟合，进而构建各晶格单元的总材料占比ρ、各晶格单元总材料中各构型材料的百分比以及设计变量组的刚度矩阵

的一一对应关系模型，所述二维参数拟合的计算模型为：

其中，w_i为权重系数，q为多项式基向量，m为保留特征基底向量的个数，a＝[a₁,a₂,...]^T为插值误差最小的系数向量，α(ρ_i)为第i个晶格单元与密度相关的单元系数。

5.根据权利要求3所述的一种多构型晶格结构多尺度拓扑优化设计方法，其特征在于，所述一一对应关系模型为：

其中，

为本征正交分解拟合系数，m为保留特征基底向量的个数。

6.根据权利要求3所述的一种多构型晶格结构多尺度拓扑优化设计方法，其特征在于，步骤S3中，所述多构型晶格结构多尺度拓扑优化模型为：

find:ρ＝[ρ₁ ρ₂…ρ_N]

:ratio＝[ratio₁,ratio₂,...,ratio_N]

min:

subject to:KU＝F

:V(ρ)/V₀＝θ

:0<ρ_min≤ρ≤1

:0≤ratio≤1

其中，K表示整体刚度矩阵，U表示整体位移矩阵，F表示整体力向量，V(ρ)表示材料体积，V₀表示设计域体积，θ表示预设的体积分数，ρ_min表示最小相对密度值，以避免刚度矩阵的奇异性，U_i为晶格单元的位移矩阵，

为U_i的转置。

7.根据权利要求3所述的一种多构型晶格结构多尺度拓扑优化设计方法，其特征在于，步骤S3中，还可采用全局收敛性GCMMA算法迭代更新宏观尺度及微观尺度中的设计变量组，从而确定多构型晶格结构多尺度拓扑优化设计的最优结果。

Images