基于IGA各向异性多材料导热拉胀耦合微结构拓扑优化方法

本发明属于计算机辅助工程中的结构优化设计领域，具体涉及一种基于等几何分析方法(isogeometric analysis，iga)各向异性多材料的导热拉胀耦合微结构拓扑优化方法。

背景技术：

1、负泊松比超材料具有独特的力学特性，如高吸能能力、高压缩韧性等，在减震、防弹、生物医学和能源等领域有着广泛的应用潜力。目前，人工制造负泊松比超材料的主要方法是通过微观结构进行设计，可以创造出具有负泊松比效应和其他特定性能的材料。另外，各向异性复合材料具有轻质、高强度、高导热性和出色的设计能力等优点，而多材料的加入可以使得在微结构设计中的材料形态和性能更为多样化，更好的满足工程实际的需要。采用结构拓扑优化对各向异性多材料微结构进行高导热拉胀耦合轻量化设计能解决单独设计负泊松比超材料时导热性能较弱的问题，拓展了负泊松比超材料的应用范围，为材料科学的研究和创新开辟了新的途径。

2、结构拓扑优化是在给定的几何设计空间内对设计人员设置的定义规则集优化材料的布局及结构的过程。目标是通过对设计范围内的外力、荷载条件、边界条件、约束以及材料属性等因素进行数学建模和优化，从而最大限度的提高零件的性能。拓扑优化理论中需要求得某些未知场的精确解，因此需要结合数值计算方法，如有限元法(finite elementmethod，fem)、无网格法(meshless methods，mms)、有限体积法(finite volume method，fvm)、有限差分法(finite difference method，fdm)等。然而这些传统数值计算方法存在一些局限性，比如单元之间的连续性低、分析模型和几何模型之间的脱节问题以及高阶单元的效率低，并且由于这些数值方法中的离散模型难以很好地拟合复杂曲面边界上的几何模型，因此很难保证曲面边界上离散模型的光滑性。这些限制通过一种称为等几何分析(isogeometric analysis，iga)的新方法成功解决。iga作为一种数值计算方法，直接使用非均匀有理b样条(non-uniform rational b-spline，nurbs)函数作为分析模型的形状函数，其精确几何表示融合了拓扑优化和形状优化中的分析和设计模型，实现了分析模型和几何模型的统一，所得到的计算结果更加精确，并且iga的样条函数不需要额外的滤波器就可以避免拓扑优化中的棋盘格现象，得到的拓扑构型更加清晰。目前基于iga的拓扑优化集中在静力学、传热学领域，而且通常仅考虑单物理场作用下的结构优化设计，在设计负泊松比超材料时往往忽略了它的导热性能，材料种类大多仅考虑了单材料，结构设计域也局限于各向同性结构，基于iga和各向异性多材料的高导热拉胀耦合微结构拓扑优化尚未见公开报道。

技术实现思路

1、负泊松比超材料独特的“拉胀”特性有利于增强结构的压缩韧性、抗震性和降噪性等，可以在某些情况下替代传统材料，从而提高结构的性能，这种特性使得负泊松比超材料的应用前景十分广阔。另外，各向异性材料的研究和使用可以提高制造工艺和性能要求，而引入多材料可以更加有效地利用材料的优势和特性，创造出性能更为卓越的新型材料以满足现代工业的需要。为解决负泊松比超材料结构的应用仍受到热导率较低的限制的问题，本发明提出基于iga各向异性多材料的导热拉胀耦合微结构拓扑优化方法，根据固体各向同性材料惩罚(solid isotropic material with penalization，simp)模型引入一种假想的相对密度在0～1之间可变的材料并选择设计域中等几何离散控制点的相对密度作为设计变量来构造相对密度场，利用均匀化理论分别预测材料等效弹性矩阵和等效传热系数矩阵，以实现负泊松比超材料“拉胀”特性和实现高导热性能的加权函数为优化目标，以结构的总体积为约束条件，结合多材料交替活跃相(alternating active phase，aap)算法，建立各向异性多材料的高导热拉胀耦合微结构拓扑优化模型，并编写算法程序针对不同的权系数和不同各向异性材料参数的各向异性多材料求出最优的高导热拉胀耦合微结构。

