小应变模式下的应变增量验证
假设一正方形截面弹性杆件（不计重量），长度为l=10m，截面尺寸为1m*1m，其弹性模量E=1e9Pa，泊松比v=0.2。在计算时，固定杆件左端，同时在右端施加一个大小为1000N、方向指向截面外的荷载。划分网格时，杆件横截面网格数量为10*10，长度方向的网格数量为100，网格见图1。
图1 模型图
2. 杆件变形理论值计算
材料力学中对于拉压杆的纵向变形、纵向线应变的计算公式分别如下：
(1)
（2）
将荷载FN=1000N、杆长l=10m、弹性模量E=1e9Pa、横截面积A=1*1=1m^2代入式（1）可计算得到纵向变形值delta_l=1e-5m，纵向线应变epsilon=1e-6。
3. 数值结果对比
3.1 FLAC3D 计算结果分析
计算命令流如下，计算结果见图2~图3。
;建立模型并计算model newzone create brick point 0(0,0,0) point 1(10,0,0) point 2(0,1,0) point 3(0,0,1) size 100 10 10zone cmodel assign elasticzone property young 1e9 poisson 0.2zone face skinzone face apply velocity-normal 0 range group 'West'zone face apply stress-normal 1e3 range group 'East'model solve;计算变形与应变的理论值fish def _dispCalE = 1e9L = 10.0A = 1.0*1.0F = 1e3delta_l = (F*L)/(E*A)strain = delta_l/Lend@_dispCal;计算网格点位移计算误差、单元应变增量极差以及单元应变增量上、下界误差fish def _errCalmax_strain_inc = 0.0min_strain_inc = 1e9max_disp_x = 0.0;获取最大x位移值loop foreach gpnt gp.listif gp.disp.x(gpnt) >= max_disp_x thenmax_disp_x = gp.disp.x(gpnt)endifendloop;获取最大/最小x应变增量loop foreach zpnt zone.listif zone.strain.inc.xx(zpnt) >= max_strain_inc thenmax_strain_inc = zone.strain.inc.xx(zpnt)endifif zone.strain.inc.xx(zpnt) <= min_strain_inc thenmin_strain_inc = zone.strain.inc.xx(zpnt)endifendloop;最大位移的计算误差err_disp = math.abs((max_disp_x-delta_l)/delta_l);单元x应变增量的极差range = max_strain_inc - min_strain_inc;单元x应变增量的上界误差err_ub = math.abs((max_strain_inc - strain)/strain);单元x应变增量的下界误差= math.abs((min_strain_inc - strain)/strain)end@_errCallist @err_displist @rangelist @err_ublist @err_lb
图2 模型x方向位移云图
图3 模型x方向应变增量云图
由图2可知，模型最大伸长量为1.0001e-5m，与理论计算值间的误差为1.248207039804774e-04，可忽略不计。
由图3可知，模型的最大、最小应变增量分别为9.9958e-7与1.0005e-6，极差为9.765875761028484e-10，因此可认为整个模型的x应变增量相等，其上、下界误差（绝对值）分别为5.512123256025262e-4、4.253752505003220e-4，其值极小，均可忽略不计。对于云图颜色不均匀的问题，可对Attributes中的云图插值间隙值进行调整，以使其颜色颜色均匀，见图4。由图4可知，在本例中系统自动设定的 contour interval为5e-11，其数量级远小于单元应变增量值的数量级，因此将其值放大为1e-8，调整后的应变增量云图见图5。
图4 云图interval调整界面
图5 调整后的x方向应变增量云图
3.2 ANSYS APDL计算结果分析
为简化计算流程，本次计算采用beam188单元进行分析，命令流如下，计算结果见图6~图8。
finish/clear/prep7et,1,beam188mp,ex,1,1e9mp,prxy,1,0.2sectype,1,beam,rect,beamsecdata,1,1,10,10k,1k,2,10l,1,2lesize,all,,,100lmesh,alldk,1,allfk,2,fx,1000/solusolve/post1/eshape,1/gline,,-1plnsol,u,xplnsol,epel,xplesol,epel,x
图6 纵向变形云图
图7 节点应变云图（Nodal strain）
图8 单元应变云图（Element strain）
由图6~图8可知，ANSYS的计算结果与理论计算值一致，因此，其与FLAC3D计算结果间的误差可忽略不计，三者的计算结果可认定为一致。
4 小结
