CN105160085B

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

Info

Publication number: CN105160085B
Application number: CN201510527831.9A
Authority: CN
Inventors: 张永存; 尚士朋; 陈文炯; 刘书田
Original assignee: Dalian University of Technology
Current assignee: Dalian University of Technology
Priority date: 2015-08-24
Filing date: 2015-08-24
Publication date: 2018-05-15
Anticipated expiration: 2035-08-24
Also published as: CN105160085A

Abstract

本发明公开了一种渐进均匀化预测周期性复合材料热膨胀系数的简易实现方法，具有以下步骤：通过商业有限元软件建立单胞有限元模型，设置单胞有限元模型中材料参数；对单胞有限元模型进行有限元网格划分，得到单胞有限元模型的节点、单元信息；由单胞有限元模型的节点坐标值确定初始位移场；分别计算周期性结构等效弹性模量和周期性结构等效热弹性常数；计算周期性复合材料热膨胀系数。本发明具有执行方便，能够以现有商业有限元软件为黑箱，预测周期性复合材料的热膨胀系数，大幅度拓展了均匀化方法的适应性，将原均匀化中求解应变场并在每个单元上积分的过程转换为商业有限元软件中可以得到的对应位移场和节点反力场的乘积，过程简化。

Description

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

技术领域

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

背景技术

热膨胀是指由于环境温度变化而导致材料/结构体积变化的现象。这种由于温度引起的形状改变经常被用于传感器、可变体结构的主动控制设计。另一方面，材料这种热胀冷缩的固有现象也对结构的功能设计产生不同程度的危害，为设计带来巨大的挑战。而实现结构膨胀设计与精确变形控制的基础是构建热膨胀系数与材料组分及其分布的函数关系，实现复合材料/结构化材料热膨胀系数的高精度预测。

复合材料热膨胀系数的预测方法一直备受关注。针对单向纤维增强复合材料的热膨胀系数预测，提出了多种解析模型。解析模型为获得一个相对简单的解析表达式，通常需要做一些不切实际的假设，这些假设会导致与实验结果产生较大误差；为与实验结果吻合较好，则需要考虑更多的因素，所建立的解析模型则变得复杂不实用。解析模型针对简单规则的结构较为实用。代表体元法能够考虑更为复杂的细观结构，且能够考虑更多的细节，是一个通用的计算方法。代表体元方法是目前最为流行的热膨胀预测方法。该方法基于热弹性能量等效原理，概念清晰，执行简单，是一种近似模型。

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

发明内容

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

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

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

S1、通过商业有限元软件建立单胞有限元模型，设置单胞有限元模型中材料参数，材料参数包括单胞有限元模型中各个材料的弹性模量，膨胀系数和泊松比；

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

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

S4、分别计算周期性结构等效弹性模量和周期性结构等效热弹性常数；

S5、计算周期性复合材料热膨胀系数。

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

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

所述步骤S4中计算周期性结构等效弹性模量具有如下步骤：

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

其中，和分别为在步骤A1条件下第i个节点在工况下x，y方向的节点反力，和分别为在步骤A1条件下第i个节点在工况下x，y方向的节点反力，和分别为在步骤A1条件下第i个节点在工况下x，y方向的节点反力，K为单胞有限元模型的总体刚度阵，所述步骤A1条件为将施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分析；

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

其中，和分别为在步骤A2条件下第i个节点在f⁽¹¹⁾工况下x，y方向的特征位移，和分别为在步骤A2条件下第i个节点在f⁽²²⁾工况下x，y方向的特征位移，和分别为在步骤A2条件下第i个节点在f⁽¹²⁾工况下x，y方向的特征位移，为施加周期性边界条件后单胞有限元模型的总体刚度阵，所述步骤A2条件为删除所有约束后，将f^(ij)施加到单胞有限元模型的节点上，并施加周期性边界条件，之后对每个工况单独进行有限元静力分析；

A3、删除所有载荷及约束后，将χ^*(ij)施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分析，由f^*(ij)＝Kχ^*(ij)求得χ^*(ij)对应的节点反力场f^*(ij)为：

其中，和分别为在步骤A3条件下第i个节点在χ^*(11)工况下x，y方向的节点反力，和分别为在步骤A3条件下第i个节点在χ^*(22)工况下x，y方向的节点反力，和分别为在步骤A3条件下第i个节点在χ^*(12)工况下x，y方向的节点反力，K为单胞有限元模型的总体刚度阵，所述步骤A3条件为删除所有载荷及约束后，将χ^*(ij)施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分析；