2、本发明解决其技术问题所采用的技术方案是：基于iga各向异性多材料的导热拉胀耦合微结构拓扑优化方法，利用iga离散各向异性多材料微结构的静力学和稳态传热控制方程，基于能量均匀化理论(homogenization method，hm)分别预测了材料的等效弹性矩阵和等效传热系数矩阵来评估结构的等效属性并对多材料微结构施加周期性边界条件，结合多材料aap算法，建立各向异性多材料的导热拉胀耦合微结构拓扑优化数学模型，通过权系数控制各向异性微结构综合性能中负泊松比与传热性能的比重，通过调整泊松比因子控制各向异性微结构的负泊松比，通过调整热导率因子控制各向异性微结构的传热性能。

3、本发明所述技术方案的具体实施步骤如下：

4、(1)根据实际工程中结构的性能要求，确定导热拉胀耦合微结构设计域、体积约束和权系数，输入各向异性多材料微结构的材料种类数m，各种材料的弹性模量e1、泊松比v12、剪切模量g12、导热系数泊松比因子bt、热导率因子ht，材料方向角θ，对多材料微结构施加周期性边界条件，根据设计域的iga参数空间节点信息、iga控制点信息和给定的边界条件，给定负泊松比目标函数最大值及最小值和传热性能目标函数最大值及最小值；

5、(2)根据所输入iga参数空间节点信息和iga控制点信息离散各向异性多材料微结构设计域，划分iga单元和相对应的iga单元控制点；

6、(3)根据多材料交替活跃相算法，设材料种类数为m，选取第一相和第二相材料，记a＝1，b＝2，形成多材料二元相拓扑优化子模型；

7、(4)构造nurbs基函数：给定一组在一维参数空间中非递减的由实数序列组成的节点矢量ξ＝{ξ1，ξ2，…，ξn+p+1}，第i个节点用ξi表示，而用于形成b样条曲线的基函数的数量表示为n，样条曲线的阶数为p，基于这些参数，可以通过cox-de boor公式递归地计算b样条基函数：

8、

9、

10、通过将单个权重ωi添加到第i个相应的iga控制点来定义nurbs基函数：

11、

12、在二维参数空间中ξ方向定义节点矢量ξ＝{ξ1，ξ2，…，ξn+p+1}，在η方向定义节点矢量ψ＝{η1，η2，…，ηm+q+1}，则二维nurbs基函数可以构造为：

13、

14、式中，p和q分别代表二维nurbs基函数在ξ方向和η方向的阶次，n和m分别表示ξ方向和η方向的基函数个数；

15、(5)基于iga理论求iga控制点密度，根据nurbs函数插值材料i(i＝1，2，3，...，m)在iga单元内任意点相对密度

16、

17、式中，和为材料i在第i个iga控制点的nurbs函数和相对密度；

18、(6)基于iga和能量均匀化理论求解材料的等效弹性矩阵：(a)对iga控制点进行分类，选择iga控制点并施加位移约束；(b)根据各向异性的力学性能求解各向异性微结构的弹性矩阵

19、

20、式中，其中，e1、e2以及v12、v21分别为材料坐标系中方向上的弹性模量及泊松比，且满足关系式g12为剪切弹性模量；(c)经过simp插值模型得到iga单元弹性矩阵

21、

22、式中，是给定材料i的弹性矩阵，pf为材料惩罚因子；(d)建立iga单元设计域ω内的iga几何矩阵并求iga单元刚度矩阵

23、

24、式中，ω为iga单元面积，bf为力学问题中的iga几何矩阵，其表达式为

25、

26、式中，rf为iga形函数矩阵，其表达式为

27、

28、(e)给定单元初始测试应变后，求解iga控制点位移值，基于iga和能量均匀化理论预测材料的等效弹性矩阵dh；

29、(7)基于iga和能量均匀化理论求解材料的等效传热系数矩阵：(a)对iga控制点进行分类，选择iga控制点并施加温度约束；(b)根据各向异性的传热性能求解各向异性微结构的导热系数矩阵

30、

31、式中，和为材料坐标系中两个方向上的导热系数，γ为转换矩阵；(c)经过simp插值模型得到iga单元传热矩阵

32、

33、式中，是给定材料i的传热矩阵；(d)建立iga单元的iga几何矩阵并求iga单元温度刚度矩阵

