2.3 算例

2.3.1 Mohr-Coulomb材料的三轴固结排水试验模拟

本例通过对三轴排水压缩试验的模拟介绍Mohr-Coulomb模型的使用，相应的cae文件名为ex2-1.cae。

1.问题描述

如图2-21所示，对一试样进行三轴固结排水（CD）剪切试验。三轴试验可分为两步，第一步是在固结应力σc的基础上增加围压增量Δσ3，使围压达到σ3；然后在竖直方向时间偏应力σ1−σ3，使得竖向应力达到σ1。若这两步中均允许排水，称为固结排水（CD）剪切试验。本例中试样的直径为 4cm，高度为 8cm。土体弹性模量E=10MPa，泊松比 v=0.3， σ3=100kPa，黏聚力c=10kPa，摩擦角ϕ=30°，求破坏时的σ1。

2.算例学习重点

Geostatic分析步的使用。

初始应力的设置。

Mohr-Coulomb材料的设置。

在ABAQUS后处理模块中处理曲线数据。

3.模型建立及求解

Step 1 建立部件。利用问题的对称性，将其简化为轴对称问题。启动 ABAQUS/CAE，在Part模块中执行【Part】/【Creat】命令，在Create Part对话框中选择Modeling space区域的Axisymmetric单选按钮，将Approximate Size 设为 0.32（为方便显示，一般设为模型最大尺寸的 2~4 倍），接受其余默认选项，单击【Continue】按钮后进入图形编辑界面，建立一个0.02m × 0.08m的矩形，完成后单击提示区中的【Done】按钮，确认退出。

Step 2 设置材料及截面特性。在Property 模块中，执行【Material】/【Creat】命令，建立名称为 soil的材料，在Edit Material对话框中执行【Mechanical】/【Elasticity】/【Elastic】命令设置弹性模型参数，继续执行【Mechanical】/【Plasticity】/【Mohr coulomb plasticity】命令定义Mohr-Coulomb模型参数，这里我们不设置剪胀角。

注意：

ABAQUS默认剪胀角最小为0.1°。

执行【Section】/【Create】命令，创建名为soil的均匀实体截面属性，对应的材料为soil，并执行【Assign】/【Section】命令赋给相应的区域。

提示：

本例中长度单位取m，应力单位取kPa，在设置材料弹性模量时的单位也为kPa。

Step 3 装配部件。在Assembly 模块中，执行【Instance】/【Create】命令，接受默认设置，建立相应的Instance。

Step 4 定义分析步。如前所述，三轴试验分两个步骤，即施加围压增量和偏应力增量。由于本例为CD试验，且σ3已知，因此，可将等压固结完成之后的应力状态作为初始状态，不考虑其在等压固结过程中的变形，在此基础上施加偏应力。为此，在Step模块中需设定两个分析步。一个是Geostatic分析步，ABAQUS在这一步中将判断用户定义的初始应力和对应荷载、边界条件之间是否平衡，如果平衡将不产生位移，以此模拟初始状态。第二个分析步是加载分析步。具体步骤如下：

执行【Step】/【Create】命令，在Create Step 对话框中将分析步名改为 Geoini，选择分析步类型为Geostatic，单击【Continue】按钮，接受Edit Step对话框中的默认选项。

再次执行【Step】/【Create】命令，在Create Step对话框中将分析步名改为Load，选择分析步类型为Static, General，单击【Continue】按钮，在Edit Step对话框的Basic选项卡中将Time period（时间总长）设为 1；在Incrementation 选项卡中将初始时间步长设为 0.01，最大时间步长设为 0.05；由于Mohr-Coulomb模型采用非相关联的塑性势面，在Other选项卡中将Matrix Storage选为Unsymmetric。

提示：

非线性分析步的起始步长取为总长的1%是一个比较保守的选择，这是为了保证收敛。这里的最大步长选得较小是为了获得足够的数据点。

