热力耦合问题数学均匀化方法的物理意义

朱晓鹏， 黄俊， 陈磊,*， 邢誉峰

(1. 安徽华电工程咨询设计有限公司， 合肥 230022； 2. 北京航空航天大学航空科学与工程学院， 北京 100083；3. 北京航空航天大学 合肥创新研究院， 合肥 230012)

复合材料具有比强度高、比刚度大等优点,广泛应用于航天、航空工业领域。众所周知，对于很多复合材料的宏观解，如低阶频率和模态，可以使用等应变模型或等应力模型[1]及其他均匀化方法[2]求解，但相对于宏观应力分析，细观结构分析要复杂很多。为了在计算精度和效率之间达到平衡，各种多尺度方法相继被提出，如数学均匀化方法(MHM)[3-4]、广义有限单元方法(GFEM)[5-6]、多尺度有限单元方法(MsFEM)[7-8]、异类多尺度方法(HMM)[9-10]及多尺度特征单元方法(MEM)[11-12]。其中，MHM是最具代表性的多尺度方法之一，因其具有严格的数学逻辑，且可以包含复合材料所有的细观结构信息,广泛应用于分析热力耦合作用下的周期复合材料结构响应问题[13-18]。

在应用MHM处理热力耦合作用下的复合材料结构响应问题时，摄动阶次的选择是不可回避的问题，然而多数文献只考虑了一阶摄动项，这对于细观结构物理和力学信息的捕捉往往是不够的[19]，在很多实际工程问题中[20-25]，二阶摄动项对于计算精度的影响不可忽略，这就需要寻找高阶MHM求解复合材料细观结构的物理和力学属性。Han等[26-27]基于二阶均匀化方法处理功能梯度材料的静态热力耦合问题，准确估计了有效弹性模量和应力应变场，并基于周期均匀化方法推导了用于求解复合材料弹性问题的二阶摄动项格式。Guan[20]和Cui[22]等给出了一种二阶多尺度方法用于预测周期复合材料的物理和力学属性，并应用于实际工程问题中。Allaire和Habibi[23]使用双尺度二阶渐进展开方法研究了具有辐射边界条件下的热传导模型，并证明了均匀化过程中的双尺度收敛。通过二阶摄动项可以更加准确地捕捉到材料内部物理和力学行为的细观信息，这一结论得到越来越多学者的认可，然而关于二阶摄动项的必要性问题基本上都是通过具体物理和力学问题的计算精度进行论证和说明，其结论往往不具有普遍性。

MHM本身是具有严密逻辑的数学方法，在处理力学领域具体工程问题时，如何从力学角度解释MHM是应用数学学者经常关心的问题，一旦能够研究清楚各阶摄动项的物理意义，无疑将有助于从力学角度确定合适的展开项数或阶次，避免一味追求高阶次或盲目舍弃高于某阶的项。Xing和Chen[28]在分析比较各种多尺度方法基础之上，针对线弹性平面问题揭示了MHM前三阶摄动项的物理意义，为该方法在数学和力学之间建立了联系，并通过物理意义的分析得到了二阶摄动项与一阶摄动项同为摄动主项的结论。除此之外，鲜有其他关于MHM物理意义研究的报道。相对于线弹性问题MHM物理意义的研究，热力耦合问题MHM各阶摄动项中同时包含弹性影响函数和热影响函数，以致出现更多、更复杂的虚拟载荷形式，其物理解释也变得复杂，然而确定热力耦合问题MHM的物理意义将有助于更好地理解和应用该方法处理具体的力学和物理问题，并且从物理机理的角度研究合适的展开阶次所得到的结论相对于通过具体力学和物理问题的研究所得到的结论更具有普遍性。

在此背景下，本文首先通过构造MHM各阶摄动项的全解耦格式，推导热力耦合问题MHM高阶摄动项的数学表达式及矩阵列式；其次将弹性影响函数和热影响函数分别比拟为弹性虚拟位移和热虚拟位移，利用热力耦合问题MHM矩阵列式实现方法、量纲分析方法和虚拟载荷自平衡性质，分别给出热特征位移和弹性特征位移的空间矢量图，再借助构造有限元位移场和温度场中广义坐标的概念，给出了各阶弹性虚拟位移、热虚拟位移和摄动项的物理意义，并通过物理意义的分析得到了二阶摄动项的必要性依据；最后通过数值算例和最小势能泛函计算结果验证了热力耦合问题MHM有限元矩阵列式推导和物理意义分析的正确性。

1 热力耦合问题MHM的数学表达式和有限元矩阵列式

