一种点阵结构优化设计的基于水平集的变密度方法与流程

1.本发明属于结构优化设计相关技术领域，更具体地，涉及一种点阵结构优化设计的基于水平集的变密度方法。


背景技术：

2.点阵结构(lattice structure)是一种新型结构材料，点阵结构由多孔单胞排列而成，因其具有高比强度、高比刚度、高能量吸收性、超轻等优良力学性能，以及吸声、吸热、减振、电磁屏蔽等特殊性能而被广泛研究。目前点阵结构已经在航天航空、机械、生物、隔热隔声材料、医疗机械等领域发挥着重要的作用。
3.传统点阵结构是等密度的，然而在实际应用中，结构各个部分承受着不同的载荷。为了使材料性能得到充分发挥，应该对点阵结构进行优化，确定材料的最佳密度分布，使得在同样重量的情况下，最大化点阵结构的刚度。
4.点阵结构传统的优化方法是“基结构法”(ground structure approach,gsas)。在该方法中，每个晶胞单元或支柱的一个或几个参数可以通过在设计领域中使用预定义的“基结构”进行数值优化，以满足设计目标。但“基结构法”设计变量数量多，优化时间长。
5.随着拓扑优化技术的发展，拓扑优化开始大量用于点阵结构的优化设计。采用拓扑优化方法求得点阵结构的最优密度分布，再利用连续映射技术将最优密度分布转化为变密度点阵结构。
6.目前，用于点阵结构的拓扑优化方法，主要存在以下三个方面的问题：(1)大多数拓扑优化方法都像“基结构法”一样，需要预先定义一些标准晶胞单元，然而，这极大地限制了设计的自由度；(2)优化后的晶胞单元之间存在不连续性问题；(3)几乎所有这些方法都没有考虑点阵结构的可制造性。壁厚不均匀的点阵结构很难用传统制造方法制造。即使增材制造(additive manufacturing，am)可以制造它们，但也有最小壁厚的限制。此外，壁厚的不均匀变化可能会导致额外的热变形和热应力。


技术实现要素：

7.针对现有技术的以上缺陷或改进需求，本发明提供了一种点阵结构优化设计的基于水平集的变密度方法，其目的在于，对点阵结构进行优化，晶胞连续性好、非均匀的网状的点阵结构，由此解决现有点阵结构晶胞不连续、可制造性差的技术问题。
8.为实现上述目的，本发明提供了以下技术方案：
9.一种点阵结构优化设计的基于水平集的变密度方法，包括如下步骤：
10.(s1)将待优化设计的点阵结构的结构设计域划分为n个单元，并在所述结构设计域中设置m个设计点q
i
即径向基函数的中心点，并为每个设计点q
i
赋予权重系数α
i
、β
i
，结合权重系数α
i
、β
i
和径向基函数得到各单元中心点处的两个水平集函数φ(x
e
)、ψ(x
e
)；n和m均为正整数且n＞m；
11.(s2)利用余弦函数对水平集函数φ(x)、ψ(x)进行变换得到函数φ(x)、ψ(x)，并
采用可微的max函数对φ(x)、ψ(x)进行组合得到ρ(x)，基于ρ(x)，通过近似heaviside函数得到单元密度
12.(s3)分别计算各单元中心点处两水平集函数梯度的模再根据计算式分别得到然后分别对进行p
‑
norm凝聚得到进而求得约束函数其中ξ1、ξ2为小于1的正数；
13.(s4)进行有限元分析：令点阵结构边界处若干层单元密度为1，通过单元密度计算得到各单元的杨氏模量e
e
，利用e
e
建立各单元的刚度矩阵k
e
，通过组装k
e
得到整体刚度矩阵k，再根据公式ku＝f求解出整体位移向量u，其中f为合外力向量；
14.(s5)定义点阵结构优化设计问题：设计变量为设计点q
i
处的权重系数α
i
、β
i
，设计目标是使点阵结构柔度最小化，设计约束为力平衡方程ku＝f和梯度约束g1、g2以及α
i
、β
i
的上下界约束，，该优化设计问题的数学表达形式如下：
15.find：α
i
，β
i
(i＝1,2,
…
,m)，
16.min
17.st.ku＝f
[0018][0019][0020]
α
min
＜α
i
＜α
max
[0021]
β
min
＜β
i
＜β
max
[0022]
其中，u
e
为单元位移列向量；k0为单位杨氏模量单元刚度矩阵；t表示矩阵转置；e表示第e个单元；α
min
、β
min
分别表示设计变量α
i
和β
i
的最小值，α
max
、β
max
分别表示设计变量α
i
和β
i
的最大值；
[0023]
(s6)计算目标函数c对设计变量α
i
、β
i
的灵敏度计算约束函数g1对设计变量α
i
的灵敏度和约束函数g2对设计变量β
i
的灵敏度
[0024]
(s7)基于灵敏度对步骤(s5)定义的点阵结构优化设计问题进行求解，更新设计变量α
i
、β
i
；
[0025]
(s8)重复步骤(s2)至步骤(s7)，直至满足预设的优化终止条件，完成点阵结构的优化设计。
[0026]
利用最终更新的α
i
、β
i
对应的单元中心点处新的水平集函数值φ(x
e
)、ψ(x
e
)以及
单元密度绘制设计域密度的灰度图，即得到优化后的点阵结构；
[0027]
优选地，步骤(s1)中，所述各单元中心点处的两个水平集函数的计算式为：
[0028][0029][0030]
其中，m为径向基函数的中心点数，为径向基函数，其表达式为：
[0031][0032]
式中，符号(
·
)

