CN116579151B - 一种基于mmc框架的非均匀点阵结构优化设计方法 - Google Patents

Abstract

本发明公开了一种基于MMC框架的非均匀点阵结构优化设计方法，采用周期性点阵结构作为初始设计，利用在设计域内不同子分区上分别定义的局部坐标摄动函数CPFs在设计域内生成不同的微结构拓扑形式，并同时确保相邻微结构之间的连接，使用B‑spline基函数作为CPFs的基函数，并在参数空间中采用控制点网络对每个分区进行局部控制，在给定有限的材料用量下，优化初始点阵结构，使结构柔度最小化，从而产生一个具有平滑连接和清晰边界的渐变点阵结构。本发明采用上述步骤的一种基于MMC框架的非均匀点阵结构优化设计方法，可以在各种载荷条件下，在更大的设计空间中获得具有复杂拓扑构型的渐变点阵结构。

Description

一种基于MMC框架的非均匀点阵结构优化设计方法

技术领域

本发明涉及结构拓扑优化技术领域，尤其是涉及一种基于MMC框架的非均匀点阵结构优化设计方法。

背景技术

点阵结构通常是由若干具有不同几何拓扑结构的单胞在空间中连续排布而组成的多孔结构。具有该结构的材料，由于其轻质和高强度的特性使其在例如航空航天、汽车、建筑、生物医学、超材料等领域中得到了十分广泛的应用。此外，点阵结构材料作为一种多功能材料，在解决结构局部抗弯曲、拓扑绝缘、能量吸收、热膨胀等方面的问题中也有着巨大的潜力，因此这引起了相关领域的学者们的极大兴趣。近年来，增材制造(AdditiveManufacturing,AM)技术的进步不仅消除了传统工艺中许多与制造技术有关的问题，也使得这类具有复杂三维几何形状的结构能够更容易地被制造。基于增材制造技术超高的灵活性和显著的效率优势，人们提出了一个新的范式，即面向增材制造技术的设计方法(Designfor Additive Manufacturing,DFAM)，以此能够更加有效地利用增材制造技术的潜力，并在各类应用中充分实现所设计产品卓越的多功能性能。

与传统的由周期性排列的单胞构成的点阵结构相比，由具有空间变化特性的非均匀微结构构成的点阵结构通常表现出更好的性能。除了受到自然界中天然的非均匀点阵结构的启发外，人们对使用先进的数学工具和基于梯度的结构拓扑优化算法来设计创新的渐变点阵结构的兴趣也越来越大。作为一种成熟的结构设计方法，结构拓扑优化方法能够使得特定数量的材料分布在规定的设计域中，以获得具有某些特殊性能的最佳设计。这些有效的计算工具扩大了结构优化的设计空间，使得具有多尺度特点的点阵结构的优化成为可能。

在现有的设计方法中，均匀化方法通常被用于联系微结构单胞中的材料分布和其宏观等效属性。但由于该方法中存在尺度分离假设，而实际结构中的单胞尺寸并不可能做到无限小，因此在通过该方法设计的结构中，相邻的单胞之间会存在材料分布的连通性问题，且大量具有各异拓扑构型单胞的均匀化计算和优化过程也将带来巨大的计算成本。此外，该方法中，由于单胞尺寸与结构的宏观尺度之间缺乏明确的长度尺寸比，其设计的结果在进行增材制造等实际制造过程中也面临困难。虽然同时优化宏观结构和一个或几个有代表性的体积胞元可以有效地减少设计变量的数量，节省计算成本，但这种做法也在一定程度上限制了该方法的设计空间。

发明内容

本发明的目的是提供一种基于MMC框架的非均匀点阵结构优化设计方法，对用于描述设计域中材料分布的拓扑描述函数(topology description functions，TDFs)进行坐标摄动，通过优化摄动基函数中的系数，对整体结构进行优化，从而得到最优的渐变结构形式。

为实现上述目的，本发明提供了一种基于MMC框架的非均匀点阵结构优化设计方法，其特征在于，步骤如下：

S1、MMC显式拓扑优化：采用基于MMC框架的周期性点阵结构作为初始设计，通过组件的显式几何参数得到拓扑描述函数TDFs，用TDFs来表示实体材料在设计域中的分布；

S2、坐标分区映射：在设计域内的不同子分区上分别定义局部坐标摄动函数CPFs，从而在设计域内生成不同的微结构拓扑形式，并确保相邻微结构之间的连接；

