一种多材料结构动刚度拓扑优化方法

1.本发明涉及一种工程结构优化设计领域的设计方法，更具体的是，本发明涉及一种多材料结构动刚度拓扑优化方法。


背景技术：

2.多相材料拓扑优化是在给定的边界条件以及目标和约束条件下，将具有不同力学特性的多种材料合理地配置在一起，以获得最优的结构性能。与传统单一材料的拓扑优化相比，多相材料的拓扑优化极大地扩大了设计空间，能够提供综合性能更优的最优解。随着多材料增材制造技术的发展，依据最优拓扑的计算模型，使通过逐点和逐层方式打印多材料设计结构成为可能。由此，多材料的拓扑优化问题引起了研究机构的广泛关注。
3.目前，多体积约束下的多材料布局优化已可拓展至热力耦合、频域动力学和空间桁架结构的最优布局问题。在基于边界方法的框架内，水平集法和移动构件法(mmc)因其能使边界光滑清晰、便于提取设计构型等独特优势，也收到了一定的关注与研究。但是，现有的基于边界的方法未能充分考虑新孔洞的创建，优化结果过度依赖于初始设计，难于获得优化问题的全局最优解。
4.拓扑优化框架要求在消耗尽可能小的计算成本下获得高分辨率的优化结果，而上述目标的实现取决于有限元求解器、自由度数量、优化建模策略、材料插值模型以及后处理等诸多因素。将多边形单元应用于拓扑优化问题，能够显著的减小棋盘格和孤岛效应等数值奇异性问题。高分辨率拓扑优化方法(mtop)，它采用多层级网格优化建模策略，即利用粗糙网格完成有限元分析，精细的重叠网格描述设计变量和密度变量空间，从而形成高分辨率的拓扑优化结果。将多边形单元替代传统单元作为分析单元，可以精确地估计复杂结构的动力学响应，高效地实现特征值、受迫振动等结构频域动力学问题的高分辨率拓扑优化。
5.目前，多相材料布局优化集中于静态优化问题，时域动响应拓扑优化仅限于单相材料设计，而多相材料的时域动力学布局优化问题研究相对较少。因此，将多边形网格与mtop的融合方法拓展至多材料结构的动态拓扑优化问题，建立融合框架内的多相材料高分辨率布局的时域动力学优化模型，为结构设计提供更多的自由度并最终获得更为优质的结构设计方案。


技术实现要素：

