基于无网格EFGM的各向异性多材料热力耦合结构拓扑优化方法

free galerkin method，efgm)、重构核粒子法(reproducing kernel particle method，rkpm)、无网格局部petrov-galerkin法(meshless local petrov-galerkin method，mlpg)等。其中，无网格伽辽金法efgm是一种比较成熟的无网格方法，在力学、热学等拓扑优化领域已经得到一些应用，且无需使用敏度过滤技术也能表现出良好的拓扑收敛性。但目前基于efgm的拓扑优化研究多数只考虑了单物理场载荷作用下的各向同性单一材料的分布优化，而忽略了实际的工程案例通常处在复杂的多物理场载荷环境中，如热力耦合等，其中应用了更加高效的新型各向异性多相复合材料来进行多目标拓扑优化的案例更是鲜有。弥补上述不足就需要工程技术人员在提出产品设计方案的初期充分考虑到参与设计的多材料的种类、比例、分布形式、以及各组分材料各向异性的强度和方向对结构优化设计造成的影响，通过寻找到相对经济和最优的多材料匹配方式，来实现最终产品的多功能和高性能目标。


技术实现要素：

5.目前，各向异性多相复合材料已在众多工程领域得到了广泛应用，成为了能替代传统各向同性单一材料的新选择。为了将各向异性多相复合材料引入到热力耦合多目标拓扑优化中，解决网格类数值方法在进行拓扑优化时产生的数值不稳定现象及最优拓扑结构可能由于过滤不当而产生大量中间密度或陷入局部最小值等问题，本发明提出一种基于无网格efgm的正交各向异性多材料热力耦合结构拓扑优化方法，其根据固体各向同性惩罚(solid isotropic material with penalization，simp)材料插值方案引入多种假想的相对密度在0～1之间可变的材料，并选择设计域中无网格efgm节点相对密度作为设计变量来构造相对密度场，以最小柔度和最小散热弱度的加权目标函数作为多目标拓扑优化的目标函数，以多材料体积分数为约束条件，建立正交各向异性多材料热力耦合结构拓扑优化问题的无网格efgm数学模型，并编写程序，输出在不同多材料种类数、多材料体积分数、次弹性模量、剪切模量、次泊松比、次导热系数、泊松比因子、热导率因子、力热材料方向角、权系数等下的最优正交各向异性多材料热力耦合结构拓扑优化设计。
6.本发明解决其技术问题所采用的技术方案是：基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法，利用无网格efgm理论离散各向异性多材料结构的静力学和稳态传热控制方程，并结合多材料simp插值方案建立多材料热力耦合结构拓扑优化的数学模型；确定引入的多材料种类数，通过调整多材料体积分数来控制多材料组分比例，通过调整多材料的次弹性模量、泊松比因子和力学材料方向角控制正交各向异性多材料结构的力学性能，通过调整多材料的次导热系数、热导率因子和传热材料方向角控制正交各向异性多材料结构的传热性能，通过调整权系数控制正交各向异性多材料结构综合性能中力学性能与传热性能的侧重比重，以便于实施不同需求下正交各向异性多材料热力耦合结构拓扑优化设计；假设选取的正交各向异性多材料种类数为q，在总数为r
max
的每一步一级迭代下，不重复执行次第a、b种材料参与的双材料竞争优化二级迭代，其中a＝1,2,...,q-1，b＝a+1,a+2,...,q；每次二级迭代下又执行t
max
次三级迭代并更新第a、b种材料的节点相对密度，t
max
根据需要被设置为1、2、3或更大整数；当满足迭代终止条件时，输出正交各向异性多材料热力耦合结构拓扑优化最优结果。
7.本发明所述技术方案包括以下步骤：
8.(1)根据实际工程结构设计要求输入正交各向异性多材料种类数q、多材料体积分数vi(i＝1,2,...,q)、次弹性模量剪切模量次泊松比次导热系数泊松比因子prf、热导率因子tcf、力热材料方向角θf、θ
p
、一级、二级、三级迭代最大步数r
max
、s
max
、t
max
、一级迭代步数初始值r＝1、最小容差ch
min
、容差初始值ch、惩罚力度p；
9.(2)确定正交各向异性多材料热力耦合结构设计域ω，以无网格efgm节点离散ω，根据节点信息布置ω内部及边界高斯点，确定施加在ω上的力载荷、位移边界、热载荷、三类传热边界条件，赋予多材料节点相对密度矩阵初始值给定权系数ω、柔度最大值c
max
及最小值c
min
、散热弱度最大值s
max
及最小值s
min
；
10.(3)设置二级迭代步数初始值s＝1、双材料竞争优化初始值a＝1、b＝a+1；
11.(4)设置三级迭代步数初始值t＝1；
12.(5)建立正交各向异性多材料结构力学性能变换矩阵以使全局坐标系与材料坐标系下的弹性矩阵相关联，求解正交各向异性多材料结构第i种材料在全局坐标系下的弹性矩阵
13.(6)建立正交各向异性多材料结构传热性能变换矩阵以使全局坐标系与材料坐标系下的导热系数矩阵相关联，求解正交各向异性多材料结构第i种材料在全局坐标系下的传热矩阵
14.(7)基于无网格efgm理论，并根据多材料simp插值方案引入q种相对密度在0～1间可变的材料，以多材料节点相对密度为设计变量构造相对密度场，建立第i种材料第k次分级迭代前后高斯点材料属性关系；
15.(8)建立正交各向异性多材料结构静力问题的无网格efgm离散控制方程，并求解正交各向异性多材料结构的节点参数位移u、节点近似位移柔度
16.(9)建立正交各向异性多材料结构稳态传热问题的无网格efgm离散控制方程，并求解各向异性多材料结构的节点参数温度t、节点近似温度散热弱度
17.(10)以多材料体积分数为约束条件，以柔度和散热弱度的加权目标函数最小为优化目标，建立基于无网格efgm和多材料simp插值方案的正交各向异性多材料热力耦合结构拓扑优化数学模型；
18.(11)求解无网格efgm正交各向异性多材料热力耦合结构拓扑优化模型中第i种材料的体积灵敏度，采用伴随法求解柔度灵敏度、散热弱度灵敏度、加权目标函数灵敏度；
19.(12)将体积灵敏度、加权目标函数灵敏度代入优化准则(optimality criteria，oc)，根据oc更新第k次分级迭代下第a种材料节点相对密度根据更新前后第a、b种材料节点相对密度之和保持不变的关系求解本次迭代下第b种材料节点相对密度组装多材料节点相对密度矩阵
20.(13)比较t与t
max
大小关系，判断三级迭代为继续或结束；比较s与s
max
大小关系，判断二级迭代为继续或结束；比较r与r
max
及ch与ch
min
大小关系，判断一级迭代为继续或结束；
21.(14)输出无网格efgm正交各向异性多材料热力耦合结构拓扑优化的节点灰度图作为最优拓扑结构，其中多材料按照力热性能从高到低的排序分别表示为由深到浅的灰度
颜色。
22.本发明相比现有技术的有益效果是：本发明能更高效灵活地采用无网格efgm节点离散设计域，易于构造高阶形函数，计算精度高，且无需使用灵敏度过滤和密度过滤技术，减少了拓扑边界上中间密度节点的数量，避免了基于网格类数值方法的拓扑优化中存在的棋盘格、网格依赖性、局部极小值等数值不稳定现象，所求解的正交各向异性多材料热力耦合结构边界清晰光滑，具有更高的可制造性；本发明选择无网格efgm节点相对密度作为设计变量，多目标拓扑优化的加权目标函数迭代过程不会出现振荡，收敛性好，避免了采用efgm高斯计算点相对密度作为设计变量时相对密度场和目标函数灵敏度的不连续性造成的数值不稳定现象；本发明通过调整多材料体积分数能控制组分材料用料比例，通过调整次弹性模量、剪切模量、次泊松比、泊松比因子、力学材料方向角能控制正交各向异性多材料结构的力学承载性能，通过调整次导热系数、热导率因子、传热材料方向角能控制正交各向异性多材料结构的散热性能，通过调整权系数能控制正交各向异性多材料热力耦合结构力学承载性能和散热性能所占比重；本发明能与工程实践紧密结合，能简单方便地根据不同工程需求实施对应的正交各向异性多材料热力耦合结构拓扑优化设计，得到力学承载和散热综合性能最优的正交各向异性多材料拓扑结构，可操作性强，在复杂工程结构创新设计中具有重要的理论研究意义和应用价值。
附图说明
23.图1是本发明正交各向异性多材料热力耦合结构拓扑优化设计流程图；
24.图2是本发明正交各向异性多材料热力耦合结构的全局坐标系x-y与材料坐标系1-2；
25.图3是本发明正交各向异性多材料热力耦合结构实施例的力载荷、位移边界、热载荷、三类传热边界示意图；
26.图4是本发明采用301
×
61个均匀分布的efgm节点离散实施例设计域的示意图；
27.图5是本发明根据节点信息在设计域每个积分背景网格内布置2
×
2高斯积分点的示意图；
28.图6～图12是本发明实施例中权系数ω分别为0、16、13、12、23、56、1，泊松比因子prf和热导率因子tcf均为14，力热材料方向角θf、θ
p
均为0
°
时基于无网格efgm的正交各向异性多材料热力耦合最优拓扑结构；
29.图13～图19是本发明实施例中权系数ω为710，泊松比因子prf和热导率因子tcf分别为14、13、12、1、2、3、4，力热材料方向角θf、θ
p
均为0
°
时基于无网格efgm的多材料热力耦合最优拓扑结构，tcf为14、13、12、2、3、4表示所应用的多材料为正交各向异性，tcf为1表示所应用的多材料为各向同性；
30.图20～图26是本发明实施例中权系数ω为710，泊松比因子prf和热导率因子tcf均为52，力热材料方向角θf、θ
p
分别为0
°
、15
°
、30
°
、45
°
、60
°
、75
°
、90
°
时基于无网格efgm的正交各向异性多材料热力耦合最优拓扑结构。
具体实施方式
31.下面结合附图和实施例，对本发明的具体实施方式作进一步详细描述，以使本发
明的目的、内容和优点更加清楚。
32.基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法主要包括以下步骤：
33.(1)根据实际工程结构设计要求输入正交各向异性多材料种类数q、多材料体积分数vi＝[v1,v2,...,vq]
t
(i＝1,2,...,q)、多材料在图2所示材料坐标系2方向上的次弹性模量剪切模量次泊松比次导热系数泊松比因子prf、热导率因子tcf、力热材料方向角θf、θ
p
、一级、二级、三级迭代最大步数r
max
、s
max
、t
max
、一级迭代步数初始值r＝1、最小容差ch
min
、容差初始值ch、惩罚力度p；
[0034]
(2)确定正交各向异性多材料热力耦合结构设计域ω，以无网格efgm节点离散ω，根据节点信息布置ω内部及边界高斯点，确定施加在ω上的力载荷、位移边界、热载荷、三类传热边界条件，赋予多材料节点相对密度矩阵初始值其中n为efgm节点总数，给定权系数ω(0≤ω≤1)、柔度最大值c
max
及最小值c
min
、散热弱度最大值s
max
及最小值s
min
；
[0035]
(3)设置二级迭代步数初始值s＝1、双材料竞争优化初始值a＝1、b＝a+1；
[0036]
(4)设置三级迭代步数初始值t＝1；
[0037]
(5)建立正交各向异性多材料热力耦合结构力学性能变换矩阵以使全局坐标系与材料坐标系下的弹性矩阵相关联，求解正交各向异性多材料热力耦合结构每种材料在全局坐标系下的弹性矩阵具体步骤如下：
[0038]
参见图2，当全局坐标系x-y与材料坐标系1-2之间存在力学材料方向角θf时，根据广义胡克定律，在线弹性范围内，正交各向异性多材料力学结构的应力和应变存在以下关系：
[0039][0040][0041][0042][0043]
式中，为正交各向异性多材料结构第i种材料在全局坐标系x-y下的弹性矩阵，qi为正交各向异性多材料结构第i种材料在材料坐标系1-2下的弹性矩阵，t-1
为将与qi相
关联的坐标转换矩阵，在qi中，中，其中分别为第i种材料在材料坐标系1-2中1和2方向上的弹性模量，分别为第i种材料在材料坐标系1-2中1和2方向上的泊松比，为第i种材料的剪切模量；定义泊松比因子来表示力学正交各向异性强度，则当prf＝1时，多材料呈现力学各向同性，通过修改泊松比因子prf和力学材料方向角θf，可以得到具有不同力学承载性能的多材料拓扑优化结构；
[0044]
(6)建立正交各向异性多材料结构传热性能变换矩阵以使全局坐标系与材料坐标系下的传热矩阵相关联，求解正交各向异性多材料结构每种材料在全局坐标系下的传热矩阵具体步骤如下：
[0045]
参见图2，当全局坐标系x-y与材料坐标系1-2之间存在传热材料方向角θ
p
时，根据傅立叶定律，正交各向异性多材料传热结构在全局坐标系下坐标轴x、y方向上的热通量与对应方向上的温度梯度存在以下关系：
[0046][0047][0048][0049][0050]
式中，为正交各向异性多材料结构第i种材料在全局坐标系x-y下的传热矩阵，λi为正交各向异性多材料结构第i种材料在材料坐标系1-2下的传热矩阵，γ为将与λi相关联的坐标转换矩阵，在λi中，分别为第i种材料在材料坐标系1-2中1和2方向上的导热系数；定义热导率因子来表示传热正交各向异性强度，则当tcf＝1时，多材料呈现传热各向同性，通过修改热导率因子tcf和传热材料方向角θ
p
，可以得到具有不同传热性能的多材料拓扑优化结构；
[0051]
(7)基于无网格efgm理论，并根据多材料simp插值模型引入q种相对密度在0～1间可变的材料，以多材料节点相对密度为设计变量构造相对密度场，建立第i种材料第k次分级迭代前后高斯点材料属性关系，k＝t
maxsmax
(r-1)+t
max
(s-1)+t为计算当前分级迭代步数的迭代索引，具体步骤如下：
[0052]
通过多材料simp插值模型建立第k次分级迭代下第i种材料的高斯点弹性模量e
i(k)
与实体材料弹性模量e
i(0)
之间的关系为
[0053][0054]
通过多材料simp插值模型建立第k次分级迭代下第i种材料的高斯点剪切模量与实体材料剪切模量之间的关系为
[0055][0056]
通过多材料simp插值模型建立第k次分级迭代下第i种材料的高斯点导热系数λ
i(k)
与实体材料导热系数λ
i(0)
之间的关系为
[0057][0058]
式中，为第i种材料惩罚后的高斯点相对密度，p为惩罚力度，通常p≥3；可由高斯点影响域内efgm节点相对密度通过移动最小二乘法(moving least squares,mls)构造的形函数近似求得，即
[0059][0060]
式中，为第i种材料高斯点影响域内第i个节点相对密度，n为高斯点影响域ω
x
内的节点数目，φi为高斯点影响域ω
x
内第i个节点的mls形函数；
[0061]
根据无网格efgm理论，mls形函数构造的具体过程如下：
[0062]
假设未知场函数为u(x)，则其在由n个节点所构成的计算域ω内的mls近似函数uh(x)可以被构造为
[0063][0064]
式中，p
t
(x)＝[p1(x),p2(x),...,pm(x)]是空间坐标的基函数pj(x)构成的向量，m是基函数的数量，a(x)＝[a1(x),a2(x),...,am(x)]
t
是未知系数aj(x)构成的向量，a(x)可以通过最小化加权l2范函j(x)求得
[0065][0066]
式中，任意计算点(高斯点)x的支撑域ω
x
覆盖了n个节点xi(i＝1,2,...n)，ui＝u(xi)为xi处的节点参数值，w(x-xi)＝wi(x)为xi处的紧支权函数，而最小化范函j(x)需要保证解此式可得
[0067]
最小化范函j(x)需要保证解此式可得
[0068]
a(x)a(x)＝b(x)u
ꢀꢀꢀꢀꢀ
(15)
[0069]
式中的a(x)、b(x)、u可分别表示如下：
[0070][0071]
b(x)＝[w1(x)p(x1),w2(x)p(x2),...,wn(x)p(xn)]
ꢀꢀꢀꢀꢀ
(17)
[0072]
u＝[u1,u2,...,un]
t
ꢀꢀꢀꢀꢀ
(18)
[0073]
将式(16)、(17)、(18)代入式(15)解得
[0074]
a(x)＝a-1
(x)b(x)u
ꢀꢀꢀꢀꢀ
(19)
[0075]
将a(x)回代入式(13)解得
[0076][0077]
式中，支撑域ω
x
内n个节点的mls形函数φ(x)可以表示为
[0078]
φ(x)＝[φ1(x),φ2(x),...,φn(x)]＝p
t
(x)a-1
(x)b(x)
ꢀꢀꢀꢀꢀ
(21)
[0079]
对于二维问题，本发明采用线性基p
t
(x)＝[1,x,y](m＝3)，采用三次样条权函数
[0080][0081]
式中，对于圆形和矩形支撑域ω
x
，wi(x)可分别表示为w(||x-xi||/di)和w(|x-xi|/d
ix
)w(|y-yi|/d
iy
)，其中，||x-xi||/di、|x-xi|/d
ix
、|y-yi|/d
iy
分别为ω
x
的正则化特征距离，||x-xi||为计算点x(x,y)到节点xi(xi,yi)的距离，di为ω
x
半径，|x-xi|、|y-yi|为x和xi在x、y方向的距离，d
ix
、d
iy
为ω
x
长度和宽度的1/2；
[0082]
为保证a(x)的可逆性和φi(x)的局部特性，di、d
ix
、d
iy
取决于xi影响域ωi的大小，且需包含足量相邻节点，由于正方形ω
x
对规则节点排布适应性良好，因此在本发明中，采用此支撑域；
[0083]
(8)建立正交各向异性多材料结构静力问题无网格efgm离散控制方程，并求解正交各向异性多材料结构的节点参数位移u、节点近似位移柔度具体步骤如下：
[0084]
(8.1)将正交各向异性多材料结构静力问题无网格efgm离散控制方程整理为
[0085][0086]
式中，为经过simp材料插值模型惩罚了弹性模量和剪切模量后的efgm整体力刚度矩阵，k
fα
为经过罚函数法处理位移边界条件后的整体力刚度矩阵惩罚项，u为参数位移列向量，f为整体力载荷列向量，f
α
为整体力载荷列向量惩罚项，k
fα
、f、f
α
的子矩阵分别为k
fαij
、fi、f
αi
，其积分式被表示为式(24)～式(27)：
[0087][0088][0089][0090][0091]
式中，为静力分析几何矩阵，α为罚函数的罚因子，b为体力向量，为给定牵引力向量，当x或y方向上施加了位移约束，s
x
或sy为1，否则，s
x
或sy为0；
[0092]
式(23)所示正交各向异性多材料结构静力问题无网格efgm离散控制方程及式(24)～式(27)所示子矩阵的具体求解方式如下：
[0093]
(8.2)遍历设计域各高斯点影响域内的节点并计算其mls形函数，引入多材料simp插值模型，对第i种材料在材料坐标系下的弹性矩阵qi中的弹性模量和剪切模量执行如式(9)和(10)所示的惩罚；每种材料对整体力刚度矩阵的贡献通过将各材料在材料坐标系下惩罚后的弹性矩阵进行求和实现，即建立节点几何矩阵b
fi
，求解节点无网格efgm力刚度矩阵组装无网格efgm整体力刚度矩阵
[0094]
(8.3)根据力载荷分布信息求解整体力载荷列向量f：当力载荷以集中力形式施加，根据集中力位置计算其影响域内节点efgm形函数及施加在每个节点上的力载荷；当力载荷以线分布力形式施加，遍寻线分布力所在线上高斯点，求各高斯点影响域内节点efgm形函数及施加在每个节点上的力载荷；将以上两种力载荷施加形式组装为整体力载荷列向量f；
[0095]
(8.4)采用罚函数法逐一处理每条力载荷和位移边界条件：输入位移边界节点信息并求解边界上的高斯点信息，根据边界上的高斯点搜索其影响域内节点并求对应节点的mls形函数，将罚因子α、mls形函数节点给定位移值相乘作为节点力载荷惩罚项f
αi
，组装成整体力载荷惩罚项列向量f
α
；将罚因子α、节点间的mls形函数φj相乘作为位移边界的力刚度矩阵惩罚项k
fαij
，组装为整体力刚度矩阵惩罚项k
fα
；
[0096]
(8.5)根据正交各向异性多材料结构无网格efgm静力离散控制方程求解efgm节点参数位移u列向量，结合整体节点形函数矩阵求解节点近似位移列向量求解柔度
[0097]
(9)建立各向异性多材料结构稳态传热问题无网格efgm离散控制方程，并求解各向异性多材料结构的节点参数温度t、节点近似温度散热弱度具体步骤如下：
[0098]
(9.1)将正交各向异性多材料结构稳态传热问题无网格efgm离散控制方程整理为
[0099][0100]
式中，为经过simp材料插值模型惩罚了导热系数后的efgm整体热刚度矩阵，k
pα
为经过罚函数法处理dirichlet边界条件和cauchy边界条件后的整体热刚度矩阵惩罚项，t为参数温度列向量，p为整体热载荷列向量，p
α
为整体热载荷列向量惩罚项，k
pα
、p、p
α
的子矩阵分别为k
pαij
、pi、p
αi
，其积分式被表示为式(29)～式(32)：
[0101][0102][0103]
[0104][0105]
式中，为稳态传热分析几何矩阵，h为对流换热系数，α为罚函数的罚因子，为单位面积热生成率，为给定dirichlet边界γ1上的恒定温度，q为给定neumann边界γ2上的热流密度，t
∞
为给定cauchy边界γ3上的外界环境温度；
[0106]
式(28)所示正交各向异性多材料结构稳态传热问题无网格efgm离散控制方程及式(29)～式(32)所示子矩阵的具体求解方式如下：
[0107]
(9.2)遍历设计域各高斯点影响域内的节点并计算其mls形函数，引入多材料simp插值方案，对第i种材料在材料坐标系下的传热矩阵λi中的导热系数执行如式(11)所示的惩罚，每种材料对整体热刚度矩阵的贡献通过将各材料在材料坐标系下惩罚后的传热矩阵进行求和实现，即建立节点几何矩阵b
pi
，求解节点无网格efgm热刚度矩阵
[0108]
(9.3)采用罚函数法逐一处理每条dirichlet恒定温度边界γ1：输入dirichlet边界γ1节点信息并求解边界γ1上的高斯点信息，根据边界γ1上的高斯点搜索其影响域内节点并求对应节点mls形函数，将罚因子α、节点给定恒定温度值mls形函数φi相乘作为边界γ1的热载荷施加量将罚因子α、节点间的mls形函数φi、φj相乘作为dirichlet边界γ1的热刚度矩阵惩罚项
[0109]
(9.4)逐一处理每条neumann热流密度边界γ2：输入neumann边界γ2上的节点信息并求解边界γ2上的高斯点信息，根据边界γ2上高斯点搜索其影响域内节点并求对应节点mls形函数，将热流密度q、节点mls形函数φi相乘作为cauchy边界γ2的热载荷施加量
[0110]
(9.5)逐一处理每条cauchy对流换热边界γ3：输入cauchy边界γ3节点信息并求解边界γ3上的高斯点信息，根据边界γ2上的高斯点搜索其影响域内节点并求对应节点的mls形函数，将对流换热系数h、周围环境温度t
∞
、节点mls形函数φi相乘作为cauchy边界γ3的热载荷施加量将节点间的mls形函数φi、φj与对流换热系数h相乘作为cauchy边界γ3的热刚度矩阵
[0111]
(9.6)根据热源分布信息求由热源引起热载荷：当热源以点热源形式施加，根据点热源位置计算其影响域内节点mls形函数及施加在每个节点上的热载荷；当设计域内热源均匀分布，遍寻均布热源所在域内高斯点，求高斯点影响域内节点mls形函数及施加在每个节点上的热载荷，将上两种热源施加形式组装总体热载荷列向量
[0112]
(9.7)根据正交各向异性多材料结构无网格efgm稳态传热离散控制方程
求解efgm节点参数温度列向量t，结合整体节点形函数矩阵求解节点近似温度列向量求解散热弱度
[0113]
(10)以多材料节点相对密度为设计变量，以多材料体积分数vi为约束条件，以柔度和散热弱度的加权目标函数最小为优化目标，建立基于无网格efgm和多材料simp插值方案的正交各向异性多材料热力耦合结构拓扑优化数学模型为：
[0114][0115]
式中，为第i种材料的节点相对密度行向量，为第i种材料的第j个节点相对密度，为多材料节点相对密度矩阵，q为多材料种类数，n为节点总数；为柔度，为散热弱度，c
max
和c
min
分别为单目标静力结构柔度最大值和最小值，s
max
和s
min
分别为单目标稳态传热结构散热弱度最大值和最小值；ω(0≤ω≤1)为加权目标函数权系数，为对和进行归一化处理后的加权目标函数，规定ω＝0表示单目标散热弱度最小化问题，ω＝1表示单目标柔度最小化问题；vi为第i种材料的体积，为第i种材料的高斯点相对密度，vi为多材料体积分数，v0为多材料总体积；为防止矩阵奇异，令第i种材料的高斯点相对密度和节点相对密度最小值为ρ
min
；
[0116]
(11)求解无网格efgm正交各向异性多材料热力耦合结构拓扑优化模型中第i种材料的体积灵敏度采用伴随法求解柔度灵敏度散热弱度灵敏度加权目标函数灵敏度具体步骤如下：
[0117]
无网格efgm正交各向异性多材料热力耦合结构拓扑优化模型中第i种材料的体积灵敏度为
[0118][0119]
式中，φi为第i种材料高斯点影响域内第i个节点的mls形函数；
[0120]
无网格efgm正交各向异性多材料热力耦合结构拓扑优化模型中加权目标函数灵敏度为
[0121][0122]
式(35)中，采用伴随法分别求解柔度灵敏度、散热弱度灵敏度为
[0123][0124][0125]
式(36)中，无网格efgm整体力刚度矩阵与其惩罚项叠加后的灵敏度为
[0126][0127]
式(37)中，无网格efgm整体热刚度矩阵与其惩罚项叠加后的灵敏度为
[0128][0129]
(12)将体积灵敏度加权目标函数灵敏度代入优化准则oc，根据oc更新第k次分级迭代下第a种材料的节点相对密度为根据更新前后第a、b种材料节点相对密度之和保持不变的关系求解本次迭代下第b种材料的节点相对密度为组装多材料节点相对密度矩阵更新的具体步骤如下：
[0130]
(12.1)若i＝a，则根据oc公式更新
[0131][0132]
式中，m0＝0.2为正移动限制，τ＝0.5为阻尼系数；为避免矩阵奇异，取节点相对密度下限为ρ
min
＝10-3
；根据oc准则，定义为
[0133][0134]
根据多材料结构第a种材料的节点相对密度求解更新后的第a种材料的总体积v
a(k)
以及和初始体积的差值v
a(k)-v
av0
，判断和0的大小关系，以设定新的二分法插值点上下限
并判断oc内部迭代为继续或终止，若不终止则将回代到式(40)中继续更新，若终止则输出
[0135]
(12.2)若i＝b，则根据更新前后第a、b种材料节点相对密度之和保持不变的关系，求得第k次分级迭代下第b种材料的节点相对密度为
[0136][0137]
(12.3)若i≠a且i≠b，则本次分级迭代不涉及第i种材料节点相对密度更新，即第k、k-1次迭代第i种材料节点相对密度应相同，即
[0138][0139]
(12.4)结合式(40)、(42)、(43)输出多材料整体节点相对密度计算多材料节点相对密度改变值的无穷范数
[0140]
(13)参见图1，步骤(12)和步骤(14)之间的循环判断语句执行的具体步骤如下：
[0141]
(13.1)比较t与t
max
大小关系，判断三级迭代为继续或结束：若t《t
max
，则t＝t+1，返回步骤(5)，继续执行三级迭代；若t《t
max
不成立，则本轮三级迭代结束；
[0142]
(13.2)比较s与大小关系，判断二级迭代为继续或结束：若s＜s
max
且b《q，则令s＝s+1,b＝b+1，返回步骤(4)；若s＜s
max
但b《m不成立，则令s＝s+1,a＝a+1,b＝a+1，返回步骤(4)；若s＜s
max
不成立，则本轮二级迭代结束；
[0143]
(13.3)比较r与r
max
及ch与ch
min
大小关系，判断一级迭代为继续或结束：若r《r
max
且ch＞ch
min
，则令r＝r+1，返回步骤(3)；若r《r
max
不成立或ch＞ch
min
不成立，则本轮一级迭代结束，执行步骤(14)；
[0144]
(14)输出无网格efgm正交各向异性多材料热力耦合结构拓扑优化的节点灰度图作为最优拓扑结构，其中多材料按照力热性能从高到低的排序分别表示为由深到浅的灰度颜色。
[0145]
下面是本发明方法应用于工程实践的一个实施例：
[0146]
参见图3，本实施例是长为300mm，高为600mm的二维梁，梁左右边界x和y方向上的位移被固定为0，在梁几何中心处施加一个竖直向下的力f＝105n，梁内施加了均匀分布的热源梁四个顶点为dirichlet恒定温度边界γ1，恒定温度为梁左、右边界中间各长度为30mm的边界为neumann热流密度边界γ2，给定的热流密度为q＝-104w/m2，梁上、下边界中间各长度为30mm的边界为cauchy对流换热边界γ3，给定的对流换热系数为h＝1.2
×
104w/(m2·
℃)，给定对流边界γ3的外界环境温度为t
∞
＝2℃，梁的上、下、左、右边界除了边界γ1、γ2、γ3以外的其他部分均为绝热边界；
[0147]
下面在给定的三类工程设计需求下进行四种材料的多材料热力耦合结构拓扑优化设计，tcf不为1表示所应用的多材料为正交各向异性，tcf为1表示所应用的多材料为各向同性；
[0148]
需求1：权系数ω分别为0、16、13、12、23、56、1，泊松比因子prf和热导率因子tcf均为14，力热材料方向角θf、θ
p
均为0
°
；
[0149]
需求2：权系数ω为710，泊松比因子prf和热导率因子tcf分别为14、13、12、1、2、3、
4，力热材料方向角θf、θ
p
均为0
°
；
[0150]
需求3：权系数ω为710，泊松比因子prf和热导率因子tcf均为52，力热材料方向角θf、θ
p
分别为0
°
、15
°
、30
°
、45
°
、60
°
、75
°
、90
°
；
[0151]
针对该实施例，根据具体实施方式的步骤(1)输入多材料种类数q＝4、多材料体积分数vi＝[0.25,0.2,0.15,0.4]
t
、多材料在图2所示材料坐标系2方向的次弹性模量次泊松比次导热系数剪切模量通过计算，根据实施例需求1、2、3输入规定的泊松比因子prf、热导率因子tcf、力热材料方向角θf、θ
p
，设置一级、二级、三级迭代的最大步数r
max
＝120、s
max
＝4
×
(4-1)2＝6、t
max
＝3、一级迭代步数初始值r＝1、最小容差ch
min
＝10-5
，容差初始值ch＝0.2、惩罚力度p＝4；
[0152]
针对该实施例，根据具体实施方式的步骤(2)考虑图3中的设计域ω尺寸信息，将ω离散为图4中的n＝301
×
61＝18361个均匀分布的efgm节点，根据节点信息布置图5中的ω内每个积分背景网格的2
×
2高斯积分点，输入施加的集中力f＝105n、边界γ1的恒定温度边界γ2的热流密度q＝-104w/m2、边界γ3的对流换热系数h＝1.2
×
104w/(m2·
℃)、边界γ3的外界环境温度t
∞
＝2℃、均布热源在位移边界、边界γ1、γ2、γ3上布置高斯积分点用于步骤(8)和步骤(9)施加力热边界条件，赋予多材料节点相对密度矩阵初始值其中，其中，给定实施例需求1、2、3中规定的权系数ω、柔度最大值c
max
及最小值c
min
、散热弱度最大值s
max
及最小值s
min
；
[0153]
针对该实施例，根据具体实施方式的步骤(3)～步骤(14)继续依次执行，直至迭代最终结束；
[0154]
图6～图26是本实施例在给定的三类工程需求下进行四种材料的多材料热力耦合结构拓扑优化设计：图6～图12是本发明实施例中权系数ω分别为0、16、13、12、23、56、1，泊松比因子prf和热导率因子tcf均为14，力热材料方向角θf、θ
p
均为0
°
时基于无网格efgm的正交各向异性多材料热力耦合最优拓扑结构；图13～图19是本发明实施例中权系数ω为710，泊松比因子prf和热导率因子tcf分别为14、13、12、1、2、3、4，力热材料方向角θf、θ
p
均为0
°
时基于无网格efgm的多材料热力耦合最优拓扑结构，tcf为14、13、12、2、3、4表示所应用的多材料为正交各向异性，tcf为1表示所应用的多材料为各向同性；图20～图26是本发明实施例中权系数ω为710，泊松比因子prf和热导率因子tcf均为52，力热材料方向角θf、θ
p
分别为0
°
、15
°
、30
°
、45
°
、60
°
、75
°
、90
°
时基于无网格efgm的正交各向异性多材料热力耦合最优拓扑结构；
[0155]
图6～图26表明多材料的竞争优化关系使得多材料节点相对密度矩阵中每个节点位置处都有一种材料的节点相对密度趋于1，而其他材料的节点相对密度趋于0，因此正交各向异性多材料热力耦合最优拓扑结构的节点最终呈现出每种材料规定的从深到浅的灰度颜色，使得本发明所得正交各向异性多材料热力耦合最优拓扑结构的多材料拓扑边
界清晰；通过修改多材料种类数q、多材料体积分数vi、次弹性模量剪切模量次泊松比次导热系数泊松比因子prf、热导率因子tcf、力热材料方向角θf、θ
p
、权系数ω等可以实现不同工程需求下的正交各向异性多材料热力耦合结构拓扑优化设计。
[0156]
虽然参考本实施例对本发明进行了详细描述，但以上所述并不限定本发明的保护范围，任何依据本发明思路下的修改和改进等，均视为本发明保护范围。