S3、基于B-spline基函数的CPFs构造：通过B-spline基函数保持各个子分区之间CPFs的高阶连续性，并在参数空间中采用控制点网络对每个分区进行局部控制，形成基于显式描述框架的渐变点阵结构；

S4、数值实现：在满足体积约束的前提下，优化初始周期性点阵结构，使结构柔度最小化，从而产生一个具有平滑连接和清晰边界的渐变点阵结构。

优选的，S1中，在MMC框架中用TDFs来表示实体材料在设计域中的分布。可移动变形组件(Moving Morphable Components,MMC)利用具有显式几何参数的组件作为结构拓扑优化的基本设计基元，通过组件的移动、变形、交叠和融合来实现结构宏观拓扑的变化。拓扑描述函数TDFs的构建具体包括：

在欧拉描述的MMC框架中，结构的TDFs可以描述如下：

其中，D代表预先给定的设计域，是由n根组件所构成的实体材料占据的空间，φ^s(x)＝max(φ₁(x),…,φ_n(x)表示整个结构的TDFs；而其中的φ_i(x),(i＝1,2,…,n)表示第i根组件的TDF，并可按下式详细写为：

且有

其中，p是描述超椭圆的参数，指定为一个较大的偶数，在公式(1.2)和公式(1.3)中，符号(x_0i,y_0i)、a_i、f_i(x′)、θ_i分别表示第i根组件的中心点位置、半长、半宽和从水平轴逆时针测量的倾斜角，如图1所示。

描述组件厚度变化曲线的函数f_i(x′)按下式给出：

其中，和/>分别为第i根组件不同位置处的半宽，函数f_i(x′)是在局部坐标下定义的，基于以上组件几何描述，与第i根组件相关的设计变量可以写为其中/>包含了与f_i(x′)相关的设计参数；将整个结构的拓扑构型通过若干设计变量来描述，并将其写为向量的形式

优选的，S2中，在设计域内的不同子分区上分别定义局部坐标摄动函数CPFs。

按照连续力学的思想，非均匀点阵结构可以看作是均匀点阵结构(即初始构型)叠加特定变形梯度场后得到的当前构型。从几何上进行这一过程的关键是变形梯度场的参数化描述。因此从这个角度来看，由特定基函数组成的非均匀坐标变换场可以用来构建具有梯度渐变效应的点阵结构，并可进一步用于优化该结构的力学性能。具体来说，就是在公式(1.3)所表示的坐标变换场中，对全局坐标施加坐标摄动，表述如下：

其中，表示摄动施加后的全局坐标，且有：

其中，f(x,y)和g(x,y)分别表示不同坐标轴方向的坐标摄动函数coordinateperturbation functions,CPFs)，可按下式给出：

其中，α_m和β_n分别为摄动基函数和ψ_n(x,y)的控制系数。

基于上述处理方式，CPFs可以应用于整个设计域进行全局摄动，如图2所示。CPFs的基函数通常可以选用三角函数、多项式或指数函数等，并在整个设计域内采用相同的形式。这种做法虽然可以保证摄动场的连续性，然而低阶的线性插值在实现设计域的局部控制方面效果较差，从而限制了设计空间；而高阶插值将会带来数值振荡现象。因此有必要进一步地将CPFs在设计域内的几个子域中分别构造。这样做地问题在于，点阵结构的梯度效应将在子域中分别施加，这会造成相邻单胞间的材料分布连续性较差，如图3所示，而使得整个结构的力学性能变差。

优选的，S3中，选用B-spline基函数作为来保持各个分区之间CPFs的高阶连续性。

在有限元分析(finite element analysis,FEA)中，用于在单元内(参数空间)构造位移场的多项式插值形函数，可以同时保持单元内位移场至少C¹的连续性和单元间位移场C⁰的连续性。受此启发，可以使用类似构造单元内位移场的方法构造基于设计域分区的CPFs，这样相邻分区之间CPFs的连续性足以维持子域边界材料分布的基本连通性。然而，数值实验表明，仅有C⁰阶连续性的CPFs将导致设计上出现尖角，如图4所示，这将对结构性能产生不利影响。而虽然提高形函数中多项式阶数可以缓解这一问题，但拉格朗日多项式固有的振荡性和昂贵的计算成本仍然是难以解决的。因此采用B-spline基函数作为来保持各个分区之间CPFs的高阶连续性。

具体包括：

