CN103984869B - 一种预测复合材料热弹性有效属性和局部场的方法 - Google Patents

Abstract

本发明提供了一种预测复合材料热弹性有效属性和局部场的方法，以非均质、具有周期性分布的复合材料细观结构为研究对象，构建存储在所有单胞中的Helmholtz自由能之和，将热弹性分析求解问题转换为约束条件下能量方程的最小化‑取驻值问题，并通过对能量方程变分分析，得到求解波动函数的Euler‑Lagrange方程组和相应的非齐次边界条件，使用数值分析技术‑有限元法将能量方程改写为离散形式，求解得到单胞的Helmholtz自由能密度，并将其作为有效介质的本构模型应用于复合材料中，通过改变复合材料的荷载和温度条件，对复合材料进行局部场分析。本发明实用性强，通用性高，可显著提高此类问题的解算速度和效率。

Description

一种预测复合材料热弹性有效属性和局部场的方法

技术领域

本发明涉及材料力学性能分析领域，具体涉及一种预测复合材料热弹性有效属性和局部场的方法。

背景技术

许多复合材料都应用在在诸如电子封装、热保护系统等温度敏感的环境中，在这种环境下，热膨胀系数和热导率等热性系数的准确预测变得十分重要。复合材料的有效导热系数受成分性质、体积分数、组份分布和方向等诸多参数的影响，为此各国学者提出了各种细观力学模型来预测其有效导热系数。最简单的模型是基于Voigt和Reuss假设的混合规则，提供有效属性的上界和下界，但两界限之间差异较大时难以在一般的复合材料中加以应用。为此，研究人员提出了各种技术来减少界限间的差异，或是找到位于上界和下界之间的近似值。典型方法包括自洽法(Self-consistent Scheme-SCS)及其衍生物、变分法、三阶界限法、有效单元细胞法、递归细胞法、数学均匀化理论(Mathematical homogenizationtheories,MHT)、代表体元有限元法(Representative volume element,RVE)等。这些方法最初是分析弹性性能引入的，其中许多方法也扩展到预测热机性能。其中均匀化理论是用均质的宏观结构和非均质的具有周期性分布的细观结构描述复合材料结构：将力学量表示成关于宏观大集体才细观坐标的函数，并用细观和宏观两种尺度之比为小参数展开，用摄动技术将原问题转化为细观均匀化问题和宏观均匀化问题。

上述方法大多基于某些特定的假设，因此都有其特定的适用范围，有待进一步完善，有的模型计算结果有一定的偏差。此外，虽然建立起了复合材料细观量与宏观量(如位移、应力、应变、波函数等)的关系，但只能给出这些量变化的大致趋势和均匀化平滑结果，并不能描述这些量因组分变化而产生的局部涨落。要作为复合材料强度理论的判断依据，这显然是不够的。

发明内容

针对现有技术中存在的上述不足，本发明提供一种高效可靠的复合材料热弹性有效属性和局部场预测方法，解决现有技术存在的分析效率低、精度差不足，尤其是局部场变量分布无法精确预测之不足。

为解决上述技术问题，实现发明目的，本发明采用的技术方案如下：

一种预测复合材料热弹性有效属性和局部场的方法，包括如下步骤：

1)以非均质、具有周期性分布的复合材料为对象，构建存储在所有单胞中的Helmholtz自由能之和，并建立单胞界面上的位移场连续性条件；

2)根据最小势能原理，将热弹性分析求解问题转换为约束条件下能量方程的最小化-取驻值问题；

3)用位移场精确解表示平均值与含波动函数的差值之和，并通过对能量方程变分分析，得到求解Euler-Lagrange方程组和相应的非齐次边界条件；

4)使用数值分析技术-有限元法将能量方程改写为离散形式，求解得到单胞的Helmholtz自由能密度，从而得到复合材料热弹性的有效属性，并将其作为有效介质的本构模型应用于复合材料中，通过改变复合材料的荷载和温度条件，对复合材料进行局部场分析，得到复合材料热弹性有效属性和局部场数据结果并输出。

作为上述方案的进一步优化，所述步骤1具体为：从虚构的、未加载的完全由三维空间R和无穷多重复单胞组成的异质材料中推导均质化细观力学模型，该异质材料的广义总势能Π等于存储在所有单胞中的Helmholtz自由能之和，即：

$\Π=\underset{n=-\mathrm{\∞}}{\overset{\mathrm{\∞}}{\Σ}}\frac{1}{2}{\∫}_{\mathrm{\Ω}}(C_{ijkl}{\∈}_{ij}{\∈}_{kl}+2{\mathrm{\β}}_{ij}{\∈}_{ij}\mathrm{\θ}+c_v\frac{{\mathrm{\θ}}^2}{T_0})d\mathrm{\Ω};$

式中：Ω表示单胞的体积域，C_ijkl为四阶弹性张量分量，i＝1,2,3；j＝1,2,3；k＝1,2,3；l＝,1,2,3，β_ij为热应力系数二阶张量，c_v为体积不变时的单位体积比热容，T₀为无应力下的参考温度，θ表示实际温度和参考温度之间的差异，∈_ij,∈_kl均为线性理论中三维应变张量分量。

作为上述方案的进一步优化，所述步骤2具体为：根据最小势能原理，将热弹性分析求解问题转换为约束条件下能量方程的最小化-取驻值问题；