A4、根据求得周期性结构等效弹性模量，其中，为周期性结构等效弹性模量，Y为所述单胞有限元模型的面积。

所述步骤S4中计算周期性结构等效热弹性常数具有如下步骤：

B1、约束单胞有限元模型的各节点的位移，施加单位负温升，进行有限元热分析，求得节点支撑反力场，并将所述节点支撑反力场作为等效的节点力场，所述等效的节点力场为：

其中，和分别为在步骤B1条件下第i个节点在x，y方向的节点反力，B为单胞有限元模型中单元的应变-位移矩阵，E为单胞有限元模型中各个材料的单元刚度阵，α是单胞有限元模型中各个材料的热膨胀系数，Y为所述单胞有限元模型的面积，所述步骤B1条件为约束单胞有限元模型的各节点的位移，施加单位负温升，进行有限元热分析，所述单位负温升指的是温度降低1℃；

B2、将R_α施加到单胞有限元模型的节点上，并施加周期性边界条件，进行有限元静力分析，由求得R_α对应的节点特征位移场Λ：

Λ＝{Λ₁,Λ₂,…,Λ_i,…Λ_N}^T，

其中，和分别为在步骤B2条件下第i个节点在x，y方向的位移，为施加周期性边界条件后单胞有限元模型的总体刚度阵，所述步骤B2条件为将R_α施加到单胞有限元模型的节点上，并施加周期性边界条件，进行有限元静力分析；

B3、删除所有载荷及约束后，将Λ施加到单胞有限元模型的节点上，并进行有限元静力分析，由R_Λ＝KΛ求得Λ对应的节点反力场R_Λ为：

其中，和分别为在步骤B3条件下第i个节点在x，y方向的节点反力，K为单胞有限元模型的总体刚度阵，所述步骤B3条件为删除所有载荷及约束后，将Λ施加到单胞有限元模型的节点上，并进行有限元静力分析；

B4、根据求周期性结构等效热弹性常数，其中，为周期性结构等效热弹性常数，Y为所述单胞有限元模型的面积。

周期性复合材料热膨胀系数通过以下公式求得：

其中，为周期性复合材料热膨胀系数。

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

其中，

ij＝{11,22,33,23,31,12}，x_i,y_i,z_i为第i个节点的坐标值，N表示单胞有限元模型的节点总个数。

所述步骤S4中计算周期性结构等效弹性模量具有如下步骤：

C1、将施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分析，由求得对应的节点反力场f^(ij)为：

其中，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，K为单胞有限元模型的总体刚度阵，所述步骤C1条件为将施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分析；

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

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

C3、删除所有载荷及约束后，将χ^*(ij)施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分析，由f^*(ij)＝Kχ^*(ij)求得χ^*(ij)对应的节点反力场f^*(ij)为：

其中，

和分别为在步骤C3条件下第i个节点在χ^*(11)工况下x，y，z方向的节点反力，和分别为在步骤C3条件下第i个节点在χ^*(22)工况下x，y，z方向的节点反力，和分别为在步骤C3条件下第i个节点在χ^*(33)工况下x，y，z方向的节点反力，和分别为在步骤C3条件下第i个节点在χ^*(23)工况下x，y，z方向的节点反力，和分别为在步骤C3条件下第i个节点在χ^*(31)工况下x，y，z方向的节点反力，和分别为在步骤C3条件下第i个节点在χ^*(12)工况下x，y，z方向的节点反力，K为单胞有限元模型的总体刚度阵，所述步骤C3条件为删除所有载荷及约束后，将χ^*(ij)施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分析；

C4、根据求得周期性结构等效弹性模量，其中，为周期性结构等效弹性模量，Y为所述单胞有限元模型的体积。

D1、约束单胞有限元模型的各节点的位移，施加单位负温升，进行有限元热分析，求得节点支撑反力场，并将所述节点支撑反力场作为等效的节点力场，所述等效的节点力场为：

其中，和分别为在步骤D1条件下第i个节点在x，y，z方向的节点反力，B为单胞有限元模型中单元的应变-位移矩阵，E为单胞有限元模型中各个材料的单元刚度阵，α是单胞有限元模型中各个材料的热膨胀系数，Y为所述单胞有限元模型的面积，所述步骤D1条件为约束单胞有限元模型的各节点的位移，施加单位负温升，进行有限元热分析，所述单位负温升指的是温度降低1℃；