Step 5 定义边界条件。在Load 模块中，执行【BC】/【Create】命令，在创建边界条件对话框中将分析步选为Geoini分析步，限定模型左侧的水平位移U1（轴对称边界）和模型底部的竖向位移U2。为模拟应变式三轴试验，再次执行【BC】/【Create】命令，在创建边界条件对话框中将分析步选为Load分析步，指定试样顶面的位移U2为−0.008（对应的竖向应变为10%）。

Step 6 定义荷载。在Load模块中，执行【Load】/【Create】命令，选择分析步Geoini作为荷载施加步，在土样的顶面和右侧施加压力荷载100kPa。这里的压力荷载和初始应力对应。

Step 7 设置初始应力。在Load模块中，执行【Predefined Field】/【Create】命令，弹出图2-22所示的创建预定义场对话框，将Step选为Initial（ABAQUS中的初始步），类型选为Mechanical，单击【Continue】按钮后按提示区中的提示选择整个土样作为施加初始应力的区域，确认后弹出编辑预定义场对话框，直接指定6个方向的初始应力，即−100、−100、−100、0、0、0，这里的负号是因为ABAQUS以拉为正，单击【Ok】按钮确认退出。

提示：

用户定义的初始应力一定要和Geostatic分析步中生效的荷载及边界条件对应。如果Geostatic分析步出现不收敛的情况，用户应首先检查这一点。

Step 8 划分网格。进入Mesh模块，将环境栏中的Object选项选为Part，意味着网格划分是在Part的层面上进行的。执行【Mesh】/【Controls】命令，在Mesh Controls对话框中选择Element shape（单元形状）为Quad（四边形），选择Technique（划分技术）为Structured。执行【Mesh】/【Element Type】命令，在Element Type对话框中，选择单元类型为轴对称应力单元CAX4。执行【Seed】/【Part】命令，将单元Approximate global size设为0.08，即本例中只划分1个单元。执行【Mesh】/【Part】命令，对网格进行划分。

提示：

由于本例中土样均匀受力，为简便起见，只划分了一个单元。

Step 9 建立任务。进入Job模块，执行【Job】/【Create】命令，建立名为ex2-1的任务并提交运算。

4.结果分析

Step 1 进入Visualization后处理模块并打开结果数据库。

Step 2 检查Geostatic分析步初始应力平衡的效果。执行【Resutl】/【Field out】命令，选择Step/Frame为Geostaic终止时的帧（本例中为Step：1，Frame：1），在Primary Variables选项卡中选择U（位移），在Component中选择Magnitude（大小）。执行【Plot】/【Contours】/【On Undeformed Shape】命令，绘出的位移大小如图2-23所示。可见位移很小，表示Geostatic中未产生位移，初始应力平衡较好。读者也可绘出Geostatic分析步结束时的应力云图，检查是否满足预期，即与设置的初始应力一致。

Step 3 执行【Resutl】/【Field out】命令，绘出Load加载分析步结束时的竖向应力S22（轴对称问题中2方向为竖向），如图2-24所示。可见最终竖向应力为 334.6kPa，这与极限平衡条件得到的σ1f =σ3 tan 2(45+ϕ/2)+2c tan(45+ϕ/2)=334.64kPa完全一致。

提示：

读者也可利用查询功能获取应力等相关信息。

Step 4 绘制竖向应力与应变的关系。执行【Tools】/【XY Data】/【Create】命令，选择Source（数据源）为ODB Field Output，单击【Continue】按钮后在图2-25所示的对话框的Variables选项卡中选择Position为Centroid，选择竖向应力S22、应变E11、E22和E33作为输出变量；用户可通过对话框右上方的Active Steps/Frames按钮控制输出结果的帧。切换到Elements/Nodes选项卡，在该选项卡中设置欲输出变量的单元，其有4种方法：

Pick from viewport：在屏幕上选择。

Element labels：输入单元编号选择。

Element sets：根据已定义的单元集合选择。

Internal sets：根据内部集合选择，内部集合是ABAQUS在用户操作时自动定义的集合。

