一种单材料多孔结构的拓扑优化设计方法

本发明属于单材料结构拓扑优化设计相关技术领域，更具体地，涉及一种单材料多孔结构的拓扑优化设计方法。

背景技术：

增材制造(am，也被称为3d打印)能够制造出结合复杂度远远超过传统制造技术的组件，这使得一些复杂的轻量化制造成为可能。自然界中有许多常见的轻型壳多孔结构，如植物茎、鸟喙和人类骨骼。多孔结构是指结构由一个固体壳与多孔的内部所组成，而不完全为固体。多孔结构因其具有高的比刚度、良好吸声性和减震性受到了广泛研究和工程应用。

目前，多孔结构设计可分为单尺度设计和多尺度设计两大类。单尺度拓扑优化设计是将壳体和填充物分别视为两种材料。然而，尽管填充体可以按照预先定义的微观结构进行设计，但单尺度的填充体结构拓扑优化限制了结构性能，使得设计与实际的最优值相差较大。因此，常采用多尺度设计方法对壳体填充结构进行设计，多尺度方法设计多孔结构是以数值均匀化方法将微观结构和宏观材料性能之间联系起来进行设计多孔结构，但从宏观和微观同时设计多孔结构，这无疑给多孔结构拓扑优化设计带来了计算上的挑战。

技术实现要素：

针对现有技术的以上缺陷或改进需求，本发明提供了一种单材料多孔结构的拓扑优化设计方法，所述设计方法是一种整体式系统结构优化设计方法，能够优化出满足一定要求的单材料多孔结构，可以满足航空航天、火箭、导弹等军工制造中的结构性能设计需求。

为实现上述目的，按照本发明的一个方面，提供了一种单材料多孔结构的拓扑优化设计方法，所述设计方法包括以下步骤：

(1)分别基于初始水平集函数及径向基函数构建单材料多孔结构，进而通过扩展系数矩阵分块法进行插值以计算得到初始扩展系数，即设计变量初值；

(2)计算单材料多孔结构中连续变化的单元伪密度，并计算结构总体积；接着计算结构的每个单元的局部伪密度，并使用p范数方法求取局部伪密度最大值；

(3)基于得到的单元伪密度对单材料多孔结构的弹性模量进行插值以得到等效单元弹性模量；

(4)基于等效单元弹性模量在结构设计域中求解得到整体结构的位移场，进而计算单材料多孔结构最小柔度拓扑优化模型的目标函数；接着，计算目标函数、结构总体积及局部伪密度对设计变量的灵敏度，并更新全局设计变量，继而得到单材料多孔结构的结构拓扑。

进一步地，分别采用两种初始化结构构建初始水平集函数及径向基函数以构建单材料多孔结构，并通过离散小波变换方法压缩初始水平集函数及径向基函数矩阵，基于压缩后的初始水平集函数和径向基函数通过扩展系数矩阵分块法计算得到初始扩展系数，即设计变量初值。

进一步地，通过隐式水平集函数构建单材料多孔结构，以插值隐式水平集函数的扩展系数为结构设计变量。

进一步地，采用全局径向基函数对水平集函数进行插值。

进一步地，采用高斯积分将原始单元的4个节点上2*2个水平集函数值插值为41*41个水平集函数值。

进一步地，基于heaviside函数获得单材料多孔结构中连续变化的单元伪密度，并计算结构总体积；接着，通过圆形过滤法计算结构的每个单元的局部伪密度，并使用p范数方法求取局部伪密度最大值。

进一步地，步骤(2)包括以下子步骤：

首先，将设计域ω均匀划分，并采用heaviside函数将水平集函数值映射成为有限元模型中的单元伪密度，单元伪密度由以下公式求得：

ρi＝∫ωh(φi)dω

式中，φ为插值后的单元中节点的水平集函数值，φi是在插值之前的单元四个节点的水平集函数值；ω为设计域；h为heaviside函数；

接着，通过单元伪密度计算结构的总体积；

之后，通过圆形过滤法求解局部伪密度，ve为任意单元e为中心的局部伪密度：

其中，n表示在该单元的所有单元的数量，

n＝{i|||ρi-ρc||2≤r}

其中，r和ρc分别表示过滤半径和圆心单元；

接着，使用p范数求解单元伪密度的最大值。

进一步地，使用单元伪密度插值得到等效单元弹性模量的计算公式如下：

e＝ρi*e0

其中，e为插值后的等效单元弹性模量；ρi为单元伪密度；e0为材料的弹性模量。