B-spline的基函数的定义有许多种，其中Cox-de-Boor递归公式因其在数值实现上的便利性而被广泛采用。采用Cox-de-Boor递归公式定义的B-spline的基函数，在构造B-spline之前首先需要在参数空间中定义一个节点向量，它由一组非减的参数组成：Ξ＝{ξ₀ξ₁,ξ₂,…,ξ_n+p+1}，i＝0,1,…,n+p)；其中，ξ_i代表第i个节点，n和p分别代表B-spline基函数的数量和阶数；节点序列中任意两个节点之间的距离可用于区分基函数的均匀性；若在同一节点序列中某一特定的节点出现了k次，则称该节点具有k次重复度；

对于零阶基函数，相应的B-spline基函数为：

而大于等于1阶的基函数可以按如下公式给出：

其中，若出现0/0的情况，则该项直接指定为0；

通常，一条p阶B-spline曲线可以按如下表达式给出：

其中，P_k，k＝0,1,…,n表示物理空间中控制多边形上的第k个控制点；若考虑分别沿ξ，η两个参数轴方向各给出一组控制点向Ξ¹＝{ξ₀,ξ₁,…,ξ_n+p+1}和Ξ²＝{η₀,η₁,…,η_m+q+1}，则二者的张量积可按以下方式用于定义B-spline曲面：

其中，P_i,j表示物理空间中相应的控制点坐标；

根据以上关于B-spline基函数的构造，基于各个子分区定义的CPFs可以表达为：

其中，α_i,j和β_i,j分别用于构造不同方向CPFs中基函数物理空间中的控制点网络；

至此，整个渐变点阵结构在MMC方法框架下的拓扑描述函数可以表达为：

其中，更新后的设计变量向中的包含了与CPFs相关的摄动系数变量。

优选的，S4中，基于上述渐变点阵结构的显式描述框架，相应的结构柔度最小化的优化问题列式可以写为：

S.t.

其中，D，f，t，和分别表示整个设计域，体力密度，面力边界Γ_t上的给定面力和位移边界Γ_u上的给定位移；/>表示一阶Sobolev空间；H＝H(x)表示Heaviside函数，且当x>0时，H(x)＝1，否则H(x)＝0；此外，点x处的四阶弹性张量被表达为v＝0on S_u}代表虚位移v的可行解集；V和/>分别表示设计域中实体材料的用量和给定的材料用量上限；

采用统一的四边形平面应力单元来离散设计域，其中单元刚度矩阵k_e可以按下式计算：

其中，B为应变-位移矩阵；Ω_e表示第e个单元所占区域；D^s表示实体材料的本构矩阵；ρ_e表示第e个单元的密度，ρ_e可以按下式计算：

其中，表示整个结构的拓扑描述函数在第i个节点位置上的值；q>0和H_∈(x)分别表示惩罚因子和正则化Heaviside函数，在具体数值实现中，q＝2，H_∈(x)按下式给出：

其中，∈和α分别为用于控制光滑化程度和避免总刚奇异的两个取正值的参数，具体来说，∈取2-4倍的最小网格尺寸，α＝10^-3；

基于以上公式，目标函数和约束函数对任意设计变量的灵敏度可以分别表示为：

其中，K和分别为总刚和单刚；u代表位移向量；NE代表设计域内单元总数，对于描述渐变点阵结构中的第k根组件的第j个设计变量/>和用于构造CPFs的第t个摄动系数α_t，上两式中的/>项可以根据链式法则解析地得到：

最大值算子可用以下K-S函数代替：

其中，参数l＝100。

因此，本发明采用上述步骤的一种基于MMC框架的非均匀点阵结构优化设计方法，采用周期性点阵结构作为初始设计，利用在设计域内不同子分区上分别定义的局部坐标摄动函数CPFs在设计域内生成不同的微结构拓扑形式，结构具有清晰明确的边界，相比以往方法扩大了设计空间。使用B-spline基函数作为CPFs的基函数，并在参数空间中采用控制点网络对每个分区进行局部控制，确保相邻微结构之间材料分布的光滑过渡。优化过程中所涉及的设计变量也仅包括组件的几何参数和CPFs的摄动系数，使用更少的设计变量，可以在各种载荷条件下，在更大的设计空间中获得具有复杂拓扑构型的渐变点阵结构。

下面通过附图和实施例，对本发明的技术方案做进一步的详细描述。

附图说明

图1为本发明一种基于MMC框架的非均匀点阵结构优化设计方法的组件几何描述示意图；

图2为本发明通过CPFs应用于整个设计域进行全局摄动，实现通过坐标变换引入梯度渐变效应的效果图；

图3为本发明采用CPFs应用于整个设计域进行全局摄动时相邻单胞间的材料分布图示意；

