一种有限变形下三维随机非均质材料的热弹性均化方法

摘要

本发明公开了一种有限变形条件下三维随机非均质材料的热弹性均化方法，包括：1)基于随机序列添加RSA方法构建颗粒随机分布的复合材料代表体积单元三维RVE数值模型，并消除了当颗粒体积分数较大时的颗粒重叠现象；2)在热弹性环境下对该复合材料三维RVE数值模型进行有限元分析计算，求得RVE数值模型有效性质的数值解；3)建立随机均化模型求解复合材料的宏观有效性质，并以该复合材料的宏观有效性质作为其真实有效性质，为复合材料的宏观有效性质提供了更为可靠的均化结果，为新型先进材料的使用和结构的优化设计提供充分的依据。

说明书

技术领域

本发明涉及复合材料计算力学，具体地说是针对一种有限变形下三维随机 非均质材料的热弹性均化方法。

背景技术

对非均质材料的研究和分析是一个经典的问题，在连续介质力学范围内， 存在着宏观力学和细观力学两种研究方法。宏观力学方法着重于研究由复合材 料组成的结构(如层合板)的性能，如刚度、强度、热应力的计算等；细观力 学方法主要在细观结构尺度上研究各相材料的相互作用，以及建立材料宏观性 能与细观结构参数之间的关系。这里的宏观性能也称为有效性质，是指能够在 宏观尺度上实验测量的材料性能，包括有效弹性常数、有效传导率等。

由于非均质材料的宏观长度尺寸显著地大于微尺度下不同组分的长度尺 度，若对所有微尺度组分均采用显式求解，无论是解析法还是数值解法，都是 极具挑战性且代价昂贵的工作。一种显著减少分析此类问题所需成本的方法就 是用有效均质材料替代原先的非均质材料并通过构建反映宏观性质的有效本构 方程进而获得材料的宏观有效性质，即均化法(如图1所示)。该方法于20世 纪70年代由法国科学家提出并应用到具有周期性结构的材料分析中，通过使用 势能将材料的有效特征量引入并进行求解。此方法使用了从材料中取出的一个 表征体积单元(RVE:RepresentativeVolumeElement)。因其可以极大地减少复合 设计参数的数量，且随着工程需求的变化，该方法得到不断发展，特别是将有 限元方法与之结合后，形成了目前国际上多尺度分析非均质材料性质的最有效 方法之一，期间我国著名学者刘书田、崔俊芝和张洪武等人也做出了重要的贡 献，并用其来预测材料的弹性常数、导热系数、热膨胀系数，以及弹塑性问题、 黏弹性问题等。

而对于非均质材料在多物理场中的均质化问题，前人已经对线性热弹性情 况进行了探索。Rosen、Hashin和Torquato等人较早得出并使用了等效比热及热 膨胀系数的边界值。然而，求得这种边界值在有限变形条件下却具有较大难度， 即使是纯力学分析的情况。目前，研究的重心在于表征非均质材料的宏观尺度 热弹性响应方面。出于此种目的，Francfort提出了线性热弹性周期介质的均化 方法。此方法后又经L’Hostis和Devries、Yu和Tang、Terada等多人引用及发 展。此外Ghosh和Liu提出了基于微极理论的均化方法。

综上可知，线性情况下非均质材料的热弹性均化问题已经有所成就。但是， 具有较大温度偏移的有限变形情况却鲜有发展。Laschet和Alzina等人导出了忽 略了变形相位的非线性热传导方程。另外，如果忽略能量平衡的均化，则可对 非线性力学分析中的热诱导效应进行分析。这方面的研究起始于Aboudi和 Arnold，此二者建立了均化的近似解法。此外，Khisaeva和Osroja-Starzewski 及Miehe等人在这方面也有所贡献，但是只有Aboudi及等人较完整地 提出了非均质材料在耦合场中有限热弹性分析的框架。

发明内容

基于传统均匀化方法的不足之处，本发明在有限变形条件下，对非均质材 料微观结构在热弹性物理场中的热力学响应进行了分析，并提出了非均质材料 的热弹性随机均化框架。基于热弹性力学，提出了用于均化法的基本尺度转换 方法。施加宏观结构温度以得到基于温度场的微观结构参量，并对材料微观结 构的参数的随机性与相关性给予了充分考虑。微观力学分析步骤可分解为纯力 学分析及纯热学分析两步，以此来获得宏观尺度的均化响应。最终得到任意颗 粒体积分数情况下的加细数值解，从而得到更为精确的宏观有效参数，为新型 先进材料和结构的优化设计提供客观、充分和真实的依据。

本发明是通过以下技术方案来解决的：

一种有限变形下三维随机非均质材料的热弹性均化方法，该方法包括下述 步骤：