$J={\&Integral;}_R\left\{\begin{array}{c}<\frac{1}{2}{C}_{ijkl}{u}_{(i,|j)}{u}_{\left(k\right|l)}+{\mathrm{\&beta;}}_{ij}{u}_{\left(i\right|j)}\mathrm{\&theta;}+\frac{1}{2}{c}_v\frac{{\mathrm{\&theta;}}^2}{T_0}>\\ +{\mathrm{\&lambda;}}_i(<u_i>-v_i)+{\&Integral;}_{S_j}{\mathrm{\&gamma;}}_{ij}(u_i^{+j}-u_i^{-j}){\mathrm{dS}}_j\end{array}\right\}dR;$

式中：R为三维空间，u_i,u_j均为真实位移值，v_i为平均位移值，尖括号表示沿厚度方向的积分，u_(k|l)为u_(i|j)中i＝k,j＝l时的取值，λ_i,γ_ij为Lagrange乘子，d_i为单胞直径，S_j为单胞表面，<u_i>为u_i沿厚度方向的积分，

作为上述方案的进一步优化，所述步骤3具体为：用位移场精确解u_i表示平均值v_i与含波动函数的差值w_i之和，即：

u_i(x；y)＝v_i(x)+w_i(x；y)；

式中：x,y为直角坐标系，x＝(x₁,x₂,x₃),y＝(y₁,y₂,y₃)；

对能量方程变分分析，得到求解w_i的Euler-Lagrange方程组和相应的非齐次边界条件；

w_i(x；y)＝y_iv_i,j+χ_i(x；y)

式中：χ_i为波动函数；

作为上述方案的进一步优化，所述步骤4具体为：使用数值分析技术-有限元法将能量方程改写为离散形式：

${\&Pi;}_{\mathrm{\&Omega;}}=\frac{1}{2\mathrm{\&Omega;}}(V^{T}EV+2V^T{D}_{h\&Element;}\stackrel{\&OverBar;}{\&Element;}+{\stackrel{\&OverBar;}{\&Element;}}^T{D}_{\&Element;\&Element;}\stackrel{\&OverBar;}{\&Element;}+2V^T{D}_{h\mathrm{\&theta;}}\mathrm{\&theta;}+2{\stackrel{\&OverBar;}{\&Element;}}^T{D}_{\&Element;\mathrm{\&theta;}}\mathrm{\&theta;}+D_{\mathrm{\&theta;}\mathrm{\&theta;}}\frac{{\mathrm{\&theta;}}^2}{T_0})$

式中：

E＝∫_Ω(Γ_hS)^TD(Γ_hS)dΩ,D_h∈＝∫_Ω(Γ_hS)^TDdΩ,D_∈∈＝∫_ΩDdΩ,

D_hθ＝∫_Ω(Γ_hS)^TβdΩ,D_∈θ＝∫_ΩβdΩ,D_θθ＝∫_Ωc_vdΩ

其中：D为6×6阶含四阶弹性张量C_ijkl的材料矩阵，β为6×1阶含β_ij列阵；Γ_h为算子矩阵；

单胞的Helmholtz自由能密度为：

${\&Pi;}_{\mathrm{\&Omega;}}=\frac{1}{2}{\stackrel{\&OverBar;}{\&Element;}}^T\stackrel{\&OverBar;}{D}\stackrel{\&OverBar;}{\&Element;}+{\stackrel{\&OverBar;}{\&Element;}}^T\stackrel{\&OverBar;}{\mathrm{\&beta;}}\mathrm{\&theta;}+\frac{1}{2}{\stackrel{\&OverBar;}{c}}_v\frac{{\mathrm{\&theta;}}^2}{T_0};$

式中：为有效弹性材料矩阵，为有效热应力系数列阵，为有效比热容，含全局应变：

有效热膨胀系数为：

$\stackrel{\&OverBar;}{\mathrm{\&alpha;}}=-{\stackrel{\&OverBar;}{D}}^{-1}\stackrel{\&OverBar;}{\mathrm{\&beta;}};$

求解得到单胞的Helmholtz自由能密度，根据宏观性能重构单胞内的局部场：

$u=v+\left[\begin{array}{ccc}{v}_{1,1}& v_{1,2}& v_{1,3}\\ v_{2,1}& v_{2,2}& v_{2,3}\\ v_{3,1}& v_{3,2}& v_{3,3}\end{array}\right]\left\{\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right\}+\stackrel{\&OverBar;}{S}\stackrel{\&OverBar;}{V}$

式中：u,v分别为u_i,v_i的列阵；为重构后的形函数，为修正后波动函数的节点值列阵；

重构的局部应变场为：

使用组成材料的三维本构关系重构局部应力场：σ＝D∈+βθ；

将作为复合材料热弹性有效属性数据结果输出，将σ作为局部场数据结果输出；V为所有活动节点波动函数的节点值列阵，Ω表示单胞的体积域，S表示形函数。相比于现有技术，本发明具有如下优点：

1、本发明提供的预测复合材料热弹性有效属性和局部场的方法，采用变分渐近法作为数学基础，仅含细观力学的最基本假设，不需任何周期性波动函数和边界条件假设，具有数学上的严密性。

