RVE周期性边界条件施加

记录一下学习周期性边界条件施加的迷思。

理论

代表性体积单元（RVE），指的是在多尺度分析中可使宏观分析与微观分析解耦的细观模型，满足宏观上充分小，微观上充分大的尺寸要求。
通过对RVE施加周期性边界条件，可以模拟具有周期结构的RVE宏观尺度载荷工况，大大减少计算成本。
对于复合材料这种一般各向异性材料，我们有如下公式进行刚度计算：
$$\overline{\sigma }=\overline{C}\overline{\epsilon }$$
$$\overline{S}={\overline{C}}^{-1}$$
在复合材料等效材料属性的测试中，一般需要在单轴应力/应变状态下进行分析。下面在单轴应变状态下进行分析。
周期性边界条件顾名思义，指以单元为周期应变应具有规律。建立沿x,y,z方向边长分别为a,b,c的立方体，可对6个面进行如下约束：
$$u{|}_{x=a}-u{|}_{x=0}=a{\epsilon }_{x}u{|}_{y=b}-u{|}_{y=0}=b{\gamma }_{yx}u{|}_{z=c}-u{|}_{z=0}=c{\gamma }_{yx}$$

$$v{|}_{x=a}-v{|}_{x=0}=a{\gamma }_{yx}v{|}_{y=b}-v{|}_{y=0}=b{\epsilon }_{y}v{|}_{z=c}-v{|}_{z=0}=c{\gamma }_{zy}\text{\hspace{1em}(1)}$$

$$w{|}_{x=a}-w{|}_{x=0}=a{\gamma }_{xz}w{|}_{y=b}-w{|}_{y=0}=b{\gamma }_{yz}w{|}_{z=a}-w{|}_{z=0}=c{\epsilon }_z$$
其中 { ε x , ε y , ε z , γ x y , γ x z , γ y z } \{\varepsilon_x,\varepsilon_y,\varepsilon_z,\gamma_{xy},\gamma_{xz},\gamma_{yz}\} {εx​,εy​,εz​,γxy​,γxz​,γyz​}即为单元的应变状态，在计算等效模量时一般设置为1即可。
打个比方，如果需要计算$E_{11}$，应该施加${\epsilon }_x=1$，其他皆为0的应变。其他模量则同理。注意工程剪应变是剪应变的两倍，所以应设置为0.5。

ABAQUS实现

在ABAQUS中有一个插件叫Micromechanics Plugin可以方便的实现周期性边界条件的施加，施加原理同上面是一样的。

Micromechanics Plugin使用

下载插件压缩包，解压后复制到D:…\ABAQUS2023\plugins，打开ABAQUS即可使用。使用教程详见Abaqus FE-RVE（fe-rve）微观力学插件介绍Micromechanics Plug-in
需要事先建立好模型，完成除边界条件、约束、载荷之外一切工作并建立好job，点击插件loading设定好选项即可计算。
成功的话就不需要任何操作了，下面可以当做插件的讲解。

约束方程施加

第一部分提到的(1)式我们通过约束方程施加。
首先需要建立每个面所有结点的set，需要注意的是，每个边的结点不要重复添加进其他的面，故有一些面会少一行的结点。
在相互作用/约束管理器/创建中创建方程，设置参数如下
与$u{|}_{x=a}-u{|}_{x=0}=a{\epsilon }_x$等效，其中80是我沿着x方向的边长，前两项自由度为1表示沿x方向的位移u，第三项自由度序号则表示向量RP_Normal { ε x , ε y , ε z } \{\varepsilon_x,\varepsilon_y,\varepsilon_z\} {εx​,εy​,εz​}的第一项。该插件在坐标原点的结点设置了两个向量来施加线应变和剪应变。
类似的，可以添加其余8个约束方程。

载荷工况施加

计算每个模量的载荷工况应该单独施加，可以通过施加载荷工况满足这个要求。
以计算$E_{11}$为例，需要设置两个位移/转角边界条件来完成应变的施加，分别是RP_Normal { ε x , ε y , ε z } = { 1 , 0 , 0 } \{\varepsilon_x,\varepsilon_y,\varepsilon_z\}=\{1,0,0\} {εx​,εy​,εz​}={1,0,0}和RP_Shear { γ x y , γ x z , γ y z } = { 0 , 0 , 0 } \{\gamma_{xy},\gamma_{xz},\gamma_{yz}\}=\{0,0,0\} {γxy​,γxz​,γyz​}={0,0,0}。在载荷工况中添加$E_{11}$，边界条件添加刚创建的这两项，无需添加载荷。其他五个工况也与此类似。
如图所示是RP_Normal { ε x , ε y , ε z } = { 1 , 0 , 0 } \{\varepsilon_x,\varepsilon_y,\varepsilon_z\}=\{1,0,0\} {εx​,εy​,εz​}={1,0,0}边界条件。

注意到最下面的边界条件施加在initial分析步，点开发现固定了立方体最中间的节点三个方向的位移为0，施加对称/反对称/固定约束，以免出现单元整体平移的情况。
至此就完成所有边界条件的设置了。计算完成后可以得到6个工况下的变形，可以查看进行检查。
最后贴一个需要创建的集。中间字打错了，是去除重合边来着。
介绍就到这里。希望对你有所帮助！