表示max{0,
·
}，参数r根据径向基函数中心点坐标q
i
＝(x
i
,y
i
)和单元中心点坐标x
e
＝(x
e
,y
e
)计算得到，计算式为：
[0033][0034]
式中，ε为预设的正数，d
min
取值大于所述结构设计域最大尺寸的一半。
[0035]
优选地，步骤(s2)中，所述利用余弦函数对所述两个水平集函数φ(x)、ψ(x)进行变换得到函数φ(x)、ψ(x)，计算式为：
[0036][0037][0038]
其中，λ1、λ2为余弦函数波长。
[0039]
优选地，步骤(s2)中，所述采用可微的max函数对φ(x)、ψ(x)进行组合得到ρ(x)，计算式为：
[0040][0041][0042]
其中，χ取常数。
[0043]
优选地，步骤(s2)中，所述近似heaviside函数为：
[0044][0045]
式中，参数γ、η分别控制heaviside函数的平滑性和阶跃点；
[0046]
所述单元密度的计算式为：
[0047][0048]
优选地，步骤(s3)中，水平集函数梯度的模的计算式为：
[0049][0050]
其中：
[0051][0052][0053]
式中，的计算式为：
[0054][0055][0056]
同理可计算水平集函数梯度的模
[0057]
所述对进行p
‑
norm凝聚得到的计算式为：
[0058][0059][0060]
式中，p
n
为不小于2的正整数。
[0061]
优选地，步骤(s4)中，所述各单元的杨氏模量e
e
的计算式为：
[0062][0063]
式中，e
min
为点阵结构中孔洞杨氏模量，e0为点阵结构中实体材料杨氏模量；所述各单元的刚度矩阵k
e
的计算式为：
[0064]
k
e
＝e
e
·
k0[0065]
式中，k0为单位杨氏模量单元刚度矩阵。
[0066]
优选地，步骤(s6)中，所述目标函数c对设计变量α
i
、β
i
的灵敏度的计算式为：
[0067][0068][0069]
其中：
[0070][0071][0072][0073][0074][0075][0076]
式中，e
e
为单元杨氏模量；为单元杨氏模量e
e
对单元密度的一阶微分；u
e
为单元位移列向量；t表示矩阵转置；k0为单位杨氏模量单元刚度矩阵；e表示第e个单元；n为结构设计域单元数量；是的简写，φ
e
是φ(x
e
)的简写，ψ
e
是ψ(x
e
)的简写；
[0077]
所述约束函数g1对设计变量α
i
的灵敏度的计算式为：
[0078][0079][0080]
所述约束函数g2对设计变量β
i
的灵敏度的计算式为：
[0081][0082][0083][0084]
优选地，步骤(s7)中，对进行灵敏度滤波得到基于滤波后的灵敏度以及灵敏度对步骤(s5)定义的点阵结构优化设计问题进行求解；
[0085]
其中，的计算式为：
[0086][0087][0088][0089]
式中，n
e
为到单元e的距离小于滤波半径r
min
的单元集合，是第i个单元的密度，e
i
是第i个单元的杨氏模量，i∈n
e
；μ是正数，用于避免零除；h
ei
是权重因子，h
ei
的计算式为：
[0090]
h
ei
＝max(0,r
min
‑
δ(i,e))
[0091]
式中，δ(i,e)表示单元i到单元e的距离。
[0092]
优选地，步骤(s8)中，所述预设的优化终止条件为：
[0093][0094]
其中，c
err
为柔度误差，k为当前迭代次数，δ为柔度误差的下限值，k
max
为迭代次数的上限值。
[0095]
总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有
益效果：
[0096]
1、本发明将点阵结构设计域划分为若干个单元，在结构设计域内均匀设置设计点，在设计点处定义两个权重系数作为设计变量，结合权重系数和径向基函数得到两个水平集函数，利用余弦函数、可微的max函数和近似heaviside函数对水平集函数进行变换和组合得到单元密度；通过有限元分析，利用单元密度计算单元杨氏模量，利用单元杨氏模量建立刚度矩阵，进而求解处整体位移向量；构建制造约束的约束函数，以柔度最小化为优化目标，以制造约束为约束条件，定义点阵结构优化设计问题；计算目标函数和约束函数对设计变量的灵敏度，基于灵敏度更新设计变量，直至满足优化终止条件，得到最优的点阵结构。
[0097]
本发明基于水平集函数的变密度方法，单元密度可变，利用水平集函数、余弦函数、可微的max函数、平滑的近似heaviside函数、p
‑
norm凝聚等，可得到晶胞连续性好、材料密度分布佳的网状的点阵结构；不需要预定义标准晶胞单元，提高了设计的灵活性；能够容易地对水平集函数的参数(λ1,λ2)进行调整从而得到不同孔隙率、晶胞大小的点阵结构；对水平集函数添加梯度约束，从而满足制造约束，极大提升点阵结构的可制造性；方法简单，便于实施，适用性较强，容易与现有有限元软件集成。
附图说明
[0098]
图1是本发明较佳实施例中点阵结构优化设计的基于水平集的变密度方法的流程图；
[0099]
图2是图1的点阵结构优化设计的基于水平集的变密度方法涉及的一个具体结构设计域示意图；
[0100]
图3是图1的优化例子点阵结构的初始设计示意图；
[0101]
图4是图1的优化例子点阵结构的优化结果示意图；
[0102]
图5是图1的优化例子中目标函数的迭代曲线。
具体实施方式
[0103]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0104]
如图1和图2所示，本发明提供一种点阵结构优化设计的基于水平集的变密度方法，且以带有集中载荷的平面简支梁点阵结构柔度最小化的优化问题为例，来具体解释该方法。在给定的1m
×
2m矩形结构设计域d内固定区域左下角，限制区域右下角竖直方向自由度，并在区域上边界中心处施加集中力f＝1n。图3为实施例点阵结构的初始设计，其为传统的均匀的点阵结构，应用本发明提供的方法对上述简支梁点阵结构进行优化，使其柔度最小。
[0105]
本发明提供的点阵结构优化设计的基于水平集的变密度方法主要包括如下步骤：
[0106]
步骤一，将待优化设计的点阵结构的结构设计域划分为n个单元，并在结构设计域中设置m个设计点q
i
即径向基函数的中心点，并为每个设计点q
i
赋予权重系数α
i
、β
i
，结合权
重系数α
i
、β
i
和径向基函数得到各单元中心点处的两个水平集函数φ(x
e
)、ψ(x
e
)；n和m均为正整数且n＞m。
[0107]
具体地，将图2给定的结构设计域d划分为100
×
200个的正方形单元；设计域d内均匀选取20
×
40个径向基函数的中心点q
i
，并在每个q
i
处定义权重系数α
i
、β
i
作为优化中的设计变量，结合权重系数α
i
、β
i
和径向基函数可得各单元中心点处两个水平集函数值φ(x
e
)、ψ(x
e
)，具体计算式为：
[0108][0109][0110]
其中，m为径向基函数的中心点q
i
的数量，本实施例中m＝20
×
40，为径向基函数，这里选用紧支径向基函数，其表达式为：
[0111][0112]
其中符号(
·
)

