CN105787167A - 聚合物基复合材料有效应力松弛系数的预测方法 - Google Patents

Abstract

本发明提供的聚合物基复合材料有效应力松弛系数的预测方法，首先基于玻尔兹曼叠加原理建立时间域内积分形式的线粘弹性材料本构方程；然后对于由周期性单胞构成的聚合物基复合材料，通过能量法定义有效应力松弛刚度；其次基于变分渐近均匀化理论建立求解聚合物基复合材料有效应力松弛刚度矩阵的细观力学模型；最后基于构建的细观力学模型预测聚合物基复合材料的有效应力松弛刚度。该方法仅需要一次求解过程即可获得不同方向上不同材料的属性，相对于在不同加载条件下重复运行的方法更简便、高效、快捷。计算的有效应力松弛系数与波动函数精度相一致，无需引入更多近似关系，计算精度更高，聚合物基复合材料的强度和寿命预测更加准确。

Description

聚合物基复合材料有效应力松弛系数的预测方法

技术领域

本发明涉及材料力学性能分析领域，具体涉及一种聚合物基复合材料有效应力松弛系数的预测方法。

背景技术

聚合物基复合材料主要由各种长短纤维嵌入有机聚合物基体构成，具有单一材料不可兼顾的优良性能，在飞机机身、机翼等主承力结构中的应用比例不断提高。由于聚合物基体具有粘弹性性质，聚合物基复合材料通常也表现出显著的粘弹性行为，即应力和应变的大小随时间和温度变化。聚合物基复合材料的粘弹性现象源自聚合物基体的长分子链。粘弹性诱发的徐变和应力松弛限制了先进复合材料结构在工程领域的长寿命工作。因此需对聚合物基复合材料的粘弹性性质和有效应力松弛现象进行深入研究。

细观力学模型是分析聚合物基复合材料粘弹性行为的主要工具，所以目前最常用的方法是根据细观力学模型来预测聚合物基复合材料的有效应力松弛系数，从而来了解聚合物基复合材料的粘弹性性质和有效应力松弛现象。Hashin首次提出了对应原理，通过相应的线粘弹性理论将粘弹性非均匀介质的有效松弛和徐变函数与非均匀介质的有效弹性模量关联起来；Park基于Prony级数发展了线粘弹性材料函数间互换的准确有效的数值方法，该方法直接适用于时间、频率和Laplace(拉普拉斯)转换域内的模量和协调函数的互换。最常用的表征聚合物基复合材料粘弹性行为的方法是Laplace转换和Laplace转置，并运用弹性-粘弹性对应原理。Brinson应用有限元法在Laplace转换域中分析了两相和三相粘弹性复合材料，并与Mori-Tanaka模型进行对比。Levin提出了复合材料有效粘弹性行为细观力学建模的分析方法，使用体积分数指数运算来描述组分的粘弹性属性。最近，Abadi使用商业有限元软件ABAQUS预测了单胞周期性边界条件下纤维增强聚合物基复合材料应力松弛响应和徐变响应。国内李丹、胡更开等基于Laplace变换和双夹杂相互作用的弹性模型，提出了一种新的预测颗粒增强聚合物材料有效粘弹性性质的细观力学模型。

上述方法大多基于某些特定的假设，因此都有特定的适用范围，分析效率低，有待进一步完善。通过Laplace转换、反演和对应原理计算得到的结果相比精确解有一定偏差，有时甚至会导致误差。要将上述模型结果应用于材料的强度和寿命预测显然是不够的。

发明内容

针对现有技术中存在的上述不足，本发明专利提供一种高效可靠的聚合物基复合材料有效应力松弛系数的预测方法，解决现有技术存在的分析效率低，精度差的不足，尤其是克服有效应力松弛刚度预测精度不足的缺陷。

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

聚合物基复合材料有效应力松弛系数的预测方法，其特征在于，包括如下步骤：

1)基于玻尔兹曼叠加原理建立时间域内积分形式的线粘弹性材料本构方程；

