CN103604729B - 一种颗粒随机分布的复合材料宏观有效性质的预测方法 - Google Patents

Abstract

本发明公开了一种颗粒随机分布的复合材料宏观有效性质的预测方法，包括：1）基于随机序列添加RSA方法使用MATLAB构建颗粒随机分布的复合材料表征体积单元RVE数值模型；2）针对传统RSA方法生成的复合材料RVE数值模型，消除当颗粒体积分数较大时的颗粒重叠现象；3）对颗粒随机分布的复合材料RVE数值模型进行有限元分析计算，求得RVE数值模型有效性质的数值解；4）建立随机均化模型求解复合材料的宏观有效性质，并以该复合材料的宏观有效性质作为其真实有效性质，为复合材料的宏观有效性质提供了更为可靠的预测结果，为新型先进材料的使用和结构的优化设计提供充分的依据。

Description

一种颗粒随机分布的复合材料宏观有效性质的预测方法

技术领域

本发明涉及复合材料细观力学，具体地说是针对一种颗粒随机分布的复合材料，构建数值模型对其进行表征并对其宏观有效性质进行预报。

背景技术

所谓复合材料，是指由不同材料组成的或由不同状态的同一材料组成的复合体。通常，复合材料中有许多相（相：指材料的不同组成成分。如混凝土可以认为是两相材料：一相是水泥，称为基体；一相是沙子，称为夹杂）。基体相一般是连续的，占材料的体积比较大；而其它相是离散的，分布于基体材料中，所占体积分数较小。复合材料物理性质的变化遍及整个微观结构，不仅其力学参数随空间位置的变化而改变，而且其细观结构的特性会对材料的宏观力学性质产生决定性的影响。这里的宏观性质也称为有效性质，是指能够在宏观尺度上实验测量的材料性能，它包括有效弹性常数、有效传导率以及有效耦合系数等。

由于复合材料的宏观长度尺寸显著地大于微尺度下不同组分（即相）的长度尺度，若对所有微尺度组分均采用显式求解，无论是解析法还是数值解法，都是极具挑战性且代价昂贵的工作。一种显著减少分析此类问题所需成本的方法就是用有效均质材料替代原先的复合材料并通过构建反映宏观性质的有效本构方程进而获得材料的宏观有效性质，即均匀化方法。该方法本身取决于从材料中取出的一个表征体积单元(Representative Volume Element，RVE)，使用有限元方法(Finite Element Method，FEM)进行分析计算。因其可以极大地减少复合设计参数的数量，且随着工程需求的变化，该方法得到不断发展。

然而，传统的均匀化方法存在着一些不足之处。首先，传统的均匀化方法都假设复合材料的夹杂具有相同的大小和形状，且假设细观结构呈均匀性和周期性分布。而事实上，复合材料的细观结构往往非常复杂，夹杂的分布一般是随机的，夹杂的大小和形状有一定的概率离散性。将细观结构视为均匀性和周期性后，通过均匀化方法获得的材料宏观有效性质具有局限性，这直接影响了材料在结构工程中的应用。其次，传统的均匀化方法大多采用随机序列添加(Random sequential Addition，RSA)方法[Jia,X.and Williams,R.A.A packingalgorithm for particles of arbitrary shapes.Powder Technology,120:175-786,2001]来生成用于均化分析的RVE。这种方法有一定的缺陷，如以生成圆形颗粒夹杂为例，理论上这种方法所能达到的颗粒最大体积分数为v_particle=π/4=0.785，并不能模拟v_particle超过π/4时的情形。这是因为，随着颗粒体积分数的增加，颗粒受RVE边界约束而发成了重叠现象，从而导致均化方法分析得到的有效性质很快收敛失真。这就给均化过程带来了困扰。

因此，如何有效生成均化技术中重要的一环——微观结构RVE，使得建立于这个微观结构之上的均化有效性质能更好的反映复合材料的真实性质，成为了迫切的要求，这正是本发明的主体工作。

发明内容

基于传统均匀化方法的不足之处，本发明针对一种颗粒增强复合材料，着重解决当颗粒体积分数增大时颗粒发生重叠而导致均化计算结果失真的问题，提供了更为可靠的计算结果，为新型先进材料的使用和结构的优化设计提供充分的依据。

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