2、本发明提供的预测复合材料热弹性有效属性和局部场的方法，采用非强制性边界条件，仅需一次求解即可获得不同方向上不同材料的属性，相对于在不同温度条件下重复运行程序更简便、高效、快捷。

4、本发明提供的预测复合材料热弹性有效属性和局部场的方法，计算的有效材料属性和局部场分布与波动函数精度相一致，无需引入更多近似关系(如平均温度梯度和热通量)的后处理计算。

5、本发明提供的预测复合材料热弹性有效属性和局部场的方法，可模拟完全各向异性复合材料，突破了有限元只能处理宏观正交各向异性材料的局限。

6、本发明提供的预测复合材料热弹性有效属性和局部场的方法，实用性强，通用性高，可显著提高此类问题的解算速度和效率。

附图说明

图1为本发明提供的非均匀复合材料的单胞结构的有效属性分析方法的流程图。

图2为本发明研究的非均匀复合材料单胞结构分析示意图。

图3为本发明研究的周期性非均匀材料和相关的单胞。

图4为本发明具体实施例中硼/铝复合材料轴向有效热膨胀系数随纤维体积率变化图。

图5为本发明具体实施例中硼/铝复合材料横向有效热膨胀系数随纤维体积率变化图。

图6为本发明具体实施例中钢/铝颗粒增强复合材料有效比热容随颗粒体积分数变化曲线。

图7为本发明具体实施例中硼/铝复合材料在温度上升100K下单胞内σ₂₂分布云图。

图8为本发明具体实施例中硼/铝复合材料在温度上升100K下单胞内σ₂₃分布云图。

具体实施方式

1、理论公式。

本发明非均匀复合材料的单胞结构的有效属性分析方法的流程图如图1所示，本发明基于如下两个基本假设：

假设1：根据对振荡函数的一般讨论可知，宏观量必定是细观相应量的体积平均值，如位移u_i的体积平均v_i定义为：

$v_i=\frac{1}{\mathrm{\&Omega;}}{\&Integral;}_{\mathrm{\&Omega;}}{u}_{i}d\mathrm{\&Omega;}=<u_i>---\left(1\right)$

式中：Ω表示单胞的体积域，尖括号表示沿厚度方向的积分，〈u_i>为u_i沿厚度方向的积分。

假设2：细观分析得到的有效材料性能与宏观结构的约束条件及加载条件无关，它仅与材料细部几何、物理参数的分布形态、大小有关。这说明材料的有效性能是材料本身固有的特性，这符合物理学的直观。

值得注意的是，这些假设并不是强制性的。假设1的数学意义是场变量的精确解是单胞体积域上的积分；假设2意味着在宏观分析中可忽略材料特性的尺寸效应，这也是常规连续介质力学分析中所做的假设。单胞的细观力学分析只有在h/l＜＜1时有效(h为单胞的特征尺度，l为宏观尺度)。

本发明的非均匀复合材料单胞结构分析示意图如图2所示，非均质、具有周期性分布的复合材料和相关单胞如图3所示，设其在三维空间R中占据区域为Ω，其细观结构可看成是非均质单胞在空间的周期性重复排列。理论推导中使用三种坐标系统：两个直角坐标x＝(x₁,x₂,x₃),y＝(y₁,y₂,y₃)和一个整数坐标n＝(n₁,n₂,n₃)。直角坐标x为笛卡尔坐标，其原点设置在区域Ω内，选择x_i作为全局坐标以指示该质点在哪一个单胞内，x₁为横坐标，x₂为竖坐标，x₃为纵坐标。选择与x_i平行的y_i作为局部坐标指示在该单胞内的何处(下标i＝1,2,3，下标重复时表示累加)。选择局部坐标y_i的原点作为单胞的几何中心。若单胞的直径为d_i，则y_i∈[-d_i/2,d_i/2]。为唯一确定单胞在异质材料中的位置，引入整数坐标n_i，整数坐标与全局坐标的关系是n_i＝x_i/d_i。

如假设2所暗示的那样，可从虚构的、无界的、未加载的复合材料结构中得到具有相同有效属性的真实的、有界的、加载的细观结构。同样，可从虚构的、未加载的完全由三维空间R和无穷多重复单胞组成的异质材料中推导均质化细观力学模型。这种假想材料的广义总势能Π等于存储在所有单胞中的Helmholtz自由能之和，即：