D2、将R_α施加到单胞有限元模型的节点上，并施加周期性边界条件，进行有限元静力分析，由求得R_α对应的节点特征位移场Λ：

Λ＝{Λ₁,Λ₂,…,Λ_i,…Λ_N}^T，

其中，和分别为在步骤D2条件下第i个节点在x，y，z方向的位移，为施加周期性边界条件后单胞有限元模型的总体刚度阵，所述步骤D2条件为将R_α施加到单胞有限元模型的节点上，并施加周期性边界条件，进行有限元静力分析；

D3、删除所有载荷及约束后，将Λ施加到单胞有限元模型的节点上，并进行有限元静力分析，由R_Λ＝KΛ求得Λ对应的节点反力场R_Λ为：

其中，和分别为在步骤D3条件下第i个节点在x，y，z方向的节点反力，K为单胞有限元模型的总体刚度阵，所述步骤D3条件为删除所有载荷及约束后，将Λ施加到单胞有限元模型的节点上，并进行有限元静力分析；

D4、根据求周期性结构等效热弹性常数，其中，为周期性结构等效热弹性常数，Y为所述单胞有限元模型的面积。

周期性复合材料热膨胀系数通过以下公式求得：

其中，为周期性复合材料热膨胀系数。

与现有技术相比，本发明具有执行方便，能够以现有商业有限元软件为黑箱，无须任何编程，利用商业有限元软件提供的单元类型以及任意组合，实现任何复杂周期性结构的性能预测，而计算的过程与代表体元方法的难易程度相似，并且在保持传统均匀化方法计算精度的基础上，能够很方便预测周期性复合材料的热膨胀系数；大幅度拓展了均匀化方法的适应性，克服了传统均匀化实现方法中仅能采用一种有限单元类离散分析单胞的缺陷，能够处理不同单元类型混合的复杂单胞结构；借鉴求解周期性结构等效弹性模量的新均匀化方法思路，将原均匀化中求解应变场并在每个单元上积分的过程转换为商业有限元软件中可以得到的对应位移场和节点反力场的乘积，过程简化，并且将原均匀化方法的通用性拓展到混合单元(杆、板、实体等)建模的单胞计算中去，使得新均匀化方法能够计算任意复杂结构的弹性模量及热膨胀系数。

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

附图说明

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

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

图2是本发明实施例1中单胞有限元模型的结构示意图。

图3是本发明实施例1中单胞有限元模型的结构示意图。

具体实施方式

实施例1

如图1所示，计算纤维增强复合材料膨胀系数，具有以下步骤：

S1、通过商业有限元软件建立单胞有限元模型，设置单胞有限元模型中材料参数，如图2所示；

S3、由单胞有限元模型的节点坐标值确定初始位移场为：

其中，

，ij＝{11,22,33,23,31,12}，x_i,y_i,z_i为第i个节点的坐标值，N表示单胞有限元模型的节点总个数；

S4、分别计算周期性结构等效弹性模量和周期性结构等效热弹性常数，

计算周期性结构等效弹性模量具有如下步骤：

其中，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，和分别为在步骤C1条件下第i个节点在工况下x，y，z方向的节点反力，K为单胞有限元模型的总体刚度阵；

其中，和分别为在步骤C2条件下第i个节点在f⁽¹¹⁾工况下x，y，z方向的特征位移，和分别为在步骤C2条件下第i个节点在f⁽²²⁾工况下x，y，z方向的特征位移，和分别为在步骤C2条件下第i个节点在f⁽³³⁾工况下x，y，z方向的特征位移，和分别为在步骤C2条件下第i个节点在f⁽²³⁾工况下x，y，z方向的特征位移，和分别为在步骤C2条件下第i个节点在f⁽³¹⁾工况下x，y，z方向的特征位移，和分别为在步骤C2条件下第i个节点在f⁽¹²⁾工况下x，y，z方向的特征位移，为施加周期性边界条件后单胞有限元模型的总体刚度阵；

其中，

和分别为在步骤C3条件下第i个节点在χ^*(11)工况下x，y，z方向的节点反力，和分别为在步骤C3条件下第i个节点在χ^*(22)工况下x，y，z方向的节点反力，和分别为在步骤C3条件下第i个节点在χ^*(33)工况下x，y，z方向的节点反力，和分别为在步骤C3条件下第i个节点在χ^*(23)工况下x，y，z方向的节点反力，和分别为在步骤C3条件下第i个节点在χ^*(31)工况下x，y，z方向的节点反力，和分别为在步骤C3条件下第i个节点在χ^*(12)工况下x，y，z方向的节点反力，K为单胞有限元模型的总体刚度阵；