1)基于随机序列添加法(RSA)使用FORTRAN语言构建颗粒随机分布 的复合材料表征体积单元(RVE)的三维数值模型；

对含不规则颗粒夹杂的复合材料，其颗粒形状通过设置不同的边界曲线来 描绘；颗粒中心点的坐标通过随机数产生；

2)对颗粒随机分布的复合材料表征体积单元(RVE)三维数值模型进行有 限元分析计算，求得表征体积单元(RVE)数值模型有效性质的数值解；

①对表征体积单元(RVE)进行网格划分；

②对基体与颗粒组成的随机分布的复合材料的界面施加热力学边界条件， 求得有限元方法下表征体积单元(RVE)有效热力学性质的数值解；

3)建立随机均化模型求解复合材料的宏观有效性质；

针对颗粒随机分布的复合材料表征体积单元(RVE)三维数值模型中的未 知参数，选取样本空间n，对n个颗粒随机分布的复合材料表征体积单元(RVE)， 经计算后得到一系列的随机的数值解，运用数理统计的方法对这些数值解进行 统计处理，并将数理统计的平均值作为颗粒随机分布的复合材料宏观有效性质 的预测值。

进一步地，所述步骤1)中，基于随机序列添加法(RSA)方法使用FORTRAN 语言构建颗粒随机分布的复合材料表征体积单元(RVE)三维数值模型，是在 三维基体材料中逐个生成球形或椭球形颗粒夹杂，且任意两个颗粒的中心距必 须大于或等于颗粒直径，即颗粒不能重叠。

进一步地，所述步骤1)中，采用随机序列添加法构建颗粒随机分布的复合 材料表征体积单元(RVE)三维数值模型，包括下述步骤：

1a)根据颗粒体积分数求v₂和颗粒数N_p求颗粒半径R：

$R=\sqrt[3]{3v_2/\left(4N_p\pi \right)}---\left(1\right)$

1b)利用随机数生成初始颗粒P₁的中心点坐标(x₁,y₁,z₁)，生成第二个颗粒P₂的中心点坐标(x₂,y₂,z₂)，计算(P₁，P₂)的中心距L₁₂：

$L_{12}=\sqrt{{(x_2-x_1)}^2+{(y_2-y_1)}^2+{(z_2-z_1)}^2}---\left(2\right)$

1c)若L₁₂≥2R，则第二个颗粒满足要求，接着生成第三个颗粒P₃且判断 (P₁，P₃)的中心距L₁₃、(P₂，P₃)的中心距L₂₃是否均大于2R；若L₁₂<2R，则重新生 成第二个颗粒P₂直至不发生颗粒重叠为止。

进一步地，所述步骤2)中，施加热力学边界条件，求得有限元方法下表征 体积单元(RVE)有效热力学性质的数值解，通过下述方式实现：

2a)对表征体积单元RVE进行网格划分，选用八节点六面体单元；

2b)对于施加于非均质材料上的普通热力学边界值问题，其在参考构型 上的线性动量平衡方程为：

$Div\left[P\right]+{\rho }_{0}f={\rho }_0\stackrel{··}{x}---\left(3\right)$

将其解与能量平衡方程耦合得：

$\stackrel{·}{e}=P·\stackrel{·}{F}-Div\left[q_0\right]+{\rho }_{0}r---\left(4\right)$

其中，Div[·]表示散度，P为第一皮奥拉-基尔霍夫应力，ρ₀为参考质量密度， x是位置矢量，e是单位体积上的内能，F是变形梯度，q₀表示热流量矢量， f(r)是独立于变形的单位质量体应力，r表示材料单位质量的热供应量；

对于施加于非均质材料上的普通热力学边界值问题，其本构方程如下：

$e=\tilde{e}(F,\eta )=\hat{e}(F,\theta ),P=P(F,\theta ),q_0=q_0(F,g_0,\theta )---\left(5\right)$

由热力学第二定律导出的本构限制条件如下：

$P=\frac{\partial \psi }{\partial F}=\frac{\partial \tilde{e}}{\partial F},D_0:=-\frac{q_0·g_0}{\theta }=-q_0·Grad\left[ln\theta \right]\ge 0---\left(6\right)$

其中，η是单位体积上的熵，θ表示温度，g₀表示温度梯度，ψ是参考构 型单位体积上的亥姆霍兹自由能，热耗散J＝det[F]，并有 且是热力学边界值问题中唯一的耗散源。

在一级均化框架中，通过使用勒让德变换得到确保宏观尺度组分 响应的热力学相容性条件：

其中，c是单位体积在恒定变形时的比热；

2c)在有限元求解方法中，表征体积单元RVE有效热力学性质的数值解微 观结构通过两步求解过程得到：力学求解阶段及热学求解阶段。

