渐进均匀化预测周期性复合材料热传导系数的简易实现方法

摘要

本发明公开了一种渐进均匀化预测周期性复合材料热传导系数的简易实现方法，通过商业有限元软件建立单胞有限元模型，设置单胞有限元模型的材料导热系数；对单胞有限元模型进行有限元网格划分，得到单胞有限元模型的节点、单元信息；由单胞有限元模型的节点坐标值确定初始温度场；计算周期性复合材料热传导系数。本发明利用商业有限元软件提供的单元类型以及任意组合，实现任何复杂周期性结构的性能预测，原均匀化中求解温度梯度场控制方程的右端载荷项时需要在每个单元上进行积分，与常规的有限元格式不同，通常需要自己编程实现，使用门槛较高。本发明使得均匀化方法能够通过现有的商业有限元软件进行直接计算，大幅度降低了使用难度。

说明书

技术领域

本发明涉及一种新的渐进均匀化预测周期性复合材料热传导系数的有限元 列式和基于商业有限元软件的实现方法。属于复合材料性能表征领域。

背景技术

复合材料因其优异的材料性能而在工程领域广泛使用，复合材料结构的设 计和使用都需要我们对复合材料性能有充分的了解。热传导系数就是其中一个 重要方面。针对复合材料的热传导系数预测，通常有解析模型、代表体元法及 均匀化等方法。解析模型为获得一个相对简单的解析表达式，通常需要做一些 不切实际的假设，这些假设会导致与实验结果产生较大误差；为与实验结果吻 合较好，则需要考虑更多的因素，所建立的解析模型则变得复杂不实用。解析 模型针对简单规则的结构较为实用。代表体元法能够考虑更为复杂的细观结构， 且能够考虑更多的细节，是一个通用的计算方法。代表体元方法是目前最为流 行的热传导预测方法。该方法基于热弹性能量等效原理，概念清晰，执行简单， 是一种近似模型。

均匀化方法以摄动理论为依据，有着严格的数学基础，对于无限周期微结 构的材料能给出精确解。均匀化方法已经被用于弹性模量、热膨胀系数、热传 导系数等周期性复合材料的性能预测。在渐近均匀化方法的传统有限元实现中， 需要在每个单元上积分以求得等效荷载和应变能，这就需要与单元相关矩阵的 所有细节，例如本构矩阵、应变－位移矩阵。对于不同的单元类型，这些矩阵 也是不同的，所以需要针对不同的单元写出相应的有限元列式和编写相应的代 码，因此针对实体结构和板壳等不同的周期性复合材料结构，需要发展不同的 均匀化求解列式。因此，针对不同的问题，需要编写相应的有限元代码。现有 通用商业软件均没有均匀化方法相关模块功能，这严重限制了均匀化方法的应 用。因此，为推广均匀化方法的应用，需要与现有成熟的商业有限元软件相结 合，一种基于商业有限元软件的渐进均匀化预测周期性复合材料热传导系数的 简易实现方法亟待研发。

发明内容

根据上述提出的渐进均匀化理论计算周期性结构热传导系数实施过程复 杂、效率低下、通用性差等一系列问题，提出了一种基于商业有限元软件的渐 进均匀化预测周期性复合材料热传导系数的简易实现方法。

本发明采用的技术手段如下：

一种渐进均匀化预测周期性复合材料热传导系数的简易实现方法，具有以 下步骤：

S1、通过商业有限元软件建立单胞有限元模型，设置单胞有限元模型的材 料导热系数；

S2、对单胞有限元模型进行有限元网格划分，得到单胞有限元模型的节点、 单元信息；

S3、由单胞有限元模型的节点坐标值确定初始温度场；

S4、计算周期性复合材料热传导系数。

当所述周期性复合材料呈二维结构时，所述单胞有限元模型为二维模型， 所述步骤S3中由单胞有限元模型的节点坐标值确定初始温度场Λ^0(m)为：

${\Lambda }^{0\left(1\right)}={\{{\Lambda }_1^{0\left(1\right)},{\Lambda }_2^{0\left(1\right)},...,{\Lambda }_i^{0\left(1\right)},...,{\Lambda }_N^{0\left(1\right)}\}}^T,$

${\Lambda }^{0\left(2\right)}={\{{\Lambda }_1^{0\left(2\right)},{\Lambda }_2^{0\left(2\right)},...,{\Lambda }_i^{0\left(2\right)},...,{\Lambda }_N^{0\left(2\right)}\}}^T,$

