复杂复合材料结构等效材料性能多尺度计算方法

摘要

本发明提出一种复杂复合材料结构等效材料性能多尺度计算方法，采用尺度分离的方法，将宏观、细观、微观三尺度结构分离，根据不同尺度模型的几何特征，分别建立各个尺度分析模型；将三尺度问题转化为两个多尺度问题：宏观‑细观多尺度问题、细观‑微观多尺度问题，依次对着两个多尺度问题进行分析，将微观多尺度问题得到的等效模量最终返回给宏观多尺度问题。克服了传统结构分析方法计算效率低、精度差的缺点，有效提升了复合材料结构性能预测的效率和精度，使其可以用于指导复合材料的生产、研发等工作。本发明可应用于航空航天领域复杂复合材料结构设计、分析，以及其他复合材料工程领域的结构设计热、力学分析问题。

说明书

技术领域

本发明涉及复合材料设计领域，是一种复杂复合材料结构分析设计方法，具体为一种复杂复合材料结构等效材料性能多尺度计算方法。

背景技术

复合材料由于质量轻、强度高、具有较强的可设计性等特点，广泛应用与航空航天的结构器件中。而由于复合材料结构复杂，为了研究复合材料的性能，提升复合材料结构件的使用效率，近百年来，国内外学者提出了大量用以预测复合材料行为的理论。其核心是通过求解控制方程，从而确定结构内部位移、温度等物理量的分布，从而完成对材料性能的预测。

目前，复合材料材料性能预测方法主要分为四类：

第一类是解析法，其代表的方法有：稀疏法、Mori-Tanaka方法、自洽方法、广义自洽法。该类方法通过求解无限大基体内单夹杂或多夹杂问题，得到远场应变与单个夹杂平均应变之间关系，从而得到材料的有效模量。这种方法理论较为简单，但由于实际复合材料存在一定的边界，边界效应会使得计算结果产生一定的误差，此外，部分解析方法只适用于结构简单，体积分数较低的复合材料，这使得这类方法在复合材料性能预测上存在一定的局限。

第二类是半解析法，其代表的方法为变换场分析法。该方法在细观采用显式的本构关系来联系宏观与细观场，该方法需要给定均匀化与局部化的规则，对于多相材料以及非线性的非均匀材料，该方法等效本构关系所需要的内部变量会非常多，限制了该方法的应用。

第三类是数值方法，其代表的方法为数值均匀化方法，该方法将复合材料转化为一个多尺度分析问题，通过局部化和均匀化的方法，建立宏观积分点与细观体积代表单元之间的联系，从而完成材料性能的预测，与解析法相比，该方法的计算量较小，并且由于在计算时考虑到材料的细观形貌，因此计算精度较高。

然而现有的多尺度方法仅考虑了两尺度的信息，由于大多数复合材料采用铺层的形式，材料的细观尺度并非纤维与基体简单的组合形式，而是多种纤维与基体的组合，纤维的铺层角，每次纤维的排布形式都会很大程度的影响宏观材料的性能。此外，大多数多尺度分析软件用户开发由国外航空航天科研机构开发，由于种种原因这些软件都没有对外公开，并且大多数用于学术研究的程序因精度、计算代价等问题，限制了其在工程领域内的应用。

发明内容

为了避免现有技术的不足之处，本发明提出了一种复杂复合材料结构等效材料性能多尺度计算方法，在方法中采用了复合材料结构分析的三尺度模型，由于该方法考虑了细观、微观结构，从而使得宏观结果计算精度得以提升；此外，该方法通过商业有限元软件ABAQUS的二次开发实现，从而增加了其通用性，使其能够更好地解决大规模的工程问题。

本发明的技术方案为：

所述一种复杂复合材料结构等效材料性能多尺度计算方法，其特征在于：包括以下步骤：