一种颗粒随机分布的复合材料宏观有效性质的预测方法，该方法包括下述步骤：

1）基于随机序列添加（Random Sequential Addition，RSA）方法使用MATLAB构建颗粒随机分布的复合材料表征体积单元（Representative VolumeElement，RVE）数值模型；

对含不规则颗粒夹杂的复合材料，其颗粒形状通过设置不同的边界曲线来描绘；颗粒中心点的坐标通过MATLAB随机数产生；本发明采用RSA方法生成圆形或椭圆形颗粒随机分布的复合材料RVE数值模型；

2）针对传统RSA方法生成的复合材料RVE数值模型，消除当颗粒体积分数较大时的颗粒重叠现象；

颗粒体积分数较大时，将颗粒与基体的体积模量、剪切模量及体积分数交换位置，使颗粒填充过程转变为基体填充过程；

3）对颗粒随机分布的复合材料RVE数值模型进行有限元分析计算，求得RVE数值模型有效性质的数值解；

①对RVE进行网格划分，并对基体与颗粒的界面相进行精细化处理，赋予这些界面单元不同于基体单元和颗粒单元的参数；

②施加不同的力学边界条件，求得有限元方法下RVE有效性质的数值解；

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

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

进一步地，所述步骤1）中，基于RSA方法使用MATLAB构建颗粒随机分布的复合材料RVE数值模型，是在基体材料中逐个生成圆形或椭圆形颗粒夹杂，且任意两个颗粒的中心距必须大于或等于颗粒直径，即颗粒不能重叠，RSA方法的包括下述步骤：

1a）根据颗粒体积分数求v₂和颗粒数N求颗粒半径

$R=\sqrt{v_2/\left(\mathrm{N\π}\right)}---\left(1\right)$

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

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

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

进一步地，所述步骤2）中，当颗粒体积分数较大时，传统RSA方法生成的RVE必然存在颗粒重叠现象；消除当颗粒体积分数较大时的颗粒重叠现象按下述方法进行改进：

在当颗粒体积分数v₂≥0.5时，将公式（1）中颗粒半径R替换为基体半径

$R_M=\sqrt{(1-v_2)/\mathrm{N\&pi;}}---\left(3\right)$

将基体的有效模量与颗粒的有效模量交换位置，将颗粒体积分数v₂变为基体体积分数1-v₂。这时，计算复合材料在颗粒体积分数为v₂时的有效体积模量k^*和有效剪切模量μ^*的公式将从

k^*=f₀(k₁,μ₁;k₂,μ₂;v₂)   （4）

μ^*=g₀(k₁,μ₁;k₂,μ₂;v₂)   （5）

转化为

k^*=f₁(k₂,μ₂;k₁,μ₁;1-v₂)   （6）

μ^*=g₁(k₂,μ₂;k₁,μ₁;1-v₂)   （7）

式中，下脚标1和2分别代表基体和颗粒；k₁和μ₁分别是基体的体积模量与剪切模量；k₂和μ₂分别是颗粒的体积模量与剪切模量；f和g表示有效模量（k^*，μ^*）是关于材料参数（k₁，μ₁，k₂，μ₂，v₂）的函数。本发明采用有限元方法求解，这是一种数值解法，故这里不必给出f和g的解析形式。

在RSA方法中，这一转换将使大体积分数的颗粒填充过程转变为小体积分数的基体填充过程，消除了RVE中的颗粒重叠现象。

进一步地，所述步骤3）中施加不同的力学边界条件，为施加了线性位移均匀应变边界条件，通过下式实现：

在线性位移边界条件中

$u_i^0\left(S\right)={\mathrm{\&epsiv;}}_{\mathrm{ij}}^0{x}_j---\left(8\right)$

式中，是常数，x_j是节点的位置坐标,是单元节点的位移，即

$\stackrel{\&OverBar;}{{\mathrm{\&epsiv;}}_{\mathrm{ij}}}={\mathrm{\&epsiv;}}_{\mathrm{ij}}^0---\left(9\right)$