在力学求解阶段，通过将宏观变形F施加于原微观结构并使温度保持为 可以得到变形后的微观结构并以此得到力学宏观参量的数值解其方程表 示为

$Div\left[P(F,\bar{\theta })\right]=0---\left(27\right)$

在热学求解阶段，将宏观温度梯度${\bar{g}}_0$施加于固定的微观结构从而得到热学宏 观参量的数值解，其方程表示为

$-Div\left[q_0(F,g_0,\bar{\theta })\right]=0---\left(29\right).$

进一步地，所述步骤3)中，建立随机均化模型求解复合材料的宏观有效性 质，通过下述方式实现：

3a)建立随机均匀化模型

用基体与颗粒的性质参数(k₁,μ₁,k₂,μ₂,α₀,k_c,γ₁,γ₂,γ₃,β₁,β₂,β₃)及颗粒体积分 数v₂来表示复合材料的宏观有效性质即有效量，即

式中，代表复合材料的宏观有效性质，k₁和μ₁分别是基体的体积模量与剪切 模量，k₂和μ₂分别是颗粒的体积模量与剪切模量；α₀是基体的热膨胀系数，k_c是 基体的导热系数；{γ₁,γ₂,γ₃,β₁,β₂,β₃}是Ogden材料的参数，满足f 表示各有效量是关于材料参数(k₁,μ₁,k₂,μ₂,α₀,k_c,γ₁,γ₂,γ₃,β₁,β₂,β₃,v₂)的函数，下 标FEM表示采用有限元方法求解；

当考虑表征体积单元(RVE)中颗粒位置分布的随机性时，首先选取样本 容量即选定n个表征体积单元(RVE)，每个RVE中颗粒的位置都不同，而颗 粒位置的不同直接决定了有效性质的差异；

3b)求单个表征体积单元(RVE)的有效性质

对第i个RVE，求出其有效性质

其中，代表第i个表征体积单元(RVE)的有效性质，i＝(1,2,…n)；

3c)对n个表征体积单元(RVE)有效性质进行数理统计

通过上述步骤得到n个RVE的有效性质，对其进行数理统计。

进一步地，所述步骤3a)中，复合材料的宏观有效性质包括应力张量有效热流量张量变形梯度张量F。

进一步地，所述步骤3c)中，对n个表征体积单元(RVE)有效性质进行 数理统计，包括下述方法：

求张量矩阵的2-范数，再求其平均值和均方差，表示为：

式中，表示各有效张量的2-范数，为n个表征体积单元(RVE)有 效张量的平均值，是n个表征体积单元(RVE)有效张量的均方差；上 述公式求得的值和即可视为颗粒随机分布时均化模型下的有效 常数。

进一步地，所述步骤3c)中，对n个表征体积单元(RVE)有效性质进行 数理统计，包括下述方法：

直接求张量矩阵中各元素的均值和均方差，则其结果也为矩阵形式，表示 为：

$E\left(a\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{a}_i---\left(35\right)$

$D\left(a\right)=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{(a_i-E(a\left)\right)}^2}---\left(36\right)$

式中，a表示各张量矩阵中的各元素，E(a)及D(a)分别表示其均值及均方差；

以此得到的均值及均方差矩阵与原张量矩阵形式相同。

本发明同现有技术相比，有以下几个优点：

本发明的研究内容具有明显的前瞻性和挑战性，研究工作属于连续介质力 学、热力学和材料科学领域中的前沿性课题，具有现实的应用前景、学术价值 和理论意义。就本发明来说，它具有如下特点：

(1)综合考虑了非均质材料组份参数的随机性及参数之间的相关性；

(2)建立三维模型，解决了复杂RVE的数值建模难题；

(3)充分运用Fortran语言强大的计算功能，获得了比传统的细观力学解 析边界值更为精确的计算结果，为非均质材料结构的优化设计奠定了基础。

附图说明

图1是本发明具体实施方案的流程图。

图2是用均匀的等效材料替代宏观结构非均质材料示意图。

图3是仅考虑各元素随机性时RVE的应力、变形梯度、热流量分布图。

图3(a)是变形、应力、能量的μ-3σ分布图。

图3(b)是变形、应力、能量的均值μ分布图。

图3(c)是变形、应力、能量的μ+3σ分布图。

图4是考虑各元素随机性及相关性时RVE的热流量分布图。

图4(a)是RVE的热流量μ-3σ分布图。

图4(b)是RVE的热流量μ分布图。

图4(c)是RVE的热流量μ+3σ分布图。

具体实施方式

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

