基于时域渐进理论的复合材料粘弹性均匀化方法与流程

本发明属于理论力学、多尺度复合材料力学领域，涉及基于时域渐进理论的复合材料等效粘弹性松弛刚度的计算分析，以及借助商业软件开发的基于时域渐进理论的复合材料等效粘弹性的有限元计算方法，具体是一种基于时域渐进理论的复合材料粘弹性均匀化方法。

背景技术：

复合材料通常指不同组分材料复合制作的增强复合材料。相较于传统材料，复合材料具有强度高、质量轻、耐腐蚀等优点。然而，由于复合材料的不均匀特性，直接对每个微观结构进行分析研究非常困难。此外，由于复合材料常常具有粘弹性，导致其制造和服役过程中常常发生蠕变变形和应力松弛，影响复合材料稳定性及可靠性。因此，研究复合材料微观结构的粘弹性均匀化方法，即有助于准确计算复合材料宏观粘弹性，又有助于评价复合材料宏观粘弹性在复合材料制造和服役过程中的力学性能。

渐进均匀化方法是一种广泛使用的计算复合材料料性能的方法。基于复合材料周期性假设和摄动原理，渐进均匀化方法可以实现宏观位移与微观位移耦合计算。当前，渐进均匀化方法已广泛应用于弹性、热弹性、热传导系数等复合材料的均匀化研究。然而，弹性、热弹性、热传导系数都属于对时间独立的材料参数，现有均匀化方法仍然难以计算时间相关的等效粘弹性刚度。仅有的基于时域渐进理论开发的粘弹性均匀化方法需要自编有限元程序，存在开发难度高、工作量大等问题。另外，已有的时域渐进均匀化方法计算等效粘弹性松弛刚度需要分别计算粘弹性的弹性、粘性以及两者转换的伴随系数三个特征位移场，计算量大，过程复杂。综上，现有基于时域渐进理论在时间域上对复合材料进行粘弹性均匀化的有限元计算存在计算过程复杂、计算量大、难以便捷计算等效粘弹性松弛刚度等问题。

技术实现要素：

为解决现有技术中存在的问题，降低渐进均匀化法的使用难度，推进时域渐进理论在复合材料均匀化方面的应用，本发明公开了一种基于粘弹性时域渐进理论计算复合材料时域特征位移场和时域等效粘弹性松弛刚度的有限元实现方法。借助商业有限元软件，该均匀化方法可有效降低了时域特征位移场和时域等效粘弹性松弛刚度的有限元计算难度。

本发明是这样实现的：

所述的基于时域渐进理论的复合材料粘弹性均匀化方法具体如下：

步骤一、基于线性粘弹性本构方程和时域渐进理论给出渐进均匀化特征位移场的控制方程；

步骤二、基于线性粘弹性本构方程和时域渐进理论给出复合材料等效粘弹性松弛刚度计算表达式；

步骤三、完成渐进均匀化特征位移场控制方程有限元计算，求解时间域特征位移场，求解时间域等效粘弹性松弛刚度。

进一步，本发明均匀化对象为周期性或统计意义上为周期性、具线性粘弹性的复合材料；均匀化复合材料对象的材料特性包括各向同性粘弹性、横观各向异性粘弹性、各向异性粘弹性；针对二维或三维微观结构的代表性体积单元，计算其时间相关的特征位移场，并计算其等效松弛刚度。

进一步，所述的步骤一中基于线性粘弹性本构方程和时域渐进理论给出的特征位移场控制方程表达式为：

式中，cijkl(y,t)为材料粘弹性松弛刚度，特征位移场，y是微观尺度坐标，t是时间，下标i,j,k,l,m,n代表了材料主轴方向，按照二维或三维模型，i,j,k,l,m,n分别为{1,2}或{1,2,3}。

进一步，所述的步骤二中基于线性粘弹性本构方程和时域渐进理论给出的等效粘弹性松弛刚度计算表达式为：

式中，为等效粘弹性松弛刚度，iklmn是6×6单元矩阵,δ(t)是狄拉克函数，<·>y为体平均算子y是单元体积，v^y是代表性体积单元总体积，φ(x,y)为关于y的周期性场函数，x是宏观尺度坐标。