$\&Pi;=\underset{n=-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{1}{2}{\&Integral;}_{\mathrm{\&Omega;}}(C_{ijkl}{\&Element;}_{ij}{\&Element;}_{kl}+2{\mathrm{\&beta;}}_{ij}{\&Element;}_{ij}\mathrm{\&theta;}+c_v\frac{{\mathrm{\&theta;}}^2}{T_0})d\mathrm{\&Omega;}---\left(2\right)$

式中：C_ijkl为四阶弹性张量分量，j＝1,2,3；k＝1,2,3；l＝,1,2,3。β_ij为热应力系数二阶张量，c_v为体积不变时的单位体积比热容，T₀为无应力下的参考温度，θ表示实际温度和参考温度之间的差异，∈_ij,∈_kl均为线性理论中三维应变张量分量。

${\&Element;}_{ij}(n;y)=\frac{1}{2}\&lsqb;u_{i;j}(n;y)+u_{j;i}(n;y)\&rsqb;---\left(3\right)$

式中：u_j为u_i中i＝j的值。

从无穷多个单胞形成一个连续的异质材料来看，相邻单胞界面上的位移场需满足连续性条件。整数坐标系(n₁,n₂,n₃)下的单胞需满足如下条件：

$\begin{array}{c}{u}_i(n_1,n_2,n_3;d_1/2,y_2,y_3)=u_i(n_1+1,n_2,n_3;-d_1/2,y_2,y_3),\\ u_i(n_1,n_2,n_3;y_1,d_2/2,y_3)=u_i(n_1,n_2+1,n_3;,y_1,-d_2/2,y_3),\\ u_i(n_1,n_2,n_3;y_1,y_2,d_3/2)=u_i(n_1,n_2,n_3+1;y_1,y_2,-d_3/2)\end{array}---\left(4\right)$

为获得复合材料有效热膨胀系数和比热容，可假设θ相对于时间和空间坐标是不变的，因此自动满足相邻单胞间温度场的连续性条件。根据最小势能原理，精确求解问题转换为式(1)、(4)约束下能量方程(2)的最小化问题。为避免与离散整数参数n_i相关的困难，可将式(2)、(3)、(4)用准连续介质中的连续函数改写为：

$\&Pi;=\frac{1}{2}{\&Integral;}_R<C_{ijkl}{\&Element;}_{ij}{\&Element;}_{kl}+2{\mathrm{\&beta;}}_{ij}{\&Element;}_{ij}\mathrm{\&theta;}+c_v\frac{{\mathrm{\&theta;}}^2}{T_0}>dR---\left(5\right)$

${\&Element;}_{ij}(x;y)=\frac{1}{2}\&lsqb;u_{i;j}(x;y)+u_{j;i}(x;y)\&rsqb;\&equiv;u_{\left(i\right|j)}---\left(6\right)$

$\begin{array}{c}{u}_i(x_1,x_2,x_3;d_1/2,y_2,y_3)=u_i(x_1+d_1,x_2,x_3;-d_1/2,y_2,y_3),\\ u_i(x_1,x_2,x_3;y_1,d_2/2,y_3)=u_i(x_1,x_2+d_2,x_3;,y_1,-d_2/2,y_3),\\ u_i(x_1,x_2,x_3;y_1,y_2,d_3/2)=u_i(x_1,x_2,x_3+d_3;y_1,y_2,-d_3/2)\end{array}---\left(7\right)$

式中，

使用Lagrange乘子技术，可将热弹性分析转换为求如下能量方程J的驻值点(最小化)问题：

$J={\&Integral;}_R\left\{\begin{array}{c}<\frac{1}{2}{C}_{ijkl}{u}_{(i,|j)}{u}_{\left(k\right|l)}+{\mathrm{\&beta;}}_{ij}{u}_{\left(i\right|j)}\mathrm{\&theta;}+\frac{1}{2}{c}_v\frac{{\mathrm{\&theta;}}^2}{T_0}>\\ +{\mathrm{\&lambda;}}_i(<u_i>-v_i)+{\&Integral;}_{S_j}{\mathrm{\&gamma;}}_{ij}(u_i^{+j}-u_i^{-j}){\mathrm{dS}}_j\end{array}\right\}dR---\left(8\right)$

式中：R为三维空间，u_(kl)为u_(i|j)中i＝k,j＝l时的取值， λ_i,γ_ij为将式(1)和(7)的约束强加于式(8)引入的Lagrange乘子；S_i为n_i＝1的表面，S_j为S_i中i＝j的取值。

细观力学的主要目标是由v_i找到真实的位移u_i，为此需求解式(8)在全局坐标系x下各点的驻值问题，这是该问题的求解十分困难。为解决这一难题，就必须找出某种无需将每个微结构都表示出来的等效材料模型，这一模型应当即能表征材料的平均力学性能，又能表征材料的非均质性，为此可将u_i精确解表示为平均值v_i和差值之和，即：

u_i(x；y)＝v_i(x)+w_i(x；y) (9)

根据式(1)可知，式(9)中的<w_i>＝0。由非均质材料同质化分析可知，w_i应渐近小于v_i，即w_i～(h/l)v_i。将式(9)代入式(8)，并使用式(6)可得到能量方程的主导项J₁为：

$J_1={\&Integral;}_R\left\{\begin{array}{c}<\frac{1}{2}{C}_{ijkl}{w}_{\left(i\right|j)}{w}_{\left(k\right|l)}+{\mathrm{\&beta;}}_{ij}{w}_{\left(i\right|j)}\mathrm{\&theta;}+\frac{1}{2}{c}_v\frac{{\mathrm{\&theta;}}^2}{T_0}>+{\mathrm{\&lambda;}}_i<w_i>\\ +{\&Integral;}_{S_j}{\mathrm{\&gamma;}}_{ij}(w_i^{+j}-w_i^{-j}-v_{i,j}{d}_j){\mathrm{ds}}_j\end{array}\right\}dR---\left(10\right)$

式中：

对能量方程J₁进行变分分析可得到求解w_i的Euler-Lagrange方程组和相应的非齐次边界条件。为方便推导，可通过变量变化使边界条件齐次化。考虑式(10)的最后一项，可用如下变量变化表征：

w_i(x；y)＝y_iv_i,j+χ_i(x；y) (11)

式中：χ_i为波动函数，

由于选择局部坐标系的原点为单胞的中点，这意味着对χ_i存在以下约束：

<χ_i>＝0 (12)

根据变分渐近法，将式(12)代入式(10)，得到定义在单胞上的求解χ驻值问题为：

$\begin{array}{c}{J}_{\mathrm{\&Omega;}}=\frac{1}{2}<C_{ijkl}\&lsqb;{\stackrel{\&OverBar;}{\&Element;}}_{ij}+{\mathrm{\&chi;}}_{\left(i\right|j)}\&rsqb;\&lsqb;{\stackrel{\&OverBar;}{\&Element;}}_{kl}+{\mathrm{\&chi;}}_{\left(k\right|l)}\&rsqb;>+\frac{1}{2}<c_v\frac{{\mathrm{\&theta;}}^2}{T_0}>\\ +<{\mathrm{\&beta;}}_{ij}\&lsqb;{\stackrel{\&OverBar;}{\&Element;}}_{ij}+{\mathrm{\&chi;}}_{\left(i\right|j)}\&rsqb;\mathrm{\&theta;}>+{\mathrm{\&lambda;}}_i<w_i>+{\&Integral;}_{S_j}{\mathrm{\&gamma;}}_{ij}({\mathrm{\&chi;}}_i^{+j}-{\mathrm{\&chi;}}_i^{-j}){\mathrm{dS}}_j\end{array}---\left(13\right)$

式中：为均质化材料的全局应变张量分量；中i＝k,j＝l时的取值，

式(13)的能量方程J_Ω包含热弹性细观力学模型的所有信息。对于二元复合材料的简单驻值问题，可通过变分解析求解。但对更一般情况，需使用数值方法(如有限元法)进行数值求解。

2、细观问题的有限元求解。

Lagrange乘子的引入会增加未知数的个数，直接求解式(13)的有限元解并不是最方便和最有效的方式。可将式(13)的变分表达式改写为如下能量方程Π_Ω的驻值问题：

${\&Pi;}_{\mathrm{\&Omega;}}=\frac{1}{2\mathrm{\&Omega;}}{\&Integral;}_{\mathrm{\&Omega;}}\left\{\begin{array}{c}{C}_{ijkl}\&lsqb;{\stackrel{\&OverBar;}{\&Element;}}_{ij}+{\mathrm{\&chi;}}_{\left(i\right|j)}\&rsqb;\&lsqb;{\stackrel{\&OverBar;}{\&Element;}}_{kl}+{\mathrm{\&chi;}}_{\left(k\right|l)}\&rsqb;+c_v\frac{{\mathrm{\&theta;}}^2}{T_0}\\ 2{\mathrm{\&beta;}}_{ij}\&lsqb;{\stackrel{\&OverBar;}{\&Element;}}_{ij}+{\mathrm{\&chi;}}_{\left(i\right|j)}\&rsqb;\mathrm{\&theta;}\end{array}\right\}d\mathrm{\&Omega;}---\left(14\right)$

并受如下6个约束：

${\mathrm{\&chi;}}_i^{+j}={\mathrm{\&chi;}}_i^{-j},(i,j=1,2,3)---\left(15\right)$

值得注意的是：式(15)代表了著名的周期性边界条件，以往模型通常是预先假设该周期性边界条件；而在目前的模型中，是来自式(13)变分推导。复合材料在y₁-y₂平面均匀，沿y₃不均匀，χ_i仅是y₃的函数。除χ_i,3外，波动函数的其它偏导χ_i,j会消失，因此式(12)中的约束不会影响Π_Ω的最小值求解，反而可帮助唯一确定χ_i。在实际操作中，可先约束任意节点上的波动函数为零，然后再使用该约束重构特定的波动函数。为了在有限元求解过程中实现式(15)约束，一种有效方法是在组装后的总刚中引入罚函数，但这样的处理会引入额外的近似，且解的鲁棒性取决于罚数的数量，会大量耗费计算机的内存，降低计算效率。这里选择正边界表面上的节点(y_i＝d_i/2)从属于负边界表面上的节点(y_i＝-d_i/2)。通过整合全部独立活动自由度，可隐性和精确地包含式(15)中的约束，并可减少线性系统内未知量的数目。

为推导离散形式的能量方程，引入如下矩阵符号：

$\stackrel{\&OverBar;}{\&Element;}={\left[\begin{array}{cccccc}{\stackrel{\&OverBar;}{\&Element;}}_{11}& 2{\stackrel{\&OverBar;}{\&Element;}}_{12}& {\stackrel{\&OverBar;}{\&Element;}}_{22}& 2{\stackrel{\&OverBar;}{\&Element;}}_{13}& 2{\stackrel{\&OverBar;}{\&Element;}}_{23}& {\stackrel{\&OverBar;}{\&Element;}}_{33}\end{array}\right]}^T---\left(16\right)$

${\&Element;}_1=\left\{\begin{array}{c}{\mathrm{\&chi;}}_{1;1}\\ {\mathrm{\&chi;}}_{1;2}+{\mathrm{\&chi;}}_{2;1}\\ {\mathrm{\&chi;}}_{2;2}\\ {\mathrm{\&chi;}}_{1;3}+{\mathrm{\&chi;}}_{3;1}\\ {\mathrm{\&chi;}}_{2;3}+{\mathrm{\&chi;}}_{3;2}\\ {\mathrm{\&chi;}}_{3;3}\end{array}\right\}=\left\{\begin{array}{ccc}{\left(\right)}_{;1}& 0& 0\\ {\left(\right)}_{;2}& {\left(\right)}_{;1}& 0\\ 0& {\left(\right)}_{;2}& 0\\ {\left(\right)}_{;3}& 0& {\left(\right)}_{;1}\\ 0& {\left(\right)}_{;3}& {\left(\right)}_{;2}\\ 0& 0& {\left(\right)}_{;3}\end{array}\right\}\left\{\begin{array}{c}{\mathrm{\&chi;}}_1\\ {\mathrm{\&chi;}}_2\\ {\mathrm{\&chi;}}_3\end{array}\right\}={\mathrm{\&Gamma;}}_h\mathrm{\&chi;}---\left(17\right)$

式中：Γ_h为算子矩阵，χ为含三个波动函数分量的列阵。

使用有限元将波动函数χ离散为：

χ(x_i；y_i)＝S(y_i)V(x_i) (18)

式中：S表示形函数(不含约束节点和从节点)，V为所有活动节点波动函数的节点值列阵。

将式(16)、(17)和(18)代入式(14)，得到能量方程的离散形式为：

${\&Pi;}_{\mathrm{\&Omega;}}=\frac{1}{2\mathrm{\&Omega;}}(V^{T}EV+2V^T{D}_{h\&Element;}\stackrel{\&OverBar;}{\&Element;}+{\stackrel{\&OverBar;}{\&Element;}}^T{D}_{\&Element;\&Element;}\stackrel{\&OverBar;}{\&Element;}+2U^T{D}_{h\mathrm{\&theta;}}\mathrm{\&theta;}+2{\stackrel{\&OverBar;}{\&Element;}}^T{D}_{\&Element;\mathrm{\&theta;}}\mathrm{\&theta;}+D_{\mathrm{\&theta;}\mathrm{\&theta;}}\frac{{\mathrm{\&theta;}}^2}{T_0})---\left(19\right)$

式中：

E＝∫_Ω(Γ_hS)^TD(Γ_hS)dΩ,D_h∈＝∫_Ω(Γ_hS)^TDdΩ

D_∈∈＝∫_ΩDdΩ,D_hθ＝∫_Ω(Γ_hS)^TβdΩ

D_∈θ＝∫_ΩβdΩ,D_θθ＝∫_Ωc_vdΩ

其中：D为6×6阶含四阶弹性张量C_ijkl的材料矩阵，β为6×1阶含β_ij列阵。

最小化式(19)的Π_Ω，得到如下线性系统：

$EV=-D_{h\&Element;}\stackrel{\&OverBar;}{\&Element;}-D_{h\mathrm{\&theta;}}\mathrm{\&theta;}---\left(20\right)$

由式(20)可知，波动函数节点值列阵V与θ呈线性正比关系，这意味着解可象征性地表示为：

$V=V_0\stackrel{\&OverBar;}{\&Element;}+V_{\mathrm{\&theta;}}\mathrm{\&theta;}---\left(21\right)$

式中，V₀为应变节点值列阵；V_θ为温差节点值列阵。

将式(21)代入式(19)，可得到单胞的Helmholtz自由能密度为：

${\&Pi;}_{\mathrm{\&Omega;}}=\frac{1}{2}{\stackrel{\&OverBar;}{\&Element;}}^T\stackrel{\&OverBar;}{D}\stackrel{\&OverBar;}{\&Element;}+{\stackrel{\&OverBar;}{\&Element;}}^T\stackrel{\&OverBar;}{\mathrm{\&beta;}}\mathrm{\&theta;}+\frac{1}{2}{\stackrel{\&OverBar;}{c}}_v\frac{{\mathrm{\&theta;}}^2}{T_0}---\left(22\right)$

式中：为有效弹性材料矩阵，为有效热应力系数列阵，为有效比热容，含全局应变，从而得到复合材料热弹性的有效属性：

$\stackrel{\&OverBar;}{D}=\frac{1}{\mathrm{\&Omega;}}(V_0^T{D}_{h\&Element;}+D_{\&Element;\&Element;}),\stackrel{\&OverBar;}{\mathrm{\&beta;}}=\frac{1}{\mathrm{\&Omega;}}\&lsqb;\frac{1}{2}(D_{h\&Element;}^T{V}_{\mathrm{\&theta;}}+V_0^T{D}_{h\mathrm{\&theta;}})+D_{\&Element;\mathrm{\&theta;}}\&rsqb;,{\stackrel{\&OverBar;}{c}}_v=\frac{1}{\mathrm{\&Omega;}}\&lsqb;V_{\mathrm{\&theta;}}^T{D}_{h\mathrm{\&theta;}}{T}_0+D_{\mathrm{\&theta;}\mathrm{\&theta;}}\&rsqb;---\left(23\right)$

有效热膨胀系数为：

$\stackrel{\&OverBar;}{\mathrm{\&alpha;}}=-{\stackrel{\&OverBar;}{D}}^{-1}\stackrel{\&OverBar;}{\mathrm{\&beta;}}---\left(24\right)$

若需求解单胞内的局部场，可根据宏观性能(如全局位移v_i，全局温度分布θ)重构局部场。要重构局部位移场，需根据相应的活动节点指定从属节点值，约束节点值指定为零，从而由V₀,V_θ构建两个新阵列显然，新构建的对应于式(15)约束下式(14)Π_Ω的驻值，可能并不满足式(12)中的约束。为得到满足式(12)的真实解，为V₀的真实解，为V_θ的真实解，可在每个节点上加上一个常数(等同于波动函数的平均值)。修正后的波动函数的节点值真实解为：

$\stackrel{\&OverBar;}{V}={\stackrel{\&OverBar;}{V}}_0\stackrel{\&OverBar;}{\&Element;}+{\stackrel{\&OverBar;}{V}}_{\mathrm{\&theta;}}\mathrm{\&theta;}---\left(25\right)$

式中：为修正后的应变节点值列阵；V_θ为修正后的温差节点值列阵。

唯一确定波动函数后，由式(9)、(11)和(25)重构局部位移场为：

$u=v+\left[\begin{array}{ccc}{v}_{1,1}& v_{1,2}& v_{1,3}\\ v_{2,1}& v_{2,2}& v_{2,3}\\ v_{3,1}& v_{3,2}& v_{3,3}\end{array}\right]\left\{\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right\}+\stackrel{\&OverBar;}{S}\stackrel{\&OverBar;}{V}---\left(26\right)$

式中：u,v分别为u_i,v_i的列阵。由于对从属节点和约束节点进行了重构，为重构后的形函数，不同于S。

局部应变场∈可由式(6)、(9)、(11)、(17)和(25)重构为：

$\&Element;=\stackrel{\&OverBar;}{\&Element;}+{\mathrm{\&Gamma;}}_h\stackrel{\&OverBar;}{S}\stackrel{\&OverBar;}{V}---\left(27\right)$

可直接使用组成材料的三维本构关系重构局部应力场σ：

σ＝D∈+βθ (28)

最后，通过改变复合材料的荷载和温度条件，对复合材料进行分析，得到复合材料热弹性有效属性和局部场数据结果并输出。其输出得到的数据结果可以作为后续环节的输入数据或输入到计算机中进行处理。将作为复合材料热弹性有效属性数据结果输出，将σ作为局部场数据结果输出；V为所有活动节点波动函数的节点值列阵，Ω表示单胞的体积域，S表示形函数。

为了验证上述方法预测复合材料有效热膨胀系数的准确性，考虑各向同性硼/铝构成的复合材料，其材料属性如表1所示。

表1硼/铝材料属性

表2列出了各细观力学模型预测的硼/铝复合材料有效热膨胀系数。由表中数据可看出：本发明解可与精确解吻合到小数点后二位；而GMC(Generalized method of cell)和Tamma结果与精确解相比并不可靠。

表2硼/铝复合材料的有效热膨胀系数

图4、5中绘出了不同纤维体积率下有效热膨胀系数。正如所期望那样，轴向热膨胀系数和横向热膨胀系数随纤维体积的减少而增加，同时本发明解与精确解吻合性很好。

为了验证上述方法预测复合材料有效比热容的准确性，考虑宏观各向同性的钢/铝颗粒增强复合材料，各成分材料属性如表3所示：

表3钢/铝复合材料属性

不同颗粒体积分数下的本发明预测的有效比热容绘于图6，并与“精确解”进行对比。由图中可看出：随着颗粒体积含量的增大，有效比热容也逐渐增大，且呈线性递增关系；本发明预测结果与精确解吻合很好。

为了验证上述方法预测复合材料因宏观温度改变产生的局部应力场的准确性，考虑考虑硼/铝复合材料的纤维体积分数为0.2，假设单胞处于无应力状态，温度增加100K。由于各成分的热膨胀系数不相同，单胞内将产生热应力。图7、8分别绘出了单胞内热应力σ₂₂,σ₂₃的分布。由图中可看出：本发明解能准确捕捉所有沿纤维-基体界面应力分布的突然变化。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (1)

1.一种预测复合材料热弹性有效属性和局部场的方法，其特征在于，包括如下步骤：

1)以非均质、具有周期性分布的复合材料为对象，构建存储在所有单胞中的Helmholtz自由能之和，并建立单胞界面上的位移场连续性条件；