步骤1：按照复合材料实际尺度建立宏观有限元分析模型，宏观有限元分析模型材料坐标系为(X₁,X₂,X₃)；通过显微CT扫描实验，得到复合材料细观结构的物理模型，根据复合材料细观结构物理模型的体积分数、增强相与基体相的几何特征以及排布形式、缺陷位置、铺层数量和铺层角信息，建立细观有限元模型，细观有限元模型材料坐标系记为(Y₁,Y₂,Y₃)；通过电子显微镜实验，得到复合材料微观单胞的物理模型，根据复合材料微观单胞物理模型增强相的体积分数、形状、以及缺陷位置，建立微观有限元模型，微观有限元模型材料坐标系记为(Z₁,Z₂,Z₃)；其中Y_i＝X_i/ξ，Z_i＝Y_i/η，i＝1,2,3，ξ，η分别为宏观-细观，细观-微观尺度间的桥接系数，且满足ξ<＜1，η<＜1；

步骤2：根据需要计算的复合材料，赋予微观有限元模型材料属性；

步骤3：将多尺度分析分为两步，首先通过细观-微观两尺度分析，得到细观尺度的等效材料属性；根据细观尺度的等效材料属性，通过宏观-细观两尺度分析，得到宏观结构的等效材料属性。

进一步的优选方案，所述一种复杂复合材料结构等效材料性能多尺度计算方法，其特征在于：计算的等效材料性能为等效刚度矩阵；步骤3的具体步骤为：

步骤3.1：在周期性假设的条件下，将微观有限元模型的位移渐进展开式带入弹性力学控制方程

$\frac{\partial {\sigma }_{ij}^{\xi }}{\partial x_j}+b_j=0$

${\sigma }_{ij}^{\xi }=C_{ijmn}^{\xi }{ϵ}_{mn}^{\xi }$

${ϵ}_{mn}^{\xi }=\frac{1}{2}(\frac{\partial u_k^{\xi }}{\partial x_l}+\frac{\partial u_l^{\xi }}{\partial x_k})$

中，得到微观等效的刚度表达式：

$C_{ijmn}^H=\frac{1}{\left|Y\right|}\underset{Y}{\int }{C}_{ijkl}[\frac{1}{2}(\frac{\partial {\chi }_k^{mn}}{\partial y_l}+\frac{\partial {\chi }_l^{mn}}{\partial y_k})+\frac{1}{2}({\delta }_{mk}{\delta }_{nl}+{\delta }_{nk}{\delta }_{ml})]dY$

其中，上角标代表微观有限元模型，下角标代表6个不同应力的方向，下角标k，l代表3个不同位移的方向，的上角标代表微观有限元模型均匀化，下角标代表刚度矩阵中6个不同的方向，Y代表单胞体积，为微观位移特征函数，与位移对应k，l代表3个不同的位移特征函数的方向，C_ijkl为单一组分材料的弹性模量，δ_mk为Kronecker张量，且满足：

${\delta }_{mk}=\left(\begin{array}{cc}1& m=k\\ 0& m\ne k\end{array}\right)$

步骤3.2：采用等效热应力加载，将步骤3.1中的微观等效的刚度表达式转化为：

$C_{ijmn}^H=\frac{1}{\left|Y\right|}\underset{Y}{\int }{C}_{ijkl}[\frac{1}{2}(\frac{\partial {\chi }_k^{mn}}{\partial y_l}+\frac{\partial {\chi }_l^{mn}}{\partial y_k})-{ϵ}_{kl}^{mn}]dY$

${ϵ}_{kl}^{mn}={\alpha }_{kl}^{mn}\Delta T$

${\alpha }_{kl}^{mn}=\left(\begin{array}{cccccc}-1& 0& 0& 0& 0& 0\\ 0& -1& 0& 0& 0& 0\\ 0& 0& -1& 0& 0& 0\\ 0& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& -1& 0\\ 0& 0& 0& 0& 0& -1\end{array}\right)$

其中，代表等效热应变大小，为单位热膨胀系数，ΔT为单位温度变化；

步骤3.3：得到步骤3.2中等效的微观有限元模型刚度矩阵后，根据细观有限元模型内每个铺层的铺层角，依据经典层合板理论，得到细观有限元模型每层铺层的等效刚度矩阵，并依此对细观有限元模型的刚度矩阵进行组装，形成总刚度矩阵：

$C_{ijkl}^t=T_t^{-1}{C}_{ijkl}^H{\left(T_t^{-1}\right)}^T$