其中，m＝{1，2}，x_i，y_i为第i个节点的坐标值，N表示单 胞有限元模型的节点总个数。

所述步骤S4中计算周期性复合材料热传导系数具体如下步骤：

A1、将Λ^0(m)施加到单胞有限元模型的节点上，并对每个工况单独进行有限 元静力分析，由f^0(m)＝KΛ^0(m)得到Λ^0(m)对应的节点热流反力场f^0(m)为：

$f^{0\left(1\right)}={\{f_1^{0\left(1\right)},f_2^{0\left(1\right)},...,f_i^{0\left(1\right)},...,f_N^{0\left(1\right)}\}}^T,$

$f^{0\left(2\right)}={\{f_1^{0\left(2\right)},f_2^{0\left(2\right)},...,f_i^{0\left(2\right)},...,f_N^{0\left(2\right)}\}}^T,$

其中，为在步骤A1条件下第i个节点在Λ⁰⁽¹⁾工况下的节点热流反力，为在步骤A1条件下第i个节点在Λ⁰⁽²⁾工况下的节点热流反力，K为单胞有限元 模型的总体导热矩阵，所述步骤A1条件为将Λ^0(m)施加到单胞有限元模型的节点 上，并对每个工况单独进行有限元静力分析；

A2、删除所有约束后，将f^0(m)施加到单胞有限元模型的节点上，并施加周 期性边界条件，之后对每个工况单独进行有限元静力分析，由求得 f^0(m)对应的特征位移场Λ^*(m)：

${\Lambda }^{*\left(1\right)}={\{{\Lambda }_1^{*\left(1\right)},{\Lambda }_2^{*\left(1\right)},...,{\Lambda }_i^{*\left(1\right)},...,{\Lambda }_N^{*\left(1\right)}\}}^T,$

${\Lambda }^{*\left(2\right)}={\{{\Lambda }_1^{*\left(2\right)},{\Lambda }_2^{*\left(2\right)},...,{\Lambda }_i^{*\left(2\right)},...,{\Lambda }_N^{*\left(2\right)}\}}^T,$

其中为在步骤A2条件下第i个节点在f⁰⁽¹⁾工况下的特征温度，为在 步骤A2条件下第i个节点在f⁰⁽²⁾工况下的特征温度，为施加周期性边界条件 后单胞有限元模型的总体导热系数矩阵，所述步骤A2条件为删除所有约束后， 将f^0(m)施加到单胞有限元模型的节点上，并施加周期性边界条件，之后对每个工 况单独进行有限元静力分析；

A3、删除载荷、约束及周期性边界条件，将Λ^*(m)施加到单胞有限元模型的 节点上，并对每个工况单独进行有限元静力分析，由f^*(m)＝KΛ^*(m)求得Λ^*(m)对应的 节点热流反力场为f^*(m)为：

$f^{*\left(1\right)}={\{f_1^{*\left(1\right)},f_2^{*\left(1\right)},...,f_i^{*\left(1\right)},...f_N^{*\left(1\right)}\}}^T,$

$f^{*\left(2\right)}={\{f_1^{*\left(2\right)},f_2^{*\left(2\right)},...,f_i^{*\left(2\right)},...f_N^{*\left(2\right)}\}}^T,$

其中为在步骤A3条件下第i个节点在Λ^*(1)工况下的节点热流反力，为 在步骤A3条件下第i个节点在Λ^*(2)工况下的节点热流反力，K为单胞有限元模 型的总体导热系数矩阵，所述步骤A3条件为删除载荷、约束及周期性边界条件， 将Λ^*(m)施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分 析；

A4、根据求得周期性复合材料热传导系数， 其中，为周期性复合材料热传导系数，Y为所述单胞有限元模型的面积， n＝{1，2}，m和n表示施加的温度场工况，m＝1和n＝1代表x方向，m＝2和n＝2 代表y方向。

当所述周期性复合材料呈三维结构时，所述单胞有限元模型为三维模型， 所述步骤S3中由单胞有限元模型的节点坐标值确定初始温度场Λ^0(m)为：

${\Lambda }^{0\left(1\right)}={\{{\Lambda }_1^{0\left(1\right)},{\Lambda }_2^{0\left(1\right)},...,{\Lambda }_i^{0\left(1\right)},...,{\Lambda }_N^{0\left(1\right)}\}}^T,$