图4为本发明中仅采用仅C⁰连续性基函数带来的“尖角”问题示意图；

图5为本发明实施例中短梁算例的边界条件及初始布局；

图6为本发明实施例中短梁算例布置32根初始组件的优化结果；

图7为本发明实施例中短梁算例布置128根初始组件的优化结果。

具体实施方式

以下通过附图和实施例对本发明的技术方案作进一步说明。

实施例

本方法的有效性和有益效果将通过经典的短梁算例来说明，通过本发明的优化方法对短梁结构设计进行优化。

在算例中涉及的量都是无量纲的，设计域的厚度为单位值。同时，所用的各向同性的实体材料的杨氏模量和泊松比分别为E^s＝1和ν^s＝0.3。此外，选择MMA算法作为优化器来更新设计变量。当每个设计变量在连续两次迭代之间的相对变化的最大值低于指定的阈值(即5％)或达到最大迭代步数时，优化过程终止。除非另有说明，统一采用3阶B-spline基函数来构建CPFs。

如图5所示，本算例中采用的是长度为4宽度为2的矩形设计域。网格规模为400×200；边界条件为左端固支，右侧边界中点施加一竖直向下单位集中力F；给定材料体积分数上限为0.4；设计域外侧有一层厚度为0.04的不可设计域；初始设计共布置8×4共32根组件；设计域被均匀地划分为8*4个分区，摄动系数的变化范围取为[-2,2]。

优化结果如图6所示，最终的目标函数值为C_opt＝88.10。

在初始设计域中布置更多的初始组件，即16×8共128根初始组件，得到的优化如图7所示，其目标函数值为C_opt＝85.84。

从优化结果中可以看出，优化后结构中单胞的大小、尺寸和方向相对初始结构都发生了较大变化，其中微结构的方向与结构在所施加边界条件下产生的主应力方向趋近，且单胞之间过度平滑，几乎看不到微结构之间连通性较差的部分。材料在靠近固支端的主要传力路径位置产生聚集，而这里也正是应力水平较高的区域。而在应力水平较低的位置分配更少的材料也能有效地提高材料的使用效率，提高结构性能。当使用的初始组件数量更多时，优化后的结果中也呈现多尺度结构的特征。相对于传统连续体的拓扑优化结果，本方法所得的渐变结构保留了更多点阵结构的复杂微结构特征。

此外，尽管使用的网格中单元数量达到80000个，但由于结构的几何模型与分析模型解耦的优势，上述优化结果所使用的设计变量数仅分别为160和448。仅通过较少地设计变量就可得到具有相当复杂程度且边界清晰的渐变点阵结构，这是传统基于像素描述的拓扑优化方法难以实现的。

最后应说明的是：以上实施例仅用以说明本发明的技术方案而非对其进行限制，尽管参照较佳实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对本发明的技术方案进行修改或者等同替换，而这些修改或者等同替换亦不能使修改后的技术方案脱离本发明技术方案的精神和范围。

Claims (4)

1.一种基于MMC框架的非均匀点阵结构优化设计方法，其特征在于，步骤如下：

S1、MMC显式拓扑优化：采用基于MMC框架的周期性点阵结构作为初始设计，通过组件的显式几何参数得到拓扑描述函数TDFs，用TDFs来表示实体材料在设计域中的分布；

S2、坐标分区映射：在设计域内的不同子分区上分别定义局部坐标摄动函数CPFs，从而在设计域内生成不同的微结构拓扑形式，并确保相邻微结构之间的连接；

S3、基于B-spline基函数的CPFs构造：通过B-spline基函数保持各个子分区之间CPFs的高阶连续性，并在参数空间中采用控制点网络对每个分区进行局部控制，形成基于显式描述框架的渐变点阵结构，采用Cox-de-Boor递归公式定义的B-spline的基函数，在构造B-spline之前首先需要在参数空间 中定义一个节点向量，它由一组非减的参数组成：Ξ＝{ξ₀ξ₁,ξ₂,…,ξ_n+p+1}，/>其中，ξ_i代表第i个节点，n和p分别代表B-spline基函数的数量和阶数；节点序列中任意两个节点之间的距离可用于区分基函数的均匀性；若在同一节点序列中某一特定的节点出现k次，则称该节点具有k次重复度；

对于零阶基函数，相应的B-spline基函数为：

而大于等于1阶的基函数按如下公式给出：

其中，若出现0/0的情况，则该项直接指定为0；

一条p阶B-spline曲线按如下表达式给出：