技术特征：
1.基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法，其特征在于，该方法包括以下步骤：(1)根据实际工程结构设计要求输入正交各向异性多材料种类数q、多材料体积分数v
i
(i＝1,2,...,q)、次弹性模量剪切模量次泊松比次导热系数泊松比因子prf、热导率因子tcf、力热材料方向角θ
f
、θ
p
、一级、二级、三级迭代最大步数r
max
、s
max
、t
max
、一级迭代步数初始值r＝1、最小容差ch
min
、容差初始值ch、惩罚力度p；(2)确定正交各向异性多材料热力耦合结构设计域ω，以无网格efgm节点离散ω，根据节点信息布置ω内部及边界高斯点，确定施加在ω上的力载荷、位移边界、热载荷、三类传热边界条件，赋予多材料节点相对密度矩阵初始值给定权系数ω、柔度最大值c
max
及最小值c
min
、散热弱度最大值s
max
及最小值s
min
；(3)设置二级迭代步数初始值s＝1、双材料竞争优化初始值a＝1、b＝a+1；(4)设置三级迭代步数初始值t＝1；(5)建立正交各向异性多材料结构力学性能变换矩阵以使全局坐标系与材料坐标系下的弹性矩阵相关联，求解正交各向异性多材料结构第i种材料在全局坐标系下的弹性矩阵(6)建立正交各向异性多材料结构传热性能变换矩阵以使全局坐标系与材料坐标系下的导热系数矩阵相关联，求解正交各向异性多材料结构第i种材料在全局坐标系下的传热矩阵(7)基于无网格efgm理论，并根据多材料simp插值方案引入q种相对密度在0～1间可变的材料，以多材料节点相对密度为设计变量构造相对密度场，建立第i种材料第k次分级迭代前后高斯点材料属性关系；(8)建立正交各向异性多材料结构静力问题的无网格efgm离散控制方程，并求解正交各向异性多材料结构的节点参数位移u、节点近似位移柔度(9)建立正交各向异性多材料结构稳态传热问题的无网格efgm离散控制方程，并求解各向异性多材料结构的节点参数温度t、节点近似温度散热弱度(10)以多材料体积分数为约束条件，以柔度和散热弱度的加权目标函数最小为优化目标，建立基于无网格efgm和多材料simp插值方案的正交各向异性多材料热力耦合结构拓扑优化数学模型；(11)求解无网格efgm正交各向异性多材料热力耦合结构拓扑优化模型中第i种材料的体积灵敏度，采用伴随法求解柔度灵敏度、散热弱度灵敏度、加权目标函数灵敏度；(12)将体积灵敏度、加权目标函数灵敏度代入优化准则(optimality criteria，oc)，根据oc更新第k次分级迭代下第a种材料节点相对密度根据更新前后第a、b种材料节点相对密度之和保持不变的关系求解本次迭代下第b种材料节点相对密度组装多材料节点相对密度矩阵(13)比较t与t
max
大小关系，判断三级迭代为继续或结束；比较s与s
max
大小关系，判断二级迭代为继续或结束；比较r与r
max
及ch与ch
min
大小关系，判断一级迭代为继续或结束；
(14)输出无网格efgm正交各向异性多材料热力耦合结构拓扑优化的节点灰度图作为最优拓扑结构，其中多材料按照力热性能从高到低的排序分别表示为由深到浅的灰度颜色。2.根据权利要求1所述基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法，其特征在于，步骤(5)包括以下具体步骤：根据各向异性的力学参数求解各向异性多材料结构第i种材料在全局坐标系x-y下的弹性矩阵式中，为各向异性多材料结构第i种材料在材料坐标系1-2下的导热系数矩阵，为将与q
i
相关联的坐标转换矩阵，且其中，分别为材料坐标系中1和2方向上的第i种材料的弹性模量，分别为材料坐标系中1和2方向上的泊松比；定义泊松比因子来表示力学各向异性强度，则当prf＝1时，多材料呈现力学各向同性；通过修改泊松比因子prf和力学材料方向角θ
f
，可以得到具有不同力学承载性能的多材料拓扑优化结构。3.根据权利要求1所述基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法，其特征在于，步骤(6)包括以下具体步骤：根据各向异性多材料结构传热参数求解各向异性多材料结构第i种材料在全局坐标系x-y下的传热矩阵式中，为各向异性多材料结构第i种材料在材料坐标系1-2下的导热系数矩阵，为将与λ
i
相关联的坐标转换矩阵，分别为材料坐标系中1和2方向上第i种材料的导热系数；定义热导率因子来表示传热各向异性强度，则当tcf＝1时，多材料呈现传热各向同性，通过修改热导率因子tcf和传热材料方向角θ
p
，可以得到具有不同传热性能的多材料拓扑优化结构。4.根据权利要求1所述基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法，其特征在于，步骤(7)包括以下具体步骤：基于无网格efgm理论，并根据simp材料插值模型引入q种假想的相对密度在0～1之间可变的材料，选择设计域中efgm离散节点的多材料节点相对密度作为设计变量来构造相对密度场，建立第i种材料第k次分级迭代时的高斯点材料力热属性与初始材料力热属性之间的关系，其中k＝t
max
s
max
(r-1)+t
max
(s-1)+t为计算当前分级迭代步数的迭代索引；第i种材料第k次分级迭代的高斯点相对密度结合高斯点影响域内的节点相对密度和移动最小二乘法(moving least squares,mls)构造的形函数近
似求得，即式中，为第i种材料第k次分级迭代高斯点影响域内第i个节点相对密度，n为高斯点影响域ω
x
内的节点数目，φ
i
为高斯点影响域ω
x
内第i个节点的mls形函数。5.根据权利要求1所述基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法，其特征在于，步骤(8)包括以下具体步骤：(8.1)建立正交各向异性多材料结构静力问题的无网格efgm离散控制方程(8.2)求解引入多材料simp插值模型后的efgm整体力刚度矩阵(8.3)根据力载荷分布信息求解整体力载荷列向量f；(8.4)采用罚函数法逐一处理每条力载荷和位移边界条件，得到efgm整体力刚度矩阵惩罚项k
fα
和efgm整体力载荷向量惩罚项f
α
；(8.5)根据正交各向异性多材料结构无网格efgm静力离散控制方程求解efgm节点参数位移u列向量，结合整体节点形函数矩阵求解节点近似位移列向量求解柔度6.根据权利要求1所述基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法，其特征在于，步骤(9)包括以下具体步骤：(9.1)建立正交各向异性多材料结构稳态传热问题的无网格efgm离散控制方程(9.2)求解引入多材料simp插值模型后efgm整体热刚度矩阵(9.3)～(9.6)对给定的三类传热边界γ1、γ2、γ3条件和热载荷分别进行处理，求解刚度矩阵和热载荷列向量各分量，其中采用罚函数法处理第一类传热边界γ1，得到efgm整体温度刚度矩阵惩罚项k
pα
和efgm整体热载荷向量惩罚项p
α
，求解efgm整体热载荷列向量p；(9.7)根据正交各向异性多材料结构无网格efgm稳态传热离散控制方程求解efgm节点参数温度列向量t，结合整体节点形函数矩阵求解节点近似温度列向量求解多材料结构散热弱度7.根据权利要求1所述基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法，其特征在于，步骤(10)以多材料节点相对密度为设计变量，以多材料体积分数为约束条件，以柔度和散热弱度的加权目标函数最小为优化目标，建立基于无网格efgm和多材料simp插值方案的正交各向异性多材料热力耦合结构拓扑优化数学模型为：
式中，为第i种材料的节点相对密度行向量，为第i种材料的第j个节点相对密度，为多材料节点相对密度矩阵，q为多材料种类数，n为节点总数；为柔度，为散热弱度，c
max
和c
min
分别为单目标静力结构柔度最大值和最小值，s
max
和s
min
分别为单目标稳态传热结构散热弱度最大值和最小值；ω(0≤ω≤1)为加权目标函数权系数，为对和进行归一化处理后的加权目标函数，规定ω＝0表示单目标散热弱度最小化问题，ω＝1表示单目标柔度最小化问题；v
i
为第i种材料的体积，为第i种材料的高斯点相对密度，v
i
为多材料体积分数，v0为多材料总体积；为防止矩阵奇异，令第i种材料的高斯点相对密度和节点相对密度最小值为ρ
min
。8.根据权利要求1所述基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法，其特征在于，步骤(11)包括以下具体步骤：无网格efgm正交各向异性多材料热力耦合结构拓扑优化模型中第i种材料的体积灵敏度为式中，φ
i
为第i种材料高斯点影响域内第i个节点的mls形函数；无网格efgm正交各向异性多材料热力耦合结构拓扑优化模型中加权目标函数灵敏度为采用伴随法分别求解柔度灵敏度、散热弱度灵敏度为采用伴随法分别求解柔度灵敏度、散热弱度灵敏度为
式中，无网格efgm整体力刚度矩阵与其惩罚项叠加后的灵敏度为式中，无网格efgm整体热刚度矩阵与其惩罚项叠加后的灵敏度为9.根据权利要求1所述基于无网格efgm的各向异性多材料热力耦合结构拓扑优化方法，其特征在于，步骤(12)将体积灵敏度加权目标函数灵敏度代入优化准则oc，根据oc更新第k次分级迭代下第a种材料的节点相对密度为根据更新前后第a,b种材料节点相对密度之和保持不变的关系求解本次迭代下第b种材料的节点相对密度为组装多材料节点相对密度矩阵更新的具体步骤如下：若i＝a，则根据oc公式更新若i＝a，则根据oc公式更新式中，m0＝0.2为正移动限制，τ＝0.5为阻尼系数；为避免矩阵奇异，取节点相对密度下限为ρ
min
＝10-3
；根据oc准则，定义为若i＝b，则根据更新前后第a、b种材料节点相对密度之和保持不变的关系，求得第k次分级迭代下第b种材料的节点相对密度为若i≠a且i≠b，则本次分级迭代不涉及第i种材料节点相对密度更新，即第k、k-1次迭代第i种材料节点相对密度应相同，即输出多材料节点相对密度计算多材料节点相对密度改变值的无穷范数

技术总结
本发明公开了基于无网格EFGM的各向异性多材料热力耦合结构拓扑优化方法，主要步骤为：输入各向异性多材料种类数、体积分数、次弹性模量、剪切模量、次泊松比、泊松比因子、次导热系数、热导率因子、力热材料方向角、迭代步数、容差、惩罚力度；以EFGM节点离散设计域，给定力热载荷和边界、权系数等参数；求解各向异性多材料弹性和传热矩阵；建立基于EFGM和多材料SIMP插值方案的各向异性多材料热力耦合结构拓扑优化数学模型；求解位移场、温度场、体积和加权目标函数灵敏度；判断多级迭代是否终止，输出最优解。本发明实施的各向异性多材料热力耦合结构拓扑优化，实用性强，消除了棋盘格等数值不稳定现象，多材料边界清晰。多材料边界清晰。多材料边界清晰。


技术研发人员：张建平 张海明 张大兵 赵磊 龚曙光 吴淑英 刘洋 吴世雄
受保护的技术使用者：湘潭大学
技术研发日：2023.05.05
技术公布日：2023/7/27