进一步，所述的步骤三中利用不同商业有限元软件和脚本程序完成计算，其中商业软件的选择中：建模，定义周期性边界是digimat，定义热应变载荷是matlab，求解与后处理是abaqus。所述的渐进均匀化特征位移场控制方程的实施步骤如下：

3.1、通过商业软件建立复合材料代表性体积单元模型，设置各相材料弹性、粘弹性参数，划分有限元网格，获得有限元模型节点信息、单元信息；

3.2、借助商业软件或脚本程序定义周期性边界；

3.3、设置分析步，设置自定义热膨胀系数，加载温度载荷；

3.4、求解控制方程，计算时间相关的特征位移场；

3.5、计算等效粘弹性松弛刚度。

进一步，所述的步骤3.3中借助材料热应变构造特征位移场控制方程右侧应力载荷完成计算，其表达式为：

式中，f(y,t)是特征位移场控制方程应力载荷，由各相材料自定义热膨胀构造。是单位热应变。由于所述的单位热应变自身较大，容易导致求解复杂微观结构复合材料的特征位移场时发生单元畸变，最终产生计算错误。本发明在权利要求3所述的特征位移场控制方程两侧引入等比例缩小系数η，以实现复杂结构复合材料特征位移场的正确求解。引入缩小系数的特征位移场控制方程左侧为：

式中，为等比例缩小的特征位移场，为等比例缩小的特征位移场控制方程应力载荷；引入缩小系数的特征位移场控制方程右侧为：

式中，是等比例缩小的热应变，其表达式为：

式中，ψ为坐标转换矩阵；二维有限元模型的坐标转换矩阵为：

三维有限元模型的坐标转换矩阵为：

其中，li,mi,ni为局部坐标系各坐标轴与整体坐标系各坐标轴之间夹角的方向余弦。

式中，α^mn为自定义热膨胀系数，二维有限元模型的α^mn为：

其中，mn＝11,22,12；三维有限元模型的α^mn为：

其中，mn＝11,22,33,12,23,31。

式中，δt(t)＝h(t)为单位温度载荷，h(t)是阶梯函数。

此外，针对局部坐标与整体坐标一致的简单结构复合材料，计算过程可以省略等比例缩小载荷以及坐标转换，所需自定义热应变载荷为：

进一步，所述的步骤3.4完成特征位移场控制方程有限元求解。根据载荷缩小比例系数η对仿真结果进行修正，其表达式为：

对仿真位移场进行等比例修正后，获得特征位移场

对于应力载荷未缩小的情况，计算结果即为特征位移场。其表达式为：

进一步，所述的步骤3.5中通过计算单位热应变条件下随时间变化的体平均应力来计算等效粘弹性松弛刚度，其表达式为：

式中，为0阶应力。

进一步，利用载荷等比例缩小系数对计算结果进行修正，可得等效粘弹性松弛刚度的有限元列式为：

其中，为第n个增量步，时间为tn时松弛刚度，v^e代表第e个单元的体积，代表第e个单元的应力。

对于应力载荷未缩小的情况，有限元计算结果即为等效粘弹性松弛刚度，其表达式为：

进一步，对于二维模型，步骤3.4及3.5需重复三次以计算完整等效粘弹性松弛刚度；对于三维模型，步骤3.4及3.5需重复六次以计算完整等效粘弹性松弛刚度。

进一步，所述的步骤三中复合材料等效粘弹性采用与其树脂材料相同的时间-温度转换方程，以计算不同温度条件下等效粘弹性松弛刚度。

本发明与现有技术的有益效果在于：

本发明首先基于线性粘弹性本构方程以及时域渐进理论给出渐进均匀化特征位移场控制方程以及复合材料等效粘弹性松弛刚度计算表达式；在此基础上，本发明结合商业有限元软件开发了基于渐进均匀化法求解时域特征位移场和等效粘弹性松弛刚度的便捷实现方法。

本发明借助材料热应变构造时间、位置相关的特征位移场控制方程应力载荷以计算粘弹性的特征位移场，有效降低了渐进均匀化法的计算难度；相较于编写有限元程序计算复合材料等效粘弹性，借助商业有限元软件实现该均匀化方法可以满足任意微观结构复合材料的建模和计算需要；本发明解决了现有技术中基于时域渐进理论计算复合材料等效材料参数时实施过程复杂、效率低等问题，开发难度更低，应用潜力更大。