基于微观结构的周期性和单胞域的一致性假设，均匀化理论可以解耦成为微观结构异质边界值的单胞问题和宏观问题，对于结构线性小变形问题，热力耦合的应力-应变关系为

(1)

总应变-位移和热应变分别表示为

(2)

(3)

式中：uε为真实位移；amn为热膨胀系数；T为结构温差。

将式(2)和式(3)代入式(1)，并考虑弹性常数张量的对称性可以得到

(4)

二维周期复合材料控制方程为一致椭圆偏微分方程，其具有Dirichlet边界条件[29]：

(5)

(6)

表1 摄动位移及影响函数控制方程Table 1 Perturbation displacements and governing equation of influence functions

3) 三阶及以上阶次影响函数控制方程具有相同的递推形式。

热力耦合问题MHM的应用分3步实施：①求解单胞问题，得到各阶热影响函数、弹性影响函数及均匀化材料参数EH和aH；②求解各阶宏观位移场导数和温度场导数；③求解MHM真实位移和应力。

表1中各阶弹性影响函数和热影响函数控制方程的矩阵列式可以归纳为

(7)

(8)

(9)

(10)

(11)

式中：

(12)

其中：s=3,4,5…；e为单胞内的子单元数量；Ye表示单胞内子单元的区域；式(8)和式(9)中Eε表示周期结构内不同材料的弹性常数矩阵，式(10)和式(11)中的弹性常数矩阵在文献[28]附录A中给出；B为应变矩阵；χ为弹性影响函数矩阵；ψ为热影响函数矩阵；N为形函数矩阵；B、χ、ψ、N的维度取决于使用有限元(FEM)计算时所选择的单元阶次。

以线性矩阵单元为例，N为2×8矩阵，形式如下：

(13)

式中：Ni(i=1,2,3,4)为单胞内子单元4个节点的形函数。此时B为3×8矩阵，与标准有限元方法形式相同；一阶影响函数χ1和ψ1分别为8×3矩阵和8×1矩阵，二阶影响函数χ2和ψ2分别为8×6矩阵和8×2矩阵，三阶影响函数χ3和ψ3分别为8×12矩阵和8×4矩阵，更高阶次影响函数矩阵的维度以此类推。

(14)

结构由于热膨胀作用只产生线应变，剪应变为0，所以这种由于热变形产生的应变可以看作结构的初始应变ε0。对于二维问题，ε0的表达式为

(15)

所以均匀化宏观问题可以表达为

(16)

式中：f为外载荷。

在实际工程应用中，一般将一个单胞或多个单胞作为一个宏观单元处理，提高计算效率。

在得到宏观位移及其各阶导数、各阶温差场导数及摄动位移后，结构真实细观位移可以使用如下渐进展开表达式得到：

(17)

式中：∂u0/∂x、∂2u0/∂x2及∂3u0/∂x3分别与χ1、χ2和χ3的上标kl、klp和klpq的组合相对应求得；∂T/∂x和∂2T/∂x2分别与ψ2、ψ3的上标klp和klpq的组合相对应求得。

从以上求解公式可以看出，热力耦合问题MHM解耦形式都可以使用有限元方法得到。由式(17)可以看出，MHM的显著优点为细观尺度的解可以用宏观尺度来表示，这意味着在得到结构单胞的影响函数后，其计算精度取决于宏观位移、宏观位移场导数和宏观温度场导数的求解精度。

2 热力耦合问题MHM的物理解释

2.1 二维周期复合材料单胞结构

对于单胞模型内部一个线性子单元，用于求解一阶弹性影响函数χ1的弹性自平衡虚拟载荷形式为

(18)

式中：ν为泊松比；l1和l2分别为矩形子单元沿着坐标轴x1和x2方向的长度。显然，如果l1=l2或c=1，即正方形子单元，F1的第1列(与kl=11相对应)中每个节点沿着坐标轴x1和x2方向的比值为±1/ν，第2列(与kl=22相对应)中每个节点沿着坐标x1和x2方向的比值为±ν，第3列(与kl=12相对应)中每个节点沿着坐标x1和x2方向的比值为±1。

同理，对于单胞模型内部一个线性子单元，用于求解一阶弹性影响函数ψ1的热自平衡虚拟载荷形式为

(1+ν)aε-(1+ν)caε(1+ν)aε

(1+ν)caε-(1+ν)aε(1+ν)caε]T

(19)