C4、根据求得周期性结构等效弹性模量，其中，为周期性结构等效弹性模量，Y为所述单胞有限元模型的体积；

计算周期性结构等效热弹性常数具有如下步骤：

其中，和分别为在步骤D1条件下第i个节点在x，y，z方向的节点反力，B为单胞有限元模型中单元的应变-位移矩阵，E为单胞有限元模型中各个材料的单元刚度阵，α是单胞有限元模型中各个材料的热膨胀系数，Y为所述单胞有限元模型的面积；

Λ＝{Λ₁,Λ₂,…,Λ_i,…Λ_N}^T，

其中，和分别为在步骤D2条件下第i个节点在x，y，z方向的位移，为施加周期性边界条件后单胞有限元模型的总体刚度阵；

其中，和分别为在步骤D3条件下第i个节点在x，y，z方向的节点反力，K为单胞有限元模型的总体刚度阵；

D4、根据求周期性结构等效热弹性常数，其中，为周期性结构等效热弹性常数，Y为所述单胞有限元模型的面积；

S5、计算纤维增强复合材料膨胀系数，纤维增强复合材料膨胀系数通过以下公式求得：

实施例2

如图1所示，计算二维零膨胀结构膨胀系数，具有以下步骤：

S1、通过商业有限元软件建立单胞有限元模型，设置单胞有限元模型中材料参数，如图3所示；

S3、由单胞有限元模型的节点坐标值确定初始位移场为：

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

计算周期性结构等效弹性模量具有如下步骤：

其中，和分别为在步骤A1条件下第i个节点在工况下x，y方向的节点反力，和分别为在步骤A1条件下第i个节点在工况下x，y方向的节点反力，和分别为在步骤A1条件下第i个节点在工况下x，y方向的节点反力，K为单胞有限元模型的总体刚度阵；

其中，和分别为在步骤A2条件下第i个节点在f⁽¹¹⁾工况下x，y方向的特征位移，和分别为在步骤A2条件下第i个节点在f⁽²²⁾工况下x，y方向的特征位移，和分别为在步骤A2条件下第i个节点在f⁽¹²⁾工况下x，y方向的特征位移，为施加周期性边界条件后单胞有限元模型的总体刚度阵；

其中，和分别为在步骤A3条件下第i个节点在χ^*(11)工况下x，y方向的节点反力，和分别为在步骤A3条件下第i个节点在χ^*(22)工况下x，y方向的节点反力，和分别为在步骤A3条件下第i个节点在χ^*(12)工况下x，y方向的节点反力，K为单胞有限元模型的总体刚度阵；

A4、根据求得周期性结构等效弹性模量，其中，为周期性结构等效弹性模量，Y为所述单胞有限元模型的面积；

计算周期性结构等效热弹性常数具有如下步骤：

其中，和分别为在步骤B1条件下第i个节点在x，y方向的节点反力，B为单胞有限元模型中单元的应变-位移矩阵，E为单胞有限元模型中各个材料的单元刚度阵，α是单胞有限元模型中各个材料的热膨胀系数，Y为所述单胞有限元模型的面积；

Λ＝{Λ₁,Λ₂,…,Λ_i,…Λ_N}^T，

其中，和分别为在步骤B2条件下第i个节点在x，y方向的位移，为施加周期性边界条件后单胞有限元模型的总体刚度阵；

其中，和分别为在步骤B3条件下第i个节点在x，y方向的节点反力，K为单胞有限元模型的总体刚度阵；

B4、根据求周期性结构等效热弹性常数，其中，为周期性结构等效热弹性常数，Y为所述单胞有限元模型的面积；

S5、计算二维零膨胀结构膨胀系数，二维零膨胀结构膨胀系数通过以下公式求得：

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

Claims

1.一种渐进均匀化预测周期性复合材料热膨胀系数的简易实现方法，其特征在于具有以下步骤：

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

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

S5、计算周期性复合材料热膨胀系数；

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

当所述周期性复合材料呈二维结构时，所述步骤S4中计算周期性结构等效热弹性常数具有如下步骤：

B1、约束单胞有限元模型的各节点的位移，施加单位负温升，进行有限元热分析，求得节点支撑反力场，并将所述节点支撑反力场作为等效的节点力场，所述等效的节点力场R_α＝∫_YB^TEαdY为：