如图1所示，本发明主要有三个主要流程。首先基于随机序列添加法(RSA)， 使用FORTRAN语言构建了一种颗粒随机分布的三维RVE(表征体积单元)模 型。同时，在有限元方法的基础上，通过蒙特卡洛模拟法获得了材料在线性边 界条件下宏观有效量的数值解，建立随机均化模型求解复合材料的宏观有效性 质，并绘制了RVE内部变形梯度场、应力场以及热流量场的三维分布图。具体 技术方案如下：

1、构建随机RVE数值模型

基于随机序列添加法即RSA方法使用FORTRAN语言构建颗粒随机分布的 复合材料表征体积单元(RVE)三维数值模型，是在三维基体材料中逐个生成 球形或椭球形颗粒夹杂，且任意两个颗粒的中心距必须大于或等于颗粒直径， 即颗粒不能重叠，包括下述步骤：

根据颗粒体积分数v₂和颗粒数N_p求颗粒半径R：

$R=\sqrt[3]{3v_2/\left(4N_p\pi \right)}---\left(1\right)$

产生随机数得到初始颗粒P₁的中心点坐标(x₁,y₁,z₁)，及第二个颗粒P₂的中心 点坐标(x₂,y₂,z₂)，计算(P_1，P₂)的中心距L₁₂：

$L_{12}=\sqrt{{(x_2-x_1)}^2+{(y_2-y_1)}^2+{(z_2-z_1)}^2}---\left(2\right)$

若L₁₂≥2R，则第二个颗粒满足要求，接着生成第三个颗粒P₃且判断(P₁，P₃)的 中心距L₁₃、(P₂，P₃)的中心距L₂₃是否均大于2R；若L₁₂<2R，则重新生成第二个 颗粒P₂直至不发生颗粒重叠为止。

当确定第n(n≥2)个颗粒的中心点坐标时，需计算第n个颗粒与所有前n-1个 颗粒的中心距，若L_1n,L_2n,...,L_n-1,n这n-1个中心距都大于2R，则第n个颗粒生成成 功，否则重新确定第n个颗粒的中心点坐标直至它与所有已生成的颗粒都不相 交。图2是用均匀的等效材料替代宏观结构非均质材料示意图。

2、使用有限元法进行热力学分析

2.1对表征体积单元RVE进行网格划分

选用八节点六面体单元对表征体积单元RVE进行网格划分，网格的密度可 根据颗粒体积分数的增长而适当增大。

2.2基于连续介质力学的均化框架

2.2.1平衡方程

非均质材料的参考(现)构型为位置矢量表示为X(x)，边界 上的单位外法矢为N(n)，无限小体积为dV(dv)，无限小表面积dA(da)， 绝对温度值θ。上的梯度算子表示为$Grad[·]\equiv \partial [·]/\partial X\left(Grad\right[·]\equiv \partial [·]/\partial x),$散 度算子表示为Div[·](Div[·])。位置变动泛函通过时间t将和联系起来： 现构型的变形梯度表示为若令J＝det[F]则有 dv＝JdV。由变形梯度可以得到右(左)柯西-格林应变C＝F^TF(b＝FF^T)并可将第 一皮奥拉-基尔霍夫应力P投影到基尔霍夫应力τ＝PF^T及柯西应力T＝J^-1τ。参考 (现)构型中得到的温度梯度和热流量矢量分别表示为$(g=\partial \theta /\partial x=F^{-T}{g}_0)$及q₀(q＝J^-1Fq₀)。

对于施加于上的普通热力学边界值问题(BVP)，其在上的线性动量 平衡方程为

$Div\left[P\right]+{\rho }_{0}f={\rho }_0\stackrel{··}{x}---\left(3\right)$

将其解与能量平衡方程耦合得

$\stackrel{·}{e}=P·\stackrel{·}{F}-Div\left[q_0\right]+{\rho }_{0}r---\left(4\right)$

其中，ρ₀(ρ＝J^-1ρ₀)为参考(现)质量密度，f(r)是独立于变形的单位质量体 应力，e是内部能量。为了简洁，之后的各热力学参量表示的皆为单位体积的 情况。

2.2.2本构限制条件

热力学BVP问题的本构方程如下：

$e=\tilde{e}(F,\eta )=\hat{e}(F,\theta ),P=P(F,\theta ),q_0=q_0(F,g_0,\theta )---\left(5\right)$

由于热弹性响应的特殊性，本构方程中没有内部变量，且因η＝η(F,θ)，故有 内能导出赫尔姆霍兹自由能，其勒让德变换形式为ψ＝ψ(F,θ)＝e-θη， 并可由此得到因此，作用于P和q₀的本构限制可由热力学第二定律 导出如下：

其中热耗散并有是现问题中唯一的耗散源。 而对函数形式的观测不变性要求中导出的其他本构限制条件使得ψ和q₀独立于 F的旋度。