用体积平均法和相应的边界值计算出应力、应变和应变能的有效量；对应变ε_ij=(u_i,j+u_j,i)/2，在公式（9）中用散度定理得

$\stackrel{\&OverBar;}{{\mathrm{\&epsiv;}}_{\mathrm{ij}}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Gamma;}}\frac{1}{2}(u_i{n}_j+u_j{n}_i)\mathrm{d\&Gamma;}---\left(10\right)$

式中，Γ是RVE的边界，n_i是Γ外法向量的第i个分量，V是RVE的体积，u_i和u_j是位移分量。

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

4a）建立随机均匀化模型

用基体与颗粒的弹性参数（k₁,μ₁;k₂,μ₂）及颗粒体积分数v₂这五个参数来表示复合材料的有效体积模量k^*和有效剪切模量μ^*，即

k^*=f_FEM(k₁,μ₁;k₂,μ₂;v₂)   （11）

μ^*=g_FEM(k₁,μ₁;k₂,μ₂;v₂)   （12）

式中，下脚标1和2分别代表基体和颗粒；k₁和μ₁分别是基体的体积模量与剪切模量；k₂和μ₂分别是颗粒的体积模量与剪切模量；f和g表示有效模量（k^*，μ^*）是关于材料参数（k₁，μ₁，k₂，μ₂，v₂）的函数，下脚标FEM表示采用有限元方法求解；

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

4b）求单个RVE的有效性质

对第i个RVE，可以求出其有效体积模量和有效剪切模量

$k_i^{*}=f_{\mathrm{FEM}}(k_{1i},k_{2i},{\mathrm{\&mu;}}_{1i},{\mathrm{\&mu;}}_{2i},v_{2i})---\left(13\right)$

${\mathrm{\&mu;}}_i^{*}=g_{\mathrm{FEM}}(k_{1i},k_{2i},{\mathrm{\&mu;}}_{1i},{\mathrm{\&mu;}}_{2i},v_{2i})---\left(14\right)$

其中，i=(1,2,…n)；

4c）对n个RVE的有效性质进行数理统计

通过上述步骤得到n个RVE的有效体积模量和有效剪切模量对其进行数理统计，其平均值和均方差（标准差）分别为

$E_{k_{\mathrm{FEM}}^{*}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_i^{*}---\left(15\right)$

$E_{{\mathrm{\&mu;}}_{\mathrm{FEM}}^{*}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_i^{*}---\left(16\right)$

${\mathrm{SD}}_{k_{\mathrm{FEM}}^{*}}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(k_i^{*}-E_{k_{\mathrm{FEM}}^{*}})}^2}---\left(17\right)$

${\mathrm{SD}}_{{\mathrm{\&mu;}}_{\mathrm{FEM}}^{*}}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&mu;}}_i^{*}-E_{{\mathrm{\&mu;}}_{\mathrm{FEM}}^{*}})}^2}---\left(18\right)$

式中，和分别n个RVE体积模量和剪切模量的平均值，和分别代表n个RVE体积模量和剪切模量的标准差；

上述公式求得的值和即可视为颗粒随机分布时均化模型下的有效弹性常数。

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

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

A、使用MATLAB构建任意形状且夹杂位置随机分布的复合材料RVE数值模型。其中，夹杂的形状通过设置不同的边界曲线来描绘，颗粒中心点的坐标通过MATLAB随机数产生。这一模型消除了传统均匀化方法中夹杂周期分布这一假设的不全面性。在接下来的步骤中，选取一种圆形颗粒随机分布的RVE作为研究对象。本发明中的数值模型不仅适用于颗粒周期分布的RVE，也适用于夹杂随机分布的复合材料RVE，其适用性更广。

B、在夹杂体积分数较大的情况下，改进了生成RVE的RSA方法，将颗粒与基体的体积模量、剪切模量及体积分数交换位置，使颗粒填充过程转变为基体填充过程，消除了颗粒重叠现象。

C、使用构建了均化模型，选取较大的样本容量，将计算结果的数理统计平均值近似地视为材料的宏观有效性质。可以获得任意颗粒体积分数情况下的任意加细数值解，为复合材料结构的优化设计奠定了理论和技术基础。

附图说明

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

