热-力耦合问题的高阶多尺度分析

Higher order multiscale analysis for thermo-mechanical coupling problems

摘要

本文针对拟周期性复合材料结构的热-力耦合问题提出了一种高阶多尺度分析方法， 基于稳态非线性热-力耦合控制方程发展了适用于拟周期材料热力学行为预测的二阶双尺度渐近展开理论，并设计了相应的有限元算法。数值算例表明，相较传统有限元方法，本文的方法能够在保持计算精度的同时大大节约计算资源，因而具有重要工程应用价值。

Abstract

In this paper higher-order multiscale analysis method for thermo-mechanical coupling problems in quasi-periodic composite structures is considered.A second-order two-scale （SOTS） asymptotic expansion method is developed based on the steady-state nonlinear thermo-mechanical coupling governing equations for effectively predicting the thermomechanical behavior of quasi-periodic materials.In the method， multiscale asymptotic expansions of the temperature and displacement fields is firstly utilized， then the first-order and second-order cell functions and homogenized coefficients are computed to establish corresponding homogenized equations， finally the second-order two-scale approximate solutions are obtained.A specialized finite element algorithm is designed to address the nonlinear nature of the problem through representative macroscopic temperature sampling， computation of temperature-dependent cell functions and homogenized coefficients， interpolation-based determination of homogenized material parameters and cell functions， and direct iteration method for solving homogenized equations.Numerical examples demonstrate that the proposed method can achieve remarkable computational efficiency while maintaining excellent accuracy in comparison with the conventional finite element approaches.This significant improvement in computational performance makes the method particularly valuable for large-scale engineering applications， especially in nuclear reactor safety assessment where both accuracy and efficiency are crucial for thermal-mechanical analysis.Furthermore， the method has great potential for optimal design applications in advanced composite materials， where repeated multiscale simulations are often required for parameter optimization and performance evaluation.

Graphical abstract

关键词

拟周期性 / 复合材料 / 热-力耦合问题 / 二阶双尺度渐进展开

Key words

quasi-periodic / composite materials / thermo-mechanical coupling / second-order two-scale asymptotic expansion（2020 MSC 65M60）

核反应堆极端环境（高温、辐照、复杂热-力载荷）对材料性能提出了严苛要求，使用传统均匀材料难以兼顾耐辐照性、热管理能力与力学可靠性。近年来，拟周期复合材料（Quasi-periodic Composites）因其独特的长程有序但非严格周期的微结构设计在核反应堆燃料元件、包壳涂层及堆内构件制造中展现出显著优势。例如，高温气冷堆中TRISO颗粒的拟周期分布可分散辐照损伤，梯度功能材料（Functionally Graded Material， FGM）的拟周期成分调制则能有效缓解热应力集中。

拟周期复合材料的多场耦合行为涉及跨尺度效应，即微观拟周期结构的局部热-力响应与宏观性能之间存在关联，因而传统均匀化方法或单尺度数值模拟难以精确刻画其物理机制。针对这一问题，本文提出了一种求解拟周期性复合材料结构稳态热-力耦合问题的二阶双尺度方法。

二阶双尺度（Second-Order Two-Scale， SOTS）方法^［1-6］是一种基于渐近均匀化方法^［7-9］的多尺度计算方法^［10-11］，常被用于分析具有周期性微结构的复合材料或多孔介质。该方法的核心思想是，通过建立宏观和微观两个尺度的耦合模型，在保证计算效率的同时，准确捕捉材料的微观结构效应。对定义在周期复合材料上的热-力耦合问题的高阶多尺度方法的研究，目前已有很多^［12-14］。

在实际工程应用中，受到材料制备工艺差异和服役过程中疲劳损伤累积等因素的影响，原本设计为周期性的复合材料，其物理参数和力学参数不再具有严格的周期性特征，而会随宏观位置的变化而改变，此即本文需要研究的拟周期复合材料。在拟周期结构材料的渐近分析研究方面，Cui等^［15-17］对于拟周期结构复合材料的稳态热传导问题中的温度场及稳态线弹性问题中的位移场进行了渐近展开分析，结合数值算例验证了理论结果。Dong等^［18］研究了直角坐标系下拟周期复合材料的动态热-力耦合问题，提出了相应的二阶双尺度数值算法。

在实际工程应用中，通常材料参数不会呈现简单的线性性质，而是具有依赖于温度变化的性质。为了精确分析复合材料结构的热-力耦合行为，必须考虑复合材料参数依赖温度变化这一性质^［19-22］。Dong等^［23］考虑了具有温度依赖性材料参数和非线性焦耳热效应的周期复合材料结构热-电耦合问题的高阶多尺度计算方法。本文对拟周期复合材料的非线性稳态热-力耦合问题进行了二阶双尺度渐近分析并给出相应的有限元算法。除非另有说明，全文统一采用Einstein求和约定：当同一项中出现重复指标时，默认对该指标在所有维度上求和。

1 模型

拟周期性复合材料结构下的非线性稳态热-力耦合模型可以用如下问题表述：

\left\{\begin{array}{l}-\frac{\partial }{\partial x_j}\left(C_{ijkl}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)\frac{\partial u_k^{\epsilon }\left(\mathit{x}\right)}{\partial x_l}-\right.\\ \left.{\beta }_{ij}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)\left(T^{\epsilon }\left(\mathit{x}\right)-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)\right)=f_k\left(\mathit{x}\right)\mathrm{i}\mathrm{n}\mathrm{\Omega },\\ -\frac{\partial }{\partial x_i}\left(k_{ij}\left(\mathit{x},\frac{\mathit{x}}{\epsilon },T^{\epsilon }\right)\frac{\partial T^{\epsilon }\left(\mathit{x}\right)}{\partial x_j}\right)=h\left(\mathit{x}\right)\mathrm{i}\mathrm{n}\mathrm{\Omega },\\ {\mathit{u}}^{\epsilon }\left(\mathit{x}\right)=\mathrm{0},T^{\epsilon }\left(\mathit{x}\right)=T_{\mathrm{r}\mathrm{e}\mathrm{f}}\mathrm{o}\mathrm{n}\partial \mathrm{\Omega }\end{array}\right.(1)

其中，\mathrm{\Omega }\subset {\mathbf{ℝ}}^n是有界区域，具有Lipschitz连续边界，\partial \mathrm{\Omega }是\mathrm{\Omega }的边界；{\mathit{u}}^{\epsilon }\left(\mathit{x}\right)={\left(u_i^{\epsilon }\left(\mathit{x}\right)\right)}_{\mathrm{1}\le i\le n}和T^{\epsilon }\left(\mathit{x}\right)分别表示待求解的位移场和温度场，\epsilon是小周期参数；\mathit{f}\left(\mathit{x}\right)=\left(f_i\left(\mathit{x}\right)\right)\in {\left(L^{\mathrm{2}}\left(\mathrm{\Omega }\right)\right)}^n为体积力，h\left(\mathit{x}\right)\in L^{\mathrm{2}}\left(\mathrm{\Omega }\right)为热源，T_{\mathrm{r}\mathrm{e}\mathrm{f}}为参考温度。此外，对于材料参数，C_{ijkl}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)表示四阶弹性张量，{\beta }_{ij}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)表示二阶热模量张量，k_{ij}\left(\mathit{x},\frac{\mathit{x}}{\epsilon },T^{\epsilon }\right)表示二阶热传导系数张量。我们假设材料参数满足以下两个条件：

（A1） 令\mathit{y}=\frac{\mathit{x}}{\epsilon }为微观尺度单胞Y={\left[\mathrm{0,1}\right]}^n局部坐标，则

C_{ijkl}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)=C_{ijkl}\left(\mathit{x},\mathit{y}\right),{\beta }_{ij}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)={\beta }_{ij}\left(\mathit{x},\mathit{y}\right),k_{ij}\left(\mathit{x},\frac{\mathit{x}}{\epsilon },T^{\epsilon }\right)=k_{ij}\left(\mathit{x},\mathit{y},T^{\epsilon }\right)

均关于微观变量\mathit{y}是1-周期函数；

（A2） C_{ijkl}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right),{\beta }_{ij}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right),k_{ij}\left(\mathit{x},\frac{\mathit{x}}{\epsilon },T^{\epsilon }\right)均为有界可测函数，满足：