$C_{ijkl}^{total}=\underset{t=1}{\overset{n}{\Sigma }}{C}_{ijkl}^t$

其中，T_t为细观模型每层铺层的转换矩阵，t＝1,2…n，为细观模型单层铺层在总体坐标系下的刚度矩阵，为细观有限元模型总刚度矩阵；

步骤3.4：将步骤3.3得到的细观有限元模型总刚度矩阵赋予宏观有限元分析模型中，并对宏观有限元分析模型施加载荷，得到宏观有限元分析模型的响应。

进一步的优选方案，所述一种复杂复合材料结构等效材料性能多尺度计算方法，其特征在于：计算的等效材料性能为等效热传导系数；步骤3的具体步骤为：

步骤3.1：在周期性假设的条件下，将微观有限元模型的温度渐进展开式带入稳态热传导控制方程

$\frac{\partial }{\partial x_i}\left(k_{ij}^{\xi }\right(x)\frac{\partial T^{\xi }}{\partial x_j})+\rho Q=0$

$T^{\xi }=\bar{T}\left(x\right)$

中，得到微观等效的热传导系数表达式：

$K_{ip}^H=\frac{1}{\left|Y\right|}\underset{Y}{\int }{k}_{ij}^{\xi }({\delta }_{jp}+\frac{\partial H^p}{\partial y_j})dY$

其中，代表微观模型热传导系数，上角标代表微观尺度，下角标代表3个不同的方向，ρ代表材料密度，Y代表单胞体积，Q代表内热流密度，代表微观模型温度边界条件；δ_jp代表Kronecker张量；

步骤3.2：采用等效热应变加载，将步骤3.1中的微观等效的的热传导系数表达式转化为：

$K_{ip}^H=\frac{1}{\left|Y\right|}\underset{Y}{\int }{k}_{ij}[\frac{\partial H_k^p}{\partial y_l}-{ϵ}_{kl}^{mn}]dY$

${ϵ}_{kl}^{mn}={\alpha }_{kl}^{mn}\Delta T$

${\alpha }_{kl}^{mn}=\left(\begin{array}{cccccc}-1& 0& 0& 0& 0& 0\\ 0& -1& 0& 0& 0& 0\\ 0& 0& -1& 0& 0& 0\\ 0& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& -1& 0\\ 0& 0& 0& 0& 0& -1\end{array}\right)$

其中，代表等效热应变大小，为单位热膨胀系数，ΔT为单位温度变化；

步骤3.3：得到步骤3.2中等效的微观有限元模型热传导系数矩阵后，根据细观有限元模型内每个铺层的铺层角，依据经典层合板理论，得到细观有限元模型每层铺层的等效热传导系数矩阵，并依此对细观有限元模型的热传导系数矩阵进行组装，形成总热传导系数矩阵：

$K_{ip}^t=T_t^{-1}{K}_{iP}^H{\left(T_t^{-1}\right)}^T$

$K_{ip}^{total}=\underset{t=1}{\overset{n}{\Sigma }}{K}_{ip}^t$

其中，T_t为细观模型每层铺层的转换矩阵，t＝1,2…n，为细观模型单层铺层在总体坐标系下的热传导系数矩阵，为细观有限元模型总热传导系数矩阵；

步骤3.4：将步骤3.3得到的细观有限元模型总热传导系数矩阵赋予宏观有限元分析模型中，并对宏观有限元分析模型施加载荷，得到宏观有限元分析模型的响应。

进一步的优选方案，所述一种复杂复合材料结构等效材料性能多尺度计算方法，其特征在于：计算的等效材料性能为等效热膨胀系数；步骤3的具体步骤为：

步骤3.1：微观等效的热弹性常数均匀化计算式为：

${\beta }_{ij}^H=\frac{1}{\left|Y\right|}{\int }_Y{C}_{ijkl}^{\zeta }({\alpha }_{kl}+{\theta }_{k,y_l})dY$

步骤3.2：采用等效热应变加载，计算等效热弹性系数为