图1为包含3块夹杂的二维单胞细分网格模型。为了更加清晰地给出虚拟载荷矩阵中每一列矢量图，在本文中选取夹杂和基体的弹性模量分别为50 GPa和60 GPa，热膨胀系数分别为1.5×10-6/℃和1×10-6/℃，泊松比均为0.2。前三阶虚拟载荷和影响函数的物理意义讨论如下：

① 只要夹杂的材料性质相同，在每一块夹杂上相同位置节点上的矢量都是相同的，这个结论普遍成立，与单胞包含夹杂的数量无关。

③ψ1和χ1反映了单胞内材料的不连续性，因为它们是由分布在夹杂和基体交界处节点的分

图1 包含3块夹杂的二维周期复合材料单胞结构Fig.1 2D periodical composites unit cell structures with 3 inclusions

段线性载荷计算得到的。另一方面也反映了夹杂和基体之间的相互作用，是单胞内材料不均匀性的最基本关系，因为在基体和夹杂的内部节点是没有载荷作用的。

④u1为一阶摄动主项，本质上是χ1矩阵的3列和ψ1的1列虚拟位移的线性组合，所以一般来说，如果基体和夹杂的弹性模量、热膨胀系数相差不大时，使用MHM处理周期结构复合材料热力耦合问题摄动到一阶的精度就足够了。

图2 一阶虚拟载荷矢量图Fig.2 First-order virtual load vector

图3 二阶虚拟载荷矢量图Fig.3 Second-order virtual load vector

图4 三阶虚拟载荷矢量图Fig.4 Third-order virtual load vector

综合前三阶虚拟载荷矢量示意图可以得到以下结论：

1) 一阶弹性影响函数χ1和一阶热影响函数ψ1反映的是单胞内基体和夹杂之间的相互作用的结果，一阶摄动项u1是χ1矩阵中3列向量和ψ1的线性组合，反映单胞细观信息的位移场；包含一阶摄动项的MHM计算精度对于基体和夹杂之间弹性模量和热膨胀系数相差不大的情况下是足够精确的。

3) 一阶影响函数χ1和ψ1反映了基体内部、夹杂内部及基体和夹杂相互作用的结果，二阶摄动项u2是χ2矩阵中6列向量线性组合与ψ2矩阵中2列向量线性组合的和，反映单胞内材料更加准确的细观信息，所以对于一般周期复合材料而言，包含二阶摄动项的MHM的计算精度是足够的，但对于单胞内的基体和夹杂材料参数相差特别大的周期复合材料，极限情况下，如具有孔洞的周期结构，将孔洞视为夹杂，则基体和夹杂的材料弹性模量比例无限大，或某些颗粒增强结构，基体和夹杂的弹性模量比值无限小，此时包含二阶摄动项的计算精度是否足够需要专门研究。

2.2 周期复合材料杆单胞结构

文献[28]中给出了线性单元杆单胞弹性虚拟载荷和弹性特征位移分析结果，本文使用二次单元分析杆单胞的弹性影响函数和热影响位移，以及与之相对应的弹性虚拟载荷和热虚拟载荷。

考虑包含一块夹杂的杆单胞，如图5所示，对于这样简单的问题可以推导出其解析解形式。使用3个二次杆单元计算该单胞模型，该情况下，m=k=l=p=q=1。E1和E2分别表示基体和夹杂的弹性模量，a1和a2分别为基体和夹杂的热膨胀系数，L为单元长度，A为横截面积。

图5 包含一块夹杂的周期复合材料杆单胞结构Fig.5 Periodical composites rod unit cell structure with one inclusion

针对图5中的一个子单元，各阶影响函数控制方程如下：

1) 一阶热影响函数控制方程

(20)

2) 一阶弹性影响移函数控制方程

(21)

(22)

(23)

分别求解方程可得到单胞结构弹性影响函数和热影响函数为

(24)

(25)

值得注意的是，当相邻2个单元的弹性模量、热膨胀系数均相同时，即均质材料，式(24)右端项为0，同时对于非均质材料，若耦合参数相同(E2a2=E1a1)时，方程右端项也为0；对比式(25)，弹性自平衡虚拟载荷只有在不同弹性模量交界处才非0，即非均质材料，而热自平衡虚拟载荷只有在不相同耦合模量交界处才非0，基体和夹杂内部没有虚拟载荷，所以一阶虚拟载荷反映的是基体和夹杂之间的相互作用，是反映单胞内部结构最基本的信息载体。

式(12)、式(24)、式(25)联立，可得得到均匀化弹性模量和均匀化热膨胀系数分别为

(26)

二阶热虚拟载荷和弹性虚拟载荷分别为

(27)