\begin{array}{l}{C}_{ijkl}\left(\mathit{x},\mathit{y}\right)=C_{ijlk}\left(\mathit{x},\mathit{y}\right)=C_{jikl}\left(\mathit{x},\mathit{y}\right)=C_{klij}\left(\mathit{x},\mathit{y}\right),\\ {\beta }_{ij}\left(\mathit{x},\mathit{y}\right)={\beta }_{ji}\left(\mathit{x},\mathit{y}\right),k_{ij}\left(\mathit{x},\mathit{y},T^{\epsilon }\right)=k_{ji}\left(\mathit{x},\mathit{y},T^{\epsilon }\right),\end{array}

并且存在常数{\alpha }_{\mathrm{1}}\mathrm{、}{\alpha }_{\mathrm{2}}\mathrm{、}{\chi }_{\mathrm{1}}\mathrm{、}{\chi }_{\mathrm{2}}\mathrm{、}{\kappa }_{\mathrm{1}}\mathrm{、}{\kappa }_{\mathrm{2}}满足：

{\alpha }_{\mathrm{1}}{\eta }_{ij}{\eta }_{ij}\le C_{ijkl}\left(\mathit{x},\mathit{y}\right){\eta }_{ij}{\eta }_{kl}\le {\alpha }_{\mathrm{2}}{\eta }_{ij}{\eta }_{ij},\mathrm{0}{\alpha }_{\mathrm{1}}\le {\alpha }_{\mathrm{2}},

\begin{array}{l}{\chi }_{\mathrm{1}}{\left|\xi \right|}^{\mathrm{2}}\le {\beta }_{ij}\left(\mathit{x},\mathit{y}\right){\xi }_i{\xi }_j\le {\chi }_{\mathrm{2}}{\left|\xi \right|}^{\mathrm{2}},\mathrm{0}{\chi }_{\mathrm{1}}\le {\chi }_{\mathrm{2}},\\ {\kappa }_{\mathrm{1}}{\left|\xi \right|}^{\mathrm{2}}\le k_{ij}\left(\mathit{x},\mathit{y},T^{\epsilon }\right){\xi }_i{\xi }_j\le {\kappa }_{\mathrm{2}}{\left|\xi \right|}^{\mathrm{2}},\mathrm{0}\\ {\kappa }_{\mathrm{1}}\le {\kappa }_{\mathrm{2}},\end{array}

其中，\left\{{\eta }_{ij}\right\}\in {\mathbf{ℝ}}^{n\times n}是任意对称矩阵，\xi ={\left({\xi }_{\mathrm{1}},\cdots ,{\xi }_n\right)}^{\mathrm{T}}\in {\mathbf{ℝ}}^n是任意实向量。

2 二阶双尺度分析

首先，假设位移场{\mathit{u}}^{\epsilon }\left(\mathit{x}\right)={\left(u_i^{\epsilon }\left(\mathit{x}\right)\right)}_{\mathrm{1}\le i\le n}和温度场T^{\epsilon }\left(\mathit{x}\right)有如下渐近展开形式：

\left\{\begin{array}{l}{u}_k^{\epsilon }\left(\mathit{x}\right)=u_k^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y}\right)+\epsilon u_k^{\left(\mathrm{1}\right)}\left(\mathit{x},\mathit{y}\right)+\\ {\epsilon }^{\mathrm{2}}{u}_k^{\left(\mathrm{2}\right)}\left(\mathit{x},\mathit{y}\right)+O\left({\epsilon }^{\mathrm{3}}\right),\\ T^{\epsilon }\left(\mathit{x}\right)=T^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y}\right)+\epsilon T^{\left(\mathrm{1}\right)}\left(\mathit{x},\mathit{y}\right)+\\ {\epsilon }^{\mathrm{2}}{T}^{\left(\mathrm{2}\right)}\left(\mathit{x},\mathit{y}\right)+O\left({\epsilon }^{\mathrm{3}}\right)\end{array}\right.（2）

根据如下的多指标记号泰勒公式：

f\left(x_{\mathrm{0}},y_{\mathrm{0}},z_{\mathrm{0}}+h\right)=f\left(x_{\mathrm{0}},y_{\mathrm{0}},z_{\mathrm{0}}\right)+D^{\left(\mathrm{0,0},\mathrm{1}\right)}f\left(x_{\mathrm{0}},y_{\mathrm{0}},z_{\mathrm{0}}\right)h+\frac{\mathrm{1}}{\mathrm{2}}{D}^{\left(\mathrm{0,0},\mathrm{2}\right)}f\left(x_{\mathrm{0}},y_{\mathrm{0}},z_{\mathrm{0}}\right)h^{\mathrm{2}}+O\left(h^{\mathrm{3}}\right)

对依赖于温度变化的材料参数进行渐近展开，可以得到如下材料参数k_{ij}\left(\mathit{x},\frac{\mathit{x}}{\epsilon },T^{\epsilon }\right)的渐近展开形式：

\begin{array}{l}{k}_{ij}\left(\mathit{x},\mathit{y},T^{\epsilon }\right)=k_{ij}(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}+\epsilon T^{\left(\mathrm{1}\right)}+{\epsilon }^{\mathrm{2}}{T}^{\left(\mathrm{2}\right)}+O({\epsilon }^{\mathrm{3}}))=k_{ij}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)+D^{\left(\mathrm{0,0},\mathrm{1}\right)}{k}_{ij}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\left[\epsilon T^{\left(\mathrm{1}\right)}+{\epsilon }^{\mathrm{2}}{T}^{\left(\mathrm{2}\right)}+O\left({\epsilon }^{\mathrm{3}}\right)\right]+\\ \frac{\mathrm{1}}{\mathrm{2}}{D}^{\left(\mathrm{0,0},\mathrm{2}\right)}{k}_{ij}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)·\\ {\left[\epsilon T^{\left(\mathrm{1}\right)}+{\epsilon }^{\mathrm{2}}{T}^{\left(\mathrm{2}\right)}+O\left({\epsilon }^{\mathrm{3}}\right)\right]}^{\mathrm{2}}+\cdots \approx \\ k_{ij}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)+\epsilon T^{\left(\mathrm{1}\right)}{D}^{\left(\mathrm{0,0},\mathrm{1}\right)}{k}_{ij}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)+\\ {\epsilon }^{\mathrm{2}}\left[T^{\left(\mathrm{2}\right)}{D}^{\left(\mathrm{0,0},\mathrm{1}\right)}{k}_{ij}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)+\right.\\ \left.\frac{\mathrm{1}}{\mathrm{2}}{\left(T^{\left(\mathrm{1}\right)}\right)}^{\mathrm{2}}{D}^{\left(\mathrm{0,0},\mathrm{2}\right)}{k}_{ij}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\right]=k_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)+\epsilon k_{ij}^{\left(\mathrm{1}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)+{\epsilon }^{\mathrm{2}}{k}_{ij}^{\left(\mathrm{2}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\end{array}（3）

其中，偏微分算子D^{\left({\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},{\alpha }_{\mathrm{3}}\right)}=\frac{{\partial }^{\left|\alpha \right|}}{{\partial }_{x_{\mathrm{1}}}^{{\alpha }_{\mathrm{1}}}{\partial }_{x_{\mathrm{2}}}^{{\alpha }_{\mathrm{2}}}{\partial }_{x_{\mathrm{3}}}^{{\alpha }_{\mathrm{3}}}},n重指标\alpha =\left({\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},{\alpha }_{\mathrm{3}}\right)\in {\mathbb{N}}^{\mathrm{3}},\left|\alpha \right|={\alpha }_{\mathrm{1}}+{\alpha }_{\mathrm{2}}+{\alpha }_{\mathrm{3}}\mathrm{。}将式（2）和式（3）代入问题（1），利用链式求导法则，

\frac{\partial }{\partial x_i}\to \frac{\partial }{\partial x_i}+{\epsilon }^{-\mathrm{1}}\frac{\partial }{\partial y_i}（4）

并分别比较\epsilon各幂次系数，对O\left({\epsilon }^{-\mathrm{2}}\right)项有

\left\{\begin{array}{l}-\frac{\partial }{\partial y_j}\left(C_{ijkl}\frac{\partial u_k^{\left(\mathrm{0}\right)}}{\partial y_l}\right)=\mathrm{0},\\ -\frac{\partial }{\partial y_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\frac{\partial T^{\left(\mathrm{0}\right)}}{\partial y_j}\right)=\mathrm{0}\end{array}\right.（5）

由式（5）可知，u_k^{\left(\mathrm{0}\right)}和T^{\left(\mathrm{0}\right)}均与微观尺度\mathit{y}无关，从而有

u_k^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y}\right)=u_k^{\left(\mathrm{0}\right)}\left(\mathit{x}\right),T^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y}\right)=T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)（6）