${\beta }_{ij}^H=\frac{1}{\left|Y\right|}\underset{Y}{\int }{C}_{ijkl}^{\xi }[\frac{\partial {\theta }_k}{\partial y_l}-{ϵ}_{kl}]dY$

得到等效热膨胀系数的计算式为

步骤3.3：得到步骤3.2中等效的微观有限元模型热膨胀系数后，根据细观有限元模型内每个铺层的铺层角，依据经典层合板理论，得到细观有限元模型每层铺层的等效热膨胀系数，并依此对细观有限元模型的热膨胀系数进行组装，形成热膨胀系数矩阵：

${\alpha }_{kl}^t=T_t^{-1}{\alpha }_{kl}^H{\left(T_t^{-1}\right)}^T$

${\alpha }_{kl}^{total}=\underset{t=1}{\overset{n}{\Sigma }}{\alpha }_{kl}^t$

其中，T_t为细观模型每层铺层的转换矩阵，t＝1,2…n，为细观模型单层铺层在总体坐标系下的热膨胀系数，为细观有限元模型热膨胀系数矩阵；

步骤3.4：将步骤3.3得到的细观有限元模型热膨胀系数矩阵赋予宏观有限元分析模型中，并对宏观有限元分析模型施加载荷，得到宏观有限元分析模型的响应。

有益效果

本发明提出的三尺度复合材料分析方法，有益效果是：

1、利用了多尺度方法，在复合材料分析过程中，充分考虑了细观、微观结构几何形貌对于宏观结构的影响，与传统的复合材料分析手段相比有更好地精度，此外，对于损伤、失效判断，可以通过观测细观结构应力分布的变化，明确损伤机理。

2、通过建立三尺度模型，使得在复合材料等效性质的计算时，考虑了纤维方向、铺层厚度因素对于整体刚度矩阵的影响，而传统的两尺度方法忽略了细观尺度铺层的厚度、纤维的方向因素。因此与传统多尺度方法相比有更好的精度。

3、三尺度方法能够通过基于ABAQUS平台的二次开发得以实现，具有较好的适用性，从而促进了多尺度方法在工程材料计算领域的应用。

本发明的附加方面和优点将在下面的描述中部分给出，部分将从下面的描述中变得明显，或通过本发明的实践了解到。

附图说明

本发明的上述和/或附加的方面和优点从结合下面附图对实施例的描述中将变得明显和容易理解，其中：

图1：本发明的计算流程图；

图2：实施例中某型压力容器几何模型；

图3：压力容器横截面图；

图4：CT扫描下预制体的微观形貌；

图5：简化后的细观有限元模型；

图6：CT扫描下无纬布的微观形貌；

图7：针刺、网胎微观计算模型；

图8：无纬布微观计算模型；

图9：压力容器边界条件。

具体实施方式

下面详细描述本发明的实施例，所述实施例是示例性的，旨在用于解释本发明，而不能理解为对本发明的限制。

本实施例以某型压力容器等效材料性能计算为例，按照本发明技术方案进行实施，给出了详细的实施过程。

步骤1：压力容器由碳/碳复合材料构成，细观结构的预制体由无纬布、±45°铺层、0°、90°铺层组合而成。按照算例实际尺寸，如图2和图3所示，圆柱体长20mm。在商用有限元软件ABAQUS中建立压力容器宏观有限元分析模型，宏观有限元分析模型材料坐标系为(X₁,X₂,X₃)。通过CT扫描和电镜扫描，分别得到压力容器结构的真实细观、微观模型。如图4和图6所示。

根据针刺碳/碳复合材料预制件微结构显微照片分析，可以确定微结构细观单胞的基本形式。细观单胞由若干不同铺层角的无纬布和复合网胎叠层铺设而成，在厚度方向通过针刺纤维束进行增强。网胎纤维在面内是杂乱分布的，因此是一种面内准各向同性材料。针刺和网胎类似，也属于各项同性材料。0°无纬布、环向90°无纬布及斜向无纬布的纤维排列较为紧凑，纤维体积分数较大。据此建立细观有限元模型，如图5。细观有限元模型材料坐标系记为(Y₁,Y₂,Y₃)。