34、

35、式中，bt为传热问题中的iga几何矩阵，其表达式为

36、

37、式中，rt为iga形函数矩阵，其表达式为

38、

39、(e)给定单元初始测试热应变后，求解iga控制点温度值，基于能量均匀化理论预测材料的等效传热系数矩阵kh；

40、(8)基于iga和各向异性多材料的导热拉胀耦合微结构拓扑优化的数学模型为

41、式中，为材料i在第i个iga控制点处的相对密度，n为iga控制点数，m为材料种类数，f(ρ)为加权目标函数，w为权系数，fp(ρ)为以等效弹性张量的松弛组合形式表示的负泊松比目标函数，fpmax和fpmin分别是单目标负泊松比微结构拓扑优化下目标函数的最大值和最小值，ft(ρ)为以等效传热系数矩阵主分量之和表示的传热性能目标函数，ftmax和ftmin分别是单目标传热微结构拓扑优化下目标函数的最大值和最小值，kf、u和f分别为各向异性多材料微结构iga全局位移刚度矩阵、全局位移向量和全局周期性力载荷，kt、t和p分别为各向异性多材料微结构iga全局温度刚度矩阵、全局温度向量和全局周期性温度载荷，ne为iga离散单元数，p和q分别为iga参数空间两个方向上nurbs基函数的阶次，为iga单元初始面积，ρmin＝0.001为iga控制点相对密度最小值，vi及v0分别表示材料i优化后的材料体积和材料初始总体积，为规定的材料i的全局体积分数约束；基于iga和各向异性多材料的导热拉胀耦合微结构拓扑优化的数学模型的详细步骤为：(a)据步骤(6)输出的等效弹性矩阵dh求解各向异性多材料微结构的负泊松比目标函数fp(ρ)；(b)根据步骤(7)输出的等效传热系数矩阵kh求解各向异性多材料微结构的传热性能目标函数ft(ρ)；(c)运用归一化理论对负泊松比目标函数和传热性能目标函数进行归一化处理，并对两者进行加权，以所得到的加权函数作为优化目标函数并进行灵敏度分析；(d)输出负泊松比目标函数fp(ρ)、传热性能目标函数ft(ρ)和加权目标函数值f(ρ)；

42、(9)利用oc法更新设计变量：输入当前iga控制点相对密度，根据oc准则更新iga控制点相对密度并求更新后的设计域的总体积，并由更新前后的总体积差来设定新的插值点以判断是否迭代终止，若不终止则采用更新后的iga控制点相对密度并依据oc准则继续迭代，若迭代终止则停止计算并输出更新的iga控制点相对密度；

43、(10)在多材料二元相拓扑优化子模型中判断b＝m是否成立，如果不成立，则设置b＝b+1，形成新的多材料二元相拓扑优化子模型，求解并更新设计变量，如果b＝m成立，则判断a＝m-1是否成立，如果不成立，则设置a＝a+1及b＝a+1，形成另一个新的多材料二元相拓扑优化子模型，求解并更新设计变量，最后，如果a＝m-1成立，判断是否满足终止迭代条件，若满足，输出拓扑优化结果，若不满足，跳到步骤(3)继续进行迭代，直到满足迭代终止条件为止。

44、本发明的有益效果是：本发明基于等几何分析方法和固体各向同性材料惩罚模型的拓扑优化技术较好地消除了基于常见拓扑优化方法获得的拓扑结构易出现锯齿、棋盘格和中间密度等数值不稳定性问题；本发明基于等几何分析方法，以nurbs函数作为分析模型的形状函数，实现分析模型和几何模型的统一，避免了传统分析方法的空间离散化误差，并大大提高了高阶单元的计算效率，求得的计算结果也更加精确；本发明得到的拓扑构型清晰，结构边界光滑，便于后续的加工制造；本发明通过调整权系数、泊松比因子和热导率因子控制各向异性多材料微结构的负泊松比与传热性能，可根据工程实际需求得到相应的各向异性多材料微结构，计算流程简单；本发明实现了负泊松比超材料和高导热性能材料之间的平衡，能根据不同设计需求得到具有高导热性能的各向异性多材料负泊松比微结构，对于多载荷工况下的材料设计具有重要理论研究意义和实际应用价值。