对于O\left({\epsilon }^{-\mathrm{1}}\right)项，有

\left\{\begin{array}{l}\begin{array}{l}-\frac{\partial }{\partial y_j}\left(C_{ijkl}\frac{\partial u_k^{\left(\mathrm{1}\right)}}{\partial y_l}\right)-\frac{\partial }{\partial x_j}\left(C_{ijkl}\frac{\partial u_k^{\left(\mathrm{0}\right)}}{\partial y_l}\right)-\\ \frac{\partial }{\partial y_j}\left(C_{ijkl}\frac{\partial u_k^{\left(\mathrm{0}\right)}}{\partial x_l}\right)+\frac{\partial }{\partial y_j}\left({\beta }_{ij}\left(T^{\left(\mathrm{0}\right)}-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)\right)=\mathrm{0},\end{array}\\ \begin{array}{l}-\frac{\partial }{\partial y_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\frac{\partial T^{\left(\mathrm{1}\right)}}{\partial y_j}\right)-\frac{\partial }{\partial x_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\frac{\partial T^{\left(\mathrm{0}\right)}}{\partial y_j}\right)-\\ \frac{\partial }{\partial y_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\frac{\partial T^{\left(\mathrm{0}\right)}}{\partial x_j}\right)=\mathrm{0}\end{array}\end{array}\right.(7)

将式（6）代入式（7），整理可得到关于u_k^{\left(\mathrm{1}\right)}\mathrm{、}{T}^{\left(\mathrm{1}\right)}的如下表达式：

\left\{\begin{array}{l}\begin{array}{l}{u}_k^{\left(\mathrm{1}\right)}\left(\mathit{x},\mathit{y}\right)=M_{pkm}\left(\mathit{x},\mathit{y}\right)\frac{\partial u_m^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}-\\ R_k\left(\mathit{x},\mathit{y}\right)\left(T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right),\end{array}\\ T^{\left(\mathrm{1}\right)}\left(\mathit{x},\mathit{y}\right)=N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}\end{array}\right.（8）

一阶单胞函数M_{pkm}\left(\mathit{x},\mathit{y}\right)\mathrm{、}{R}_k\left(\mathit{x},\mathit{y}\right)\mathrm{、}N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)分别满足如下问题：

\left\{\begin{array}{l}-\frac{\partial }{\partial y_j}\left(C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial M_{pkm}\left(\mathit{x},\mathit{y}\right)}{\partial y_l}\right)=\\ \frac{\partial C_{ijmp}\left(\mathit{x},\mathit{y}\right)}{\partial y_j},\mathrm{i}\mathrm{n}Y,\\ M_{pkm}\left(\mathit{x},\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（9）

\left\{\begin{array}{l}-\frac{\partial }{\partial y_j}\left(C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial R_k\left(\mathit{x},\mathit{y}\right)}{\partial y_l}\right)=\\ \frac{\partial {\beta }_{ij}\left(\mathit{x},\mathit{y}\right)}{\partial y_j},\mathrm{i}\mathrm{n}Y,\\ R_k\left(\mathit{x},\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（10）

\left\{\begin{array}{l}-\frac{\partial }{\partial y_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_j}\right)=\frac{\partial k_{ip}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_i},\mathrm{i}\mathrm{n}Y,\\ N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（11）

接下来，比较O\left({\epsilon }^{\mathrm{0}}\right)项：

\left\{\begin{array}{l}\begin{array}{l}-\frac{\partial }{\partial y_j}\left(C_{ijkl}\frac{\partial u_k^{\left(\mathrm{2}\right)}}{\partial y_l}\right)-\frac{\partial }{\partial x_j}\left(C_{ijkl}\frac{\partial u_k^{\left(\mathrm{1}\right)}}{\partial y_l}\right)-\frac{\partial }{\partial y_j}\left(C_{ijkl}\frac{\partial u_k^{\left(\mathrm{1}\right)}}{\partial x_l}\right)-\frac{\partial }{\partial x_j}\left(C_{ijkl}\frac{\partial u_k^{\left(\mathrm{0}\right)}}{\partial x_l}\right)+\frac{\partial }{\partial x_j}\left({\beta }_{ij}\left(T^{\left(\mathrm{0}\right)}-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)\right)+\\ \frac{\partial }{\partial y_j}\left({\beta }_{ij}{T}^{\left(\mathrm{1}\right)}\right)=f_k,\end{array}\\ \begin{array}{l}-\frac{\partial }{\partial y_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\frac{\partial T^{\left(\mathrm{2}\right)}}{\partial y_j}\right)-\frac{\partial }{\partial x_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\frac{\partial T^{\left(\mathrm{1}\right)}}{\partial y_j}\right)-\frac{\partial }{\partial y_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\frac{\partial T^{\left(\mathrm{1}\right)}}{\partial x_j}\right)-\frac{\partial }{\partial x_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\frac{\partial T^{\left(\mathrm{0}\right)}}{\partial x_j}\right)-\frac{\partial }{\partial x_i}\left(k_{ij}^{\left(\mathrm{1}\right)}\frac{\partial T^{\left(\mathrm{0}\right)}}{\partial y_j}\right)-\\ \frac{\partial }{\partial y_i}\left(k_{ij}^{\left(\mathrm{1}\right)}\frac{\partial T^{\left(\mathrm{0}\right)}}{\partial x_j}\right)=h\end{array}\end{array}\right.（12）

将（6）式和（8）式代入（12）式，并对（12）式在单胞Y上做积分平均，计算并整理得问题（1）对应的均匀化问题如下：

\left\{\begin{array}{l}-\frac{\partial }{\partial x_j}\left({\overline{C}}_{ijkl}\left(\mathit{x}\right)\frac{\partial u_k^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_l}-{\overline{\beta }}_{ij}\left(\mathit{x}\right)\left(T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)\right)=f_k\left(\mathit{x}\right),\mathrm{i}\mathrm{n}\mathrm{\Omega },\\ -\frac{\partial }{\partial x_i}\left({\overline{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},T^{\left(\mathrm{0}\right)}\right)\frac{\partial T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_j}\right)=h\left(\mathit{x}\right),\mathrm{i}\mathrm{n}\mathrm{\Omega },\\ {\mathit{u}}^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)=\mathrm{0},T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)=T_{\mathrm{r}\mathrm{e}\mathrm{f}},\mathrm{o}\mathrm{n}\partial \mathrm{\Omega }\end{array}\right.（13）

其中均匀化材料的参数定义如下：

\left\{\begin{array}{l}{\overline{C}}_{ijkl}\left(\mathit{x}\right)=\frac{\mathrm{1}}{\left|Y\right|}{\int }_Y{C}_{ijkl}\left(\mathit{x},\mathit{y}\right)+C_{ijpq}\left(\mathit{x},\mathit{y}\right)\frac{\partial M_{lpk}\left(\mathit{x},\mathit{y}\right)}{\partial y_q}\mathrm{d}\mathit{y},\\ {\overline{\beta }}_{ij}\left(\mathit{x}\right)=\frac{\mathrm{1}}{\left|Y\right|}{\int }_Y{\beta }_{ij}\left(\mathit{x},\mathit{y}\right)+C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial R_k\left(\mathit{x},\mathit{y}\right)}{\partial y_l}\mathrm{d}\mathit{y},\\ {\overline{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},T^{\left(\mathrm{0}\right)}\right)=\frac{\mathrm{1}}{\left|Y\right|}{\int }_Y{k}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)+k_{ik}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial N_j\left(\mathit{x},\mathit{y}\right)}{\partial y_k}\mathrm{d}\mathit{y}\end{array}\right.（14）

利用式（12）减去均匀化问题（13），并将（6）式和（8）式代入可得关于u_k^{\left(\mathrm{2}\right)}\mathrm{、}{T}^{\left(\mathrm{2}\right)}的如下表达式：