${\Lambda }^{0\left(2\right)}={\{{\Lambda }_1^{0\left(2\right)},{\Lambda }_2^{0\left(2\right)},...,{\Lambda }_i^{0\left(2\right)},...,{\Lambda }_N^{0\left(2\right)}\}}^T,$

${\Lambda }^{0\left(3\right)}={\{{\Lambda }_1^{0\left(3\right)},{\Lambda }_2^{0\left(3\right)},...,{\Lambda }_i^{0\left(3\right)},...,{\Lambda }_N^{0\left(3\right)}\}}^T,$

其中，

m＝{1,2,3}，x_i,y_i,z_i为第i个节点的坐标值，N表示单 胞有限元模型的节点总个数。

所述步骤S4中计算周期性复合材料热传导系数具有如下步骤：

B1、将Λ^0(m)施加到单胞有限元模型的节点上，并对每个工况单独进行有限 元静力分析，由f^0(m)＝KΛ^0(m)得到Λ^0(m)对应的节点热流反力场f^0(m)为：

$f^{0\left(1\right)}={\{f_1^{0\left(1\right)},f_2^{0\left(1\right)},...,f_i^{0\left(1\right)},...f_N^{0\left(1\right)}\}}^T,$

$f^{0\left(2\right)}={\{f_1^{0\left(2\right)},f_2^{0\left(2\right)},...,f_i^{0\left(2\right)},...f_N^{0\left(2\right)}\}}^T,$

$f^{0\left(3\right)}={\{f_1^{0\left(3\right)},f_2^{0\left(3\right)},...,f_i^{0\left(3\right)},...f_N^{0\left(3\right)}\}}^T,$

其中为在步骤B1条件下第i个节点在Λ⁰⁽¹⁾工况下的节点热流反力，为 在步骤B1条件下第i个节点在Λ⁰⁽²⁾工况下的节点热流反力，为在步骤B1条 件下第i个节点在Λ⁰⁽³⁾工况下的节点热流反力，K为单胞有限元模型的总体导热 刚度阵，所述步骤B1条件为将Λ^0(m)施加到单胞有限元模型的节点上，并对每个 工况单独进行有限元静力分析；

B2、删除所有约束后，将f^0(m)施加到单胞有限元模型的节点上，并施加周期 性边界条件，之后对每个工况单独进行有限元静力分析，由求得f^0(m)对应的特征位移场Λ^*(m)：

${\Lambda }^{*\left(1\right)}={\{{\Lambda }_1^{*\left(1\right)},{\Lambda }_2^{*\left(1\right)},...,{\Lambda }_i^{*\left(1\right)},...{\Lambda }_N^{*\left(1\right)}\}}^T,$

${\Lambda }^{*\left(2\right)}={\{{\Lambda }_1^{*\left(2\right)},{\Lambda }_2^{*\left(2\right)},...,{\Lambda }_i^{*\left(2\right)},...{\Lambda }_N^{*\left(2\right)}\}}^T,$

${\Lambda }^{*\left(3\right)}={\{{\Lambda }_1^{*\left(3\right)},{\Lambda }_2^{*\left(3\right)},...,{\Lambda }_i^{*\left(3\right)},...{\Lambda }_N^{*\left(3\right)}\}}^T,$

其中为在步骤C2条件下第i个节点在f⁰⁽¹⁾工况下的特征温度，为在 步骤B2条件下第i个节点在f⁰⁽²⁾工况下的特征温度，为在步骤B2条件下第 i个节点在f⁰⁽³⁾工况下的特征温度，为施加周期性边界条件后单胞有限元模型 的总体导热刚度阵，所述步骤B2条件为删除所有约束后，将f^0(m)施加到单胞有 限元模型的节点上，并施加周期性边界条件，之后对每个工况单独进行有限元 静力分析；

B3、删除载荷、约束及周期性边界条件，将Λ^*(m)施加到单胞有限元模型的 节点上，并对每个工况单独进行有限元静力分析，由f^*(m)＝KΛ^*(m)求得Λ^*(m)对应的 节点热流反力场为f^*(m)为：

$f^{*\left(1\right)}={\{f_1^{*\left(1\right)},f_2^{*\left(1\right)},...,f_i^{*\left(1\right)},...f_N^{*\left(1\right)}\}}^T,$

$f^{*\left(2\right)}={\{f_1^{*\left(2\right)},f_2^{*\left(2\right)},...,f_i^{*\left(2\right)},...f_N^{*\left(2\right)}\}}^T,$