2)根据最小势能原理，将热弹性分析求解问题转换为约束条件下能量方程的最小化-取驻值问题；

3)用位移场精确解表示平均值与含波动函数的差值之和，并通过对能量方程变分分析，得到求解Euler-Lagrange方程组和相应的非齐次边界条件；

4)使用数值分析技术-有限元法将能量方程改写为离散形式，求解得到单胞的Helmholtz自由能密度，从而得到复合材料热弹性的有效属性，并将其作为有效介质的本构模型应用于复合材料中，通过改变复合材料的荷载和温度条件，对复合材料进行分析，得到复合材料热弹性有效属性和局部场数据结果并输出；

所述步骤1具体为：从虚构的、未加载的完全由三维空间R和无穷多重复单胞组成的异质材料中推导均质化细观力学模型，该异质材料的广义总势能Π等于存储在所有单胞中的Helmholtz自由能之和，即：

$\&Pi;=\underset{n=-\mathrm{\&infin;}}{\overset{\mathrm{\&infin;}}{\&Sigma;}}\frac{1}{2}{\&Integral;}_{\mathrm{\&Omega;}}(C_{ijkl}{\&Element;}_{ij}{\&Element;}_{kl}+2{\mathrm{\&beta;}}_{ij}{\&Element;}_{ij}\mathrm{\&theta;}+c_v\frac{{\mathrm{\&theta;}}^2}{T_0})d\mathrm{\&Omega;};$