进一步地，单材料多孔结构最小柔度拓扑优化模型的表达式为：

find:α＝[α1α2...αn]

αi,min≤αi≤αi,max

式中，vmax表示为整体结构体积约束；ve和vp分别表示结构局部伪密度及局部伪密度约束；u和v分别表示允许的位移空间u中的实位移场和虚位移场；u0表示狄利克雷边界上的位移；h(φ)是heaviside函数；α是设计变量，用来表示径向基函数插值后的扩展系数，且αi,max和αi,min表示设计变量的上下界限；a(u,v)＝l(v)表示为弹性平衡方程弱形式。

进一步地，使用离散小波变换方法压缩优化过程中的整体刚度矩阵、扩展系数矩阵及水平集函数矩阵。

总体而言，通过本发明所构思的以上技术方案与现有技术相比，本发明提供的单材料多孔结构的拓扑优化设计方法主要具有以下有益效果：

1.采用基于水平集方法的局部伪密度约束生成单材料多孔结构，这样不仅可以从宏观角度生成多孔结构，而且优化后的多孔结构具有光滑边界且不存在中间密度。

2.采用矩阵分块法来降低优化过程中所需的内存，可以在有限内存下实现大规模设计域的多孔结构拓扑优化设计工作。

3.引入全局径向基函数gsrbf对水平集函数进行插值，克服其在优化过程中需要重新初始化、速度扩展、迭代步长需满足逆风差分格式cfl条件及无法与许多成熟优化算法结合等数值缺陷，同时采用离散小波变换dwt方法压缩插值系数矩阵，从而形成一种极其稀疏的矩阵插值系统，在保证足够插值精度的前提下，进一步提高优化求解效率。

4.本发明提供的单材料多孔结构拓扑优化设计方法，与现有的技术方案相比，所研究结构是基于参数化水平集方法的单材料多孔结构，具有更加清晰的结构边界，可以满足航空航天、火箭、导弹等军工制造中的结构性能设计需求。

5.所述优化设计方法开发了参数化水平集方法与多孔结构问题相结合，可以得到具有较高稳定性且具有最佳材料相分布的多孔结构。

6.通过矩阵分块、矩阵范数及局部伪密度法的应用，可以实现仅在宏观尺度上对多孔结构的设计，降低优化复杂性及计算量。

附图说明

图1是本发明提供的单材料多孔结构的拓扑优化设计方法的流程示意图；

图2中的(a)、(b)及(c)分别是本发明实施例长悬臂梁初始设计域及初始孔洞位置的示意图；

图3中的(a)、(b)、(c)及(d)分别是本发明实施例均匀初始化下悬臂梁的优化结果示意图；(a)对应的过滤半径为10，体积约束为0.5；(b)对应的过滤半径为10，体积约束为0.6；(c)对应的过滤半径为15，体积约束为0.5；(d)对应的过滤半径为15，体积约束为0.6；

图4中的(a)、(b)、(c)及(d)分别是本发明实施例非均匀初始化下悬臂梁的优化结果示意图；(a)对应的过滤半径为10，体积约束为0.4；(b)对应的过滤半径为10，体积约束为0.5；(c)对应的过滤半径为15，体积约束为0.4；(d)对应的过滤半径为15，体积约束为0.5。

具体实施方式

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

请参阅图1，本发明提供的单材料多孔结构的拓扑优化设计方法，所述设计方法包括以下步骤：

步骤一，分别采用两种初始化结构构建初始水平集函数及径向基函数以构建单材料多孔结构，并通过离散小波变换方法压缩初始水平集函数及径向基函数矩阵，基于压缩后的初始水平集函数和径向基函数通过扩展系数矩阵分块法计算得到初始扩展系数，即设计变量初值。

具体地，通过隐式水平集函数构建单材料多孔结构，以插值隐式水平集函数的扩展系数为结构设计变量，在n个固定水平集节点处的结构隐式水平集函数为：

式中，x＝x1,x2,...,xn表示所有插值节点坐标，即水平集节点；n表示节点总数；αn表示水平集函数在节点n处的扩展系数；φ表示结构中表示单材料多孔结构的水平集函数；t为时间。该水平集函数由高斯径向基函数φn(x)插值；φn(x)表示高斯径向基函数，公式为：

其中，c是形状参数，等于水平集网格面积或者体积的倒数；xn表示水平集函数的第n个节点的坐标；||x-xn||是用于计算当前采样点x到xn节点距离的欧几里得范数。