二阶弹性影响函数和热影响函数分别为

(28)

三阶热虚拟载荷和弹性虚拟载荷分别为

(29)

三阶热影响函数和弹性影响函数分别为

(30)

从二阶虚拟载荷和三阶虚拟载荷形式可以看出，当弹性模量相同时，弹性虚拟载荷为0，当耦合参数相同时，热虚拟载荷为0，即在非均质材料交界处，弹性虚拟载荷一定非0，而热虚拟载荷可能为0，所以弹性问题只考虑不同材料之间弹性参数的异同，而热问题需要同时考虑弹性参数和热膨胀系数；另一方面，弹性虚拟载荷和热虚拟载荷在单元内部节点不为0，所以二阶和三阶虚拟载荷除了反映基体和夹杂之间的相互作用，也刻画了基体和夹杂内部更加详细的细观信息。

3 数值分析

3.1 二维周期复合材料单胞结构

结构大小为15 mm×15 mm，单胞内含有4块夹杂，如图6所示。基体的弹性模量和热膨胀系数分别为E1=2×109Pa,a1=3×10-6/℃，夹杂的弹性模量和热膨胀系数分别为E2=6×1010Pa,a2=1×10-6/℃，基体和夹杂的泊松比均为ν=0.2，结构温差T=16(x2+y2)。

分别使用式(12)计算得到均匀化弹性模量和均匀化热膨胀系数分别为EH=2.67×109Pa,aH=2.5×10-6/℃。图7给出了结构沿A、B、C、D上节点在x1方向的位移曲线。表2给出了结构总势能泛函Π计算结果，其中FEM为有限元解，作为标准检验MHM的求解精度，MHM1表示包含一阶摄动项的计算结果，MHM2表示包含二阶摄动项的计算结果，MHM3表示包含三阶摄动项的计算结果。结合图7和表2可以得到：①包含一阶摄动位移的计算结果MHM1误差较大；②MHM2和MHM3的精度明显高于MHM1；③MHM3精度高于MHM2。

图6 二维周期复合材料单胞结构Fig.6 Unit cell of 2D periodical composite structure

3.2 二维周期复合材料多胞结构

结构大小为45 mm×45 mm，包含5×5个单胞，每个单胞内部含有一块夹杂，如图8所示，材料参数及温差和3.1节相同。影响函数和宏观位移导数分别使用超单胞周期边界条件和微分求积有限单元法[30]计算得到，均匀化弹性模量和均匀化热膨胀系数分别为EH=2.42×109Pa,aH=2.63×10-6/℃。表3给出了结构的总势能泛函计算结果，FEM、MHM1、MHM2、MHM3的含义与3.1节相同。图9给出了结构沿A′、B′、C′、D′上节点x1方向的位移曲线。结合图9和表3可以得到：①包含一阶摄动位移的计算结果MHM1误差较大；②MHM2和MHM3的计算精度明显好于MHM1；③MHM2和MHM3的计算精度较高，且两者相差不大；④二阶摄动项对计算精度的影响较大，而三阶摄动项对计算精度的影响可以忽略。

图7 沿纵线A、B、C、D上节点位移曲线Fig.7 Nodal displacement curves along longitudinal lines A,B,C and D

表2 二维周期复合材料单胞结构的势能泛函
Table 2 Potential energy function for 2D periodical
composites unit cell

摄动位移Π/(10-6J)MHMFEMΠFEM-ΠMHMΠFEM×100%MHM1-2.47515.8MHM2-2.916-2.9390.78MHM3-2.9370.068

图8 二维周期复合材料多胞结构Fig.8 2D periodical composite multi-cell structure

表3 二维周期复合材料多胞结构的势能泛函Table 3 Potential energy function for 2D periodical composite multi-cell structure

图9 沿纵线A′、B′、C′、D′上节点位移曲线Fig.9 Nodal displacement curves along longitudinal lines A′，B′，C′ and D′

4 结 论

本文通过构造摄动项的全解耦格式推导了热力耦合问题高阶MHM的数学表达式和有限元矩阵列式，将弹性影响函数和热影响函数分别比拟为弹性虚拟位移和热虚拟位移，利用虚拟载荷的自平衡性质、量纲分析和空间矢量分布揭示了各阶摄动项的物理意义；数值算例分析证明了公式推导结果和物理意义分析的正确性。通过热力耦合问题MHM各阶摄动项物理意义的分析和数值计算结果，给出了二阶摄动项对于周期复合材料计算精度的必要性依据，为热力耦合问题MHM的应用奠定了力学基础。