图2是传统均匀化方法中的周期性RVE模型。

图3是生成随机RVE的RSA算法流程图。

图4是两种颗粒位置随机分布的RVE：（a）圆形颗粒RVE；（b）椭圆颗粒RVE。

图5是改进RSA算法后获得的RVE模型。

图6是两种不同的有限元网格划分示意图。

图7是对基体与颗粒界面相的网格精细化处理过程。

图8是传统均匀化方法求得的材料宏观有效性质随颗粒体积分数变化的曲线图。

图9是改进后的均匀化方法求得的材料宏观有效性质随颗粒体积分数变化的曲线图。

具体实施方式

下面结合附图对本发明的具体实施方式作进一步的说明。如图1所示，本发明共有5个主要流程：

1.构建随机RVE数值模型

对于一种由颗粒和基体组成的二维两相复合材料，传统的均匀化方法通常假设颗粒在RVE内部呈周期性分布，从而将周期性的微观结构进一步简化为单颗粒RVE，如图2所示。这种简化方法虽然大大降低了计算难度，但却忽略了颗粒分布上客观存在的随机性。

若采用RSA方法，则可生成颗粒随机分布的RVE，其流程图如图3所示。传统RSA方法的核心思想及操作流程：在基体材料中逐个生成颗粒夹杂。生成颗粒的实质是在微观结构中随机确定颗粒中心点的坐标值。其流程的详细步骤如下：

首先选取微观结构的边长为单位1，任意给定需要生成的颗粒数N及颗粒体积分数V₂，即可求得每个颗粒的半径

$R=\sqrt{v_2/\left(\mathrm{N\&pi;}\right)}---\left(1\right)$

然后在微观结构内第一个颗粒P₁的中心点坐标(x₁,y₁)，接着随机确定第二个颗粒P₂的中心点坐标(x₂,y₂)。接下来计算颗粒P₂与颗粒P₁的中心距

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

若L₁₂>2R（即是说颗粒P₂与颗粒P₁不相交），则颗粒P₂生成成功，否则重新生成颗粒P₂直至颗粒P₂与颗粒P₁不相交。

当确定第n(n≥2)个颗粒的中心点坐标时，需计算第n个颗粒与所有前n-1个颗粒的中心距，若L_1n,L_2n,...,L_n-1,n这n-1个中心距都大于2R，则第n个颗粒生成成功，否则重新确定第n个颗粒的中心点坐标直至它与所有已生成的颗粒都不相交。图4给出了传统RSA方法在体积比v₂=0.4时产生的两种颗粒随机分布的微观结构，图4-a为圆形颗粒，图4-b为椭圆颗粒。

2.改进的RSA技术方案

当颗粒在复合材料中所占的体积分数逐渐增大如达到0.9时，根据图3中的流程生成的微观结构中4个颗粒将发生严重重叠，如图5-a所示。为了消除这种现象，可以在v₂增大到一定程度如当v₂≥0.5时，采用基体填充代替颗粒填充，假设体积比是v₂的颗粒已经存在于微观结构中，这时需将1-v₂的基体材料填充进去，在使用MATLAB程序实现该技术方案时只需将原随机序列添加法中基体与颗粒的弹性参数互换即可，这时生成的RVE如图5-b所示，

将公式（1）中颗粒半径R替换为基体半径

$R_M=\sqrt{(1-v_2)/\mathrm{N\&pi;}}---\left(3\right)$

颗粒体积分数v₂变为基体体积分数1-v₂；将1-v₂的黑色基体材料随机填充进去之后，微观结构中剩余的白色圆形部分就是微观结构中体积比是v₂的颗粒夹杂,很明显四个随机分布的颗粒未发生重叠。也就是说，颗粒体积分数为v₂的颗粒生成过程（图5-a）等价于基体体积分数为1-v₂的基体生成过程（图5-b）。

这时，计算复合材料在颗粒体积分数为v₂时的有效体积模量k^*和有效剪切模量μ^*的公式将从

k^*(v₂)=f(k₁,μ₁;k₂,μ₂;v₂)   （4）

μ^*(v₂)=g(k₁,μ₁;k₂,μ₂;v₂)   （5）

转化为