\left\{\begin{array}{l}\begin{array}{l}{u}_k^{\left(\mathrm{2}\right)}\left(\mathit{x},\mathit{y}\right)=M_{pqkm}\left(\mathit{x},\mathit{y}\right)\frac{{\partial }^{\mathrm{2}}{u}_m^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p\partial x_q}+F_{pkm}\left(\mathit{x},\mathit{y}\right)\frac{\partial u_m^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}-S_k\left(\mathit{x},\mathit{y}\right)\left(T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)-\\ G_{pk}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial }{\partial x_p}\left(T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right),\end{array}\\ T^{\left(\mathrm{2}\right)}\left(\mathit{x},\mathit{y}\right)=N_{pq}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{{\partial }^{\mathrm{2}}{T}^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p\partial x_q}+I_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}\end{array}\right.（15）

对于u_k^{\left(\mathrm{2}\right)}，单胞函数M_{pqkm}\left(\mathit{x},\mathit{y}\right)\mathrm{、}{F}_{pkm}\left(\mathit{x},\mathit{y}\right)\mathrm{、}{S}_k\left(\mathit{x},\mathit{y}\right)\mathrm{、}{G}_{pk}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)分别满足如下问题：

\left\{\begin{array}{l}\begin{array}{l}-\frac{\partial }{\partial y_j}\left(C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial M_{pqkm}\left(\mathit{x},\mathit{y}\right)}{\partial y_l}\right)=C_{iqkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial M_{pkm}\left(\mathit{x},\mathit{y}\right)}{\partial y_l}+\frac{\partial }{\partial y_j}\left(C_{ijkq}\left(\mathit{x},\mathit{y}\right)M_{pkm}\left(\mathit{x},\mathit{y}\right)\right)+\\ C_{iqmp}\left(\mathit{x},\mathit{y}\right)-{\overline{C}}_{iqmp}\left(\mathit{x}\right),\mathrm{i}\mathrm{n}Y,\end{array}\\ M_{pqkm}\left(\mathit{x},\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（16）

\left\{\begin{array}{l}\begin{array}{l}-\frac{\partial }{\partial y_j}\left(C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial F_{pkm}\left(\mathit{x},\mathit{y}\right)}{\partial y_l}\right)=\frac{\partial }{\partial x_j}\left(C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial M_{pkm}\left(\mathit{x},\mathit{y}\right)}{\partial y_l}\right)+\frac{\partial }{\partial y_j}\left(C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial M_{pkm}\left(\mathit{x},\mathit{y}\right)}{\partial x_l}\right)+\\ \frac{\partial C_{ijmp}\left(\mathit{x},\mathit{y}\right)}{\partial x_j}-\frac{\partial {\overline{C}}_{ijmp}\left(\mathit{x}\right)}{\partial x_j},\mathrm{i}\mathrm{n}Y,\end{array}\\ F_{pkm}\left(\mathit{x},\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（17）

\left\{\begin{array}{l}\begin{array}{l}-\frac{\partial }{\partial y_j}\left(C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial S_k\left(\mathit{x},\mathit{y}\right)}{\partial y_l}\right)=\frac{\partial }{\partial x_j}\left(C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial R_k\left(\mathit{x},\mathit{y}\right)}{\partial y_l}\right)+\frac{\partial }{\partial y_j}\left(C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial R_k\left(\mathit{x},\mathit{y}\right)}{\partial x_l}\right)+\\ \frac{\partial {\beta }_{ij}\left(\mathit{x},\mathit{y}\right)}{\partial x_j}-\frac{\partial {\overline{\beta }}_{ij}\left(\mathit{x}\right)}{\partial x_j},\mathrm{i}\mathrm{n}Y,\end{array}\\ S_k\left(\mathit{x},\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（18）

\left\{\begin{array}{l}\begin{array}{l}-\frac{\partial }{\partial y_j}\left(C_{ijkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial G_{pk}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_l}\right)=C_{ipkl}\left(\mathit{x},\mathit{y}\right)\frac{\partial R_k\left(\mathit{x},\mathit{y}\right)}{\partial y_l}+\frac{\partial }{\partial y_j}\left(C_{ijkp}\left(\mathit{x},\mathit{y}\right)R_k\left(\mathit{x},\mathit{y}\right)\right)+\\ \frac{\partial }{\partial y_j}\left({\beta }_{ij}\left(\mathit{x},\mathit{y}\right)N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\right)+{\beta }_{ip}\left(\mathit{x},\mathit{y}\right)-{\overline{\beta }}_{ip}\left(\mathit{x}\right),\mathrm{i}\mathrm{n}Y,\end{array}\\ G_{pk}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（19）

对于T^{\left(\mathrm{2}\right)}，单胞函数N_{pq}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\mathrm{、}{I}_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)分别满足如下问题：

\left\{\begin{array}{l}\begin{array}{l}-\frac{\partial }{\partial y_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial N_{pq}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_j}\right)=\frac{\partial }{\partial y_k}\left(k_{kq}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\right)+k_{pq}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)+\\ k_{qk}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_k}-{\overline{k}}_{pq}^{\left(\mathrm{0}\right)}\left(\mathit{x},T^{\left(\mathrm{0}\right)}\right),\mathrm{i}\mathrm{n}Y,\end{array}\\ N_{pq}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（20）

\left\{\begin{array}{l}\begin{array}{l}-\frac{\partial }{\partial y_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial I_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_j}\right)=\frac{\partial k_{ip}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial x_i}+\frac{\partial }{\partial y_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial x_j}\right)+\\ \frac{\partial k_{ip}^{\left(\mathrm{1}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_i}+\frac{\partial }{\partial x_i}\left(k_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_j}\right)-\frac{\partial {\overline{k}}_{ip}^{\left(\mathrm{0}\right)}\left(\mathit{x},T^{\left(\mathrm{0}\right)}\right)}{\partial x_i},\mathrm{i}\mathrm{n}Y,\end{array}\\ I_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（21）

所以，拟周期性复合材料结构非线性热-力耦合问题的位移场和温度场具有如下的二阶双尺度近似解：

\left\{\begin{array}{l}\begin{array}{c}\begin{array}{l}{u}_k^{\epsilon }\left(\mathit{x}\right)=u_k^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)+\epsilon \left[M_{pkm}\left(\mathit{x},\mathit{y}\right)\frac{\partial u_m^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}-R_k\left(\mathit{x},\mathit{y}\right)\left(T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)\right]+{\epsilon }^{\mathrm{2}}\left[M_{pqkm}\left(\mathit{x},\mathit{y}\right)\frac{{\partial }^{\mathrm{2}}{u}_m^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p\partial x_q}\right.+\\ F_{pkm}\left(\mathit{x},\mathit{y}\right)\frac{\partial u_m^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}-\left.S_k\left(\mathit{x},\mathit{y}\right)\left(T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)-T_{ref}\right)-G_{pk}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}\right]+O\left({\epsilon }^{\mathrm{3}}\right),\end{array}\end{array}\\ \begin{array}{c}\begin{array}{l}{T}^{\epsilon }\left(\mathit{x}\right)=T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)+\epsilon N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}+{\epsilon }^{\mathrm{2}}\left[N_{pq}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{{\partial }^{\mathrm{2}}{T}^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p\partial x_q}+I_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}\right]+\\ O\left({\epsilon }^{\mathrm{3}}\right)\end{array}\end{array}\end{array}\right.（22）

在分析拟周期材料时，宏观变量\mathit{x}与微观变量\mathit{y}之间存在耦合关系。因此，每个单胞问题的求解需将宏观坐标\mathit{x}作为参数输入。为了降低计算难度，我们假设拟周期复合材料结构问题具有如下的尺度分离系数：

