【FCA-Chp-7】第7章 粘弹性（2）
7.6 粘弹性复合材料微观力学
7.6.1 一维情况
回顾Carson域内的本构方程式7.51。根据对应原理，对于弹性材料的所有微观力学方程在Carson域内对线性粘弹性材料有效。例如，Reuss微观力学模型假设在基体和纤维中的应变是均匀且完全相同的（6.1.1节的讨论）。因此，复合材料的刚度C是各成分（基体和纤维）刚度根据各自的体积分数Vm、Vf的加权组合
由此，刚度张量在Laplace域内是（见式7.47）
最后，通过进行Laplace逆变换式7.48可得时域内的刚度张量如下
解7.3 纤维的弹性行为和基体的粘弹性行为定义如下：
纤维（弹性）：
然后，Carson变换为
使用Reuss模型计算复合材料行为
回到Laplace域
回到时域（Laplace逆变换）
7.6.2 三维情况
1张量用加粗形式表示，或用带索引的分量表示。
使用卷积定理（表7.1），式7.65的Laplace变换为
或用Carson变换为
假设纤维是弹性的，基体是粘弹性的，后者用Maxwell模型表示
其Carson变换为
使用对应原理得到
使用式1.75并假设泊松比为常数，在Carson域内基体的拉梅常数为
基体的剪切模量为
式中
系数
用于解释微观结构的几何，包括微观结构的包含物及其几何排列[36]。对于按照方阵排列的圆柱形纤维[37]，则有
注意，式7.73产生6个独立的弛豫张量元素。这是因为式7.73表示的复合材料的微观结构按照方阵排列。如果微观结构是随机的（图1.12），复合材料是横观各向同性的（见1.12.4节），则弛豫张量仅有5个独立的元素。当轴x1为复合材料横观各向同性的轴时，平均化过程（式6.7）产生具有横观各向同性的弛豫张量如下
基于[50]的MATLAB程序可以从[5]获得，用于数值上的Laplace逆变换。另一个算法在[8，附录D]中提供。
i. 从实现式7.73~7.76的部分代码得到的输出是弛豫模量在Carson域内关于s的方程。注意，必须声明作为符号的变量s。
ii. 用得到的方程除以s，回到Laplace域。
iii. 使用函数invlapFEAcomp变回到时域，该函数来源于[50]。
iv. 最后，使用粘弹性模型方程拟合E2(t)的数值。通常，基体弛豫使用的模型方程可以方便地用于复合材料弛豫；在本例中，使用Maxwell模型。该步操作由fitfunFEAcomp.m实现。
MATLAB代码PMMViscoMatrix.m、invlapFEAcomp.m、fitfunFEAcomp.m可从[5]获得。结果在图7.4中给出。对于复合材料的Maxwell参数的完整集 合在例7.7中计算。
7.7 粘弹性复合材料宏观力学
7.7.1 均衡对称层压板
均衡对称层压板的面内粘弹性行为可以使用1.15节（表观层压板特性）的过程得到，但是在Carson域内。在Carson域内，从单层板坐标系下的单层板的刚度开始。将每个单层板的矩阵转动到层压板坐标系下。然后，使用式1.102对它们求平均值。使用式1.05，得到在Carson域内的层压板工程特性参数，除以s回到Laplace域。最后，进行Laplace逆变换，得到时域内的层压板刚度。然后，使用一个模型方程拟合它们，如同例7.4。
7.7.2 一般层压板
由于对应原理，从经典层压板理论（CLT，见第3章）而来的应力-合力与应变-曲率方程在Carson域内对于线性粘弹性层压复合材料仍然有效。在Carson域内的层压板的A、B、D、H矩阵可以使用从一阶剪切变形理论（FSDT，3.1.1节）而来的方程进行计算。这种方法 论在[54]中使用过。
7.8 粘弹性复合材料的有限元分析
大多数商业软件已经实现了各向同性材料的粘弹性（蠕变）。这严重限制了对聚合物基体复合材料的粘弹性行为分析感兴趣的用户。
然而，可以利用商用软件中用户可编程的特征，实现本章中出现的公式。这样相对简单，因为在本章中使用的方法不是依赖于应力的，但是线性粘弹性方法及其实现并不复杂。在例7.7中使用UMAT子程序实现粘弹性公式。
解7.5 模型的大部分可以通过修改例5.2的模型得到，具体过程如下伪代码：
i. 找到例5.2的cae文件
菜单：File，Set Work Directory，[C:\SIMULIA\User\Ex_7.5]，OK
菜单：File，Open，[C:\SIMULIA\User\Ex_5.2\Ex_5.2.cae]，OK
菜单：File，Save As，[C:\SIMULIA\User\Ex_7.5\Ex_7.5.cae]，OK
ii. 修改部件
模块：Part
# 编辑实体拉伸，改变90°铺层厚度
# 在模型树，# 展开：Parts(1)，Part-1，Features
# 右单击：Solid extrude，Edit，Depth [11.25]，OK
# 移动分割面
# 右单击：Datum plane-1，Edit，Offset [10]，OK
iii. 修改材料和截面特性。在例5.2中定义为Material-1的材料必须删除，因为在Abaqus中粘性分析只接受各向同性材料。
模块：Property
菜单：Material，Manager
# 选择：Material-1，Delete，Yes
# 创建各向同性粘弹性材料
Create，Name [iso-visco]，Mechanical，Elasticity，Elastic
Type：Isotropic，[11975, 0.1866]
Moduli time scale：Instantaneous
# 粘性部分
Mechanical，Elasticity，Viscoelastic，Domain：Time，Time：Prony
g_1 [0.999]，k_1 [0.999]，tau_1 [58.242]，OK
# 关闭Material Manager对话框
菜单：Section，Edit，Section-1，Material：iso-visco，OK
iv. 修改分析步
模块：Step
菜单：Step，Manager
# 选择：Step-1，Delete，Yes，# 这样也会删掉已创建的边界条件
Create，Name [Step-1]，Procedure type：General，Visco，Cont
标签页：Basic，Time period [150]
标签页：Incrementation，Type：Fixed，Max number of increments [200]
Increment size [1]，OK，# 关闭Step Manager对话框
v. 定义约束和相互作用
模块：Interaction
# 将X=40的面设置为刚性面，用于施加力
菜单：Tools，Reference Point，[40,10,11.25]，X # RP-1
菜单：Constraint，Create
Name [Constraint-1]，Type：Coupling，Cont
# 选取：RP-1 # 常控制点
Surface，# 选取X=40处的面，Done
Constrained degrees of freedom：U1 # 不勾选：U1以外的其他所有自由度，OK
vi. 创建一个集 合用于历史输出
模块：Step，Model：Model-1，Step：Step-1 # Step-1 must be selected
菜单：Tools，Set，Create
Name [Set-1]，Cont，# 选取：RP-1，Done
菜单：Output，History Output Requests，Edit，H-Output-1
Domain：Set ：Set-1，Output Variables：选取下面列出的变量
# 展开：Stresses，# 勾选：S # 确认选择所有项
# 展开：Strains，# 勾选：E # 确认选择所有项
# 展开：Displacement/Velocity/Acceleration，# 勾选：U，OK
vii. 添加载荷和边界条件
模块：Load
# 约束参考点的转动UR2，UR3
菜单：BC，Create
Name [BC-RP-1]，Step：Initial，Type：Disp/Rota，Cont
# 选取：RP-1，Done，# 勾选：UR2，UR3，OK
# 在参考点施加载荷
菜单：Load，Create
Step：Step-1，Mechanical，Concentrated force，Cont
# 选取：RP-1，Done，CF1 [1000.0]，OK
viii. 对模型重划分网格
模块：Mesh
菜单：Mesh，Instance，Yes
ix. 求解和可视化结果
模块：Job
菜单：Job，Manager
Create，Cont，OK，Submit，# 当完成后，Results
模块：Visualization
# 绘制位移-时间曲线
菜单：Result，History Output，U1，Plot，# 关闭对话框
菜单：Plot，Contours，On Deformed Shape
工具条：Field Output，Primary，U，U1
菜单：Result，Step/Frame，Index：0，Apply
菜单：Result，Step/Frame，Index：1，Apply
菜单：Result，Step/Frame，Index：150，Apply
# 类似地显示E11和S11结果
结果在表7.3中给出。注意，在例5.2中应力和应变云图显示的是局部坐标系下的结果，因此表7.3中所示的结果是从外侧单层板提取的。这样做的目的是强调每个单层板的云图显示的是在局部坐标系CSYS下的结果。
表7.3 结果（例7.5）
解7.6 模型的大部分可以通过修改例7.5的模型而得到，具体操作如下面的伪代码所示。注意，初始位移在短时间内（1min）使用自动增量的方式加载。然后，第二个分析步设置为计算150min的弛豫，使用固定间隔以便数据解读。
i. 找到例7.5的cae文件
菜单：Set Work Directory，[C:\SIMULIA\User\Ex_7.6]，OK
菜单：File，Open，[C:\SIMULIA\User\Ex_7.5\Ex_7.5.cae]，OK
菜单：File，Save As，[C:\SIMULIA\User\Ex_7.6\Ex_7.6.cae]，OK
ii. 修改分析步
模块：Step
菜单：Step，Edit，Step-1
标签页：Basic，Time period [1]
标签页：Incrementation，Type：Automatic
Creep/swelling/viscoelastic strain error tolerance：[1E-6]，OK
菜单：Step，Create
Name [Step-2]，Procedure type：General，Visco，Cont
标签页：Basic，Time period [150]
标签页：Incrementation，Type：Fixed，Max number of increments [200]
Increment Size [1]，OK
iii. 修改施加的载荷和边界条件
模块：Load
# 删除RP-1上的载荷
菜单：Load，Delete，Load-1，Yes
# 在RP-1上定义一个位移
菜单：BC，Create，Step：Step-1，Disp/Rota，Cont
# 选取：RP-1，Done，# 勾选：U1 [4.0]，OK
iv. 求解和可视化结果
模块：Job
# 创建新的分析作业。不要覆盖蠕变算例的结果。
菜单：Job，Manager
Create，Cont，OK，Submit，# 当完成后，Results
模块：Visualization
# 绘制时间-位移曲线
菜单：Result，History Output，U1，Plot，# 关闭对话框
# 可视化不同时刻的结果，与例7.5中的操作相同。
结果在表7.4中给出。
表7.4 结果（例7.6）
表7.5 单层板粘弹性特性（例7.7用）
解7.7 单层板的Maxwell参数在表7.5中给出。在Abaqus中使用UMAT子程序和下列时间相关的特性，可以实现正交各向异性材料的本构方程。
子程序可从[5，umat3dvisco.for]获得。模型其他部分的建立可以使用与例7.6相似的过程。
i. 找到例7.6的cae文件，并设置工作路径
菜单：Set Work Directory，[C:\SIMULIA\User\Ex_7.7]，OK
菜单：File，Open，[C:\SIMULIA\User\Ex_7.6\Ex_7.6.cae]，OK
菜单：File，Save As，[C:\SIMULIA\User\Ex_7.7\Ex_7.7.cae]，OK
ii. 创建使用UMAT的材料
模块：Property
菜单：Material，Create
Name [user]，General，User Material，Type：Mechanical
[102417 11975 0.401 0.1886 5553.8 5037.3 16551 58.424 44.379 54.445]
# E1o，E2o，nu12o，nu23o，G12o，G23o，tau1，tau2，tau12，tau23
# 作为一列向量输入
General，Depvar，
Number of solution-dependent state variables [6]，OK
iii. 修改分析步
模块：Step
菜单：Step，Manager
# 选取：Step-1，Edit，标签页：Basic，Time period [0.001]
标签页：Incrementation，Type：Automatic，Error tolerance [1E-6]，OK
# 保持Step-2不变，# 关闭Step Manager对话框
分析步：Step-1，# 在工作区上方的环境栏选择
菜单：Output，Field Output Requests，Edit，F-Output-1
Output Variables：Edit variables [S,E,U,RF,SDV]，OK
菜单：Output，History Output Requests，Edit，H-Output-1
Output Variables：Edit variables [S11,E11,U1,RF1,SDV]，OK
iv. 将单元类型修改为二阶单元，不使用沙漏刚度
模块：Mesh
菜单：Mesh，Element Type，# select all，Done
Family：3D Stress，Geometric Order：Quadratic
# 勾选：Reduced integration，OK
菜单：Mesh，Instance，Yes
v. 求解和可视化结果
模块：Job
# 创建第3个分析作业，但不覆盖前面的结果
菜单：Job，Manager
Create，Cont，标签页：General，User subroutine file [umat3dvisco.for]，OK
Submit，# 当完成后，Results
模块：Visualization
# 可视化支反力，用于计算层压板刚度
工具条：Field Output，Symbol，RF，RESULTANT
# 保存数据到Excel文件
# 在左侧树
# 展开：Output Databases，# 展开：Job-4.odb (使用的Job名称)
# 展开：History Output，# 右单击：Reaction Force
Save As，Name [XYData-RF1]，OK
# 展开：XYData，# 右单击：XYData-RF1，Plot
# 右单击：XYData-RF1，Edit，# 选择所有数据
# 右单击：Copy，# 粘贴到Excel文件
# 关闭对话框
结果在表7.6中给出。
表7.6 随时间变化的平均应力σ4
在Abaqus中，分析被划分为分析步s = 1,...,ns。进一步，每个分析步通常又划分为增量步i = 1,...,ni。Abaqus记录从分析开始经过的（总）时间，包括重启动分析步。另外，Abaqus还记录从分析步开始经过的时间，称为分析步时间ts。
关于图7.5，当前增量步（在当前分析步内）开始时的分析步时间ts传递到用户子程序中的变量TIME(1)。当前增量步开始时的（总）时间是TIME(2)。当前时间增量ΔT是DTIME。
图7.5 分析步时间ts、总时间t和时间增量ΔT的定义
这样，当前分析步时间可以计算如下
当前（总）时间可以计算如下
图7.6 在例7.7中计算得到的层压板刚度
解7.8 在Abaqus中，可以使用弹性特性表示纤维，使用粘弹性特性表示基体。模型定义与例6.4很相似。有两种方法求解本例：一种是修改Python脚本Ex_6.4-CE.py；另一种是修改Ex_6.4-CE.cae。下面介绍第一种方法。
首先，将PBC_2D.py和srecover2D.py复制到本地路径。然后，将Ex_6.4-CE.py复制为Ex_7.8.py，并使用Notepad 打开。通过修改Ex_7.8.py中相应的行，将工作路径更新到当前路径C:\SIMULIA\User\Ex_7.8，如下：
os.chdir(r’C:\SIMULIA\User\Ex_7.8’)
然后，修改材料特性，使其包含基体所有的粘弹性参数。弹性模量的单位为TPa（即1000GPa），因为模型尺寸使用微米为单位。
Ef, nuf = 168.4E-3, 0.443 # TPa
Em, num, tau, g_1, k_1 = 4.082E-3, 0.311, 39.15, 0.999, 0.999 # TPa,,min,,,
strain = [0.000, 0.0, 0.01] # epsilon_11，epslion_22，gamma_12
下一步，需要新定义粘弹性基体材料。纤维保持为弹性。为了得到粘弹性材料的Python脚本，打开Abaqus/CAE，（空的模型）保存为Ex_7.8.cae。然后，按照例7.5的操作定义粘弹性材料。在Abaqus/CAE中详细操作的伪代码如下：
模块：Property
菜单：Material，Create
Name [matrix]，Mechanical，Elasticity，Elastic，Type：Isotropic
Moduli time scale：Instantaneous，Data [4.082E-3, 0.311]
Mechanical，Elasticity，Viscoelastic，Domain：Time，Time：Prony
g_1，k_1，tau_1 [0.999, 0.999, 39.15]，OK
保存模型之前，从Ex_7.8.rec中复制下面的命令行到Ex_7.8.py中，替换在# Materials程序段中定义基体材料特性的两行。
mdb.models[’Model-1’].Material(name=’matrix’)
mdb.models[’Model-1’].materials[’matrix’].Elastic(moduli=INSTANTANEOUS,
table=((0.004082, 0.311), ))
mdb.models[’Model-1’].materials[’matrix’].Viscoelastic(domain=TIME,
table=((0.999, 0.999, 39.15), ), time=PRONY)
并且，立即使用变量名替换数值（用于脚本参数化），如下：
mdb.models[’Model-1’].Material(name=’matrix’)
mdb.models[’Model-1’].materials[’matrix’].Elastic(moduli=INSTANTANEOUS,
table=((Em, num), ))
mdb.models[’Model-1’].materials[’matrix’].Viscoelastic(domain=TIME,
table=((g_1, k_1, tau), ), time=PRONY)
现在，保存模型以重置rec文件。创建一个粘弹性分析步Step-1，用于施加应变，创建分析步Step-2，用于追踪接着发生的弛豫，该操作与例7.6相同。在Abaqus中详细操作的伪代码如下：
模块：Step
菜单：Step，Create
Name [Step-1]，Procedure type：General，Visco，Cont
标签页：Basic，Time period [0.001]
标签页：Incrementation，Type：Automatic，viscoelastic tolerance [1E-6]，OK
菜单：Step，Create
Name [Step-2]，Procedure type：General，Visco，Cont
标签页：Basic，Time period [100]
标签页：Incrementation，Type：Fixed，Maximum number of increments：[200]
Increment Size [1]，OK
保存模型之前，从Ex_7.8.rec中复制以下命令行到Ex_7.8.py，替换# Step程序段中定义Step-1的命令行。
mdb.models[’Model-1’].ViscoStep(cetol=1e-06, initialInc=0.001, maxInc=0.001,
minInc=1e-08, name=’Step-1’, previous=’Initial’, timePeriod=0.001)
mdb.models[’Model-1’].ViscoStep(cetol=0.0, initialInc=1.0, maxNumInc=200,
name=’Step-2’, previous=’Step-1’, timeIncrementationMethod=FIXED,
timePeriod=100.0)
下一步，更新文件名称以便添加边界条件到模型。
execfile(’C:/SIMULIA/User/Ex_7.8/PBC_2D.py’)
下一步，更新下面一行中的文件名称。
mdb.saveAs(pathName=’C:/SIMULIA/User/Ex_7.8/Ex_7.8.cae’)
下一步，添加下面注释行。
# Calculate Stresses and Strains
# execfile(’srecover2D.py’)
# visualize
# o3 = session.openOdb(name=’Job-1.odb’)
# session.viewports[’Viewport：1’].setValues(displayedObject=o3)
此时，更新后的脚本应该看起来像从[5，Ex_7.8.py]获得的脚本。下一步，在Abaqus/CAE新建一个模型，然后运行脚本，操作如下：
菜单：File，New Model Database，With Standard/Explicit Model
菜单：File，Run Script [Ex_7.8.py]，OK
可视化结果，操作如下：
模块：Job
菜单：Job，Results，Job-1
模块：Visualization
菜单：Plot，Contours，On Deformed Shape
菜单：Result，Step/Frame
注意，如何打开图7.7所示的窗口，单击菜单：Result，Step/Frame，打开Step/Frame对话框。现在，可以选择分析步和增量步。在本例中，Step-1是加载分析步，Step-2是弛豫分析步。因为在Step-2中使用了固定增量步，增量步编号与时间（分钟）对应。对于时间相关的分析，例如本例，帧（frame）和增量步（increment）含义相同的。
图7.7 用于在输出数据库中选择分析步和帧（增量步）的Step/Frame窗口
i. 选择分析步。在本例中，Step-1是加载分析步，Step-2是弛豫分析步。在脚本中选择分析步如下：
# In Ex. 7.8, Step-2 is the relaxation step
frameRepository = odb.steps[’Step-2’].frames;
ii. 脚本参数化选择存储增量步的帧：
# Get the results for frame [i], where i is the increment number
i = 0
frameS = frameRepository[i].fieldOutputs[’S’].values;
frameE = frameRepository[i].fieldOutputs[’E’].values;
frameIVOL = frameRepository[i].fieldOutputs[’IVOL’].values;
剩余的脚本保存不变。完整的脚本如下（也可从[5，srecover2D.py]获取）。
# srecover2D.py modified for Ex. 7.8
from visualization import *
# Open the output data base for the current Job
odb = openOdb(path=’Job-1.odb’);
myAssembly = odb.rootAssembly;
# Temporary variable to hold the frame repository speeds up the process
# In Ex. 7.8, Step-2 is the relaxation step
frameRepository = odb.steps[’Step-2’].frames;
# Get the results for frame [i], where i is the increment number
i = 80
frameS = frameRepository[i].fieldOutputs[’S’].values;
frameE = frameRepository[i].fieldOutputs[’E’].values;
frameIVOL = frameRepository[i].fieldOutputs[’IVOL’].values;
Tot_Vol=0.; # Total Volume
Tot_Stress=0.; # Stress Sum
Tot_Strain = 0.; # Strain Sum
# Calculate Average
for II in range(0,len(frameS)):
Tot_Vol =frameIVOL[II].data;
Tot_Stress =frameS[II].data * frameIVOL[II].data;
Tot_Strain =frameE[II].data * frameIVOL[II].data;
Avg_Stress = Tot_Stress/Tot_Vol;
Avg_Strain = Tot_Strain/Tot_Vol;
# from Abaqus Analysis User’s Manual - 1.2.2 Conventions -
# Convention used for stress and strain components
print ’2D Abaqus/Standard Stress Tensor Order：11-22-33-12’
print ’Average stresses Global CSYS：11-22-33-12’;
print Avg_Stress;
print ’Average strain Global CSYS：11-22-33-12’;
print Avg_Strain;
odb.close()
图7.8 时间相关的剪切模量
练习题