k^*(v₂)=f(k₂,μ₂;k₁,μ₁;1-v₂)   （6）

μ^*(v₂)=g(k₂,μ₂;k₁,μ₁;1-v₂)   （7）

式中，下脚标1和2分别代表基体和颗粒；k₁和μ₁分别是基体的体积模量与剪切模量；k₂和μ₂分别是颗粒的体积模量与剪切模量；f和g表示有效模量（k^*，μ^*）是关于材料参数（k₁，μ₁，k₂，μ₂，v₂）的函数。

本发明采用有限元方法求解，这是一种数值解法，故这里不必给出f和g的解析形式。

3.有限元分析

3.1网格划分及精细化处理

如前所述，当前最有效地均化分析法是有限元方法，而有限元网格的划分方法直接决定了材料宏观弹性常数的计算精度。理论上，使用图6-b所示的划分方法比图6-a所示的方法效果更好。在图6-b中，可以对每个单元指定基体属性或颗粒属性，这种方法组集起来的平衡方程将是一个超大维数的线性方程组，给求解带来困难；当颗粒形状不规则时，这种方法不再适用。采用图6-a所示的正方形方法可以降低网格划分的难度，但必须对基体与颗粒的界面单元进行后处理。图7-a和图7-b为网格的精细化处理示意图，中间空白部分是颗粒单元。

3.2施加边界条件

求解材料的宏观有效响应如弹性常数，必须对RVE网格施加边界条件。常用的边界条件有四种，分别是线性位移边界条件、均匀应力边界条件、周期边界条件、混合边界条件。本发明将对从微观结构中取出的RVE施加线性位移边界条件来得出材料的有效响应。

在线性位移边界条件中

$u_i^0\left(S\right)={\mathrm{\&epsiv;}}_{\mathrm{ij}}^0{x}_j---\left(8\right)$

式中，是常数，x_j是节点的位置坐标,事单元节点的位移。可以证明，平均应变等于边界上的常应变，即

$\stackrel{\&OverBar;}{{\mathrm{\&epsiv;}}_{\mathrm{ij}}}={\mathrm{\&epsiv;}}_{\mathrm{ij}}^0---\left(9\right)$

用体积平均法和相应的边界值可计算出应力、应变和应变能的有效量。对应变ε_ij=(u_i,j+u_j,i)/2，在公式（6）中用散度定理得

$\stackrel{\&OverBar;}{{\mathrm{\&epsiv;}}_{\mathrm{ij}}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Gamma;}}\frac{1}{2}(u_i{n}_j+u_j{n}_i)\mathrm{d\&Gamma;}---\left(10\right)$

式中，Γ是RVE的边界，n_i是Γ外法向量的第i个分量，V是RVE的体积，u_i和u_j是位移分量。

4.构建随机均化模型

4.1均化方法与材料的有效性能

均匀化方法是基于复合材料的细观场和细观性能来估计宏观有效场和宏观有效性能的一种处理方法。在传统的均匀化方法中，局部或者细观的应力σ_ij、应变ε_ij的体积平均可以定义为

$\stackrel{\&OverBar;}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&sigma;}}_{\mathrm{ij}}\mathrm{d\&Omega;}---(4-1)$

$\stackrel{\&OverBar;}{{\mathrm{\&epsiv;}}_{\mathrm{ij}}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}{\mathrm{\&epsiv;}}_{\mathrm{ij}}\mathrm{d\&Omega;}---(4-2)$

式中，Ω表示RVE所占的区域，V是RVE的体积，横杠表示有效量。对一个弹性体，应变能密度的体积平均可表示为

$\stackrel{\&OverBar;}{\mathrm{\&omega;}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}\mathrm{\&omega;d\&Omega;}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{\mathrm{\&sigma;}}_{\mathrm{ij}}{\mathrm{\&epsiv;}}_{\mathrm{ij}}\mathrm{d\&Omega;}$

$=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{c}_{\mathrm{ijkl}}{\mathrm{\&epsiv;}}_{\mathrm{ij}}{\mathrm{\&epsiv;}}_{\mathrm{kl}}\mathrm{d\&Omega;}$