这里通过直接选择的方法确定。确认后单击【Save】按钮保存数据，单击【Dismiss】按钮退出。

Step 5 已保存的XY数据可通过【Tools】/【XY Data】/【Manger】命令或者工具箱区中的按钮调出XY Data管理器进行重命名、绘制或编辑等操作。用户也可通过【Report】/【XY】命令，将数据导出，利用Excel等软件进行处理。相应工作也可在ABAQUS中进行。执行【Tools】/【XY Data】/【Create】命令，选择Source（数据源）为Operation on XY Data，单击【Continue】按钮后弹出图2-26所示的对话框。对话框的右侧提供了一系列运算函数，拖动活动条，找到Combine函数，单击后自动填入对话框上方的公式区。对话框下方的是已定义的XY Data名称，双击后自动填入公式区相应位置。按照图2-26设置，绘出竖向应力-S22与竖向应变-E22（负号是将对应数值取正，方便表达）的关系曲线，如图2-27所示。由图可见，在加载初期，材料处于弹性区，竖向应力随着竖向应变的增加而增加，直至屈服。由于没有指定材料的硬化，达到屈服之后，竖向应力不再随着竖向应变而变化。

提示：

在Operation on XY Data对话框中定义的数据可以利用【Save As】按钮保存。

Step 6 参照Step 5，利用Sum函数，得到体积应变（E11+E22+E33）随时间的变化关系。利用Combine函数，绘出体积应变与轴线应变的关系曲线，如图2-28所示（图中应变已取负号）。由图可见，在土体屈服前，竖向压缩应变比侧向膨胀大，体积呈压缩趋势。达到剪切屈服后，由于剪胀角为 0.1°，体积有膨胀的趋势。读者可自行改变剪胀角的大小，观察计算结果的异同。

2.3.2 修正剑桥模型材料的三轴固结排水试验模拟

本例的cae文件名为ex2-2.cae。

1.问题描述

有一正常固结土试样，在σ3=100kPa下固结稳定，即初始平均应力p0 =100kPa；土样的泊松比v=0.3， λ=0.15，κ=0.05，M =1.0，e1=2，由于其是正常固结土，则初始孔隙比e0 =e1−λln(p0)=1.309。求排水条件下破坏时的σ1。

提示：

如果没有等向压缩试验数据，可由侧限压缩试验结果估计参数，即λ=Cc/2.3，κ=Cs/2.3，其中Cc、Cs分别是压缩指数和回弹指数。但需注意，侧限压缩试验的压缩曲线并不是等向压缩曲线，其在e坐标轴上的截距并不是e1，读者可参考相关土力学书籍。另外，M =6sinϕ′/(3−sinϕ′)，ϕ′是临界状态摩擦角。

2.算例学习重点

修正剑桥模型的设置。

初始应力及初始孔隙比的设置。

超固结比对土体性状的影响。

3.模型建立及求解

Step 1模型的建立与算例2-1类似，这里直接在其基础上修改。将算例ex2-1.cae另存为ex2-2.cae。

Step 2 修改材料。进入Property模块，执行【Material】/【Edit】/【Soil】命令，或通过材料属性管理器对材料Soil进行编辑，在编辑材料对话框中利用按钮，删去原有的材料定义。然后分别按照图2-3和图2-20的设置，定义孔隙介质弹性模型和修正剑桥塑性模型。本例中通过定义Intercept的e1及当前初始孔隙比来定义初始屈服面位置，无须在Data数据列表中设置a0。或者，读者也可直接输入a0 =pc/2=p0/2=50kPa （正常固结土）。

Step 3 定义初始孔隙比。进入Load模块，执行【Predefined Field】/【Create】命令，弹出图2-29所示的创建预定义场对话框，将Step选为Initial（ABAQUS中的初始步），类型选为Other，单击【Continue】按钮后按提示区中的提示选择整个土样作为施加初始孔隙比的区域，确认后弹出编辑预定义场对话框，输入孔隙比为1.309，单击【OK】按钮确认退出。