2)对于由周期性单胞构成的聚合物基复合材料，通过能量法定义有效应力松弛刚度；

3)基于变分渐近均匀化理论建立求解聚合物基复合材料有效应力松弛刚度矩阵的细观力学模型，并求得单胞能量密度；

4)基于构建的细观力学模型预测聚合物基复合材料的有效应力松弛系数。

进一步，所述步骤1具体为：

基于玻尔兹曼叠加原理建立时间域内积分形式的线粘弹性材料本构方程，即：

${\mathrm{\σ}}_{ij}\left(t\right)={\∫}_{-\mathrm{\∞}}^t\[C_{ijkl}\left(t-\mathrm{\τ}\right){\stackrel{\·}{\mathrm{\ϵ}}}_{kl}\left(\mathrm{\τ}\right)\]d\mathrm{\τ};$

式中：C_ijkl(t)为应力松弛刚度；为应变率，σ_ij(t)为应力张量，t为时间。

进一步，所述步骤2具体为：

对于由周期性单胞构成的聚合物基复合材料，通过能量法定义有效应力松弛刚度，即：

${\stackrel{\‾}{\mathrm{\σ}}}_{ij}\left(t\right)=C_{ijkl}^{*}\left(t\right){\stackrel{\‾}{\mathrm{\ϵ}}}_{ij}^0;$

$\frac{1}{\mathrm{\Ω}}{\∫}_{\mathrm{\Ω}}\frac{1}{2}{C}_{ijkl}^{*}\left(t\right){\mathrm{\ϵ}}_{ij}\left(t\right){\mathrm{\ϵ}}_{kl}\left(t\right)d\mathrm{\Ω}=\frac{1}{2}{C}_{ijkl}^{*}\left(t\right){\stackrel{\‾}{\mathrm{\ϵ}}}_{ij}^0{\stackrel{\‾}{\mathrm{\ϵ}}}_{kl}^0;$