$f^{*\left(3\right)}={\{f_1^{*\left(3\right)},f_2^{*\left(3\right)},...,f_i^{*\left(3\right)},...f_N^{*\left(3\right)}\}}^T,$

其中为在步骤B3条件下第i个节点在Λ^*(1)工况下x，y，z方向的节点热 流反力，为在步骤B3条件下第i个节点在Λ^*(2)工况下的节点热流反力，为 在步骤B3条件下第i个节点在Λ^*(3)工况下的节点热流反力，K为单胞有限元模 型的总体导热刚度阵，所述步骤B3条件为删除载荷、约束及周期性边界条件， 将Λ^*(m)施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分 析；

B4、根据求得周期性复合材料热传导系数， 其中，为周期性复合材料热传导系数，Y为所述单胞有限元模型的体积， n＝{1,2,3}，m和n表示施加的温度场工况，m＝1和n＝1代表x方向，m＝2和n＝2 代表y方向，m＝3和n＝3代表z方向。

原均匀化中求解温度梯度场控制方程的右端载荷项时需要在每个单元上进 行积分，与常规的有限元格式不同，通常需要自己编程实现，使得该方法的使 用门槛较高。与现有技术相比，本发明能够在保持传统均匀化方法计算精度的 基础上，执行非常方便，能够以现有商业有限元软件为黑箱，实现任何复杂周 期性结构的性能预测，计算的过程与代表体元方法的难易程度相似，大幅度降 低了该方法的使用难度。另外，该发明在计算过程中，能够利用商业有限元软 件提供的单元类型或任意单元组合(如杆、板、实体等)，计算任意复杂结构的 热传导系数，拓展了该方法的适用性。

基于上述理由本发明可在复合材料性能表征等领域广泛推广。

附图说明

下面结合附图和具体实施方式对本发明作进一步详细的说明。

图1是本发明具体实施方式中计算周期性复合材料热传导系数的流程图。

图2是本发明实施例1中含方形纤维增强复合材料单胞的结构示意图。

图3是本发明实施例2中空心球结构单胞的结构示意图。

具体实施方式

实施例1

如图1和图2所示，含方形纤维增强复合材料单胞的内部为增强纤维项1， 外部为树脂基体2，计算含方形纤维增强复合材料单胞的热传导系数，具有以下 步骤：

S1、通过商业有限元软件建立单胞有限元模型，设置单胞有限元模型的材 料导热系数；

S2、对单胞有限元模型进行有限元网格划分，得到单胞有限元模型的节点、 单元信息；

S3、由单胞有限元模型的节点坐标值确定初始温度场Λ^0(m)为：

${\Lambda }^{0\left(1\right)}={\{{\Lambda }_1^{0\left(1\right)},{\Lambda }_2^{0\left(1\right)},...,{\Lambda }_i^{0\left(1\right)},...,{\Lambda }_N^{0\left(1\right)}\}}^T,$

${\Lambda }^{0\left(2\right)}={\{{\Lambda }_1^{0\left(2\right)},{\Lambda }_2^{0\left(2\right)},...,{\Lambda }_i^{0\left(2\right)},...,{\Lambda }_N^{0\left(2\right)}\}}^T,$

其中，m＝{1，2}，x_i，y_i为第i个节点的坐标值，N表示单 胞有限元模型的节点总个数；

S4、计算周期性复合材料热传导系数具体如下步骤：

A1、将Λ^0(m)施加到单胞有限元模型的节点上，并对每个工况单独进行有限 元静力分析，由f^0(m)＝KΛ^0(m)得到Λ^0(m)对应的节点热流反力场f^0(m)为：

$f^{0\left(1\right)}={\{f_1^{0\left(1\right)},f_2^{0\left(1\right)},...,f_i^{0\left(1\right)},...,f_N^{0\left(1\right)}\}}^T,$

$f^{0\left(2\right)}={\{f_1^{0\left(2\right)},f_2^{0\left(2\right)},...,f_i^{0\left(2\right)},...,f_N^{0\left(2\right)}\}}^T,$

其中，为在步骤A1条件下第i个节点在Λ⁰⁽¹⁾工况下的节点热流反力，为在步骤A1条件下第i个节点在Λ⁰⁽²⁾工况下的节点热流反力，K为单胞有限元 模型的总体导热矩阵，所述步骤A1条件为将Λ^0(m)施加到单胞有限元模型的节点 上，并对每个工况单独进行有限元静力分析；