根据显微扫描得到材料的两种体积分数，在ABAQUS有限元软件中建立两种微观单胞，单胞1纤维体积分数为50％，用来模拟纤维含量较少的网胎以及针刺微观有限元模型，如图7。单胞2体积分数为81％，用来模拟纤维较为紧凑的无纬布微观有限元模型，如图8。微观有限元模型材料坐标系记为(Z₁,Z₂,Z₃)。

Y_i＝X_i/ξ，Z_i＝Y_i/η，i＝1,2,3，ξ，η分别为宏观-细观，细观-微观尺度间的桥接系数，且满足ξ<＜1，η<＜1。

步骤2：根据需要计算的复合材料，赋予微观有限元模型材料属性。

步骤3：将多尺度分析分为两步，首先通过细观-微观两尺度分析，得到细观尺度的等效材料属性；根据细观尺度的等效材料属性，通过宏观-细观两尺度分析，得到宏观结构的等效材料属性。

当计算的等效材料性能为等效刚度矩阵，步骤3的具体步骤为：

步骤3.1：在ABAQUS中施加Tie约束，从而实现周期性边界条件的施加，使得对应面位移相同。在周期性假设的条件下，将微观有限元模型的位移渐进展开式带入弹性力学控制方程

$\frac{\partial {\sigma }_{ij}^{\xi }}{\partial x_j}+b_j=0$

${\sigma }_{ij}^{\xi }=C_{ijmn}^{\xi }{ϵ}_{mn}^{\xi }$

${ϵ}_{mn}^{\xi }=\frac{1}{2}(\frac{\partial u_k^{\xi }}{\partial x_l}+\frac{\partial u_l^{\xi }}{\partial x_k})$

中，得到微观等效的刚度表达式：

$C_{ijmn}^H=\frac{1}{\left|Y\right|}\underset{Y}{\int }{C}_{ijkl}[\frac{1}{2}(\frac{\partial {\chi }_k^{mn}}{\partial y_l}+\frac{\partial {\chi }_l^{mn}}{\partial y_k})+\frac{1}{2}({\delta }_{mk}{\delta }_{nl}+{\delta }_{nk}{\delta }_{ml})]dY$

其中，上角标代表微观有限元模型，下角标代表6个不同应力的方向，下角标k，l代表3个不同位移的方向，的上角标代表微观有限元模型均匀化，下角标代表刚度矩阵中6个不同的方向，Y代表单胞体积，为微观位移特征函数，与位移对应k，l代表3个不同的位移特征函数的方向，C_ijkl为单一组分材料的弹性模量，δ_mk为Kronecker张量，且满足：

${\delta }_{mk}=\left(\begin{array}{cc}1& m=k\\ 0& m\ne k\end{array}\right).$

这里实施例中在ABAQUS中设置6个线性扰动分析步，从而完成热载荷在不同方向(11、22、33、12、13、23)的加载。

步骤3.2：采用等效热应力加载，将步骤3.1中的微观等效的刚度表达式转化为：

$C_{ijmn}^H=\frac{1}{\left|Y\right|}\underset{Y}{\int }{C}_{ijkl}[\frac{1}{2}(\frac{\partial {\chi }_k^{mn}}{\partial y_l}+\frac{\partial {\chi }_l^{mn}}{\partial y_k})-{ϵ}_{kl}^{mn}]dY$

${ϵ}_{kl}^{mn}={\alpha }_{kl}^{mn}\Delta T$

${\alpha }_{kl}^{mn}=\left(\begin{array}{cccccc}-1& 0& 0& 0& 0& 0\\ 0& -1& 0& 0& 0& 0\\ 0& 0& -1& 0& 0& 0\\ 0& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& -1& 0\\ 0& 0& 0& 0& 0& -1\end{array}\right)$

其中，代表等效热应变大小，为单位热膨胀系数，ΔT为单位温度变化。

将微观单胞的各个方向分析结果均匀化后，由于各个等效热载荷均为单位1载荷，因此均匀化后，得到微观结构等效材料属性，如表1，表2，至此微观分析结束。

表1：多尺度方法得到的针刺和网胎微观模型等效刚度矩阵

表2：微观多尺度方法得到的无纬布的等效刚度矩阵