式中：Ω表示单胞的体积域，C_ijkl为四阶弹性张量分量，i＝1,2,3；j＝1,2,3；k＝1,2,3；l＝,1,2,3，β_ij为热应力系数二阶张量，c_v为体积不变时的单位体积比热容，T₀为无应力下的参考温度，θ表示实际温度和参考温度之间的差异，∈_ij,∈_kl均为线性理论中三维应变张量分量；

所述步骤2具体为：根据最小势能原理，将热弹性分析求解问题转换为约束条件下能量方程的最小化-取驻值问题；

式中：R为三维空间，u_i,u_j均为真实位移值，v_i为平均位移值，尖括号表示沿厚度方向的积分，u_(k|l)为u_(i|j)中i＝k,j＝l时的取值，λ_i,γ_ij为Lagrange乘子，d_i为单胞直径，S_j为单胞表面，<u_i>为u_i沿厚度方向的积分，

所述步骤3具体为：用位移场精确解u_i表示平均值v_i与含波动函数的差值w_i之和，即：

u_i(x；y)＝v_i(x)+w_i(x；y)；

式中：x,y为直角坐标系，x＝(x₁,x₂,x₃),y＝(y₁,y₂,y₃)；

对能量方程变分分析，得到求解w_i的Euler-Lagrange方程组和相应的非齐次边界条件；