A2、删除所有约束后，将f^0(m)施加到单胞有限元模型的节点上，并施加周 期性边界条件，之后对每个工况单独进行有限元静力分析，由求得 f^0(m)对应的特征位移场Λ^*(m)：

${\Lambda }^{*\left(1\right)}={\{{\Lambda }_1^{*\left(1\right)},{\Lambda }_2^{*\left(1\right)},...,{\Lambda }_i^{*\left(1\right)},...,{\Lambda }_N^{*\left(1\right)}\}}^T,$

${\Lambda }^{*\left(2\right)}={\{{\Lambda }_1^{*\left(2\right)},{\Lambda }_2^{*\left(2\right)},...,{\Lambda }_i^{*\left(2\right)},...,{\Lambda }_N^{*\left(2\right)}\}}^T,$

其中为在步骤A2条件下第i个节点在f⁰⁽¹⁾工况下的特征温度，为在 步骤A2条件下第i个节点在f⁰⁽²⁾工况下的特征温度，为施加周期性边界条件 后单胞有限元模型的总体导热系数矩阵，所述步骤A2条件为删除所有约束后， 将f^0(m)施加到单胞有限元模型的节点上，并施加周期性边界条件，之后对每个工 况单独进行有限元静力分析；

A3、删除载荷、约束及周期性边界条件，将Λ^*(m)施加到单胞有限元模型的 节点上，并对每个工况单独进行有限元静力分析，由f^*(m)＝KΛ^*(m)求得Λ^*(m)对应的 节点热流反力场为f^*(m)为：

$f^{*\left(1\right)}={\{f_1^{*\left(1\right)},f_2^{*\left(1\right)},...,f_i^{*\left(1\right)},...f_N^{*\left(1\right)}\}}^T,$

$f^{*\left(2\right)}={\{f_1^{*\left(2\right)},f_2^{*\left(2\right)},...,f_i^{*\left(2\right)},...f_N^{*\left(2\right)}\}}^T,$

其中为在步骤A3条件下第i个节点在Λ^*(1)工况下的节点热流反力，为 在步骤A3条件下第i个节点在Λ^*(2)工况下的节点热流反力，K为单胞有限元模 型的总体导热系数矩阵，所述步骤A3条件为删除载荷、约束及周期性边界条件， 将Λ^*(m)施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分 析；

A4、根据求得周期性复合材料热传导系数， 其中，为周期性复合材料热传导系数，Y为所述单胞有限元模型的面积， n＝{1，2}，m和n表示施加的温度场工况m＝1和n＝1代表x方向，m＝2和n＝2代 表y方向。

实施例2

如图1和图3所示，空心球结构单胞的中间为一定厚度的金属壳3(Al或 者新型材料St)，外部为粘结剂材料4(树脂)，计算空心球结构单胞的热传导系 数，具有以下步骤：

S1、通过商业有限元软件建立单胞有限元模型，设置单胞有限元模型的材 料导热系数；

S2、对单胞有限元模型进行有限元网格划分，得到单胞有限元模型的节点、 单元信息；

S3、由单胞有限元模型的节点坐标值确定初始温度场Λ^0(m)为：

${\Lambda }^{0\left(1\right)}={\{{\Lambda }_1^{0\left(1\right)},{\Lambda }_2^{0\left(1\right)},...,{\Lambda }_i^{0\left(1\right)},...,{\Lambda }_N^{0\left(1\right)}\}}^T,$

${\Lambda }^{0\left(2\right)}={\{{\Lambda }_1^{0\left(2\right)},{\Lambda }_2^{0\left(2\right)},...,{\Lambda }_i^{0\left(2\right)},...,{\Lambda }_N^{0\left(2\right)}\}}^T,$

${\Lambda }^{0\left(3\right)}={\{{\Lambda }_1^{0\left(3\right)},{\Lambda }_2^{0\left(3\right)},...,{\Lambda }_i^{0\left(3\right)},...,{\Lambda }_N^{0\left(3\right)}\}}^T,$

其中，

m＝{1,2,3}，x_i,y_i,z_i为第i个节点的坐标值，N表 示单胞有限元模型的节点总个数。；

S4、计算周期性复合材料热传导系数具有如下步骤：

B1、将Λ^0(m)施加到单胞有限元模型的节点上，并对每个工况单独进行有限 元静力分析，由f^0(m)＝KΛ^0(m)得到Λ^0(m)对应的节点热流反力场f^0(m)为：

$f^{0\left(1\right)}={\{f_1^{0\left(1\right)},f_2^{0\left(1\right)},...,f_i^{0\left(1\right)},...f_N^{0\left(1\right)}\}}^T,$