其中，P_k，k＝0,1,…,n表示物理空间中控制多边形上的第k个控制点；若考虑分别沿ξ，η两个参数轴方向各给出一组控制点向量Ξ¹＝{ξ₀,ξ₁,…,ξ_n+p+1}和Ξ²＝{η₀,η₁,…,η_m+q+1}，则二者的张量积按以下方式用于定义B-spline曲面：

其中，P_i,j表示物理空间中相应的控制点坐标；

根据以上关于B-spline基函数的构造，基于各个子分区定义的CPFs表达为：

其中，α_i,j和β_i,j分别用于构造不同方向CPFs中基函数物理空间中的控制点网络；

至此，整个渐变点阵结构在MMC方法框架下的拓扑描述函数可以表达为：

其中，更新后的设计变量向中的包含了与CPFs相关的摄动系数变量；

S4、数值实现：在满足体积约束的前提下，优化初始周期性点阵结构，使结构柔度最小化，从而产生一个具有平滑连接和清晰边界的渐变点阵结构。

2.根据权利要求1所述的一种基于MMC框架的非均匀点阵结构优化设计方法，其特征在于：S1具体包括：

在欧拉描述的MMC框架中，结构的TDFs描述如下：

其中，D代表预先给定的设计域，是由n根组件所构成的实体材料占据的空间，φ^s(x)＝max(φ₁(x),…,φ_n(x)表示整个结构的TDFs；而其中的φ_i(x),(i＝1,2,…,n)表示第i根组件的TDF，并可按下式详细写为：

且有

其中，p是描述超椭圆的参数，指定为偶数，在公式(1.2)和公式(1.3)中，符号(x_0i,y_0i)、a_i、f_i(x')、θ_i分别表示从水平轴逆时针测量的第i根组件的中心点位置、半长、半宽和从水平轴逆时针测量的倾斜角，描述组件厚度变化曲线的函数f_i(x')按下式给出：

其中，和/>分别为第i根组件不同位置处的半宽，函数f_i(x')是在局部坐标下定义的，基于以上组件几何描述，与第i根组件相关的设计变量写为其中/>包含了与f_i(x')相关的设计参数；将整个结构的拓扑构型通过若干设计变量来描述，并将其写为向量的形式

3.根据权利要求2所述的一种基于MMC框架的非均匀点阵结构优化设计方法，其特征在于：S2具体包括：

在公式(1.3)所表示的坐标变换场中，对全局坐标施加坐标摄动，表述如下：

其中，表示摄动施加后的全局坐标，且有：

其中，f(x,y)和g(x,y)分别表示不同坐标轴方向的坐标摄动函数CPFs，可按下式给出：

其中，α_m和β_n分别为摄动基函数和ψ_n(x,y)的控制系数。

4.根据权利要求3所述的一种基于MMC框架的非均匀点阵结构优化设计方法，其特征在于：S4中，基于上述渐变点阵结构的显式描述框架，相应的结构柔度最小化的优化问题列式写为：

S.t.

其中，D，f，t，和分别表示整个设计域，体力密度，面力边界Γ_t上的给定面力和位移边界Γ_u上的给定位移；/>表示一阶Sobolev空间；H＝H(x)表示Heaviside函数，且当x>0时，H(x)＝1，否则H(x)＝0；此外，点x处的四阶弹性张量被表达为/> v＝0 on S_u}代表虚位移v的可行解集；V和/>分别表示设计域中实体材料的用量和给定的材料用量上限；

采用统一的四边形平面应力单元来离散设计域，其中单元刚度矩阵k_e按下式计算：

其中，B为应变-位移矩阵；Ω_e表示第e个单元所占区域；D^s表示实体材料的本构矩阵；ρ_e表示第e个单元的密度，ρ_e按下式计算：

其中，表示整个结构的拓扑描述函数在第i个节点位置上的值；q>0和H_∈(x)分别表示惩罚因子和正则化Heaviside函数，在具体数值实现中，q＝2，H_∈(x)按下式给出：

其中，∈和α分别为用于控制光滑化程度和避免总刚奇异的两个取正值的参数，具体来说，∈取2-4倍的最小网格尺寸，α＝10^-3；

基于以上公式，目标函数和约束函数对任意设计变量的灵敏度分别表示为：

其中，K和分别为总刚和单刚；u代表位移向量；NE代表设计域内单元总数，对于描述渐变点阵结构中的第k根组件的第j个设计变量/>和用于构造CPFs的第t个摄动系数α_t，上两式中的/>项可以根据链式法则解析地得到：

最大值算子可用以下K-S函数代替：

其中，参数l＝100。