$=\frac{1}{V}{\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{f}_{\mathrm{ijkl}}{\mathrm{\&sigma;}}_{\mathrm{ij}}{\mathrm{\&sigma;}}_{\mathrm{kl}}\mathrm{d\&Omega;}---(4-3)$

式中，σ_ijε_ij/2=ω是应变能密度，c_ijkl是材料的刚度系数，f_ijkl(f=c^-1)是相应的柔度系数，其取值因相的变化而不同。另外，宏观应变能密度也可表示为

$\stackrel{\&OverBar;}{\mathrm{\&omega;}}=\frac{1}{2}\stackrel{\&OverBar;}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}\stackrel{\&OverBar;}{{\mathrm{\&epsiv;}}_{\mathrm{ij}}}---(4-4)$

由能量守恒定律及公式（4-1）-（4-3）三式可知，复合材料的有效刚度系数可定义为

$\stackrel{\&OverBar;}{{\mathrm{\&sigma;}}_{\mathrm{ij}}}=\stackrel{\&OverBar;}{c_{\mathrm{ijkl}}}\stackrel{\&OverBar;}{{\mathrm{\&epsiv;}}_{\mathrm{kl}}}---(4-5)$

4.2有效性质进行的均化处理

复合材料的有效应力场、应变场、能量场是通过有限元方法（FEM）计算出来的，在求得有效场之后，即可求得材料的宏观有效响应。用基体与颗粒的弹性参数（k₁,μ₁;k₂,μ₂）及颗粒体积分数v₂这五个参数来表示复合材料的有效体积模量k^*和有效剪切模量μ^*，即

k^*=f_FEM(k₁,μ₁;k₂,μ₂;v₂)   （11）

μ^*=g_FEM(k₁,μ₁;k₂,μ₂;v₂)   （12）

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

对第i个RVE，可以求出材料的有效性能

$k_i^{*}=f_{\mathrm{FEM}}(k_{1i},k_{2i},{\mathrm{\&mu;}}_{1i},{\mathrm{\&mu;}}_{2i},v_{2i})---\left(13\right)$

${\mathrm{\&mu;}}_i^{*}=g_{\mathrm{FEM}}(k_{1i},k_{2i},{\mathrm{\&mu;}}_{1i},{\mathrm{\&mu;}}_{2i},v_{2i})---\left(14\right)$

其中，i=(1,2,…n)，最后对这n组有效性能数据进行数理统计，其平均值和均方差分别为

$E_{k_{\mathrm{FEM}}^{*}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_i^{*}---\left(15\right)$

$E_{{\mathrm{\&mu;}}_{\mathrm{FEM}}^{*}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_i^{*}---\left(16\right)$

${\mathrm{SD}}_{k_{\mathrm{FEM}}^{*}}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(k_i^{*}-E_{k_{\mathrm{FEM}}^{*}})}^2}---\left(17\right)$

${\mathrm{SD}}_{{\mathrm{\&mu;}}_{\mathrm{FEM}}^{*}}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&mu;}}_i^{*}-E_{{\mathrm{\&mu;}}_{\mathrm{FEM}}^{*}})}^2}---\left(18\right)$

式中，和分别n个RVE体积模量和剪切模量的平均值，和分别代表n个RVE体积模量和剪切模量的标准差。

上述公式求得的值和即可视为颗粒随机分布时均化模型下的有效弹性常数。

5.模拟结果

第4步在改进随机序列添加法生成颗粒随机分布的RVE和有限元方法的基础上，依据数理统计的方法求解了一种复合材料的宏观有效弹性常数。设基体的体积模量和剪切模量分别为k₁=4，μ₁=1；颗粒的体积模量和剪切模量分别为k₂=40，μ₂=10；选取模拟次数n=1000次。施加线性位移边界条件后,获得了复合材料的有效体积模量及剪切模量，并绘制了该有效性能随颗粒体积分数变化的曲线图。（注：有效模量的单位为GPa）

5.1改进前和改进后获得的有效弹性常数对比

当颗粒体积分数v₂从0.05逐步递增至0.95时，对每一个v₂，均按照上述随机均化模型进行n=1000次模拟，求得材料的宏观有效体积模量有效剪切模量及其均方差的计算值见表1所示，线性位移边界条件下及其3σ边界曲线图如图8所示。