Step 4 进入Job模块，执行【Job】/【Manger】命令，在编辑任务对话框中将任务名改为ex2-2，重新提交运算。

4.结果分析

Step 1 进入Visualization后处理模块并打开结果数据库。执行【Tools】/【XY Data】/【Create】命令，绘出a0（PEEQ）随时间的变化曲线，如图2-30所示。计算结果表明，在剪切过程中产生塑性体积应变，屈服面不断扩大，体现了材料的硬化。

Step 2 按照算例2-1中介绍的步骤，绘出轴向应力-轴向应变和体积应变-轴向应变的关系，分别如图2-31和图2-32所示。由图可见，其反映了典型的正常固结土的剪切性状，应力应变曲线呈应变硬化型，剪切过程中呈剪缩。同时可以发现，本算例中土样尚未达到临界状态，读者可自行改变土样顶面位移大小，观察计算结果的异同。

5.算例拓展

假设土体的超固结比OCR=4，其余参数不变。

Step 1 将ex2-2.cae另存为ex2-3.cae。

Step 2 重新定义孔隙比。此时土体为超固结土，则意味着土体处于回填曲线上，初始孔隙比发生了变化，由于OCR=4，初始平均应力 p 0 =100kPa，则前期固结应力 p c =400kPa，对应的初始孔隙比e0 =e1−λln(pc)+κln(OCR)=1.171。进入Load模块，执行【Predefined Field】/【Manger】命令，调出管理器，将初始孔隙比改为1.171。

Step 3 进入Job模块，执行【Job】/【Manger】命令，在编辑任务对话框中将任务名改为ex2-3，重新提交运算。

Step 4 进入Visualization后处理模块并打开结果数据库，绘制屈服面大小、轴向应力-轴向应变、体积应变-轴向应变关系曲线，分别如图2-30、图2-31和图2-32所示。由图可见，计算结果较好地反映了超固结土的剪切性状，即应力应变关系曲线为应变软化型，应力路径达屈服面后，产生剪胀，屈服面随着收缩。

2.3.3 黏弹性材料的循环剪切试验

本例cae文件名为ex2-4.cae。

1.问题描述

某材料的剪切行为符合图2-33所示的黏弹性模型，其中k1=1MPa，k2 =10MPa，η=10MPa⋅s；材料体积模量为K =100MPa，与时间无关。求该材料在周期剪应变γ=γ*sin(ωt)作用下的响应，本例中剪应变幅值γ*=0.05，圆频率ω=1，时间总长为10.0s。

2.算例学习重点

准静态分析步Visco的使用。

线黏弹性模型的设置。

施加随时间变化的荷载或边界条件。

3.模型建立及求解

Step 1 建立部件。为简便起见，本例对一单元体进行分析。启动 ABAQUS/CAE，在Part 模块中执行【Part】/【Creat】命令，基于拉伸建立一个1×1×1的三维变形实体。

Step 2 设置材料线黏弹性参数。在Property模块中，执行【Material】/【Creat】命令，建立名称为Soil的材料，在Edit Material对话框中执行【Mechanical】/【Elasticity】/【Viscoelastic】命令，在Domain下拉列表中选择Time后，选项卡上出现Time下拉列表，选择Prony（见图2-4）。本例中的线黏弹性模型，其剪切模量实为GR(t)=k1+k2e−t/τR ，τR =η/k2 =1s。对比式（2-16），这里在Prony级数中采用1项模拟，即。这里G0 =k1+k2 =11MPa，则=10/11=0.9091， =1s。在Data区的相应区域输入以上参数。注意因为体积模量参数与时间无关，k_i Prongy一栏保持空白。

Step 3 设置材料弹性参数。在Edit Material对话框中执行【Mechanical】/【Elasticity】/【Elastic】命令，这里定义的模量可以是Instantaneous（瞬时）或者Long Term（长期）模量（见图2-1）。本例设置瞬时模量。弹性模量E0 =9KG0/(3K+G0)=31.83MPa，泊松比v0 =(3K−2G0)/2(3K+G0)=0.45。

