基于VB的ANSYS二次开发之Duncan-Chang本构模型算法介绍
所谓的本构模型是指材料的应力应变关系，比如遵循虎克定律的线弹性应力应变关系就是一种常用的本构模型。ANSYS为用户提供了许多本构模型，但在某些特殊领域，现有的本构模型却很少，完全不能满足分析要求。为了解决这个问题，ANSYS为用户提供了usermat等UPFs用户子程序，这些用户子程序拥有强大的二次开发功能，可以实现各种复杂本构模型的开发。但是，对于一些简单的本构模型，用户也可以利用APDL语言进行开发，比如Duncan-Chang本构模型。
Duncan-Chang本构模型介绍
对于应力应变关系复杂的材料而言，变形主要是由颗粒间位置的变化引起的，不同应力水平下相同的应力增量引起的应变增量并不相同，从而表现出复杂的非线性特征。为了反映材料特性，人们提出了许多本构模型，以期更好的反映材料的应力应变规律。
Duncan-Chang本构模型主要优点是可以利用常规三轴固结排水剪试验确定模型参数，因此在工程实践中得到了广泛应用。Duncan-Chang本构模型是非线性弹性模型，弹性矩阵中的弹性系数不再是常量，而是随应力状态而改变。由于不考虑塑性变形，所以一般只适用于载荷不大的情况（即不太接近破坏的情况）。Duncan-Chang模型有E-V和E-B模型两类。
当然，该模型也存在一些缺陷，如无法反应不同应力路径的影响、加卸载判断不明确等，不可避免造成了计算分析误差，长期以来许多学者试图来对其进行修正。有限元软件通常采用传统塑性理论的应力符号，以拉应力为正，下面对Duncan-Chang模型采用有限元软件的应力形式进行说明。
Duncan-Chang本构模型算法
该模型是Duncan和Chang根据Konder关于岩土材料的三轴试验的偏应力与轴向应变近似呈双曲线的假定而提出的。根据Konder的建对应变硬化的模型，其试验的应力应变可用双曲线关系来近似地描述，如图3-1所示，即在不变时：
式中a和b是试验常数。
上式也可以写成
其斜率为b,截距为a。
Duncan和Chang利用上述关系推导出了弹性模量公式：
可见增量虎克定律中所用的弹性模量实际是常规试验曲线的切线斜率。这样的模量叫切线弹性模量，可用来表示，
详细介绍见附件文档！！！！
生成并调用宏文件
在编程实现本构模型的过程中，需要重复执行某一部分，用户可以将该部分独立编写后放入指定位置，由APDL命令来调用，也可以通过*CREATE命令创建宏文件，并用*END命令结束宏的创建。利用*USE命令调用宏文件，并向宏文件传递参数：
*USE,Name,ARG1,ARG2,ARG3,ARG4,ARG5,ARG6,ARG7,ARG8,ARG9,AR10,AR11,AR12,AR13,AR14,AR15,AR16,AR17,AR18
其中，Name是宏文件名，ARGI到AR18是宏文件用到的参数值。我们将用*CREATE命令创建名为Duncan-Chang的宏文件，其中包含9个参数，使用*USE命令对模型内的每个单元反复调用Duncan-Chang宏文件，不断计算得到新的切线模量。
APDL实现过程
Duncan-Chang E-v模型是一种建立在增量广义虎克定律基础上的非线性变弹性模型，是通过不断改变其切线弹性模量来实现非线性的，完全可以通过ANSYS APDL进行编程分析。计算过程中主要通过如下方式来实现：取初始材料参数，施加第一步载荷，计算并读取单元应力，根据单元的当前应力调用Duncan-Chang模型宏命令计算新的材料参数（主要是材料的弹性模量和泊松比），代替初始材料参数；施加第二步载荷，计算并读取应力增量，根据单元的当前应力调用Duncan-Chang模型宏命令计算新的材料参数，以此类推。下面给出了 Duncan-Chang E-v模型创建宏文件，并用来模拟压缩试验的APDL详细代码(\chp6\Duncan-E-v.inp):
！（1）定义单元类型：
/PREP7
ET,1,185
！（2）设置材料属性
MPTEMP,,,,,,,,
MPTEMP,1,0
MPDATA,EX,1,,2e11
MPDATA,PRXY,1,,0.3
MPTEMP,,,,,,,,
MPTEMP,1,0
MPDATA,DENS,1,,7800
!3.建立几何模型
！（1）显示工作平面
WPSYYLE,,,,,,,,1
!(2)建立模型
BLOCK,0,10,0,20,0,2,
/REPLO
wpoff,3,5,0
/REPLO
CYLIND,1, ,0,10,0,360,
/REPLO
wpoff,0,10,0
CYLIND,1, ,0,10,0,360,
wpoff,0,-5,0
CYLIND,1, ,0,10,0,360,
wpoff,4,0,0
CYLIND,1, ,0,10,0,360,
wpoff,0,5,0
CYLIND,1, ,0,10,0,360,
wpoff,0,-10,0
CYLIND,1, ,0,10,0,360,
！（3）体布尔操作
FLST,3,6,6,ORDE,2
FITEM,3,2
FITEM,3,-7
VSBV, 1,P51X
！4.生成有限元模型
/REPLO
FINISH
/AUX12
FINISH
/PREP7
MSHAPE,1,3D
MSHKEY,0
CM,_Y,VOLU
VSEL, , , , 8
CM,_Y1,VOLU
CHKMSH,'VOLU'
CMSEL,S,_Y
VMESH,_Y1
CMDELE,_Y
CMDELE,_Y1
CMDELE,_Y2
！5.加载以及求解
！（1）加载
//REPLO
VPLOT
FINISH
/SOL
FLST,2,1,5,ORDE,1
FITEM,2,3
/GO
DA,P51X,ALL,0
/REPLO
FLST,2,1,4,ORDE,1
FITEM,2,7
/GO
FLST,2,1,5,ORDE,1
FITEM,2,4
/GO
!*
SFA,P51X,1,PRES,1000
!*
/STATUS,SOLU
SOLVE
(2) 创建Duncan-Chang计算
！
！***********
*dim,Str_max,array,120
*dim,S_max,ARRAY,120
*create,Duncan-Chang.mac !创建Duncan-Chang计算的宏命令
*afun，deg ！角度单位用度，不是弧度
*set，Pa，le5
*set，sigm1，-ArrS3(i) ！ArrS3为第三主应力，拉负压正
*set，sigm3，-ArrS3(i) ！ArrS3为第一主应力，拉负压正
*if，sigm3，LT,0.1*Pa，then
sigm3=0.1*Pa
*endif
str=2*(c*cos（Fai）+sigm3*sin（Fai））/(1-sin(Fai)) !破坏时的主应力差
S=(sigml-sigm3)/Slr !应力水平
*if,S.GT,0.95,then
S=0.95 !应力水平最大取值
*endif
!判断加卸绫荷，如果(sigml-sigm3)、S均小于历史最大值时视为卸载一再加载过程
*if,Str_max(i),gl,sigml-sigm3,and,S_max(i),gt,S,thcn Et=Kur*Pa*(sigm3/Pa)**n !卸载模置
*elseif，Str_max(i),gt,sigm1-sigm3，and，S_max(i),gt,S,then Ei=k*Pa*(sigm3/Pa)**n
Et=Ei*(l-Rf*S)**2 !加载情况的切线模量
S_max(i)=S ！保存历史最大应力水平
*elseif,Str_max(i),le,sigml-sigm3,and,S_max(i),gt,S,then Ei=k*Pa*(sigm3/Pa)**n
Et=Ei*(1-Rf*S)**2 !加载情况的切线模量
Str_max(i)=sigm1-sigm3 ！保存历史最大应力
*elseif，Str_max(i),gt,sigm1-sigm3，and，S_max(i),gt,S,then Ei=k*Pa*(sigm3/Pa)**n
Et=Ei*(l-Rf*S)**2 !加载情况的切线模量
Str_max(i)=sigm1-sigm3 !保存历史最大主应力差
S_max(i)=S ！保存历史最大应力水平
*endif
Ei=k*Pa*(sigm3/Pa)**n
A=(sigm 1 -sigm3)*D/(Ei*( 1 -Rf*S))
Mu=(G-F*loglO(sigm3/Pa)/(l-A)**2 ！计算泊松比
*if,Mu,GE,0.49,then
Mu=0.49
*endif
Mp,ex,I,Et !修改单元i的Et
Mp,nuxy,I,Mu
Mpchg,i,i !修改i单元的材料属性，mpchg,mat,elem
!
*vwrite,Et !格式化输出Et，用于观察中间计算结果
（E13.5）
*end !结束创建宏文件
！
！*****************************************
*dim ArrSl,array,120
*dim ArrS3,array,120
*do,j,2,9
!取出计算结果，修改弹性模*
/post1
etable,etabs 1,s,1 !取各单元第一主应力
etable,etabs 3,s,3 !取各单元第一主应力
*do,num,1,120 !num为单元编号
*get，Arrs1(num),elem,num,etab,etabs1 !将单元结果存入数组
*get，Arrs3(num),elem,num,etab,etabs3
*enddo
/prep7
*do,i,1,120,1 !各单元依次计算
c=130E3
Fai=31.2
Rf=0.9
k=538
n=0.35
Kur=645
G-0.33
F=0.033
D=3.2
!c一粘聚力，Fai—内摩擦角，Rf—破坏比，k—初始切线模最数，
!n—指数，Mu-泊松比， Kur—考虑卸载再加载时的模量数,G、F、D一参数
*use,Duncan-Chang.mac,c,Fai,Rf,k,n,Kur,G,F,D!调用Duncan-Chang宏文件
*enddo
！*******************************************
/solu
Time,j !定义第j载荷步
Nsubst,20,100,5 !指定子步数
Sfl,3,pres,3e5+le5*j !施加载荷，增量100kPa，最终上表面压力为1.2MPa
Outres,all,all
Solve
*enddo
参考文献
[1] 师访编.ANSYS二次开发及应用实例讲解[M].中国水利水电出版社.2012.1
[2] 隋丽娜，迟剑，郭立峰编. Visual Basic范例开发大全[M].清华大学出版社.
[3] 廖孟柯编. 基于VB的ANSYS二次开发与应用[J].
[4] 张艳.Visual Basic程序设计教程[M].徐州:中国矿业大学出版社,2001.
[5] 张晋西.用VB增强ANSYS前处理能力[J].计算机应用,2002,22(3).
[6] 龚曙光,谢桂兰. ANSYS操作命令与参数化编程[M].北京:机械工业出版社,2004.
![ANSYS/ABAQUS使用(带孔平板拉伸实例)[初识有限元CAE分析]](https://img.jishulink.com/cimage/0c683776c596fd2f7dde46fd773a666b_cdn.jpg?image_process=resize,fw_576,fh_320,)