<mrow> <msub> <mi>R</mi> <mi>&alpha;</mi> </msub> <mo>=</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>R</mi> <mi>&alpha;</mi> <mn>1</mn> </msubsup> <mo>,</mo> <msubsup> <mi>R</mi> <mi>&alpha;</mi> <mn>2</mn> </msubsup> <mo>,</mo> <mo>...</mo> <mo>,</mo> <msubsup> <mi>R</mi> <mi>&alpha;</mi> <mi>i</mi> </msubsup> <mo>,</mo> <mo>...</mo> <msubsup> <mi>R</mi> <mi>&alpha;</mi> <mi>N</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> </mrow>

Λ＝{Λ₁,Λ₂,…,Λ_i,…Λ_N}^T，

<mrow> <msub> <mi>R</mi> <mi>&Lambda;</mi> </msub> <mo>=</mo> <msup> <mrow> <mo>{</mo> <msubsup> <mi>R</mi> <mi>&Lambda;</mi> <mn>1</mn> </msubsup> <mo>,</mo> <msubsup> <mi>R</mi> <mi>&Lambda;</mi> <mn>2</mn> </msubsup> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msubsup> <mi>R</mi> <mi>&Lambda;</mi> <mi>i</mi> </msubsup> <mo>,</mo> <mn>...</mn> <msubsup> <mi>R</mi> <mi>&Lambda;</mi> <mi>N</mi> </msubsup> <mo>}</mo> </mrow> <mi>T</mi> </msup> <mo>,</mo> </mrow>

其中，

ij＝{11,22,33,23,31,12}，x_i,y_i,z_i为第i个节点的坐标值，N表示单胞有限元模型的节点总个数；

当所述周期性复合材料呈三维结构时，所述步骤S4中计算周期性结构等效热弹性常数具有如下步骤：

D1、约束单胞有限元模型的各节点的位移，施加单位负温升，进行有限元热分析，求得节点支撑反力场，并将所述节点支撑反力场作为等效的节点力场，所述等效的节点力场R_α＝∫_YB^TEαdY为：

其中，和分别为在步骤D1条件下第i个节点在x，y，z方向的节点反力，B为单胞有限元模型中单元的应变-位移矩阵，E为单胞有限元模型中各个材料的单元刚度阵，α是单胞有限元模型中各个材料的热膨胀系数，Y为所述单胞有限元模型的体积；

Λ＝{Λ₁,Λ₂,…,Λ_i,…Λ_N}^T，

D4、根据求周期性结构等效热弹性常数，其中，为周期性结构等效热弹性常数，Y为所述单胞有限元模型的体积。

2.根据权利要求1所述的一种渐进均匀化预测周期性复合材料热膨胀系数的简易实现方法，其特征在于：当所述周期性复合材料呈二维结构时，所述步骤S4中计算周期性结构等效弹性模量具有如下步骤：

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

A3、删除所有载荷及约束后，将χ^*(ij)施加到单胞有限元模型的节点上，并对每个工况单独进行有限元静力分析，由f^*(ij)＝Kχ^*(ij)求得χ^*(ij)对应的节点反力场为f^*(ij)为：

3.根据权利要求1所述的一种渐进均匀化预测周期性复合材料热膨胀系数的简易实现方法，其特征在于：当所述周期性复合材料呈三维结构时，所述步骤S4中计算周期性结构等效弹性模量具有如下步骤：

其中，

4.根据权利要求1-3任一权利要求所述的一种渐进均匀化预测周期性复合材料热膨胀系数的简易实现方法，其特征在于：周期性复合材料热膨胀系数通过以下公式求得：

<mrow> <msubsup> <mi>&alpha;</mi> <mrow> <mi>k</mi> <mi>l</mi> </mrow> <mi>H</mi> </msubsup> <mo>=</mo> <msup> <mrow> <mo>&lsqb;</mo> <msubsup> <mi>E</mi> <mrow> <mi>i</mi> <mi>j</mi> <mi>k</mi> <mi>l</mi> </mrow> <mi>H</mi> </msubsup> <mo>&rsqb;</mo> </mrow> <mrow> <mo>-</mo> <mn>1</mn> </mrow> </msup> <msubsup> <mi>&beta;</mi> <mrow> <mi>i</mi> <mi>j</mi> </mrow> <mi>H</mi> </msubsup> <mo>,</mo> </mrow>

其中，为周期性复合材料热膨胀系数。