使用ψ可将能量平衡方程(4)表示为以下形式：

$c\stackrel{·}{\theta }=\theta \frac{\partial P}{\partial \theta }·\stackrel{·}{F}-Div\left[q_0\right]+{\rho }_{0}r---\left(7\right)$

其中$c=\tilde{c}(F,\theta )=\partial \hat{e}/\partial \theta =\theta (\partial \eta /\partial \theta )=-\theta ({\partial }^2\psi /\partial {\theta }^2),$表示中单位体积在恒定变形时的 比热，表示高夫-焦耳效应。通常认为此效应非常微弱并将其忽略。

然而，此处为了彰显此发明的普适性而保留此量。另外，在很多热力学分析中 比热都被视为常量。

2.2.3热弹性分析中的热力势与熵

在温度区间θ₀与θ之间对进行积分，得到下列表达式

$\hat{e}(F,\theta )=e_0\left(F\right)+{\int }_{{\theta }_0}^{\theta }c(F,{\theta }^{\prime })d\theta ---\left(8\right)$

此处，记号[·]₀(F)表示e在θ₀处的显式估值。同样，在相同的温度区间对 进行积分可以得到

$\eta (F,\theta )={\eta }_0\left(F\right)+{\int }_{{\theta }_0}^{\theta }\frac{1}{{\theta }^{\prime }}c(F,{\theta }^{\prime }){\mathrm{d\theta }}^{\prime }---\left(9\right)$

最后，将(8)、(9)式代入勒让德变换ψ＝e-θη中可以得到

其中，因此可知，求解热弹性下的热力势须已知参考温度下的等温特 性及c值。本发明在刚性媒介中考察其热传导性能时，忽略了F对其的影响。

2.3均化边界值问题

对于宏观结构的BVP问题，由于缺少一系列必要的微观参数，求解起来将 会非常繁琐。解决的办法是引入体积适当小的表征体积单元RVE利用 微观尺度特性和宏观尺度响应之间的尺度转换关系来确定非均质材料的有效宏 观物理性质。表征体积单元RVE的参数以的形式描述，其中代表微观结构 普通参数张量。将表征体积单元RVE的参数代入前述均化框架得到所谓一级均 化框架。其各方程皆与前述框架形式相同，值得一提的是，通过使用宏观勒让 德变换可得到下述保证宏观尺度组分响应的所谓热力学相容性条件：

其中，c是单位体积在恒定变形时的比热。

2.3.1尺度转换

宏观尺度到微观尺度的变换：为在上施加边界条件，将宏观结构特定点 上的基础变量用于微观结构分析中，并由此得到中的解{x,θ}。

微观尺度到宏观尺度的变换：将微观尺度解{x,θ}代入组分本构方程，可以 得到总体的张量场变量A，此变量将用于区域Ω上的如下体积平均过程以求解宏 观BVP问题：

$<A{>}_{\Omega }:=\frac{1}{\left|\Omega \right|}{\int }_{\Omega }Ad\Omega ---\left(12\right)$

宏观尺度基本热力学变量可表示为：

2.3.2施加宏观温度

假设宏观结构的参考构型无应力或热流量约束，且非均质材料的参考 温度始终为θ₀。均化后的材料温度便可较便捷地选取为此处，分别将(8) 及(9)式代入和的定义式(12)，便可得到

${\bar{e}}_0\left(\bar{F}\right)+{\int }_{{\bar{\theta }}_0}^{\bar{\theta }}\bar{c}(\bar{F},{\bar{\theta }}^{\prime })d{\bar{\theta }}^{\prime }=<e_0\left(F\right)>+<{\int }_{{\bar{\theta }}_0}^{\theta }c(F,{\theta }^{\prime })d{\bar{\theta }}^{\prime }>---\left(14\right)$

${\bar{\eta }}_0\left(\bar{F}\right)+{\int }_{{\bar{\theta }}_0}^{\bar{\theta }}\frac{1}{{\bar{\theta }}^{\prime }}\bar{c}(\bar{F},{\bar{\theta }}^{\prime })d{\bar{\theta }}^{\prime }=<{\eta }_0\left(F\right)>+<{\int }_{{\bar{\theta }}_0}^{\theta }\frac{1}{{\bar{\theta }}^{\prime }}c(F,{\theta }^{\prime })d{\bar{\theta }}^{\prime }>---\left(15\right)$

在此处，将与其对应的微观部分联系起来的关键在于设一种服 从一级均化框架条件的材料为严格纯热材料，则其在中所有点上及热动力方程 (14)、(15)中的c＝c(θ)可简化为