表示max{0,
·
},参数r可根据径向基函数中心点坐标q
i
＝(x
i
,y
i
)和单元中心点坐标x
e
＝(x
e
,y
e
)计算得到，具体计算式为：
[0113][0114]
其中ε为一个小的正数，这里取为10
‑8；d
min
推荐选取的值大于设计域最大尺寸的一半，这里d
min
取1.2；
[0115]
根据初始水平集函数φ0(x
e
)、ψ0(x
e
)求解得到初始权重系数具体地，在本实施例中，令水平集函数φ0(x
e
)、ψ0(x
e
)的值分别等于中心点坐标(x
e
,y
e
)的x
e
和y
e
。
[0116]
步骤二，利用余弦函数对水平集函数φ(x)、ψ(x)进行变换得到函数φ(x)、ψ(x)，并采用可微的max函数对φ(x)、ψ(x)进行组合得到ρ(x)，基于ρ(x)，通过近似heaviside函数得到单元密度
[0117]
具体地，利用余弦函数对两个水平集函数φ(x)、ψ(x)进行变换得到函数φ(x)、ψ(x)，具体计算式为：
[0118][0119][0120]
其中参数λ1、λ2为余弦函数波长，分别控制点阵结构单胞两个方向的尺寸，这里选取λ1＝λ2＝0.16；
[0121]
采用可微的max函数对φ(x)、ψ(x)进行组合：
[0122]
ρ(x)＝max{φ(x),ψ(x)}
[0123]
可微的max函数的表达式为：
[0124][0125]
其中函数abs
χ
(x
‑
y)通过下式计算：χ取较小的常数，可取χ＝0.01；基于得到的ρ(x)，进一步可得到单元密度计算式为：
[0126][0127]
其中h(x)为近似heaviside函数，具体计算式为：
[0128][0129]
式中参数γ、η分别控制heaviside函数的平滑性和阶跃点，分别取10、0.8。
[0130]
步骤三，分别计算各单元中心点处两水平集函数梯度的模再根据计算式再根据计算式分别得到然后分别对进行p
‑
norm凝聚得到进而求得约束函数其中ξ1、ξ2为小于1的正数。
[0131]
具体地，分别计算各单元中心点处水平集函数梯度和的模，具体计算式为：
[0132][0133][0134]
其中：
[0135][0136][0137]
式中的具体计算式为：
[0138][0139][0140]
进一步地，根据计算式分别得到然后分别对进行p
‑
norm凝聚得到具体计算式为：
[0141][0142][0143]
式中参数p
n
一般为不小于2的正整数，且不宜取值过大，这里取12。从而可求得约束函数其中ξ1、ξ2为取值通常小于1的正数，这里取0.3。
[0144]
步骤四，进行有限元分析：令点阵结构边界处若干层单元密度为1，通过单元密度计算得到各单元的杨氏模量e
e
，利用e
e
建立各单元的刚度矩阵k
e
，通过组装k
e
得到整体刚度矩阵k，再根据公式ku＝f求解出整体位移向量u，其中f为合外力向量。
[0145]
具体地，为保证结构的稳定性，在步骤二得到的单元密度的基础上，令结构边界处4层单元密度为1；
[0146]
单元的刚度矩阵k
e
的计算式为：
[0147]
k
e
＝e
e
·
k0[0148]
其中k0为单位杨氏模量单元刚度矩阵；是的简写，e
e
为单元杨氏模量，与单元密度有关，具体计算式为：
[0149][0150]
式中e
min
和e0分别为孔洞和实体材料杨氏模量，这里分别取为10
‑3、1；参数p为惩罚因子，这里取为3；通过组装单元刚度矩阵k
e
得到整体刚度矩阵k,然后根据公式ku＝f求解出整体位移向量u，其中f为合外力向量，以此完成结构有限元分析。
[0151]
步骤五，定义点阵结构优化设计问题：设计变量为设计点q
i
处的权重系数α
i
、β
i
，设计目标是使点阵结构柔度最小化，设计约束为力平衡方程ku＝f和g1、g2、α
i
、β
i
的制造约束，该优化设计问题的数学表达形式如下：
[0152]
find：α
i
，β
i
(i＝1,2,
…
,m)，
[0153]
min
[0154]
st.ku＝f
[0155][0156][0157]
α
min
＜α
i
＜α
max
[0158]
β
min
＜β
i
＜β
max
[0159]
其中，α
i
、β
i
为设计变量；m为径向基函数中心点数；u为整体位移向量；k为整体刚度矩阵；f为合外力向量；e
e
为单元杨氏模量；u
e
为单元位移列向量；k0为单位杨氏模量单元刚度矩阵；t表示矩阵转置；e表示第e个单元；n为结构设计域单元数量；分别表示对进行p
‑
norm凝聚之后的值；ξ1、ξ2为一个取值通常小于1的正数；α
min
、β
min
分别表示设计变量α
i
和β
i
的最小值，α
max
、β
max
分别表示设计变量α
i
和β
i
的最大值。
[0160]
步骤六，计算目标函数c对设计变量α
i
、β
i
的灵敏度计算约束函数g1对设计变量α
i
的灵敏度和约束函数g2对设计变量β
i
的灵敏度
[0161]
具体地，利用k
e
和e
e
的偏导数关系、e
e
和α
i
、β
i
的偏导数关系以及整体位移向量u，计算目标函数c对设计变量α
i
、β
i
的灵敏度，具体计算式为：
[0162][0163][0164]
其中：
[0165][0166][0167][0168]
[0169][0170][0171][0172]
式中，e
e
为单元杨氏模量；为单元杨氏模量e
e
对单元密度的一阶微分；u
e
为单元位移列向量；t表示矩阵转置；k0为单位杨氏模量单元刚度矩阵；e表示第e个单元；n为结构设计域单元数量；是的简写，φ
e
是φ(x
e
)的简写，ψ
e
是ψ(x
e
)的简写。
[0173]
约束函数g1对设计变量α
i
的灵敏度的计算式为：
[0174][0175][0176]
约束函数g2对设计变量β
i
的灵敏度的计算式为：
[0177][0178][0179]
利用如下公式对目标函数c对设计变量α
i
、β
i
的灵敏度进行过滤，得到过滤后的灵敏度
[0180][0181]
[0182]
其中，可由下式计算得到：
[0183][0184]
式中:n
e
为到单元e的距离小于滤波半径r
min
的单元集合，这里r
min
取1.5；是第i个单元的密度，e
i
是第i个单元的杨氏模量，i∈n
e
；μ是一个小的正数，用来避免零除；h
ei
是权重因子，可由下式计算：
[0185]
h
ei
＝max(0,r
min
‑
δ(i,e))
[0186]
其中δ(i,e)表示单元i到单元e的距离。
[0187]
步骤七，基于灵敏度对步骤五定义的点阵结构优化设计问题进行求解，更新设计变量α
i
、β
i
。
[0188]
具体地，根据目标函数c、约束函数g1，g2以及它们关于设计变量的灵敏度结合灵敏度滤波，通过移动渐进线方法(method of moving asymptotes，通常简称为mma)，更新设计变量α
i
、β
i
。
[0189]
利用更新后的α
i
、β
i
求得单元中心点处新的水平集函数值φ(x
e
)、ψ(x
e
)以及单元密度绘制设计域密度的灰度图，即得到新的点阵结构。
[0190]
步骤八，重复步骤二至步骤七，每一次重复称为一次迭代，直至满足预设的优化终止条件，得到点阵结构的最优设计，其中优化终止条件：
[0191][0192]
其中c
err
为柔度误差，k为当前迭代次数，δ为柔度误差的下限值，这里取0.05％，k
max
为迭代次数的上限值，这里取500。
[0193]
本发明较佳实施例的优化结果如下：优化后的点阵结构如图4所示，柔度值约为24，相对初始设计，柔度值减少41％；优化后的点阵结构边界光滑、清晰，能满足制造约束；图5是点阵结构最小柔度优化模型的目标函数的迭代曲线，从图中可以看出，前20步迭代中，目标函数减少迅速，20步迭代之后，因为约束函数的限制，目标函数收敛减慢，但最终会收敛到一个较稳定的值。
[0194]
本领域的技术人员容易理解，以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含
在本发明的保护范围之内。