表1线性位移边界下的宏观有效弹性常数（单位：GPa）

从表1及图8可以明显看出：在v₂∈[0，0.5]，和的均值以及均方差的变化还是正常的，没有突变和拐点。但随着v₂的进一步增大，在同一个v₂处和的取值及其变化范围明显出现了收敛和失真的现象。

当v₂>0.5时，采用改进的技术方案求得材料的宏观有效体积模量有效剪切模量及其均方差的值如表2所示。将采用改进的RSA方案的计算值绘制到图8中，并用下脚标modified区别开来，得到如图9所示的线性位移边界条件下改进的及其3σ边界曲线图。在三条不完全曲线中，带“*”的线条代表改进的方法求得的有效体积模量和有效剪切模量，带“+”的线条及不带符号的线条分别表示有效模量的3σ上界及下界(注:±3σ原则是工程中针对随机变量取值范围给出的一种原则,通常认为该随机变量的实际取值是在均值E附近上下波动,波动范围的上界是E+3SD,波动的下界是E-3SD)。可以看出随着v₂的增加，改进技术方案后获得的有效模量的均值和均方差的变化平稳，当v₂接近1时，意味着此时复合材料已经接近单质材料即有效模量的取值应该非常接近颗粒模量，图7和表1中的数据严重背离了这一常识性的原理，而表2中的数据和图形9中在v₂≥0.5时按照改进的技术方案获得的三条新曲线所表现出的趋势与此原理相吻合，明显优于表1和图8的结果。

表2改进RSA方法后的宏观有效弹性常数（单位：GPa）

通过上述对比可知，基于改进随机序列添加法产生颗粒随机分布的表征体积单元，采用均化计算方法，求出了有效弹性系数及其统计特征值，可以形成对实验观测结果的有效补充和预测，为后续复合材料结构的优化设计提供了客观、充分和真实的依据。

Claims (1)

1.一种颗粒随机分布的复合材料宏观有效性质的预测方法，其特征在于，该方法包括下述步骤：

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

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

2)针对传统RSA方法生成的复合材料RVE数值模型，消除当颗粒体积分数较大时的颗粒重叠现象；

颗粒体积分数较大时，将颗粒与基体的体积模量、剪切模量及体积分数交换位置，使颗粒填充过程转变为基体填充过程；

3)对颗粒随机分布的复合材料RVE数值模型进行有限元分析计算，求得RVE数值模型有效性质的数值解；

①对RVE进行网格划分，并对基体与颗粒的界面相进行精细化处理，赋予这些界面单元不同于基体单元和颗粒单元的参数；

②施加不同的力学边界条件，求得有限元方法下RVE有效性质的数值解；

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

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

所述步骤1)中，基于RSA方法使用MATLAB构建颗粒随机分布的复合材料RVE数值模型，是在基体材料中逐个生成圆形或椭圆形颗粒夹杂，且任意两个颗粒的中心距必须大于或等于颗粒直径，即颗粒不能重叠，RSA方法包括下述步骤：

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

$R=\sqrt{v_2/\left(\mathrm{N\&pi;}\right)}---\left(1\right)$

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

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

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

所述步骤2)中，消除当颗粒体积分数较大时的颗粒重叠现象按下述方法进行改进：

在当颗粒体积分数v₂≥0.5时，将公式(1)中颗粒半径R替换为基体半径

$R_M=\sqrt{(1-v_2)/\mathrm{N\&pi;}}---\left(3\right)$

将基体的有效模量与颗粒的有效模量交换位置，将颗粒体积分数v₂变为基体体积分数1-v₂；这时，计算复合材料在颗粒体积分数为v₂时的有效体积模量k^*和有效剪切模量μ^*的公式将从

k^*＝f₀(k₁,μ₁；k₂,μ₂；v₂)   (4)

μ^*＝g₀(k₁,μ₁；k₂,μ₂；v₂)   (5)

转化为

k^*＝f₁(k₂,μ₂；k₁,μ₁；1-v₂)   (6)

μ^*＝g₁(k₂,μ₂；k₁,μ₁；1-v₂)   (7)