${\bar{e}}_0+{\int }_{{\bar{\theta }}_0}^{\bar{\theta }}\bar{c}\left({\bar{\theta }}^{\prime }\right)d{\bar{\theta }}^{\prime }=e_0+<{\int }_{{\bar{\theta }}_0}^{\theta }c\left({\theta }^{\prime }\right)d{\bar{\theta }}^{\prime }>---\left(16\right)$

${\bar{\eta }}_0+{\int }_{{\bar{\theta }}_0}^{\bar{\theta }}\frac{1}{{\bar{\theta }}^{\prime }}\bar{c}\left({\bar{\theta }}^{\prime }\right)d{\bar{\theta }}^{\prime }={\eta }_0+<{\int }_{{\bar{\theta }}_0}^{\theta }\frac{1}{{\bar{\theta }}^{\prime }}c\left({\theta }^{\prime }\right){\mathrm{d\theta }}^{\prime }>---\left(17\right)$

其中以及皆为常量。

在当前的均化法中，既要使方程(16)得到满足，又要得到恒定的一阶均 化框架的必要条件是确保在微观力学求解步骤中所有包含温度的函数形式在宏 观温度下进行求值。这些形式包含变量

2.3.3宏观尺度下的比热

为了得到的表达式，将宏观尺度下的温度代入到表达式(14)中，可以得 到

${\bar{e}}_0\left(\bar{F}\right)+{\int }_{{\bar{\theta }}_0}^{\bar{\theta }}\bar{c}(\bar{F},{\bar{\theta }}^{\prime })d{\bar{\theta }}^{\prime }=\bar{e}(\bar{F},\bar{\theta })=<e(F,\bar{\theta })>=<e_0\left(F\right)+{\int }_{{\bar{\theta }}_0}^{\theta }\int c(F,{\bar{\theta }}^{\prime })d{\bar{\theta }}^{\prime }>---\left(18\right)$

其中，微观场量F与其宏观状态下的相关联。为了使满足式(11) 中所述的相容性条件，将的所有衍生形式皆代入得到

$\bar{c}(\bar{F},\bar{\theta })=\frac{\partial \bar{e}}{\partial \bar{\theta }}|{}_{\bar{F}}=\frac{\partial <e(F,\bar{\theta })>}{\partial \bar{\theta }}=<\frac{de}{d\bar{\theta }}>=<\frac{\partial e}{\partial F}{|}_{\bar{\theta }}·\frac{\partial F}{\partial \bar{\theta }}+\frac{\partial e}{\partial \bar{\theta }}{|}_F>=<[\frac{\partial e_0}{\partial F}+{\int }_{{\bar{\theta }}_0}^{\bar{\theta }}\frac{\partial c}{\partial F}d{\bar{\theta }}^{\prime }]·\frac{\partial F}{\partial \bar{\theta }}+c(F,\bar{\theta })>---\left(19\right)$

其中，表示在宏观变形状态为时，微观变形场对宏观温度的梯度。 该变量与非均质材料单元在恒定变形状态时温度的改变引起的能量变化有 关。当为定值时，微观结构的变形可以有所变动，相应的内能改变也会导致的变化。以下是一些特殊情形：

微观结构中，当c独立于变形时，有由此可看出， 仍然取决于和

只有在纯非弹性热材料的情况下，才有成立。

2.4微观力学求解步骤

形变场的求解仅仅取决于而跟宏观温度梯度所造成的温度场无关。因此 微观结构求解过程可分为先后两步：力学求解阶段及热学求解阶段。在力学求 解阶段，通过将宏观变形施加于并使温度保持为可以得到变形后的微观 结构并以此得到与温度梯度无关的各宏观参量。在热学求解阶段，将宏观温 度梯度施加于固定的微观结构从而得到基于温度梯度的宏观参量。

2.4.1力学求解阶段

相容性条件

对于施加于微结构体积样本的宏观求解方案：

$\frac{\partial \bar{e}}{\partial \bar{F}}{|}_{\bar{\eta }}·\stackrel{·}{\stackrel{‾}{F}}+\frac{\partial \bar{e}}{\partial \bar{\eta }}{|}_{\bar{F}}\stackrel{·}{\stackrel{‾}{\eta }}=<\frac{\partial e}{\partial F}{|}_{\eta }\stackrel{·}{F}+\frac{\partial e}{\partial \eta }{|}_F\stackrel{·}{\eta }>=<P(F,\bar{\theta })·\stackrel{·}{F}+\bar{\theta }\stackrel{·}{\eta }>=<P(F,\bar{\theta })·\stackrel{·}{F}>+\bar{\theta }<\stackrel{·}{\eta }>---\left(20\right)$

由于的两端仍然保持平衡，上述方程即可分解为下述两个方程