6.针对现有技术的上述缺点和/或改进需求，本发明利用多边形有限单元的高精度求解优势，融合多分辨率拓扑优化方法，实现粗糙位移网格条件下的高分辨率构型设计，提供了一种多材料结构动刚度拓扑优化方法。
7.为实现上述目的，本发明提出了一种多材料结构动刚度拓扑优化方法，该方法包括如下步骤：
8.(1)将多边形单元(位移场求解单元)劈分为精细的小单元，构造设计变量与密度变量的重叠网格，形成多分辨率-多边形单元的优化建模策略。
9.(2)以平均动柔度最小化为目标和多材料的体积占比为约束，建立多材料结构的动力学拓扑优化模型。
10.(3)通过hht-α方法求解结构动响应，采用伴随变量法推导目标函数和约束的灵敏度表达式。
11.(4)利用基于敏度分离技术的zpr方法构建多区域体积约束问题的优化迭代格式。
12.进一步地，步骤(1)中在mtop框架内，为了精确计算多边形位移单元的单元刚度矩阵，将形函数及其梯度的积分点设置在密度变量所在位置(也是设计变量网格的中心)。对于多材料问题，单元刚度矩阵和质量矩阵分别表示为：
[0013][0014][0015]
式中，nn为多边形位移单元积分点的个数,y
l,ij
为密度变量在积分点处的估计值,b
l
为应变矩阵,d0为线弹性材料本构矩阵,a
l,i
为设计变量网格或密度变量网格的面积,ρ0为材料密度,和为m种材料的刚度与体积差值函数。
[0016]
基于门槛投影函数，体积插值函数表示为：
[0017][0018][0019]
式中,是门槛密度,为投影控制参数,ε＜＜1为ersatz数,用于阻止引起的数值失稳.
[0020]
刚度插值函数是由密度惩罚函数mw(y
l,ij
)与多材料插值函数mm(y
l,ij
)合成形成.首先,通过惩罚函数将单元密度变量推向0或1，基于ramp材料插值模型构造密度惩罚函数为:
[0021][0022]
式中,p0＞0为惩罚参数。
[0023]
然后，构造多材料刚度插值函数：
[0024][0025]
式中,为第i种材料的弹性模。
[0026]
由此，根据上述单元刚度与质量矩阵，多边形位移单元的总体矩阵表示为：
[0027][0028]
为了抑制棋盘格与孤岛现象，利用线性滤波方法获得网格无关的优化结果，则有：
[0029]
[0030]
式中,si为相应于密度变量单元i所占的子域,xn为设计变量dn的中心坐标。
[0031]
线性权重函数定义为：
[0032][0033]
式中,r
ni
为单元密度单元i和n的中心距,r
min
为指定的过滤半径。
[0034]
进一步地，在步骤(2)中以平均动柔度最小化为目标和多材料的体积占比为约束，建立多材料结构的动力学拓扑优化模型。假设设计域中包含m种材料，按弹性模量高低线性排列，在有限材料约束下，其动刚度优化模型为：
[0035][0036]
式中,fi为t＝ti时的动载荷向量，ui、、为相应的结构位移、速度、加速度响应，c＝αrm βrk为阻尼矩阵(αr、βr为瑞利阻尼参数)，εj、ηj、χj分别为单元指标集、设计变量指标集和多材料相数指标集，为第j种材料占设计域总体积的体积分数。
[0037]
进一步地，在步骤(3)中通过hht-α方法将步骤(2)优化模型中的半离散形式的有限元方程修改为：
[0038][0039]
通过newmark-β有限差分关系，位移、速度场的更新格式为：
[0040][0041][0042]
式中，β＝(1 α)2/4，γ＝(1 2α)/2为算法参数，合理选择参数α保证算法具有至少二阶精度和无条件稳定。
[0043]
离散形式的控制方程的残差表示为：
[0044][0045][0046]
根据步骤(2)中的优化模型，在mtop框架目标函数对设计变量的灵敏度表示为：
[0047][0048]
使用伴随变量完成目标函数的敏度分析，将newmark-β有限差分关系写成残差形式：
[0049][0050][0051]
引入伴随变量ξi、μi、νi，i＝1,2,...,n
t
，拉格朗日函数其灵敏度表示为：
[0052][0053]
然后，求解伴随变量：
[0054][0055][0056]
计算拉格朗日函数与约束函数的灵敏度，公式为：
[0057][0058][0059]
进一步地，在步骤(4)中，利用凸近似方法，确定步骤(2)中优化模型的近似子问题，表示为：
[0060][0061]
式中，和是第k次迭代时的上界与下界，move为移动限，a＝{a
l,nj
}≥0和b＝{b
l,nj
}≥0为系数矩阵。
[0062]
将目标函数灵敏度分解为正、负两项，结合则有：
[0063]
[0064][0065]
式中，为目标函数的hessian矩阵，可通过psb准牛顿法估计。
[0066]
上述子问题的拉格朗日函数为：
[0067][0068][0069]
上式表明，每一个设计变量仅与一个约束函数相关，因而拉格朗日函数l可采用分离变量的方法求极值。鉴于li的可分离特性，其一阶最优条件为：
[0070][0071][0072]
式中，为子问题设计变量的最优点。由此，满足盒子约束作用下的设计变量的更新格式为：
[0073][0074]
li(d(λi),λi)的驻点条件表示为：
[0075]
附图说明
[0076]
图1为本发明实施例提供的多材料结构动刚度拓扑优化方法流程图。
[0077]
图2为本发明多分辨率-多边形单元建模策略构型设计图。
[0078]
图3为本发明实施例施加的半正弦动载荷图。
[0079]
图4为本发明悬臂梁加载工况示意图。
[0080]
图5为本发明悬臂梁两种材料的拓扑优化结果。
[0081]
图6为本发明悬臂梁六种材料的拓扑优化结果。
[0082]
图7为本发明悬臂梁十种材料的拓扑优化结果。
[0083]
图8为本发明l型梁加载工况示意图。
[0084]
图9为本发明l型梁两种材料的拓扑优化结果。
[0085]
图10为本发明l型梁六种材料的拓扑优化结果。
[0086]
图11为本发明l型梁十种材料的拓扑优化结果。
具体实施方式
[0087]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0088]
本发明提出了一种多材料结构动刚度拓扑优化方法，该方法包括如下步骤：
[0089]
(1)将多边形单元(位移场求解单元)劈分为精细的小单元，构造设计变量与密度变量的重叠网格，形成多分辨率-多边形单元的优化建模策略。
[0090]
在mtop框架内，为了精确计算多边形位移单元的单元刚度矩阵，将形函数及其梯度的积分点设置在密度变量所在位置(也是设计变量网格的中心)。对于多材料问题，单元刚度矩阵和质量矩阵分别表示为：
[0091][0092][0093]
式中，nn为多边形位移单元积分点的个数,y
l,ij
为密度变量在积分点处的估计值,b
l
为应变矩阵,d0为线弹性材料本构矩阵,a
l,i
为设计变量网格或密度变量网格的面积,ρ0为材料密度,和为m种材料的刚度与体积差值函数。
[0094]
基于门槛投影函数，体积插值函数表示为：
[0095][0096][0097]
式中,是门槛密度,为投影控制参数,ε＜＜1为ersatz数,用于阻止引起的数值失稳.
[0098]
刚度插值函数是由密度惩罚函数mw(y
l,ij
)与多材料插值函数mm(y
l,ij
)合成形成.首先,通过惩罚函数将单元密度变量推向0或1，基于ramp材料插值模型构造密度惩罚函数为:
[0099][0100]
式中,p0＞0为惩罚参数。
[0101]
然后，构造多材料刚度插值函数：
[0102][0103]
式中,为第i种材料的弹性模。
[0104]
由此，根据上述单元刚度与质量矩阵，多边形位移单元的总体矩阵表示为：
[0105][0106]
为了抑制棋盘格与孤岛现象，利用线性滤波方法获得网格无关的优化结果，则有：
[0107][0108]
式中,si为相应于密度变量单元i所占的子域,xn为设计变量dn的中心坐标。
[0109]
线性权重函数定义为：
[0110][0111]
式中,r
ni
为单元密度单元i和n的中心距,r
min
为指定的过滤半径。
[0112]
(2)以平均动柔度最小化为目标和多材料的体积占比为约束，建立多材料结构的动力学拓扑优化模型。
[0113]
以平均动柔度最小化为目标和多材料的体积占比为约束，建立多材料结构的动力学拓扑优化模型。假设设计域中包含m种材料，按弹性模量高低线性排列，在有限材料约束下，其动刚度优化模型为：
[0114][0115]
式中,fi为t＝ti时的动载荷向量，ui、、为相应的结构位移、速度、加速度响应，c＝αrm βrk为阻尼矩阵(αr、βr为瑞利阻尼参数)，εj、ηj、χj分别为单元指标集、设计变量指标集和多材料相数指标集，为第j种材料占设计域总体积的体积分数。
[0116]
(3)通过hht-α方法求解结构动响应，采用伴随变量法推导目标函数和约束的灵敏度表达式。
[0117]
通过hht-α方法将步骤(2)优化模型中的半离散形式的有限元方程修改为：
[0118][0119]
通过newmark-β有限差分关系，位移、速度场的更新格式为：
[0120][0121][0122]
式中，β＝(1 α)2/4，γ＝(1 2α)/2为算法参数，合理选择参数α保证算法具有至少二阶精度和无条件稳定。
[0123]
离散形式的控制方程的残差表示为：
[0124][0125][0126]
根据步骤(2)中的优化模型，在mtop框架目标函数对设计变量的灵敏度表示为：
[0127][0128]
使用伴随变量完成目标函数的敏度分析，将newmark-β有限差分关系写成残差形式：
[0129][0130][0131]
引入伴随变量ξi、μi、νi，i＝1,2,...,n
t
，拉格朗日函数其灵敏度表示为：
[0132][0133]
然后，求解伴随变量：
[0134][0135][0136]
计算拉格朗日函数与约束函数的灵敏度，公式为：
[0137][0138]
[0139]
(4)利用基于敏度分离技术的zpr方法构建多区域体积约束问题的优化迭代格式。
[0140]
利用凸近似方法，确定步骤(2)中优化模型的近似子问题，表示为：
[0141][0142]
式中，和是第k次迭代时的上界与下界，move为移动限，a＝{a
l,nj
}≥0和b＝{b
l,nj
}≥0为系数矩阵。
[0143]
将目标函数灵敏度分解为正、负两项，结合则有：
[0144][0145][0146]
式中，为目标函数的hessian矩阵，可通过psb准牛顿法估计。
[0147]
上述子问题的拉格朗日函数为：
[0148][0149][0150]
上式表明，每一个设计变量仅与一个约束函数相关，因而拉格朗日函数l可采用分离变量的方法求极值。鉴于li的可分离特性，其一阶最优条件为：
[0151][0152][0153]
式中，为子问题设计变量的最优点。由此，满足盒子约束作用下的设计变量的更新格式为：
[0154][0155]
li(d(λi),λi)的驻点条件表示为：
[0156]