为了提高优化效率，离散小波变换方法是减轻全插值矩阵造成的计算成本的关键，采用离散小波变换方法压缩插值系数矩阵，从而形成一种极其稀疏的矩阵插值系统，在保证足够插值精度的前提下，进一步提高优化求解效率。具体地，将原始插值矩阵a转换为其相同大小的小波形式采用小波基的矩阵可以很容易地区分其重要元素和冗余元素；因此，采用一个阈值方法来清除中适当数量的无用元素，并且重新构造一个更稀疏的插值矩阵最后，利用稀疏矩阵可以有效的计算水平集函数。

由于对水平集函数采用径向基函数插值后，在优化求解过程中，径向基函数矩阵a会随着设计域的增大而增大，这将会占用大部分的内存，导致在优化求解时出现矩阵大小超出计算机内存的问题。为了计算大规模的多孔结构，采用矩阵分块的方法对径向基函数矩阵进行分块计算，这种矩阵分块的方法要求aii(i＝1,2,…)的逆存在，以划分为2*2的矩阵块为例，求解扩展系数矩阵c＝a^-1·f。

扩展系数矩阵c为：

本实施方式中，初始化孔洞采用了两种形式对设计域进行初始化，一种初始化孔洞位置采用均匀密排圆孔结构，另一种为非均匀随机密排圆孔结构，通过这两种形式的初始化方法可以获得两种不同的多孔结构。

步骤二，基于heaviside函数获得单材料多孔结构中连续变化的单元伪密度，并计算结构总体积；接着，通过圆形过滤法计算结构的每个单元的局部伪密度，并使用p范数方法求取局部伪密度最大值。

具体地，首先将设计域ω均匀划分，并采用heaviside函数将水平集函数值映射成为有限元模型中的单元伪密度，对应的公式为：

式中，ζ是一个非常小的正数，本实施方式取ζ＝0.001，δ约等于heaviside函数带宽的一半。

单元伪密度可以由以下公式求得：

ρi＝∫ωh(φi)dω。

其中，为了计算的更为准确，使用高斯积分将原始单元的4个节点上2*2个水平集函数值插值为41*41个水平集函数值，插值公式如下：

ni＝(1+ξ0)·(1+η0)/4

在上式中，

ξ0＝ξi·ξ,η0＝ηi·η

ξi＝{-1,1,1-1},ηi＝{-1，-1,1,1}

ξ∈[-1:0.05:1],η∈[-1:0.05:1]；

φ为插值后的单元中节点的水平集函数值；φi是在插值之前的单元四个节点的水平集函数值。

接着，通过单元伪密度计算结构的总体积，公式为：v＝∫ωρidω。

之后，通过圆形过滤法求解局部伪密度，公式为：

ve为任意单元e为中心的局部伪密度：

其中，n表示在该单元的所有单元的数量，

n＝{i|||ρi-ρc||2≤r}

其中，r和ρc分别表示过滤半径和圆心单元。

接着，使用p范数求解单元伪密度的最大值，公式如下：

当p趋于无穷大时，有：

maxve＝‖ve‖p。

步骤三，基于连续变化单材料的单元伪密度，对单材料多孔结构弹性模量进行插值，以得到等效单元弹性模量。

具体地，定义材料的弹性模量为e0，泊松比为v，将单材料多孔结构的节点i的位移置为0，即

uix＝0,uiy＝0

其中，uix和uiy分别为节点i的x方向位移和y方向位移。

使用单元伪密度插值得到等效单元弹性模量的计算公式如下：

e＝ρi*e0

其中，e为插值后的等效单元弹性模量。

步骤四，基于参数化水平集理论建立单材料多孔结构最小柔度拓扑优化模型，并基于得到的等效弹性模量在结构设计域中通过有限元分析求解整体结构的位移场，根据得到的位移场计算单材料多孔结构最小柔度拓扑优化模型的目标函数；接着，基于自伴随方法计算目标函数、结构总体积和局部伪密度对设计变量的灵敏度，并采用mma移动渐近线算法更新全局设计变量，继而确定单材料多孔结构中材料的最优分布。

具体地，单材料多孔结构最小柔度拓扑优化模型的表达式为：

find:α＝[α1α2...αn]

αi,min≤αi≤αi,max