$\frac{\partial \bar{e}}{\partial \bar{\eta }}{|}_{\bar{\eta }}·\stackrel{·}{\stackrel{‾}{F}}=<P(F,\bar{\theta })·\stackrel{·}{F}>,\frac{\partial e}{\partial \bar{\eta }}{|}_{\bar{F}}\stackrel{·}{\stackrel{‾}{\eta }}=\bar{\theta }<\stackrel{·}{\eta }>---\left(21\right)$

并且因为有成立，且独立于

$\frac{\partial \bar{e}}{\partial \bar{\eta }}{|}_{\bar{F}}=\bar{\theta }---\left(22\right)$

便可得到式(11)的相容性条件。然而，尽管在微观尺度下独立的过程变 量，F却仍取决于另一方面，如果下式

$<P(F,\bar{\theta })·\stackrel{·}{F}>=<P(F,\bar{\theta })>·<\stackrel{·}{F}>---\left(23\right)$

成立，施加便可以确保相容性条件(22)成立，因以及皆趋近于 零。并且，因为独立于F，由式(21)便可得到

$\frac{\partial \bar{e}}{\partial \bar{F}}{|}_{\bar{\eta }}<P>---\left(24\right)$

以及隐含的关系式

边界值问题

前述相容性条件部分的内容是基于假设(23)式的。以下式子

$\bar{F}=<F>,\bar{P}=<P>---\left(25\right)$

可重新表述为

$\bar{P}·\stackrel{·}{\stackrel{‾}{F}}=<P·\stackrel{·}{F}>---\left(26\right)$

此表达式一般称为微观-宏观做功标准，表示宏观尺度下的机械功与微观尺度下 的功等值。在此处，它作为保证热力学相容性条件的一个必要条件而 存在。此表达式可以看做力学求解阶段对边界等条件的约束。基尔霍夫应力的 体积平均也可以用于基本宏观应力的计算：

为了满足(26)式，此处避免了体应力及动态条件的使用，因其会增大宏 观应力，而此处只需要求出其本征响应。如此可确保中的P不会发散。另一方 面，位移的不连续性会导致材料的非均匀性很难严格适用于热弹性均化框架。 因此，x在中不存在跃变。故力学求解阶段可以表示为

$Div\left[P(F,\bar{\theta })\right]=0---\left(27\right)$

并选择线性位移边界条件(ε-LN-BCs),即施加并得到从当前 情况的设定中，可以得到合适的F并施加理想的<P>，并可以对<P>和<F>的敏度 进行量化。

2.4.2热学求解阶段

从式(13)中的定义以及式(11)中相应的相容性条件，将施 加于微观尺度耗散函数可以得到

$-\frac{{\bar{q}}_0(\bar{F},{\bar{g}}_0,\bar{\theta })·{\bar{g}}_0}{\bar{\theta }}=<-\frac{q_0(F,g_0,\bar{\theta })·g_0}{\bar{\theta }}>---\left(28\right)$

这条微观到宏观的热耗散标准保证了微观尺度的热耗散同宏观尺度下的热耗散 等值。

中BVP问题在热学求解阶段的方程可表述为

$-Div\left[q_0(F,g_0,\bar{\theta })\right]=0---\left(29\right)$

其中F由力学求解阶段得到。本方程可在施加了的情况下，在线性位移 边界条件下求解。又因下式

$<-\frac{q_0(F,g_0,\bar{\theta })·g_0}{\bar{\theta }}>=-\frac{<q_0(F,g_0,\bar{\theta })>·<g_0>}{\bar{\theta }}---\left(30\right)$

成立，在此情况下由<q₀>到的赋值可确保式(28)得以满足。需要强 调的是中的θ分布仅仅是为了导出g₀，所以其大小与耗散无关且不影响其值。

值得一提的是，如果在微观尺度下由热传导系数张量Κ＝K(F,θ)得到的q₀是 不变量，则由q₀＝-Kg₀可得到因K与g₀的模值或者方向皆无关，故给 定的力学求解阶段之后的热学求解与经典的线性纯热问题求解有着相同的架 构。故当有可得出也就是说，宏观尺度下的热传导系数张量 并不决定于温度梯度

3、构建随机均化模型

3.1变量的随机性与互相相关性

考虑以下变量的随机性：基体的体积模量k₁与剪切模量μ₁、颗粒的体积模量 k₂和剪切模量μ₂、基体的热膨胀系数α₀与导热系数k_c以及所使用的Ogden材料 的参数{γ₁,γ₂,γ₃,β₁,β₂,β₃}。其中{γ₁,γ₂,γ₃,β₁,β₂,β₃}必须满足式子假 设以上变量皆服从正态分布，即其中代表以上各参数，μ和 msd分别表示各变量的均值和均方差。对于每一组确定的参数，可以用蒙特卡洛 法来产生服从正态分布的n组样本。将这n组样本代入均化步骤，得到n组不 同的结果。最后对这n组有效性能数据进行数理统计，可得到这些结果的均值、 均方差等数字特征值。