\left\{\begin{array}{l}{C}_{ijkl}\left(\mathit{x},\mathit{y}\right)={\omega }_{\mathrm{1}}\left(\mathit{x}\right){\tilde{C}}_{ijkl}\left(\mathit{y}\right),{\beta }_{ij}\left(\mathit{x},\mathit{y}\right)={\omega }_{\mathrm{2}}\left(\mathit{x}\right){\tilde{\beta }}_{ij}\left(\mathit{y}\right),\\ \begin{array}{c}{k}_{ij}\left(\mathit{x},\mathit{y},T^{\epsilon }\right)={\omega }_{\mathrm{3}}\left(\mathit{x}\right){\tilde{k}}_{ij}\left(\mathit{y},T^{\epsilon }\right)={\omega }_{\mathrm{3}}\left(\mathit{x}\right)\left({\tilde{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)+\epsilon {\tilde{k}}_{ij}^{\left(\mathrm{1}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\right.\left.+{\epsilon }^{\mathrm{2}}{\tilde{k}}_{ij}^{\left(\mathrm{2}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)+O\left({\epsilon }^{\mathrm{3}}\right)\right)\end{array}\end{array}\right.（23）

将尺度分离系数代入一阶单胞函数方程（9）~（11），整理得到与宏观变量\mathit{x}无关的新的一阶单胞问题：

\begin{array}{l}\left\{\begin{array}{l}-\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkl}\left(\mathit{y}\right)\frac{\partial {\tilde{M}}_{pkm}\left(\mathit{y}\right)}{\partial y_l}\right)=\frac{\partial {\tilde{C}}_{ijmp}\left(\mathit{y}\right)}{\partial y_j},\mathrm{i}\mathrm{n}Y,\\ {\tilde{M}}_{pkm}\left(\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y,\end{array}\right.\\ \left\{\begin{array}{l}-\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkl}\left(\mathit{y}\right)\frac{\partial {\tilde{R}}_k\left(\mathit{y}\right)}{\partial y_l}\right)=\frac{\partial {\tilde{\beta }}_{ij}\left(\mathit{y}\right)}{\partial y_j},\mathrm{i}\mathrm{n}Y,\\ {\tilde{R}}_k\left(\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y,\end{array}\right.\\ \left\{\begin{array}{l}-\frac{\partial }{\partial y_i}\left({\tilde{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial {\tilde{N}}_p\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_j}\right)=\frac{\partial {\tilde{k}}_{ip}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_j},\mathrm{i}\mathrm{n}Y,\\ {\tilde{N}}_p\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.\end{array}（24）

新的均匀化材料系数定义为：

\left\{\begin{array}{c}{\overline{C}}_{ijkl}\left(\mathit{x}\right)=\frac{{\omega }_{\mathrm{1}}\left(\mathit{x}\right)}{\left|Y\right|}{\int }_Y{\tilde{C}}_{ijkl}\left(\mathit{y}\right)+{\tilde{C}}_{ijpq}\left(\mathit{y}\right)\frac{\partial {\tilde{M}}_{lpk}\left(\mathit{y}\right)}{\partial y_q}\mathrm{d}\mathit{y}={\omega }_{\mathrm{1}}\left(\mathit{x}\right){\overline{C}}_{ijkl}^{{}^{\mathrm{*}}},\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& \end{array}& \end{array}& \end{array}\\ {\overline{\beta }}_{ij}\left(\mathit{x}\right)=\frac{{\omega }_{\mathrm{2}}\left(\mathit{x}\right)}{\left|Y\right|}{\int }_Y{\tilde{\beta }}_{ij}\left(\mathit{y}\right)+{\tilde{C}}_{ijkl}\left(\mathit{y}\right)\frac{\partial {\tilde{R}}_k\left(\mathit{y}\right)}{\partial y_l}d\mathit{y}={\omega }_{\mathrm{2}}\left(\mathit{x}\right){\overline{\beta }}_{ij}^{\mathrm{*}},\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}\begin{array}{cc}& \end{array}& \end{array}& \end{array}& \end{array}& \end{array}& \end{array}\\ {\overline{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},T^{\left(\mathrm{0}\right)}\right)=\frac{{\omega }_{\mathrm{3}}\left(\mathit{x}\right)}{\left|Y\right|}{\int }_Y{\tilde{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)+{\tilde{k}}_{ik}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial {\tilde{N}}_j\left(\mathit{y}\right)}{\partial y_k}\mathrm{d}\mathit{y}={\omega }_{\mathrm{3}}\left(\mathit{x}\right){\overline{k}}^{\mathrm{*}}{}_{ij}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right)\end{array}\right.（25）

同样，将（23）式依次代入单胞问题（16）（19）和（20）可以得到新的二阶单胞函数{\tilde{M}}_{pqkm}\left(\mathit{y}\right)\mathrm{、}{\tilde{G}}_{pk}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\mathrm{、}{\tilde{N}}_{pq}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)依次满足单胞问题：

\left\{\begin{array}{l}-\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkl}\left(\mathit{y}\right)\frac{\partial {\tilde{M}}_{pqkm}\left(\mathit{y}\right)}{\partial y_l}\right)={\tilde{C}}_{iqkl}\left(\mathit{y}\right)\frac{\partial {\tilde{M}}_{pkm}\left(\mathit{y}\right)}{\partial y_l}+\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkq}\left(\mathit{y}\right){\tilde{M}}_{pkm}\left(\mathit{y}\right)\right)+{\tilde{C}}_{iqmp}\left(\mathit{y}\right)-{\overline{C}}_{iqmp}^{\mathrm{*}},\mathrm{i}\mathrm{n}Y,\\ {\tilde{M}}_{pqkm}\left(\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（26）

\left\{\begin{array}{l}\begin{array}{l}-\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkl}\left(\mathit{y}\right)\frac{\partial {\tilde{G}}_{pk}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_l}\right)={\tilde{C}}_{ipkl}\left(\mathit{y}\right)\frac{\partial {\tilde{R}}_k\left(\mathit{y}\right)}{\partial y_l}+\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkp}\left(\mathit{y}\right){\tilde{R}}_k\left(\mathit{y}\right)\right)+{\tilde{\beta }}_{ip}\left(\mathit{y}\right)+\\ \frac{\partial }{\partial y_j}\left({\tilde{\beta }}_{ij}\left(\mathit{y}\right){\tilde{N}}_p\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\right)-{\overline{\beta }}_{ip}^{\mathrm{*}},\mathrm{i}\mathrm{n}Y,\end{array}\\ {\tilde{G}}_{pk}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（27）

\left\{\begin{array}{l}\begin{array}{l}-\frac{\partial }{\partial y_i}\left({\tilde{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial {\tilde{N}}_{pq}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_j}\right)={\tilde{k}}_{pq}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right){\tilde{k}}_{qk}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial {\tilde{N}}_p\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_k}-{\overline{k}}^{\mathrm{*}}{}_{pq}^{\left(\mathrm{0}\right)}\left(T^{\left(\mathrm{0}\right)}\right)+\\ \frac{\partial }{\partial y_k}\left({\tilde{k}}_{kq}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right){\tilde{N}}_p\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\right),\mathrm{i}\mathrm{n}Y,\end{array}\\ {\tilde{N}}_{pq}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.(28)

注意到单胞问题（17）（18）和（21）在计算时涉及到\mathit{x}的导数，则将尺度分离系数（23）依次代入单胞问题后F_{pkm}\left(\mathit{x},\mathit{y}\right)\mathrm{、}{S}_k\left(\mathit{x},\mathit{y}\right)\mathrm{、}{I}_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)也能实现变量分离，从而得到新的单胞函数如下：