w_i(x；y)＝y_iv_i,j+χ_i(x；y)

式中：χ_i为波动函数；

所述步骤4具体为：使用数值分析技术-有限元法将能量方程改写为离散形式：

${\&Pi;}_{\mathrm{\&Omega;}}=\frac{1}{2\mathrm{\&Omega;}}(V^{T}EV+2V^T{D}_{h\&Element;}\stackrel{\&OverBar;}{\&Element;}+{\stackrel{\&OverBar;}{\&Element;}}^T{D}_{\&Element;\&Element;}\stackrel{\&OverBar;}{\&Element;}+2V^T{D}_{h\mathrm{\&theta;}}\mathrm{\&theta;}+2{\stackrel{\&OverBar;}{\&Element;}}^T{D}_{\&Element;\mathrm{\&theta;}}\mathrm{\&theta;}+D_{\mathrm{\&theta;}\mathrm{\&theta;}}\frac{{\mathrm{\&theta;}}^2}{T_0})$

式中：

E＝∫_Ω(Γ_hS)^TD(Γ_hS)dΩ,D_h∈＝∫_Ω(Γ_hS)^TDdΩ,D_∈∈＝∫_ΩDdΩ,

D_hθ＝∫_Ω(Γ_hS)^TβdΩ,D∈_θ＝∫_ΩβdΩ,D_θθ＝∫_Ωc_vdΩ