式中，下脚标1和2分别代表基体和颗粒；k₁和μ₁分别是基体的体积模量与剪切模量；k₂和μ₂分别是颗粒的体积模量与剪切模量；f和g表示有效模量(k^*，μ^*)是关于材料参数(k₁，μ₁，k₂，μ₂，v₂)的函数；

所述步骤3)中施加不同的力学边界条件，求得有限元方法下RVE有效性质的数值解，通过下式实现：

在线性位移边界条件中

$u_i^0\left(S\right)={\mathrm{\&epsiv;}}_{\mathrm{ij}}^0{x}_j---\left(8\right)$

式中，是常数，x_j是节点的位置坐标,是单元节点的位移，即

$\stackrel{\&OverBar;}{{\mathrm{\&epsiv;}}_{\mathrm{ij}}}={\mathrm{\&epsiv;}}_{\mathrm{ij}}^0---\left(9\right)$

用体积平均法和相应的边界值计算出应力、应变和应变能的有效量；对应变ε_ij＝(u_i,j+u_j,i)/2，在公式(9)中用散度定理得

$\stackrel{\&OverBar;}{{\mathrm{\&epsiv;}}_{\mathrm{ij}}}=\frac{1}{V}{\&Integral;}_{\mathrm{\&Gamma;}}\frac{1}{2}(u_i{n}_j+u_j{n}_i)\mathrm{d\&Gamma;}---\left(10\right)$

式中，Γ是RVE的边界，n_i是Γ外法向量的第i个分量，V是RVE的体积，u_i和u_j是位移分量；

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

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

用基体与颗粒的弹性参数(k₁,μ₁；k₂,μ₂)及颗粒体积分数v₂这五个参数来表示复合材料的有效体积模量k^*和有效剪切模量μ^*，即

k^*＝f_FEM(k₁,μ₁；k₂,μ₂；v₂)   (11)

μ^*＝g_FEM(k₁,μ₁；k₂,μ₂；v₂)   (12)

式中，下脚标1和2分别代表基体和颗粒；k₁和μ₁分别是基体的体积模量与剪切模量；k₂和μ₂分别是颗粒的体积模量与剪切模量；f和g表示有效模量(k^*，μ^*)是关于材料参数(k₁，μ₁，k₂，μ₂，v₂)的函数，下脚标FEM表示采用有限元方法求解；

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

4b)求单个RVE的有效性质

对第i个RVE，求出其有效体积模量和有效剪切模量

$k_i^{*}=f_{\mathrm{FEM}}(k_{1i},k_{2i},{\mathrm{\&mu;}}_{1i},{\mathrm{\&mu;}}_{2i},v_{2i})---\left(13\right)$

${\mathrm{\&mu;}}_i^{*}=g_{\mathrm{FEM}}(k_{1i},k_{2i},{\mathrm{\&mu;}}_{1i},{\mathrm{\&mu;}}_{2i},v_{2i})---\left(14\right)$

其中，i＝(1,2,...n)；

4c)对n个RVE有效性质进行数理统计

通过上述步骤得到n个RVE的有效体积模量和有效剪切模量对其进行数理统计，其平均值和均方差(标准差)分别为

$E_{k_{\mathrm{FEM}}^{*}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{k}_i^{*}---\left(15\right)$

$E_{{\mathrm{\&mu;}}_{\mathrm{FEM}}^{*}}=\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{\&mu;}}_i^{*}---\left(16\right)$

${\mathrm{SD}}_{k_{\mathrm{FEM}}^{*}}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{(k_i^{*}-E_{k_{\mathrm{FEm}}^{*}})}^2}---\left(17\right)$

${\mathrm{SD}}_{{\mathrm{\&mu;}}_{\mathrm{FEM}}^{*}}=\sqrt{\frac{1}{n}\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{({\mathrm{\&mu;}}_i^{*}-E_{{\mathrm{\&mu;}}_{\mathrm{FEm}}^{*}})}^2}---\left(18\right)$

式中，和分别为n个RVE体积模量和剪切模量的平均值，和分别代表n个RVE体积模量和剪切模量的标准差；

上述公式求得的值和即可视为颗粒随机分布时均化模型下的有效弹性常数。