\left\{\begin{array}{l}-\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkl}\left(\mathit{y}\right)\frac{\partial {\tilde{F}}_{pqkm}\left(\mathit{y}\right)}{\partial y_l}\right)={\tilde{C}}_{iqkl}\left(\mathit{y}\right)\frac{\partial M_{pkm}\left(\mathit{y}\right)}{\partial y_l}+C_{iqmp}\left(\mathit{y}\right)-{\overline{C}}_{iqmp}^{\mathrm{*}},\mathrm{i}\mathrm{n}Y,\\ {\tilde{F}}_{pqkm}\left(\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（29）

\left\{\begin{array}{l}\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkl}\left(\mathit{y}\right)\frac{\partial {\tilde{S}}_{pk}^{\mathrm{1}}\left(\mathit{y}\right)}{\partial y_l}\right)=\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkp}\left(\mathit{y}\right){\tilde{R}}_k\left(\mathit{y}\right)\right),\mathrm{i}\mathrm{n}Y,\\ {\tilde{S}}_{pk}^{\mathrm{1}}\left(\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（30）

\left\{\begin{array}{l}-\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkl}\left(\mathit{y}\right)\frac{\partial {\tilde{S}}_{pk}^{\mathrm{2}}\left(\mathit{y}\right)}{\partial y_l}\right)={\tilde{C}}_{ipkl}\left(\mathit{y}\right)\frac{\partial {\tilde{R}}_k\left(\mathit{y}\right)}{\partial y_l}+\frac{\partial }{\partial y_j}\left({\tilde{C}}_{ijkp}\left(\mathit{y}\right){\tilde{R}}_k\left(\mathit{y}\right)\right)+{\tilde{\beta }}_{ip}\left(\mathit{y}\right)-{\overline{\beta }}_{ip}^{\mathrm{*}},\mathrm{i}\mathrm{n}Y,\\ {\tilde{S}}_{pk}^{\mathrm{2}}\left(\mathit{y}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（31）

\left\{\begin{array}{l}-\frac{\partial }{\partial y_i}\left({\tilde{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial {\tilde{I}}_{pq}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_j}\right)={\tilde{k}}_{qp}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)+{\tilde{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial {\tilde{N}}_p\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_j}-{\overline{k}}_{qp}^{\mathrm{*}(\mathrm{0})}\left(T^{\left(\mathrm{0}\right)}\right),\mathrm{i}\mathrm{n}Y,\\ {\tilde{I}}_{pq}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（32）

\left\{\begin{array}{l}-\frac{\partial }{\partial y_i}\left({\tilde{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial {\tilde{I}}_p\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_j}\right)=\frac{\partial {\tilde{k}}_{ip}^{\left(\mathrm{1}\right)}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)}{\partial y_i},\mathrm{i}\mathrm{n}Y,\\ {\tilde{I}}_p\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)=\mathrm{0},\mathrm{o}\mathrm{n}\partial Y\end{array}\right.（33）

最终，本文得到位移场和温度场的二阶双尺度近似解为：

\left\{\begin{array}{l}\begin{array}{c}\begin{array}{l}{u}_k^{\epsilon }\left(\mathit{x}\right)\approx u_k^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)+\epsilon \left[{\tilde{M}}_{pkm}\left(\mathit{y}\right)\frac{\partial u_m^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}-\frac{{\omega }_{\mathrm{2}}\left(\mathit{x}\right)}{{\omega }_{\mathrm{1}}\left(\mathit{x}\right)}{\tilde{R}}_k\left(\mathit{y}\right)\left(T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)\right]+{\epsilon }^{\mathrm{2}}\left[{\tilde{M}}_{pqkm}\left(\mathit{y}\right)\frac{{\partial }^{\mathrm{2}}{u}_m^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p\partial x_q}\right.-\\ \frac{{\omega }_{\mathrm{2}}\left(\mathit{x}\right)}{{\omega }_{\mathrm{1}}\left(\mathit{x}\right)}{\tilde{G}}_{pk}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}+\frac{\mathrm{1}}{{\omega }_{\mathrm{1}}\left(\mathit{x}\right)}\frac{\partial {\omega }_{\mathrm{1}}\left(\mathit{x}\right)}{\partial x_q}{\tilde{F}}_{pqkm}\left(\mathit{y}\right)\frac{\partial u_m^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}-\left(\frac{{\omega }_{\mathrm{2}}\left(\mathit{x}\right)}{{\left({\omega }_{\mathrm{1}}\left(\mathit{x}\right)\right)}^{\mathrm{2}}}\frac{\partial {\omega }_{\mathrm{1}}\left(\mathit{x}\right)}{\partial x_p}{\tilde{S}}_{pk}^{\mathrm{1}}\left(\mathit{y}\right)\right.+\\ \left.\frac{\mathrm{1}}{{\omega }_{\mathrm{1}}\left(\mathit{x}\right)}\frac{\partial {\omega }_{\mathrm{2}}\left(\mathit{x}\right)}{\partial x_p}{\tilde{S}}_{pk}^{\mathrm{2}}\left(\mathit{y}\right)\right)\left.\left(T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)-T_{\mathrm{r}\mathrm{e}\mathrm{f}}\right)\right],\end{array}\end{array}\\ \begin{array}{c}\begin{array}{l}{T}^{\epsilon }\left(\mathit{x}\right)\approx T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)+\epsilon {\tilde{N}}_p\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{\partial T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}+{\epsilon }^{\mathrm{2}}\left[{\tilde{N}}_{pq}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\frac{{\partial }^{\mathrm{2}}{T}^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p\partial x_q}+\left(\frac{\mathrm{1}}{{\omega }_{\mathrm{3}}\left(\mathit{x}\right)}\frac{\partial {\omega }_{\mathrm{3}}\left(\mathit{x}\right)}{\partial x_q}{\tilde{I}}_{pq}\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\right.\right.+\\ \left.\left.{\tilde{I}}_p\left(\mathit{y},T^{\left(\mathrm{0}\right)}\right)\right)\frac{\partial T^{\left(\mathrm{0}\right)}\left(\mathit{x}\right)}{\partial x_p}\right]\end{array}\end{array}\end{array}\right.（34）

3 有限元算法

为方便分析，分别定义温度场和位移场的均匀化解、一阶双尺度近似解和二阶双尺度近似解为：

\begin{array}{l}{T}^{\left(\mathrm{0}\epsilon \right)}=T^{\left(\mathrm{0}\right)},T^{\left(\mathrm{1}\epsilon \right)}=T^{\left(\mathrm{0}\epsilon \right)}+\epsilon T^{\left(\mathrm{1}\right)},T^{\left(\mathrm{2}\epsilon \right)}=\\ T^{\left(\mathrm{1}\epsilon \right)}+{\epsilon }^{\mathrm{2}}{T}^{\left(\mathrm{2}\right)},\end{array}

\begin{array}{l}{\mathit{u}}^{\left(\mathrm{0}\epsilon \right)}={\mathit{u}}^{\left(\mathrm{0}\right)},{\mathit{u}}^{\left(\mathrm{1}\epsilon \right)}={\mathit{u}}^{\left(\mathrm{0}\epsilon \right)}+\epsilon {\mathit{u}}^{\left(\mathrm{1}\right)},{\mathit{u}}^{\left(\mathrm{2}\epsilon \right)}=\\ {\mathit{u}}^{\left(\mathrm{1}\epsilon \right)}+{\epsilon }^{\mathrm{2}}{\mathit{u}}^{\left(\mathrm{2}\right)}\mathrm{。}\end{array}

用{‖•‖}_{L^{\mathrm{2}}}和{\left|•\right|}_{H^{\mathrm{1}}}分别表示L^{\mathrm{2}}范数和H^{\mathrm{1}}半范数，则温度场和位移场的均匀化解、一阶双尺度解和二阶双尺度解的相对误差的定义为：

\begin{array}{l}{e}_T^{\mathrm{0}}=\frac{{‖T^{\left(\epsilon \right)}-T^{\left(\mathrm{0}\epsilon \right)}‖}_{L^{\mathrm{2}}}}{{‖T^{\left(\epsilon \right)}‖}_{L^{\mathrm{2}}}},e_T^{\mathrm{1}}=\frac{{‖T^{\left(\epsilon \right)}-T^{\left(\mathrm{1}\epsilon \right)}‖}_{L^{\mathrm{2}}}}{{‖T^{\left(\epsilon \right)}‖}_{L^{\mathrm{2}}}},\\ e_T^{\mathrm{2}}=\frac{{‖T^{\left(\epsilon \right)}-T^{\left(\mathrm{2}\epsilon \right)}‖}_{L^{\mathrm{2}}}}{{‖T^{\left(\epsilon \right)}‖}_{L^{\mathrm{2}}}},\end{array}

\begin{array}{l}{E}_T^{\mathrm{0}}=\frac{{\left|T^{\left(\epsilon \right)}-T^{\left(\mathrm{0}\epsilon \right)}\right|}_{H^{\mathrm{1}}}}{{\left|T^{\left(\epsilon \right)}\right|}_{H^{\mathrm{1}}}},E_T^{\mathrm{1}}=\frac{{\left|T^{\left(\epsilon \right)}-T^{\left(\mathrm{1}\epsilon \right)}\right|}_{H^{\mathrm{1}}}}{{\left|T^{\left(\epsilon \right)}\right|}_{H^{\mathrm{1}}}},\\ E_T^{\mathrm{2}}=\frac{{\left|T^{\left(\epsilon \right)}-T^{\left(\mathrm{2}\epsilon \right)}\right|}_{H^{\mathrm{1}}}}{{\left|T^{\left(\epsilon \right)}\right|}_{H^{\mathrm{1}}}},\end{array}