附图说明

图1为本发明基于时域渐进理论计算复合材料等效粘弹性松弛刚度有限元实现流程；

图2为本发明实施例1中编织结构复合材料代表性体积单元模型图；

图3为本发明实施例1中代表性体积单元时域上的特征位移场(法向变形)；

图4为本发明实施例1中代表性体积单元时域上的特征位移场(剪切变形)；

图5为本发明实施例1中编织结构复合材料等效粘弹性松弛刚度。

具体实施方式

为使本发明的目的、技术方案及效果更加清楚，明确，以下列举实例对本发明进一步详细说明。应当指出此处所描述的具体实施仅用以解释本发明，并不用于限定本发明。

如图1所示，图1为本发明基于渐进均匀化方法计算复合材料等效粘弹性松弛刚度有限元模拟的具体实施流程；本实施例的目的是计算编织复合材料(三维模型)的等效粘弹性松弛刚度；商业软件的选择包括：建模，定义周期性边界使用digimat，定义热应变载荷是matlab，求解与后处理使用abaqus；具体实施步骤如下：

步骤1、根据纤维束尺寸、纤维束间距等数据建立便编织结构复合材料代表性体积单元的三维模型，具体如图2所示，图示中1为树脂基体，2和3分别为经纱和纬纱；建模后，分别定义纤维束的粘弹性模量、树脂的粘弹性模量；划分有限元网格，获取有限元模型节点信息、单元信息。

步骤2-1、加载周期性边界条件；

对于三维编织结构复合材料的代表性体积单元，其周期性边界的表达式为：

其中，ui，yi代表位移和坐标，li是六面体代表性体积单元边界的长度，下标i＝{1,2,3}，上标mn＝11,22,33,12,23,31。

步骤2-2、加载周期性边界条件；

由于编织复合材料结构复杂，使用abaqus中‘equation’功能进行定义；该定义过程可由商业软件digimat自动定义实现。

步骤3、根据单元位置施加热应变构造求解特征位移场控制方程所需热应力载荷，其表达式为：

式中，为构造热应变，η＝0.001是载荷等比例缩小系数，ψ为局部坐标与整体坐标间的坐标系转换矩阵，对于三维编织复合材料，坐标转换矩阵为：

其中，li,mi,ni为局部坐标系各坐标轴与整体坐标系各坐标轴之间夹角的方向余弦。α^mn为自定义热膨胀系数，根据mn的变化，α^mn有：

式中，mn＝11,22,33,12,23,31。δt(t)＝h(t)是单位温度载荷。

步骤4、完成特征位移场控制方程有限元仿真，根据载荷比例系数计算编织复合材料时域上的特征位移场，表达式为：

步骤5、计算时域上编制结构复合材料等效松弛刚度，其计算表达式为；

其中，为第n个增量步，时间为tn时松弛刚度，v^e代表第e个单元的体积，代表第e个单元的应力。

进一步，根据参数mn的变化，重复步骤4和步骤5的计算共6次，以求解编织复合材料等效松弛刚度。

进一步，不同温度下的复合材料等效粘弹性刚度可通过利用树脂材料的时间-温度转换方程结合步骤4、步骤5的结果求解。

图3～4所示的分别为本发明实施例1中代表性体积单元时域上的特征位移场(法向变形)以及代表性体积单元时域上的特征位移场(剪切变形)。图5中，圆圈标志的曲线是基于体积平均(volumeaverage,va)方法的商业软件计算的编织结构复合材料等效粘弹性结果，三角标志的曲线是本发明实施例1中利用渐进均匀化(asymptotichomogenization,ah)方法计算的编织结构复合材料等效粘弹性松弛刚度。由于该编织结构复合材料为各向异性粘弹性材料，图中(a-i)包含了粘弹性松弛刚度中共九项数据：c11，c12，c13，c22，c23，c33，c44，c55，c66。对比图中各项结果，渐进均匀化法计算的复合材料等效粘弹性松弛刚度与体积平均方法计算的等效粘弹性结果一致。

以上所述仅是本发明的优选实施方式。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，这些改进也应视为本发明的保护范围。