式中：Ω为单胞的体积域；带星号项表示细观力学模型求解得到的有效属性；上画线项表示均匀化后材料宏观分析使用的变量，ε_ij(t)为单胞内的局部应变，ε_kl(t)为应力松弛系数， ${\mathrm{\&epsiv;}}_{kl}\left(t\right)=\left\{\begin{array}{cc}0& t<0\\ {\mathrm{\&epsiv;}}_{kl}^0& t\&GreaterEqual;0\end{array}\right..$

进一步，所述步骤3具体为：基于变分渐近均匀化理论建立求解聚合物基复合材料有效应力松弛刚度矩阵的细观力学模型，即建立时间域内的单胞能量泛函：

${\mathrm{\&Pi;}}_{\mathrm{\&Omega;}}=\frac{1}{2\mathrm{\&Omega;}}{\&Integral;}_{\mathrm{\&Omega;}}{C}_{ijkl}\left(t\right)\&lsqb;{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{ij}^0+{\mathrm{\&chi;}}_{\left(i\right|j)}\&rsqb;\&lsqb;{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{kl}^0+{\mathrm{\&chi;}}_{\left(k\right|l)}\&rsqb;d\mathrm{\&Omega;};$

并受如下约束：

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

式中：χ_i称为波动函数；为全局应变；上标“+j”和“-j”表示单胞正、负边界表面上相应的量；

根据波动函数离散化得到离散形式泛函；根据变分原理，最小化离散泛函，将得到波动函数节点值列阵代入单胞泛函，得到单胞能量密度：

${\mathrm{\&Pi;}}_{\mathrm{\&Omega;}}=\frac{1}{2}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}^T{D}^{*}\stackrel{\&OverBar;}{\mathrm{\&epsiv;}};$

式中：D^*即为6×6阶有效应力松弛刚度矩阵，为由i＝1,2,3，j＝1,2,3构成的矩阵。

相比于现有技术，本发明具有如下优点：

1、本发明提供的聚合物基复合材料有效应力松弛系数的预测方法，采用变分渐近均匀化作为数学基础，仅含细观力学最基本理论，不需要任何周期性和边界假设，具有数学上的严密性。

2、本发明提供的聚合物基复合材料有效应力松弛系数的预测方法，采用非强制性边界条件，仅需要一次求解过程即可获得不同方向上不同材料的属性，相对于在不同加载条件下重复运行的方法更简便、高效、快捷。

3、本发明提供的聚合物基复合材料有效应力松弛系数的预测方法，计算的有效应力松弛系数与波动函数精度相一致，无需引入更多近似关系，计算精度更高，聚合物基复合材料的强度和寿命预测更加准确。

4、本发明提供的聚合物基复合材料有效应力松弛系数的预测方法，可模拟完全各向异性复合材料，突破了传统有限元只能处理宏观正交各向异性材料的局限性。

5、本发明提供的聚合物基复合材料有效应力松弛系数的预测方法，所有计算都在时间域内完成，不再需要传统线粘弹性复合材料计算所需的Laplace转换、反演和对应原理，可提高计算效率。

6、本发明提供的聚合物基复合材料有效应力松弛系数的预测方法，实用性强，通用性高，可显著提高聚合物基复合材料有效应力松弛问题的解算速度和效率。

附图说明

图1为周期性聚合物基复合材料略图，其中，1为聚合物基体，2为可识别单胞，a₁，a₂，a₃分别为沿局部坐标y₁，y₂，y₃轴的单胞半边长值。

图2为计算和的1/4单胞。

图3为计算横向应力松弛模量的二维模型。

图4为计算纵向应力松弛模量的三维有限元模型。

图5为不同纤维体积分数下有效应力松弛系数C₁₁(t)随时间变化图。

图6为不同纤维体积分数下有效应力松弛系数C₂₂(t)随时间变化图。

图7为不同纤维体积分数下有效应力松弛系数C₁₂(t)随时间变化图。

图8为不同纤维体积分数下有效应力松弛系数C₂₃(t)随时间变化图。

图9为不同纤维体积分数下有效应力松弛系数C₄₄(t)随时间变化图。

图10为不同纤维体积分数下有效应力松弛系数C₅₅(t)随时间变化图。

具体实施方式

下面结合实施例对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例：

1、线粘弹性材料应力松弛刚度的理论方程。

基于Bolzman(玻尔兹曼)叠加原理建立时间域内积分形式的线粘弹性材料本构方程

${\mathrm{\&sigma;}}_{ij}\left(t\right)={\&Integral;}_{-\mathrm{\&infin;}}^t\&lsqb;C_{ijkl}\left(t-\mathrm{\&tau;}\right){\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{kl}\left(\mathrm{\&tau;}\right)\&rsqb;d\mathrm{\&tau;}---\left(1\right)$

式中：C_ijkl(t)为应力松弛刚度；为应变率，σ_ij(t)为应力张量。

应力松弛在恒定应变下进行，即：

${\mathrm{\&epsiv;}}_{kl}\left(t\right)=\left\{\begin{array}{cc}0& t<0\\ {\mathrm{\&epsiv;}}_{kl}^0& t\&GreaterEqual;0\end{array}\right.---\left(2\right)$

式中：ε_kl(t)为应力松弛系数；上标“0”表示不随时间变化，但可随位置变化的常量。式(2)意味lim_t→-∞ε_kl(t)＝0。

对式(1)部分积分，得到

${\mathrm{\&sigma;}}_{ij}\left(t\right)=\left(C_{ijkl}\right(0)+{\&Integral;}_0^t\frac{\&part;C_{ijkl}\left(t-\mathrm{\&tau;}\right)}{\&part;\left(t-\mathrm{\&tau;}\right)}d\mathrm{\&tau;}){\mathrm{\&epsiv;}}_{kl}^0=C_{ijkl}\left(t\right){\mathrm{\&epsiv;}}_{kl}^0---\left(3\right)$

式(3)意味着瞬时应力值取决于应力松弛系数的瞬时值，而不是恒定应变载荷下线性粘弹性材料的历程效应。

2、有效应力松弛刚度的细观力学模型。

图1为周期性聚合物基复合材料略图。考虑由周期性单胞构成的多相粘弹性复合材料，多相复合材料的微观结构，引入两组坐标系x＝(x₁,x₂,x₃)和y＝(y₁,y₂,y₃)以方便细观力学模型推导。x_i为全局坐标系以描述宏观结构，与x_i平行的y_i为局部坐标系描述单胞，选择局部坐标系y_i的原点为单胞的几何中心。若单胞的边长为2a_i，则y_i∈[-a_i,a_i]。

2.1、线性粘弹性复合材料的有效应力松弛刚度。

式(3)可从如下瞬时势能密度泛函相对求偏导得出

$U\left(t\right)=\frac{1}{2}{C}_{ijkl}\left(t\right){\mathrm{\&epsiv;}}_{ij}^0{\mathrm{\&epsiv;}}_{kl}^0---\left(4\right)$

线性粘弹性复合材料有效应力松弛可定义为

${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{ij}\left(t\right)=C_{ijkl}^{*}\left(t\right){\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{ij}^0---\left(5\right)$

$\frac{1}{\mathrm{\&Omega;}}{\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{C}_{ijkl}^{*}\left(t\right){\mathrm{\&epsiv;}}_{ij}\left(t\right){\mathrm{\&epsiv;}}_{kl}\left(t\right)d\mathrm{\&Omega;}=\frac{1}{2}{C}_{ijkl}^{*}\left(t\right){\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{ij}^0{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{kl}^0---\left(6\right)$

式中：Ω为单胞的体积域；带星号项表示细观力学模型求解得到的有效属性；上画线项表示均匀化后材料宏观分析使用的变量，即

${\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{ij}^0=<{\mathrm{\&epsiv;}}_{kl}\left(t\right)>,{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{ij}=<{\mathrm{\&sigma;}}_{ij}\left(t\right)>---\left(7\right)$

式中： $<\&CenterDot;>=\frac{1}{\mathrm{\&Omega;}}{\&Integral;}_{\mathrm{\&Omega;}}<\&CenterDot;>d\mathrm{\&Omega;}.$

粘弹性聚合物基复合材料承受的外加应变为常量，但由于线粘弹性聚合物的应力松弛，单胞内的局部应变ε_ij(t)和局部应力σ_ij(t)随时间变化。

2.2变分渐近均匀化细观力学模型。

通过对以上泛函的最小化可得到求解聚合物基复合材料有效应力松弛刚度矩阵的变分渐近均匀化细观力学模型：

${\mathrm{\&Pi;}}_{\mathrm{\&Omega;}}=\frac{1}{2\mathrm{\&Omega;}}{\&Integral;}_{\mathrm{\&Omega;}}{C}_{ijkl}\left(t\right)\&lsqb;{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{ij}^0+{\mathrm{\&chi;}}_{\left(i\right|j)}\&rsqb;\&lsqb;{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{kl}^0+{\mathrm{\&chi;}}_{\left(k\right|l)}\&rsqb;d\mathrm{\&Omega;}---\left(8\right)$

并受如下约束

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

式中：χ_i称为波动函数；为全局应变；上标“+j”和“-j”表示单胞正、负边界表面上相应的量， ${\mathrm{\&chi;}}_i^{+j}={\mathrm{\&chi;}}_i{|}_{y_j=a_j},{\mathrm{\&chi;}}_i^{-j}={\mathrm{\&chi;}}_i{|}_{y_j=-a_j}.$

引入如下矩阵符号：

$\begin{array}{c}\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}=\&lsqb;{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{11}^0,2{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{12}^0,{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{22}^0,2{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{13}^0,2{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{23}^0,{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{33}^0\&rsqb;\\ \widehat{\mathrm{\&epsiv;}}=\&lsqb;{\widehat{\mathrm{\&epsiv;}}}_{11}^0,2{\widehat{\mathrm{\&epsiv;}}}_{12}^0,{\widehat{\mathrm{\&epsiv;}}}_{22}^0,2{\widehat{\mathrm{\&epsiv;}}}_{13}^0,2{\widehat{\mathrm{\&epsiv;}}}_{23}^0,{\widehat{\mathrm{\&epsiv;}}}_{33}^0\&rsqb;\end{array}---\left(10\right)$

式中：

$\begin{array}{c}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{ij}^0=\&lsqb;v_{i,j}\left(x\right)+v_{j,i}\left(x\right)\&rsqb;/2\\ {\widehat{\mathrm{\&epsiv;}}}_{ij}=\&lsqb;{\mathrm{\&chi;}}_{i,j}\left(t;x;y\right)+{\mathrm{\&chi;}}_{j,i}\left(t;x;y\right)\&rsqb;/2={\mathrm{\&chi;}}_{\left(i|j\right)}\end{array}---\left(11\right)$

式中：v_i(x)＝〈u_i(t；x；y)〉，u_i(t；x；y)为局部位移矢量，可表示为u_i(t；x；y)＝v_i(x)+y_jv_i,j+χ_i(t；x；y)。

式(11)中的矩阵可改写为

$\widehat{\mathrm{\&epsiv;}}=\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)}_{;2}& {\left(\right)}_{;3}\\ 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(12\right)$

式中：Γ_h为算子矩阵，χ为含位移波动函数χ_i的列阵，使用有限元将χ离散化为

χ(x_i；y_i)＝S(y_i)N(x_i)(13)

式中：S为形函数(不包括约束节点和从节点)；N为所有活动节点波动函数的节点值列阵。

将式(13)-(10)代入式(8)，得到泛函的离散形式为

${\mathrm{\&Pi;}}_{\mathrm{\&Omega;}}=\frac{1}{2\mathrm{\&Omega;}}(N^T{D}_{hh}N+2N^T{D}_{h\mathrm{\&epsiv;}}\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}+{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}^T{D}_{\mathrm{\&epsiv;}\mathrm{\&epsiv;}}\stackrel{\&OverBar;}{\mathrm{\&epsiv;}})---\left(14\right)$

式中：

$\begin{array}{c}{D}_{hh}={\&Integral;}_{\mathrm{\&Omega;}}{\left({\mathrm{\&Gamma;}}_{h}S\right)}^{T}D\left({\mathrm{\&Gamma;}}_{h}S\right)d\mathrm{\&Omega;},D_{h\mathrm{\&epsiv;}}={\&Integral;}_{\mathrm{\&Omega;}}{\left({\mathrm{\&Gamma;}}_{h}S\right)}^{T}Dd\mathrm{\&Omega;},\\ D_{\mathrm{\&epsiv;}\mathrm{\&epsiv;}}={\&Integral;}_{\mathrm{\&Omega;}}Dd\mathrm{\&Omega;},\end{array}---\left(15\right)$

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

$D_{hh}N=-D_{h\mathrm{\&epsiv;}}\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}---\left(16\right)$

由式(17)可知：波动函数节点值列阵N与线性相关，这意味着解可象征性地表示为

$N=N_0\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}---\left(17\right)$

将式(17)代入式(14)，得到单胞能量密度为

${\mathrm{\&Pi;}}_{\mathrm{\&Omega;}}=\frac{1}{2}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}^T(N_0^T{D}_{h\mathrm{\&epsiv;}}+D_{\mathrm{\&epsiv;}\mathrm{\&epsiv;}})\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}=\frac{1}{2}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}^T{D}^{*}\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}---\left(18\right)$

式中：D^*是6×6阶有效应力松弛刚度矩阵，N₀表示系数矩阵。

3、算例

基于构建的方法预测玻璃纤维增强聚合物基复合材料的有效应力松弛刚度圆形玻璃纤维在基体内呈矩形阵列(图1)。为验证构建模型的有效性和准确性，同时使用ABAQUS建立有限元单胞模型预测有效应力松弛刚度。

3.1组分材料属性

玻璃纤维：各向同性和线弹性材料，材料属性列于表1。

表1玻璃纤维的材料属性

聚合物：各向同性、线性粘弹性材料，其弹性松弛模量可使用Prony级数表示为

$E\left(t\right)=E_0\&lsqb;1-\underset{k=1}{\overset{n}{\&Sigma;}}{g}_k\left(1-e^{-t/{\mathrm{\&tau;}}_k}\right)\&rsqb;---\left(19\right)$

式中：E₀为瞬时杨氏模量；g_k为无量纲模量，τ_k为材料随时间松弛系数。为简化起见，考虑特殊情况，n＝1,g₁＝0.5,τ₁＝30，式(19)简化为

E(t)＝0.5E₀(1+e^-t/ρ)＝A+Be^-t/ρ(20)

式中：若E₀＝8000MPa，ρ＝30，则A＝B＝4000MPa。聚合物的泊松比为v＝0.4。

3.2、ABAQUS有限元单胞模型。

矩形阵列纤维增强复合材料的有效刚度矩阵具有对称性，有效应力松弛刚度矩阵可表示为：

$\left\{\begin{array}{c}{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{11}\\ {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{22}\\ {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{33}\\ {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{23}\\ {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{12}\\ {\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{13}\end{array}\right\}=\left[\begin{array}{cccccc}{C}_{11}^{*}\left(t\right)& C_{12}^{*}\left(t\right)& C_{12}^{*}\left(t\right)& 0& 0& 0\\ C_{12}^{*}\left(t\right)& C_{22}^{*}\left(t\right)& C_{23}^{*}\left(t\right)& 0& 0& 0\\ C_{12}^{*}\left(t\right)& C_{23}^{*}\left(t\right)& C_{22}^{*}\left(t\right)& 0& 0& 0\\ 0& 0& 0& C_{44}^{*}\left(t\right)& 0& 0\\ 0& 0& 0& 0& C_{55}^{*}\left(t\right)& 0\\ 0& 0& 0& 0& 0& C_{55}^{*}\left(t\right)\end{array}\right]\left\{\begin{array}{c}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{11}^0\\ {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{22}^0\\ {\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{33}^0\\ {\stackrel{\&OverBar;}{\mathrm{\&gamma;}}}_{23}^0\\ {\stackrel{\&OverBar;}{\mathrm{\&gamma;}}}_{12}^0\\ {\stackrel{\&OverBar;}{\mathrm{\&gamma;}}}_{13}^0\end{array}\right\}---\left(21\right)$

式中：表示有效应力松弛刚度系数。

3.2.1、计算

如图2所示，法向载荷和双对称轴下，仅需原单胞模型的1/4，即可计算作用在有限元模型上计算的位移约束为

$\begin{array}{cc}{u}_1{|}_{y_1=-a_1}=0,& u_1{|}_{y_1=a_1}=\mathrm{\&delta;}\\ u_2{|}_{y_2=-a_2}=0,& u_2{|}_{y_2=a_2}=0\\ u_3{|}_{y_3=-a_3}=0,& u_3{|}_{y_3=a_3}=0\end{array}---\left(22\right)$

式中：u_i为沿x_i方向的位移。

的计算式为：

$C_{11}^{*}\left(t\right)=\frac{{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{11}\left(t\right)}{{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{11}^0},C_{12}^{*}\left(t\right)=\frac{{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{22}\left(t\right)}{{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{11}^0}---\left(23\right)$

式中：由式(22)的边界条件产生。

3.2.2计算

同样使用图2的单胞模型计算作用在有限元模型上的位移约束为

$\begin{array}{cc}{u}_1{|}_{y_1=-a_1}=0,& u_1{|}_{y_1=a_1}=0\\ u_2{|}_{y_2=-a_2}=0,& u_2{|}_{y_2=a_2}=\mathrm{\&delta;}\\ u_3{|}_{y_3=-a_3}=0,& u_3{|}_{y_3=a_3}=0\end{array}---\left(24\right)$

的计算式为：

$C_{22}^{*}\left(t\right)=\frac{{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{22}\left(t\right)}{{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{22}^0},C_{23}^{*}\left(t\right)=\frac{{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{33}\left(t\right)}{{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{22}^0}---\left(25\right)$

式中：由式(24)的边界条件产生。

3.2.3、计算

如图3所示，使用二维平面应变模型计算横向应力松驰模量所需位移约束为

$\begin{array}{cc}{u}_2{|}_{y_2=-a_2}=u_2{|}_{y_2=a_2},& u_3{|}_{y_2=-a_2}=u_3{|}_{y_2=a_2},\\ u_2{|}_{y_3=-a_3}=u_2{|}_{y_3=a_3},& u_3{|}_{y_3=-a_3}=u_3{|}_{y_3=a_3},\end{array}---\left(26\right)$

进一步约束为

$\begin{array}{cc}{\mathrm{\&epsiv;}}_{33}{|}_{y_2=\&PlusMinus;a_2}=0,& {\mathrm{\&epsiv;}}_{33}{|}_{y_3=\&PlusMinus;a_3}=0\end{array}---\left(27\right)$

底角放置在滚轴上，以消除刚体位移。水平位移常量δ作用在单胞各角点上。的计算式为

$C_{44}^{*}\left(t\right)=\frac{{\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{23}\left(t\right)}{{\stackrel{\&OverBar;}{\mathrm{\&gamma;}}}_{23}^0}---\left(28\right)$

式中：全局横向剪切应变量

3.2.4、计算

如图4所示，使用三维单胞计算平均纵向应力松弛模量所需位移约束为：

$\begin{array}{c}\begin{array}{ccc}{u}_1{|}_{y_1=-a_1}=u_1{|}_{y_1=a_1},& u_3{|}_{y_1=-a_1}=u_3{|}_{y_1=a_1},& u_2{|}_{y_1=-a_1}=u_2{|}_{y_1=a_1},\end{array}\\ u_1{|}_{y_3=0}=u_2{|}_{y_3=0}=u_3{|}_{y_3=0}=0\\ \begin{array}{cc}{u}_1{|}_{y_3=a_3}=\mathrm{\&delta;}& u_2{|}_{y_3=a_3}=0\end{array}\end{array}---\left(29\right)$

的计算式为

$C_{55}^{*}\left(t\right)={\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{12}\left(t\right)/{\stackrel{\&OverBar;}{\mathrm{\&gamma;}}}_{12}^0---\left(30\right)$

式中：全局横向剪切应变常量

在所有的有限元模拟中，应变常量从t＝0开始作用，并一直持续到模拟结束，如式(2)所示。

3.3、有效应力松驰刚度预测

本发明和ABAQUS有限元预测的有效应力松驰系数如图5-10所示。由图中可看出：本发明预测结果与ABAQUS结果相一致；有效应力松驰刚度随vof(纤维体积分数)的增大而增大，而且当时间t≥50s时，有效应力松驰刚度基本保持不变；C₁₁(t)下降最小(约为5％)，而C₂₂(t),C₄₄(t),C₅₅(t)下降多达40％。本发明仅在一次分析中就可得到完整的有效瞬时系数，而ABAQUS需在各种加载和边界条件下多次运行。计算纤维增强复合材料有效属性时，本发明仅需如图3所示的二维网格划分信息(如节点和单元数目，节点坐标系等)，而ABAQUS有限元建模技术需不同的几何模型(图2～图4所示)。因此，本发明比ABAQUS有限元模型效率更高。

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

Claims (4)

1.聚合物基复合材料有效应力松弛系数的预测方法，其特征在于，包括如下步骤：

1)基于玻尔兹曼叠加原理建立时间域内积分形式的线粘弹性材料本构方程；

2)对于由周期性单胞构成的聚合物基复合材料，通过能量法定义有效应力松弛刚度；

3)基于变分渐近均匀化理论建立求解聚合物基复合材料有效应力松弛刚度矩阵的细观力学模型，并求得单胞能量密度；

4)基于构建的细观力学模型预测聚合物基复合材料的有效应力松弛系数。

2.如权利要求1所述的聚合物基复合材料有效应力松弛系数的预测方法，其特征在于，所述步骤1具体为：

基于玻尔兹曼叠加原理建立时间域内积分形式的线粘弹性材料本构方程，即：

${\mathrm{\&sigma;}}_{ij}\left(t\right)={\&Integral;}_{-\mathrm{\&infin;}}^t\&lsqb;C_{ijki}\left(t-\mathrm{\&tau;}\right){\stackrel{\&CenterDot;}{\mathrm{\&epsiv;}}}_{kl}\left(\mathrm{\&tau;}\right)\&rsqb;d\mathrm{\&tau;};$

式中：C_ijkl(t)为应力松弛刚度；为应变率，σ_ij(t)为应力张量，t为时间。

3.如权利要求2所述的聚合物基复合材料有效应力松弛系数的预测方法，其特征在于，所述步骤2具体为：

对于由周期性单胞构成的聚合物基复合材料，通过能量法定义有效应力松弛刚度，即：

${\stackrel{\&OverBar;}{\mathrm{\&sigma;}}}_{ij}\left(t\right)=C_{ijkl}^{*}\left(t\right){\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{ij}^0;$

$\frac{1}{\mathrm{\&Omega;}}{\&Integral;}_{\mathrm{\&Omega;}}\frac{1}{2}{C}_{ijkl}^{*}\left(t\right){\mathrm{\&epsiv;}}_{ij}\left(t\right){\mathrm{\&epsiv;}}_{kl}\left(t\right)d\mathrm{\&Omega;}=\frac{1}{2}{C}_{ijkl}^{*}\left(t\right){\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{ij}^0{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{kl}^0;$

式中：Ω为单胞的体积域；带星号项表示细观力学模型求解得到的有效属性；上画线项表示均匀化后材料宏观分析使用的变量，ε_ij(t)为单胞内的局部应变，ε_kl(t)为应力松弛系数， ${\mathrm{\&epsiv;}}_{kl}\left(t\right)=\left\{\begin{array}{cc}0& t<0\\ {\mathrm{\&epsiv;}}_{kl}^0& t\&GreaterEqual;0\end{array}\right..$

4.如权利要求3所述的聚合物基复合材料有效应力松弛系数的预测方法，其特征在于，所述步骤3具体为：基于变分渐近均匀化理论建立求解聚合物基复合材料有效应力松弛刚度矩阵的细观力学模型，即建立时间域内的单胞能量泛函：

${\mathrm{\&Pi;}}_{\mathrm{\&Omega;}}=\frac{1}{2\mathrm{\&Omega;}}{\&Integral;}_{\mathrm{\&Omega;}}{C}_{ijkl}\left(t\right)\&lsqb;{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{ij}^0+{\mathrm{\&chi;}}_{\left(i\right|j)}\&rsqb;\&lsqb;{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}_{kl}^0+{\mathrm{\&chi;}}_{\left(k\right|l)}\&rsqb;d\mathrm{\&Omega;};$

并受如下约束：

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

式中：χ_i称为波动函数；为全局应变；上标“+j”和“-j”表示单胞正、负边界表面上相应的量；

根据波动函数离散化得到离散形式泛函；根据变分原理，最小化离散泛函，将得到波动函数节点值列阵代入单胞泛函，得到单胞能量密度：

${\mathrm{\&Pi;}}_{\mathrm{\&Omega;}}=\frac{1}{2}{\stackrel{\&OverBar;}{\mathrm{\&epsiv;}}}^T{D}^{*}\stackrel{\&OverBar;}{\mathrm{\&epsiv;}};$

式中：D^*即为6×6阶有效应力松弛刚度矩阵，为由构成的矩阵。