\begin{array}{l}{e}_{\mathit{u}}^{\mathrm{0}}=\frac{{‖{\mathit{u}}^{\left(\epsilon \right)}-{\mathit{u}}^{\left(\mathrm{0}\epsilon \right)}‖}_{L^{\mathrm{2}}}}{{‖{\mathit{u}}^{\left(\epsilon \right)}‖}_{L^{\mathrm{2}}}},e_{\mathit{u}}^{\mathrm{1}}=\frac{{‖{\mathit{u}}^{\left(\epsilon \right)}-{\mathit{u}}^{\left(\mathrm{1}\epsilon \right)}‖}_{L^{\mathrm{2}}}}{{‖{\mathit{u}}^{\left(\epsilon \right)}‖}_{L^{\mathrm{2}}}},\\ e_{\mathit{u}}^{\mathrm{2}}=\frac{{‖{\mathit{u}}^{\left(\epsilon \right)}-{\mathit{u}}^{\left(\mathrm{2}\epsilon \right)}‖}_{L^{\mathrm{2}}}}{{‖{\mathit{u}}^{\left(\epsilon \right)}‖}_{L^{\mathrm{2}}}},\end{array}

\begin{array}{l}{E}_{\mathit{u}}^{\mathrm{0}}=\frac{{\left|{\mathit{u}}^{\left(\epsilon \right)}-{\mathit{u}}^{\left(\mathrm{0}\epsilon \right)}\right|}_{H^{\mathrm{1}}}}{{\left|{\mathit{u}}^{\left(\epsilon \right)}\right|}_{H^{\mathrm{1}}}},E_{\mathit{u}}^{\mathrm{1}}=\frac{{\left|{\mathit{u}}^{\left(\epsilon \right)}-{\mathit{u}}^{\left(\mathrm{1}\epsilon \right)}\right|}_{H^{\mathrm{1}}}}{{\left|{\mathit{u}}^{\left(\epsilon \right)}\right|}_{H^{\mathrm{1}}}},\\ E_{\mathit{u}}^{\mathrm{2}}=\frac{{\left|{\mathit{u}}^{\left(\epsilon \right)}-{\mathit{u}}^{\left(\mathrm{2}\epsilon \right)}\right|}_{H^{\mathrm{1}}}}{{\left|{\mathit{u}}^{\left(\epsilon \right)}\right|}_{H^{\mathrm{1}}}}\mathrm{。}\end{array}

基于对多尺度问题（1）的二阶双尺度分析，我们将给出求解拟周期性复合材料结构稳态热-力耦合问题的二阶双尺度算法流程。

预处理旨在确定参考单胞Y={\left[\mathrm{0,1}\right]}^N和均匀化区域\mathrm{\Omega }\subset {\mathbf{ℝ}}^N,N=\mathrm{2}\mathrm{、}\mathrm{3}，以及复合材料中各种材料的材料参数C_{ijkl}\left(\mathit{x},\mathit{y}\right)\mathrm{、}{\beta }_{ij}\left(\mathit{x},\mathit{y}\right)\mathrm{、}{k}_{ij}\left(\mathit{x},\mathit{y},T^{\epsilon }\right)\mathrm{。}

3.1 计算单胞函数

1） 设定计算温度区间\left[T_{\mathrm{m}\mathrm{i}\mathrm{n}}^{\left(\mathrm{0}\right)},T_{\mathrm{m}\mathrm{a}\mathrm{x}}^{\left(\mathrm{0}\right)}\right]，并在该区间内选取一定数量的代表点T_s^{\left(\mathrm{0}\right)},s=\mathrm{1}\mathrm{、}\mathrm{2}\mathrm{、}\cdots \mathrm{、}l\left(l\in {\mathbf{Z}}^{+}\right)，根据单胞问题（9）~（13）利用有限元方法求解代表点上一阶单胞函数M_{pkm}\left(\mathit{x},\mathit{y}\right)\mathrm{、}{R}_k\left(\mathit{x},\mathit{y}\right)\mathrm{、}N_p\left(\mathit{x},\mathit{y},T_s^{\left(\mathrm{0}\right)}\right)。

2） 利用式（14）分别求解得到均匀化材料参数{\overline{C}}_{ijkl}\left(\mathit{x}\right)\mathrm{、}{\overline{\beta }}_{ij}\left(\mathit{x}\right)\mathrm{、}{\overline{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},T_s^{\left(\mathrm{0}\right)}\right),然后通过插值方法计算具有宏观参数T^{\left(\mathrm{0}\right)}的均匀化参数{\overline{k}}_{ij}^{\left(\mathrm{0}\right)}\left(\mathit{x},T^{\left(\mathrm{0}\right)}\right)\mathrm{。}

3） 结合一阶单胞函数和均匀化参数，根据单胞问题（16）~（21）得到二阶单胞函数\begin{array}{l}{M}_{pqkm}\left(\mathit{x},\mathit{y}\right)\mathrm{、}{F}_{pkm}\left(\mathit{x},\mathit{y}\right)\mathrm{、}{S}_k\left(\mathit{x},\mathit{y}\right)\mathrm{、}{G}_{pk}\left(\mathit{x},\mathit{y},T_s^{\left(\mathrm{0}\right)}\right)\mathrm{、}\\ N_{pq}\left(\mathit{x},\mathit{y},T_s^{\left(\mathrm{0}\right)}\right)\mathrm{、}{I}_p\left(\mathit{x},\mathit{y},T_s^{\left(\mathrm{0}\right)}\right)\mathrm{。}\end{array}

4） 使用插值方法得到单胞函数N_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\mathrm{、}{G}_{pk}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)\mathrm{、}{N}_{pq}\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)及I_p\left(\mathit{x},\mathit{y},T^{\left(\mathrm{0}\right)}\right)。

3.2 计算二阶双尺度解

1） 基于均匀化方程（13），使用直接迭代法求解温度场和位移场的均匀化解u_k^{\left(\mathrm{0}\epsilon \right)}和T^{\left(\mathrm{0}\epsilon \right)}。

2） 组装一阶双尺度近似解u_k^{\left(\mathrm{1}\epsilon \right)}和T^{\left(\mathrm{1}\epsilon \right)}以及二阶双尺度近似解u_k^{\left(\mathrm{2}\epsilon \right)}和T^{\left(\mathrm{2}\epsilon \right)}。

4 数值算例

考虑如图1所示的宏观区域\mathrm{\Omega }=\left(x_{\mathrm{1}},x_{\mathrm{2}}\right)={\left[\mathrm{0,1}\right]}^{\mathrm{2}}和参考单胞Y=\left(y_{\mathrm{1}},y_{\mathrm{2}}\right)={\left[\mathrm{0,1}\right]}^{\mathrm{2}}，其中在图1b中单位正方形的单胞区域Y内有一个以（0.5， 0.5）为中心，半径为0.3的圆形夹杂，此处\epsilon =\mathrm{1}/\mathrm{8}。

假设问题的材料参数具有尺度分离性质，拟周期系数如下：

{\omega }_{\mathrm{1}}\left(\mathit{x}\right)=\mathrm{5}+\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{4}\mathrm{\pi }{x}_{\mathrm{1}}\right)+\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{4}\mathrm{\pi }{x}_{\mathrm{2}}\right),

\begin{array}{l}{\omega }_{\mathrm{2}}\left(\mathit{x}\right)=\mathrm{10}+x_{\mathrm{1}}\left(\mathrm{1}-x_{\mathrm{1}}\right)+x_{\mathrm{2}}\left(\mathrm{1}-x_{\mathrm{2}}\right),\\ {\omega }_{\mathrm{3}}\left(\mathit{x}\right)={\mathrm{e}}^{x_{\mathrm{1}}{x}_{\mathrm{2}}}+x_{\mathrm{1}}+x_{\mathrm{2}},\end{array}

复合区域由两种材料构成，具体材料参数如表1所示。其中， {\delta }_{ij}是Kronecker delta函数。

在平面应力问题中，四阶弹性张量为：

{\tilde{C}}_{ijkl}=\frac{E}{\mathrm{1}-{\nu }^{\mathrm{2}}}\left[\frac{\mathrm{1}}{\mathrm{2}}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right)+\frac{\nu }{\mathrm{1}-\nu }{\delta }_{ij}{\delta }_{kl}\right]\mathrm{。}

该问题的体积力、热源及参考温度分别设定为：

