【FCA-Chp-6】第6章 计算微观（细观）力学（2）
6.3 局部-全局分析
在局部-全局分析（图6.8）中，对于全局模型中每个Gauss积分点上给定的应变，使用一个RVE计算应力。全局模型用于计算位移和产生的应变，假设材料是均质的。局部模型考虑不均质性，并通过一个RVE来对非均质材料建模。
图6.8 使用RVE的局部-全局分析
式6.13-6.15在6.2节中用于一次施加一个应变分量，目的是得到材料的等效弹性特性。式6.13-6.15对于施加在RVE的应变的一般状态仍然有效，但是必须注意在边或顶点处的周期性边界的说明。式6.13-6.15为9个约束方程，可以施加在RVE的面上所有成对的周期性的点之间，除了边和顶点处的点。
因为每条边属于两个面，在每条边上，可以看到2个位移分量都被用于施加2个约束方程，每个约束方程来自一个面，如同式6.42-6.44给出。然而，对于一个位移分量，只能写一个约束方程。因此，这些边必须分别处理。类似地，因为3个面交汇于一个顶点，来自3个面的3个周期性约束方程的施加则需要每次使用单个位移分量。根据类似于前面给出的推导，可得下面的结果。
四对顶点需要一次对一个分析。每对顶点关于RVE的中心（坐标为(0, 0, 0)）对称。得到的约束方程如下
几何（模块：Part）、材料和截面（模块：Property）、装配（模块：Assembly）、分析步（模块：Step）和网格（模块：Mesh）可以使用例6.2中使用的方法创建，但是要创建2D模型，RVE有1个完整的纤维和4个1/4纤维，如图6.4所示。在本例中没有对称边界条件可以使用。相反，周期性边界条件（PBC）式6.13-6.15必须在Load模块中施加。这些周期性边界条件通过使用约束方程实施。本例中没有载荷本身，而是由式6.13-6.15中等式右侧的常数项（由施加的应变
产生）对模型产生加载。
为了写出约束方程，需要定义一个参考点，以及主从节点。然后，主从节点的连接通过精确表示式6.13-6.15的约束方程实现。如前面注意的，顶点的方程与面的方程要分别处理。该工作由[5，PBC_2D.py]自动完成。
求解本例的完整脚本在[5，Ex_6.4-CE.py]可得。首先，该脚本创建网格和节点集 合，稍后用于约束方程。脚本可以在Abaqus/CAE中运行，具体操作为：菜单→File，Run script，[Ex_6.4-CE.py]，OK。该脚本调用[5，PBC_2D.py]写约束方程，调用[5，srecover2D.py]计算平均应力。结果在表6.4中给出。
表6.4 在RVE中的平均应力和基体上的最大应力（例6.4）
图6.9 应变的云图（在模型坐标系X-Y-Z中的E11，例6.4）
6.4 层压的RVE
在微尺度获得RVE所用的过程类似地可以用于在中尺度分析层压板。本例中RVE表示层压板。因此，厚度方向应该保持沿厚度自由膨胀。例如，平行于x-y平面的单层板，并且没有施加式6.15，因此，厚度坐标不受约束（见图6.10）。总之，RVE必须包含整个厚度。对于受面内载荷的对称层压板，RVE可以使用对称边界条件用一半厚度定义（见例6.5）。
图6.10 层压的RVE
和
例6.5 计算[0/90/-45/45]S层压板的Gxy，其中，单层板的材料参数为E1 = 139GPa、E2 = 14.5GPa、G12 = G13 = 5.86GPa、G23 = 2.93GPa、v12 = v13 = 0.21、E1 = 139GPa、v23 = 0.38，单层板厚度tk = 1.25mm。
解6.5 作为层压板对称的一个结果，可以使用层压板厚度一半的代表性体积单元RVE，并在z = 0处使用对称边界条件。RVE的模型使用3D实体单元离散。可以使用Abaqus CAE GUI按照下列伪代码建模。
i. 设置工作路径
菜单：File, Set Work Directory，[C:\Simulia\User\Ex_6.5]，OK
菜单：File，Save As，[C:\Simulia\User\Ex_6.5\Ex_6.5.cae]，OK
ii. 创建部件
模块：Part
菜单：Part，Create
3D，Deformable，Solid，Extrusion，Approx. size [10]，Cont
菜单：Add，Line，Rectangle，[0,0]，[1,1]，X，Done，Depth [5]，OK
菜单：Tools，Datum
Type：Plane，Method：Offset from principal plane
XY Plane，Offset [1.25]
XY Plane，Offset [2.5]
XY Plane，Offset [3.75]，X，# 关闭Create Datum对话框
菜单：Tools，Partition
Cell，Use datum plane，# 选取：z=1.25的基准面，
Create Partition
# 选取：较大实体，Done，# 选取：z=2.5的基准面，
Create Partition
# 选取：较大实体，Done，# 选取：z=3.75的基准面，
Create Partition
Done，# 关闭Create Partition对话框
iii. 定义材料、截面并赋截面给部件
模块：Property
菜单：Material，Create
Mechanical，Elasticity，Elastic，Type：Engineering Constants
[139000, 14500, 14500, 0.21, 0.21, 0.38, 5860, 5860, 2930]，OK
菜单：Assign，Material Orientation
# 选取：Z=4*1.25（自由表面）处的第1层，Done
Use Default Orientation or Other Method
Definition：Coordinate system,
Additional Rotation Direction：Axis 3
Additional Rotation：Angle [0]，OK
# 选取：Z=3*1.25处的第2层，Done
Use Default Orientation or Other Method
Definition：Coordinate system,
Additional Rotation Direction：Axis 3
Additional Rotation：Angle [90]，OK
# 选取：Z=2*1.25处的第3层，Done
Use Default Orientation or Other Method
Definition：Coordinate system,
Additional Rotation Direction：Axis 3
Additional Rotation：Angle [-45]，OK
# 选取：Z=1*1.25处的第4层，Done
Use Default Orientation or Other Method
Definition：Coordinate system,
Additional Rotation Direction：Axis 3
Additional Rotation：Angle [45]，OK，Done
菜单：Section，Create
Solid，Homogeneous，Cont，Material：Material-1，OK
菜单：Assign，Section，# 选取所有，Done，OK，Done
iv. 创建装配
模块：Assembly
菜单：Instance，Create，Independent，OK
v. 定义分析步
模块：Step
菜单：Step，Create
Procedure type：Linear perturbation，Static，Cont，OK
菜单：Output，Field Output Requests，Edit，F-Output-1
Output Variables：[S,E,U,IVOL]，OK
vi. 添加边界条件和载荷
# 施加简单剪切
模块：Load
菜单：BC，Manager
Create，Name [yfix]，Step：Initial，Disp/Rota，Cont
# 选取：y=0处的表面，Done，# 勾选：U1和U2，OK
Create，Name [disp]，Step：Step-1，Disp/Rota，Cont
# 选取：Y=1处的表面，Done
# 勾选：U1 [1]，# 勾选：U2 [0]，OK
Create，Name [zsymm]，Step：Initial，Symm/Anti/Enca，Cont
# 选取：Z=0处的表面，Done，ZSYMM，OK
# 关闭BC Manager对话框
vii. 定义约束
应用Tie约束x = 常数表面，以便得到。
模块：Interaction
菜单：Constraint，Create，Tie，Cont
Choose the master type：Surface，# 选取：表面x=0，Done
Choose the slave type：Surface，# 选取：表面x=1，Done
Position Tolerance：Specify distance [1]
# 不勾选：Adjust slave surface initial position # !!!
# 不勾选：Tie rotational DOFs if applicable，OK
viii. 对模型划分网格
模块：Mesh
菜单：Seed，Instance
Approx global size [1.25]，OK
# 对于纯剪，每层厚度一个单元足够
菜单：Mesh，Instance，Yes
ix. 求解并可视化结果
模块：Job
菜单：Job，Manager
Create，Cont，OK，Submit，# 当完成后，Results
菜单：File，Save
菜单：File，Run Script，[Ex_6.5-post.py]，OK
脚本[5，Ex_6.5-post.py]可以用于计算平均剪应力，以及剪切模量。因为施加的应变等于单位应变，计算得到的平均应力σ12等于C66。因此，在脚本最后一行得到G12 = 21,441MPa。
""" Post-processing Ex. 6.5 """
# 确认工作路径正确
import os
os.chdir(r’C:\Simulia\User\Chapter_6\Ex_6.5’)
# 打开当前作业的输出数据库（ODB）
from visualization import *
odb = openOdb(path=’Job-1.odb’);
myAssembly = odb.rootAssembly;
# 创建临时变量指向frame库，提高读取速度
frameRepository = odb.steps[’Step-1’].frames;
frameS=[];
frameIVOL=[];
# 在层压板方向（全局）创建坐标系
coordSys = odb.rootAssembly.DatumCsysByThreePoints(name=’CSYSLAMINATE’,
coordSysType=CARTESIAN, origin=(0,0,0),
point1=(1.0, 0.0, 0), point2=(0.0, 1.0, 0.0) )
# 将应力从单层板坐标系转换到定义的CSYSLAMINATE层压板坐标系
# stressTrans=odb.steps[’Step-1’].frames[-1].fieldOutputs[’S’]\
# .getTransformedField(datumCsys=coordSys)
stressTrans=frameRepository[-1].fieldOutputs[’S’]\
.getTransformedField(datumCsys=coordSys)
# 将转换后的应力添加到frameS
frameS.insert(0,stressTrans.getSubset(position=INTEGRATION_POINT));
frameIVOL.insert(0,frameRepository[-1].fieldOutputs[’IVOL’]\
.getSubset(position=INTEGRATION_POINT));
Tot_Vol=0; # 总体积
Tot_Stress=0; # 应力和
for II in range(0,len(frameS[-1].values)):
Tot_Vol=Tot_Vol frameIVOL[0].values[II].data;
Tot_Stress=Tot_Stress frameS[0].values[II].data * frameIVOL[0].values[II].data;
# 计算平均值
Avg_Stress=Tot_Stress/Tot_Vol;
print ’Abaqus/Standard Stress Tensor Order:’
# 根据Abaqus Analysis User’s Manual - 1.2.2 Conventions
print ’ 11-22-33-12-13-23’;
print Avg_Stress;
print ’G12=’,Avg_Stress[3]
练习题
问题6.3 使用问题6.1中计算得到的刚度矩阵，计算E1、E2、v12、v23、G12和G23。
问题6.4 使用问题6.2中计算得到的刚度矩阵，计算E1、E2、v12、v23、G12和G23。
问题6.5 展示例6.3中去掉Tie约束后导致的非零法 正应变和G13的错误结果。
问题6.6 用例6.3中的方法计算G12，但是使用对称边界条件，仅离散RVE的1/4。
问题6.7 通过修改例6.3的边界条件和Tie约束，计算G23。
问题6.8 使用例6.4作为指导，计算G23。