Step 4 执行【Section】/【Create】命令，创建名为soil的均匀实体截面属性，对应的材料为soil，并执行【Assign】/【Section】命令赋给相应的区域。

Step 5 装配部件。在Assembly 模块中，执行【Instance】/【Create】命令，接受默认设置，建立相应的Instance。

Step 6 定义分析步。执行【Step】/【Create】命令，在Create Step对话框中将分析步名改为Cyclic，选择分析步类型为Visco，单击【Continue】按钮，在Edit Step对话框的Basic选项卡中将Time period（时间总长）设为10；在Incrementation选项卡中将Type选为Fixed（固定时间步长），将时间分析步长（Increment size）设为0.05。由于时间总长为10，需要200步，因此需将允许的最大增量数（Maximum number of increments调大至200，确认后退出。

注意：

由于采用了线黏弹性材料，材料性能随时间变化，这里不能采用Static, general分析步，而需要采用Visco分析步。ABAQUS中的Visco实际上是拟静力分析步，其可反映材料性能随时间的变化，但忽略了惯性力效应。

Step 7 定义边界条件。在Load模块中，执行【BC】/【Create】命令，在Cyclic分析步中约束模型底部（z=0平面）的所有位移U1、U2和U3。

Step 8 定义循环加载幅值曲线。ABAQUS中荷载或边界条件可以为时间的函数，其通过Amplitude幅值曲线来实现。执行【Tools】/【Amplitude】/【Create】 命令，弹出图2-34所示的创建幅值对话框。该对话框提供了一系列创建幅值曲线的方法。如 Tabular 是通过表格的形式给出时间及对应的荷载（边界条件）大小，只要有足够的数据点，该方法就可给出任意形式的幅值曲线；Periodic则通过周期函数设置幅值大小a，其周期函数为傅立叶表达形式：

本例中施加的动位移条件正弦函数，取级数的一项即可，有A0 =0，N=1，ω=1，t0 =0，A1=0，B1=1。在创建幅值对话框中选择 Periodic 后单击【Continue】按钮，弹出编辑幅值对话框。该对话框的 Time Span下拉列表中有Step time（分析步时间）和Total time（总的时间），这里选择分析步时间。Circular frequency为圆频率ω，Starting time为t0，Initial amplitude为A0，A和B对应式（2-64）中的系数A和B，该区域有多少行，则n为多少。

Step 9 在Load模块中，执行【BC】/【Create】命令，选择Category为Mechanical，Type为Displacment后继续，设置模型顶面（z=1平面）的位移边界条件如图2-35所示，即指定U1为0.05，U2和U3为0，在Amplitude下拉列表中选择之前定义的幅值曲线Amp-1。由于模型为边长为1的正方体，按以上数据设置边界条件后，可使得剪应变为0.05sin(ωt)。

提示：

ABAQUS中实际施加的荷载（边界条件）大小为指定的绝对值与幅值曲线中定义值的乘积。

Step 10划分网格。进入Mesh模块，将环境栏中的Object选项选为Part，意味着网格划分是在Part的层面上进行的。执行【Mesh】/【Element Type】命令，在Element Type对话框中，选择单元类型为三维八节点（C3D8）。执行【Seed】/【Part】命令，将单元Approximate global size设为1，即本例中只划分1个单元。执行【Mesh】/【Part】命令，对网格进行划分。

Step 11建立任务。进入Job模块，执行【Job】/【Create】命令，建立名为ex2-4的任务并提交运算。

4.结果分析

Step 1 进入Visualization后处理模块并打开结果数据库。

Step 2 按照算例2-1中介绍的方法，绘出单元中心点S13和对应的应变E13的关系曲线，如图2-36所示。由图可见，其反映了典型的黏弹性材料在循环荷载下的响应，应力应变曲线为一椭圆，反映了能量的耗散。