步骤3.3：得到步骤3.2中等效的微观有限元模型刚度矩阵后，根据细观有限元模型内每个铺层的铺层角，依据经典层合板理论，得到细观有限元模型每层铺层的等效刚度矩阵，并依此对细观有限元模型的刚度矩阵进行组装，形成总刚度矩阵：

$C_{ijkl}^t=T_t^{-1}{C}_{ijkl}^H{\left(T_t^{-1}\right)}^T$

$C_{ijkl}^{total}=\underset{t=1}{\overset{n}{\Sigma }}{C}_{ijkl}^t$

其中，T_t为细观模型每层铺层的转换矩阵，t＝1,2…n，为细观模型单层铺层在总体坐标系下的刚度矩阵，为细观有限元模型总刚度矩阵。

本实施例中，45°度铺层材料属性如表4所示，-45°度铺层材料属性如表5所示，90°铺层材料属性如表3所示，细观喷管压力容器等效材料属性如表6所示。

表3：经过坐标转换后得到的细观90°无纬布等效模量

表4：经过坐标转换后得到的细观45°无纬布等效模量

表5：经过坐标转换后得到的细观-45°无纬布等效模量

表6：多尺度方法得到的细观模型等效刚度矩阵

步骤3.4：将步骤3.3得到的细观有限元模型总刚度矩阵赋予宏观有限元分析模型中，并对宏观有限元分析模型施加载荷，本实施例中，在宏观模型左端面施加2方向载荷50MPa，约束右端面11、22、33三方向位移以及转角如图9所示，得到宏观有限元分析模型的响应。

当计算的等效材料性能为等效热传导系数，步骤3的具体步骤为：

步骤3.1：在周期性假设的条件下，将微观有限元模型的温度渐进展开式带入稳态热传导控制方程

$\frac{\partial }{\partial x_i}\left(k_{ij}^{\xi }\right(x\left)\frac{\partial T^{\xi }}{\partial x_j}\right)+\rho Q=0$

$T^{\xi }=\bar{T}\left(x\right)$

中，得到微观等效的热传导系数表达式：

$K_{ip}^H=\frac{1}{\left|Y\right|}\underset{Y}{\int }{k}_{ij}^{\xi }({\delta }_{ip}+\frac{\partial H^p}{\partial y_j})dY$

其中，代表微观模型热传导系数，上角标代表微观尺度，下角标代表3个不同的方向，ρ代表材料密度，Y代表单胞体积，Q代表内热流密度，代表微观模型温度边界条件；δ_jp代表Kronecker张量；

步骤3.2：采用等效热应变加载，将步骤3.1中的微观等效的的热传导系数表达式转化为：

$K_{ip}^H=\frac{1}{\left|Y\right|}\underset{Y}{\int }{k}_{ij}[\frac{\partial H_k^p}{\partial y_l}-{ϵ}_{kl}^{mn}]dY$

${ϵ}_{kl}^{mn}={\alpha }_{kl}^{mn}\Delta T$

${\alpha }_{kl}^{mn}=\left(\begin{array}{cccccc}-1& 0& 0& 0& 0& 0\\ 0& -1& 0& 0& 0& 0\\ 0& 0& -1& 0& 0& 0\\ 0& 0& 0& -1& 0& 0\\ 0& 0& 0& 0& -1& 0\\ 0& 0& 0& 0& 0& -1\end{array}\right)$

其中，代表等效热应变大小，为单位热膨胀系数，ΔT为单位温度变化；

步骤3.3：得到步骤3.2中等效的微观有限元模型热传导系数矩阵后，根据细观有限元模型内每个铺层的铺层角，依据经典层合板理论，得到细观有限元模型每层铺层的等效热传导系数矩阵，并依此对细观有限元模型的热传导系数矩阵进行组装，形成总热传导系数矩阵：

$K_{ip}^t=T_t^{-1}{K}_{iP}^H{\left(T_t^{-1}\right)}^T$

$K_{ip}^{total}=\underset{t=1}{\overset{n}{\Sigma }}{K}_{ip}^t$