$f^{0\left(2\right)}={\{f_1^{0\left(2\right)},f_2^{0\left(2\right)},...,f_i^{0\left(2\right)},...f_N^{0\left(2\right)}\}}^T,$

$f^{0\left(3\right)}={\{f_1^{0\left(3\right)},f_2^{0\left(3\right)},...,f_i^{0\left(3\right)},...f_N^{0\left(3\right)}\}}^T,$

其中为在步骤B1条件下第i个节点在Λ⁰⁽¹⁾工况下的节点热流反力，为 在步骤B1条件下第i个节点在Λ⁰⁽²⁾工况下的节点热流反力，为在步骤B1条 件下第i个节点在Λ⁰⁽³⁾工况下的节点热流反力，K为单胞有限元模型的总体导热 刚度阵，所述步骤B1条件为将Λ^0(m)施加到单胞有限元模型的节点上，并对每个 工况单独进行有限元静力分析；

B2、删除所有约束后，将f^0(m)施加到单胞有限元模型的节点上，并施加周期 性边界条件，之后对每个工况单独进行有限元静力分析，由求得f^0(m)对应的特征位移场Λ^*(m)：

${\Lambda }^{*\left(1\right)}={\{{\Lambda }_1^{*\left(1\right)},{\Lambda }_2^{*\left(1\right)},...,{\Lambda }_i^{*\left(1\right)},...{\Lambda }_N^{*\left(1\right)}\}}^T,$

${\Lambda }^{*\left(2\right)}={\{{\Lambda }_1^{*\left(2\right)},{\Lambda }_2^{*\left(2\right)},...,{\Lambda }_i^{*\left(2\right)},...{\Lambda }_N^{*\left(2\right)}\}}^T,$

${\Lambda }^{*\left(3\right)}={\{{\Lambda }_1^{*\left(3\right)},{\Lambda }_2^{*\left(3\right)},...,{\Lambda }_i^{*\left(3\right)},...{\Lambda }_N^{*\left(3\right)}\}}^T,$

其中为在步骤C2条件下第i个节点在f⁰⁽¹⁾工况下的特征温度，为在 步骤B2条件下第i个节点在f⁰⁽²⁾工况下的特征温度，为在步骤B2条件下第 i个节点在f⁰⁽³⁾工况下的特征温度，为施加周期性边界条件后单胞有限元模型 的总体导热刚度阵，所述步骤B2条件为删除所有约束后，将f^0(m)施加到单胞有 限元模型的节点上，并施加周期性边界条件，之后对每个工况单独进行有限元 静力分析；

B3、删除载荷、约束及周期性边界条件，将Λ^*(m)施加到单胞有限元模型的 节点上，并对每个工况单独进行有限元静力分析，由f^*(m)＝KΛ^*(m)求得Λ^*(m)对应的 节点热流反力场为f^*(m)为：

$f^{*\left(1\right)}={\{f_1^{*\left(1\right)},f_2^{*\left(1\right)},...,f_i^{*\left(1\right)},...f_N^{*\left(1\right)}\}}^T,$

$f^{*\left(2\right)}={\{f_1^{*\left(2\right)},f_2^{*\left(2\right)},...,f_i^{*\left(2\right)},...f_N^{*\left(2\right)}\}}^T,$

$f^{*\left(3\right)}={\{f_1^{*\left(3\right)},f_2^{*\left(3\right)},...,f_i^{*\left(3\right)},...f_N^{*\left(3\right)}\}}^T,$

其中为在步骤B3条件下第i个节点在Λ^*(1)工况下x，y，z方向的节点热 流反力，为在步骤B3条件下第i个节点在Λ^*(2)工况下的节点热流反力，为 在步骤B3条件下第i个节点在Λ^*(3)工况下的节点热流反力，K为单胞有限元模 型的总体导热刚度阵，所述步骤B3条件为删除载荷、约束及周期性边界条件， 将Λ^*(m)施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分 析；

B4、根据求得周期性复合材料热传导系数， 其中，为周期性复合材料热传导系数，Y为所述单胞有限元模型的体积， n＝{1,2,3}，m和n表示施加的温度场工况，m＝1和n＝1代表x方向，m＝2和n＝2 代表y方向，m＝3和n＝3代表z方向。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局 限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，根据本 发明的技术方案及其发明构思加以等同替换或改变，都应涵盖在本发明的保护 范围之内。