式中，vmax表示为整体结构体积约束，ve和vp分别表示结构局部伪密度及局部伪密度约束，u和v分别表示允许的位移空间u中的实位移场和虚位移场，u0表示狄利克雷边界上的位移，h(φ)是heaviside函数，α是设计变量，用来表示径向基函数插值后的扩展系数，且αi,max和αi,min表示设计变量的上下界限。a(u,v)＝l(v)表示为弹性平衡方程弱形式。

基于虚功原理，针对有限元平衡方程弱形式进行计算，对应的弱形式如下：

其中，a表示双线性能量式；l表示单线性荷载形式；dω为结构设计域的积分算子；h表示heaviside函数，用于表征结构形式的特征函数；ε为应变场；t表示矩阵的转置；u表示结构场的位移；v表示在动力学上允许的位移空间u中的一个虚拟位移；τ表示应用在边界的部分边界上的牵引力；p表示结构设计域的体积力；δ表示dirac函数，为heaviside函数的一阶微分；表示差分算子。

基于水平集方法构建单材料多孔结构最小柔度拓扑优化模型包括以下步骤：

(4.1)初始化设计参数、整体体积约束vmax和局部伪密度约束vup。

(4.2)基于单元伪密度求得的等效弹性模量e计算结构单元刚度ke，并组装得到整体刚度矩阵k；接着进行有限元分析求解结构位移场u。

(4.3)基于步骤(4.2)中求解的结构位移场u计算单材料多孔结构最小柔度拓扑优化模型的目标函数j：

其中，ε为应变场；t表示矩阵的转置；u表示结构场的位移；v表示在动力学上允许的位移空间u中的一个虚拟位移；e是通过单元伪密度插值后的等效弹性模量；dω为结构设计域的积分算子。

(4.4)基于自伴随方法求解目标函数与约束函数针对设计变量进行灵敏度分析，并采用mma移动渐近线算法更新全局设计变量，判断模型是否满足收敛条件，若否，则返回到步骤(4.2)，若是，则输出单材料多孔结构的最优拓扑结构。其中，通过计算本次与前一步的目标函数差值是否大于1x10^-6来决定是否进行下一次优化求解，即判断两次扩展系数的差值是否小于阈值；若是，则优化结束，输出最新的多孔结构拓扑；否则，转入下一次优化求解。

具体地，根据自伴随方法计算目标函数与约束函数针对优化设计变量的一阶微分，计算如下：

其中，

其中，α表示单材料多孔结构的设计变量，是高斯径向基函数插值时的扩展系数，仅与时间变量相关；ε为应变场；u表示结构场的位移；e是单材料多孔结构的等效弹性模量；φ表示多孔结构的水平集函数；τ表示应用在边界的部分边界上的牵引力；f表示结构设计域的体积力；h表示heaviside函数；δ表示dirac函数，为heaviside函数的一阶微分；dω为结构设计域的积分算子；为散度算子；p为p范数的指数；∑(·)为求和符号；ve为局部伪密度；表示第i个高斯径向基函数。

实施例

本实施例首先给出一个长悬臂梁结构作为算例说明本研究方法的有效性，如图2中的(a)所示为一个厚度为1的矩形区域，其结构左端被完全约束，在最右端中点处施加一个垂直向下的单位力，图2中的(b)为结构采用均匀孔洞进行初始化的结果，图2中的(c)为结构采用非均匀孔洞进行初始化的结果，比较两种初始化结果可以看出，均匀初始化可以明显看出孔洞位置及大小，而非均匀初始化由于孔洞的随机性，导致会有孔洞重叠而呈现出一种散点图的形式。在优化过程中，材料的弹性模量和泊松比分别为1和0.3，且结构被均匀离散成200×100的四节点方形单元。无量纲材料特性的设置是为了方便不同设计之间的比较，以下算例中均采用这种策略。

如图3所示给出了在不同的过滤半径r和体积约束下的均匀初始孔洞位置的优化结果，从优化结果中可以看出，通过调节过滤半径r的大小可以控制生成的多孔结构中分支结构的密度，从优化结果中还可以观察出，当初始孔洞为均匀分布时，生成的多孔结构在外边缘处会生成一圈类似壳状的填充，而分支结构是存在于壳内的。

如图4所示给出了在不同的过滤半径r和体积约束下的非均匀初始孔洞位置的优化结果。从优化结果中可以得出与均匀初始化优化结果相同的结论，即通过调节过滤半径r的大小可以控制生成的多孔结构中分支结构的密度。但与均匀初始化所得结论不同的是，当初始孔洞为非均匀分布时，没有像均匀初始孔洞优化时生成一圈壳体包罗，生成的分支结构是存在于整个结构中的。

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