其中，T_t为细观模型每层铺层的转换矩阵，t＝1,2…n，为细观模型单层铺层在总体坐标系下的热传导系数矩阵，为细观有限元模型总热传导系数矩阵；

步骤3.4：将步骤3.3得到的细观有限元模型总热传导系数矩阵赋予宏观有限元分析模型中，并对宏观有限元分析模型施加载荷，得到宏观有限元分析模型的响应。

当计算的等效材料性能为等效热膨胀系数，步骤3的具体步骤为：

步骤3.1：考虑到在周期性假设的热弹性力学边界值问题中，应变的表达式为：

${ϵ}_{kl}^T={ϵ}_{kl}^{\zeta }+{\alpha }_{kl}^{\zeta }\Delta T=\frac{1}{2}\left(\frac{\partial u_k^{\zeta }}{\partial x_l}+\frac{\partial u_l^{\zeta }}{\partial x_k}\right)$

其中，为微观总应变，为微观机械应变，为微观热膨胀系数。因此热弹性力学边界值问题的本构方程为：

${\sigma }_{ij}^{\zeta }=C_{ijkl}^{\zeta }{ϵ}_{kl}^{\zeta }=C_{ijkl}^{\zeta }({ϵ}_{kl}^T-{\alpha }_{kl}^{\zeta }\Delta T)=C_{ijkl}^{\zeta }\left[\frac{1}{2}\left(\frac{\partial u_k^{\zeta }}{\partial x_l}+\frac{\partial u_l^{\zeta }}{\partial x_k}\right)-{\alpha }_{kl}^{\zeta }\Delta T\right]$

记微观热弹性常数张量为因此本构方程可以写成除应变和本构方程外，热弹性力学边界值问题的控制方程与弹性力学边界值问题的控制方程是完全相同的。因此，基于求解等效刚度矩阵以及等效热传导系数的推导过程，得到热弹性系数微观边界值问题

其中，θ代表热膨胀问题中待求解的微观热弹性特征函数，并且与求解微观等效刚度矩阵中的特征函数χ之间有如下关系：

$\frac{\partial {\theta }_k}{\partial y_l}={\alpha }_{mn}^{\xi }\frac{\partial {\chi }_k^{mn}}{\partial y_l}$

最终得到微观等效热弹性常数的均匀化计算式：

${\beta }_{ij}^H=\frac{1}{\left|Y\right|}{\int }_Y{C}_{ijkl}^{\zeta }({\alpha }_{kl}+{\theta }_{k,y_l})dY;$

步骤3.2：采用等效热应变加载，计算等效热弹性系数为

${\beta }_{ij}^H=\frac{1}{\left|Y\right|}\underset{Y}{\int }{C}_{ijkl}^{\xi }[\frac{\partial {\theta }_k}{\partial y_l}-{ϵ}_{kl}]dY$

得到等效热膨胀系数的计算式为

步骤3.3：得到步骤3.2中等效的微观有限元模型热膨胀系数后，根据细观有限元模型内每个铺层的铺层角，依据经典层合板理论，得到细观有限元模型每层铺层的等效热膨胀系数，并依此对细观有限元模型的热膨胀系数进行组装，形成热膨胀系数矩阵：

${\alpha }_{kl}^t=T_t^{-1}{\alpha }_{kl}^H{\left(T_t^{-1}\right)}^T$

${\alpha }_{kl}^{total}=\underset{t=1}{\overset{n}{\Sigma }}{\alpha }_{kl}^t$

其中，T_t为细观模型每层铺层的转换矩阵，t＝1,2…n，为细观模型单层铺层在总体坐标系下的热膨胀系数，为细观有限元模型热膨胀系数矩阵；

步骤3.4：将步骤3.3得到的细观有限元模型热膨胀系数矩阵赋予宏观有限元分析模型中，并对宏观有限元分析模型施加载荷，得到宏观有限元分析模型的响应。

尽管上面已经示出和描述了本发明的实施例，可以理解的是，上述实施例是示例性的，不能理解为对本发明的限制，本领域的普通技术人员在不脱离本发明的原理和宗旨的情况下在本发明的范围内可以对上述实施例进行变化、修改、替换和变型。