考虑以下变量之间的线性相关性：基体的体积模量k₁与剪切模量μ₁之间、颗 粒的体积模量k₂和剪切模量μ₂之间、基体的热膨胀系数α₀与导热系数k_c之间以 及Ogden材料参数{γ₁,γ₂,γ₃,β₁,β₂,β₃}之间。其中，{γ₁,γ₂,γ₃,β₁,β₂,β₃}各参量之间 的相关性应有如下关系：ρ_γ1γ2＝ρ_γ1γ3＝ρ_γ2γ3,ρ_β1β2＝ρ_β1β3＝ρ_β2β3。

3.2进行统计处理

用基体与颗粒的性质参数(k₁,μ₁,k₂,μ₂,α₀,k_c,γ₁,γ₂,γ₃,β₁,β₂,β₃)及颗粒体积分 数v₂来表示复合材料的宏观有效性质即有效量，如应力张量P、有效热流量张量 变形梯度张量F等，即

式中，代表复合材料的宏观有效性质，如应力张量有效热流量张量变形梯度张量F等，k₁和μ₁分别是基体的体积模量与剪切模量，k₂和μ₂分别是颗 粒的体积模量与剪切模量；α₀是基体的热膨胀系数，k_c是基体的导热系数； {γ₁,γ₂,γ₃,β₁,β₂,β₃}是Ogden材料的参数，满足f表示各有效量是关 于材料参数(k₁,μ₁,k₂,μ₂,α₀,k_c,γ₁,γ₂,γ₃,β₁,β₂,β₃,v₂)的函数，下标FEM表示采用有 限元方法求解。

当考虑RVE中颗粒位置分布的随机性时，首先选取样本容量即选定n个 RVE，每个RVE中颗粒的位置都不同，而颗粒位置的不同直接决定了有效性质 的差异。

3.2.1求单个RVE的有效性质

对第i个RVE，求出其有效性质

其中，代表第i个RVE的有效性质，i＝(1,2,…n)；

3.2.2对n个RVE有效性质进行数理统计

通过上述步骤得到n个RVE的有效性质，对其进行数理统计，有两种形式：

求张量矩阵的2-范数，再求其平均值和均方差，表示为：

式中，表示各有效张量的2-范数，为n个RVE有效张量的平均值， 是n个RVE有效张量的均方差。上述公式求得的值和即 可视为颗粒随机分布时均化模型下的有效常数。

直接求张量矩阵中各元素的均值和均方差，则其结果也为矩阵形式，表 示为：

$E\left(a\right)=\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{a}_i---\left(35\right)$

$D\left(a\right)=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\Sigma }}{(a_i-E(a\left)\right)}^2}---\left(36\right)$

式中，a表示各张量矩阵中的各元素，E(a)及D(a)分别表示其均值及均方差。

以此得到的均值及均方差矩阵与原张量矩阵形式相同。

4、数值模拟结果

在有限元方法的基础上，依据蒙特卡洛模拟法构建了一种求解非均质材料 在热弹性情况下宏观有效响应的随机均化模型。在线性位移边界条件下，设基 体的体积模量和剪切模量分别为k₁＝4，μ₁＝1；颗粒的体积模量和剪切模量分别 为k₁＝40，μ₁＝10。基体的热膨胀系数及热传导系数分别是α₀＝0.001，k_c＝0.1。Ogden 材料的各参数分别为：γ₁＝0.660,γ₂＝-0.231,γ₃＝0.050；β₁＝1.8,β₂＝-2.0,β₃＝7.0。各 参量间的相关系数取为0.8。

4.1有效量计算值

当颗粒体积分数V₂从0.05逐步递增至0.25时，对每一个V₂，均按照上述均 化模型进行n＝1000次模拟，求得材料的宏观有效应力张量有效热流量张量等参数。现将颗粒体积分数为0.25、考虑所有参数随机性时，考虑各参数间相 关性与不考虑其相关性两种情况下和的均值和均方差列于下表。

表一体积分数为0.25时，相关与不相关情况下和的均值与方差

4.2物理场

此外，还求得了RVE内部应力场、变形梯度场、热流量场等物理场的均值 μ与3σ边界，分布如图3(a)、(b)、(c)、图4(a)、(b)、(c)。

本发明的有益效果：

本发明针对有限变形下的随机非均质材料，采用新的随机分析方法构建随 机均化模型，发展了一种有效的热弹性化计算方法，求出了各随机宏观响应量 及其统计特征值，以达到对实验观察的有效补充和预测，为新型先进材料和结 构的优化设计提供客观、充分和真实的依据的目的。