其中：D为6×6阶含四阶弹性张量C_ijkl的材料矩阵，β为6×1阶含β_ij列阵；Γ_h为算子矩阵；

单胞的Helmholtz自由能密度为：

${\&Pi;}_{\mathrm{\&Omega;}}=\frac{1}{2}{\stackrel{\&OverBar;}{\&Element;}}^T\stackrel{\&OverBar;}{D}\stackrel{\&OverBar;}{\&Element;}+{\stackrel{\&OverBar;}{\&Element;}}^T\stackrel{\&OverBar;}{\mathrm{\&beta;}}\mathrm{\&theta;}+\frac{1}{2}{\stackrel{\&OverBar;}{c}}_c\frac{{\mathrm{\&theta;}}^2}{T_0};$

式中：为有效弹性材料矩阵，为有效热应力系数列阵，为有效比热容，含全局应变：

有效热膨胀系数为：

$\stackrel{\&OverBar;}{\mathrm{\&alpha;}}=-{\stackrel{\&OverBar;}{D}}^{-1}\stackrel{\&OverBar;}{\mathrm{\&beta;}};$

求解得到单胞的Helmholtz自由能密度，根据宏观性能重构单胞内的局部场：

$u=v+\left[\begin{array}{ccc}{v}_{1,1}& v_{1,2}& v_{1,3}\\ v_{2,1}& v_{2,2}& v_{2,3}\\ v_{3,1}& v_{3,2}& v_{3,3}\end{array}\right]\left\{\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right\}+\stackrel{\&OverBar;}{S}\stackrel{\&OverBar;}{V}$

式中：u,v分别为u_i,v_i的列阵；为重构后的形函数，为修正后波动函数的节点值列阵；

重构的局部应变场为：

使用组成材料的三维本构关系重构局部应力场：σ＝D∈+βθ；

将作为复合材料热弹性有效属性数据结果输出，将σ作为局部场数据结果输出；V为所有活动节点波动函数的节点值列阵，Ω表示单胞的体积域，S表示形函数。