\begin{array}{l}{f}_{\mathrm{1}}\left(\mathit{x}\right)=\mathrm{0},f_{\mathrm{2}}\left(\mathit{x}\right)=\mathrm{0},h\left(\mathit{x}\right)=\mathrm{1000}\mathrm{J}/\left(\mathrm{c}\mathrm{m}·\mathrm{s}\right),\\ T_{\mathrm{r}\mathrm{e}\mathrm{f}}=\mathrm{300}\mathrm{K}\mathrm{。}\end{array}

由于无法解析求解问题（1），我们采用细网格上的传统有限元解u_k^{\left(\epsilon \right)}和T^{\left(\epsilon \right)}作为参考解与双尺度计算的各阶近似解进行比较。表2给出了对宏观区域、单胞区域及均匀化区域进行三角剖分后的单元和节点的信息。对比其中的网格数据可知，二阶双尺度方法较传统有限元法显著降低了网格数量需求。此外，为了得到可靠的数值解，传统有限元方法求解时宏观区域的单元和网格个数需要随着周期的增加而增加，而使用二阶双尺度方法网格规模不受周期的影响，从而为大规模工程计算提供了更经济的解决方案。

图2给出在x_{\mathrm{1}}=\mathrm{0.55}处的温度场和x_{\mathrm{2}}方向位移场的计算结果。从图2a可以明显看到，对于温度场来说，均匀化结果十分光滑，能够大致描绘出解的轮廓。进而，引入一阶修正后，解产生振荡，但与参考解相比仍有很大差别，只有二阶双尺度解能够准确捕捉到微观振荡。图2b为位移场的截面图。从图中可以看到，一阶渐近解能够有效逼近参考解，而二阶渐近解能够在一阶修正基础上进一步修正得到更好的逼近效果。

表3和表4分别为温度场和位移场的L^{\mathrm{2}}范数和H^{\mathrm{1}}半范数相对误差。可以看到，无论采用L^{\mathrm{2}}范数还是H^{\mathrm{1}}半范数作为评价标准，二阶双尺度解在位移场和温度场计算中的精度都明显优于均匀化解和一阶双尺度解，特别是温度场的精度提升尤为显著。此外，对于位移场，虽然E_{\mathit{u}}^{\mathrm{0}}随周期增加有所增加，但经加入一阶和二阶修正后，周期数更多则相对误差减小幅度更大、修正效果更好，说明引入一阶和二阶修正项是有必要的.图3为拟周期情况下各渐近解在不同ε下的L²范数相对误差。

5 结论

为求解拟周期复合材料结构的非线性稳态热-力耦合问题，本文发展了二阶双尺度方法，通过渐近展开方法推导了问题的二阶双尺度近似解，并建立了求解此类问题的二阶双尺度有限元数值算法。数值算例表明，相较传统的均匀化方法和一阶双尺度渐近展开法，二阶双尺度方法在捕捉拟周期复合材料微观尺度局部振荡特性方面有显著优势，说明在多尺度渐近展开中引入二阶修正项十分必要。

本研究为核反应堆关键部件（如燃料包壳和慢化剂材料）的多尺度力学建模提供了新工具。从计算效率看，采用有限元法直接求解拟周期复合材料问题需进行非常精细的网格划分，而二阶双尺度算法仅需极低计算成本，这对实现核反应堆结构完整性分析和高效高精度仿真十分重要。

参考文献

[1] Bensoussan A， Lions J L， Papanicolaou G.Asymptotic analysis for periodic structures ［M］.Providence： AMS， 2011.

[2] Oleïnik O A， Shamaev A S， Yosifian G A.Mathematical problems in elasticity and homogenization ［M］.Berlin： Elsevier， 1992.

[3] Cioranescu D， Donato P.An introduction to homogenization ［M］.Oxford： Oxford University Press， 1999.

[4] Lions J L.Some methods in the mathematical analysis of systems and their control ［M］.Beijing： Science Press， 1981.

[5] Jikov V V， Kozlov S M， Oleinik O A.Homogenization of differential operators and integral functionals ［M］.Berlin： Springer-Verlag， 1994.

[6] Cao L Q， Cui J Z， Zhu D C， et al.Spectral analysis and numerical simulation for second order elliptic operator with highly oscillating coefficients in perforated domains with a periodic structure ［J］.Sci China Ser A， 2002， 45（12）： 1588-1602.

[7] Marchenko V A， Khruslov E Y.Boundary-value problems with fine-grained boundary ［J］.Math Sb， 1964， 107（3）： 458-472.

[8] Nguetseng G.A general convergence result for a functional related to the theory of homogenization ［J］.SIAM J Math Anal， 1989， 20（3）： 608-623.

[9] Allaire G.Homogenization and two-scale convergence ［J］.SIAM J Math Anal， 1992， 23（6）： 1482-1518.

[10] Ye S Y， Ma Q， Hu B.Second order two scale calculation for eigenvalue problem of porous material model in curvilinear coordinates ［J］.J Sichuan Univ （Nat Sci Ed）， 2021， 58（3）： 031004.

[11] 叶舒愉， 马强， 胡兵.曲线坐标下多孔材料模型特征值问题的二阶双尺度分析与计算［J］.四川大学学报（自然科学版）， 2021， 58（3）： 031004.

[12] Tan L Q， Ma Q， Hu B.Multiscale asymptotic analysis for Steklov elastic eigenvalue problem in periodically perforated domain ［J］.J Sichuan Univ （Nat Sci Ed）， 2022， 59（1）： 011002.

[13] 谭理琴， 马强， 胡兵.周期多孔结构的Steklov弹性特征值问题的多尺度渐近分析［J］.四川大学学报（自然科学版）， 2022， 59（1）： 011002.

[14] Yang Z， Cui J， Wu Y， et al.Second-order two-scale analysis method for dynamic thermo-mechanical problems in periodic structure ［J］.Int J Numer Anal Model， 2015， 12（1）： 190-211.

[15] Wang X， Cao L， Wong Y.Multiscale computation and convergence for coupled thermoelastic system in composite materials ［J］.Multiscale Model Simul， 2015， 13（2）： 661-690.

[16] Dong H， Nie Y， Cui J， et al.Second-order two-scale analysis method for dynamic thermo-mechanical problems of composite structures with cylindrical periodicity ［J］.Int J Numer Anal Model， 2018， 15（5）： 834-863.

[17] Cao L Q， Cui J Z.Homogenization method for quasi-periodic structures of composite materials ［J］.J Comput Math， 1999， 21（3）： 331-344.

[18] 曹礼群， 崔俊芝.复合材料拟周期结构的均匀化方法［J］.计算数学， 1999， 21（3）： 331-344.

[19] Fang S， Zhi C J， Zhan X， et al.A second-order and two-scale computation method for the quasi-periodic structures of composite materials ［J］.Finite Elem Anal Des， 2010， 46（4）： 320-327.

[20] Su F， Xu Z， Cui J Z， et al.Multi-scale method for the quasi-periodic structures of composite materials ［J］.Appl Math Comput， 2011， 217（12）： 5847-5852.

[21] Dong H， Cui J， Nie Y， et al.Second-order two-scale computational method for damped dynamic thermo-mechanical problems of quasi-periodic composite materials ［J］.J Comput Appl Math， 2018， 343（15）： 575-601.

[22] Chen L M， Dai Z， Gu Y， et al.Integrated optimization design of lightweight thermal protection structures ［J］.Chin J Theor Appl Mech， 2011， 43（2）： 289-295.

[23] 陈立明， 戴政， 谷宇，等.轻质热防护结构的一体化优化设计研究［J］.力学学报， 2011， 43（2）： 289-295.

[24] Matsumoto T， Guzik A， Tanaka M.A boundary element method for analysis of thermoelastic deformations in materials with temperature dependent properties ［J］.Int J Numer Methods Eng， 2005， 64（11）： 1432-1458.

[25] Goupee A J， Vel S S.Transient multiscale thermoelastic analysis of functionally graded materials ［J］.Compos Struct， 2010， 92（6）： 1372-1390.

[26] Zhang H W， Yang D S， Zhang S， et al.Multiscale nonlinear thermoelastic analysis of heterogeneous multiphase materials with temperature-dependent properties ［J］.Finite Elem Anal Des， 2014， 88（9）： 97-117.

[27] Dong H， Yang Z， Nie Y.Higher-order multiscale method and its convergence analysis for nonlinear thermo-electric coupling problems of composite structures ［PP/OL］.V1.arxiv（2025-01-23）［2025-12-30］.

基金资助

国家重点研发计划(2024YFA1012803)

四川省自然科学基金(2024NSFSC0438)