一个FLAC的本构模型

ITASCA(武汉)咨询有限公司依泰斯卡(武汉)内容必要性试验总结的本构模型特定条件下的本构模型交叉学科的本构模型二次开发环境自定义本构模型的功能自定义本构模型的基本方法常用模型信息传递指针变量模型源程序分解二次开发环境FLAC3D采用面向对象的语言标准C++编写本构模型都是以动态连接库文件(.DLL文件)的形式提供

VC++6.0(SP4)或更高版本的开发环境优点自定义的本构模型和软件自带的本构模型的执行效率处在同一个水平自定义本构模型(.DLL文件)适用于高版本的FLAC(2D)、3DEC、UDEC等其他Itasca软件中返回目录自定义本构模型的功能

主要功能：对给出的应变增量得到新的应力辅助功能：模型名称、版本读写操作

模型文件的编写基类(classConstitutiveModel)的描述成员函数的描述模型的注册模型与FLAC3D之间的信息交换模型状态指示器的描述

返回目录自定义本构模型的基本方法头文件(usermodel.h)中进行新的本构模型派生类的声明修改模型的ID(>100)、名称和版本修改派生类的私有成员：私有变量和成员函数C++文件(usermodel.cpp)中修改模型结构(UserModel::UserModel(boolbRegister):ConstitutiveModel)constchar**UserModel::Properties()函数模型的参数名称字符串：用于在模型中显示constchar**UserModel::States()函数计算过程中的状态指示器：单元屈服状态返回目录doubleUserModel::GetProperty()和voidUserModel::SetProperty()函数constchar*UserModel::Initialize()函数参数和状态指示器的初始化，并对派生类声明中定义的私有变量进行赋值constchar*UserModel::Run(unsignednDin,States*ps)函数由应变增量计算得到应力增量，从而获得新的应力constchar*UserModel::SaveRestore()函数对计算结果进行保存。程序的调试在VC++的工程设置中将FLAC3D软件中的EXE文件路径加入到程序的调试范围中，并将FLAC3D自带的DLL文件加入到附加动态链接库(AdditionalDLLs)中，然后在Initialize()或Run()函数中设置断点，进行调试；在程序文件中加入return()语句，这样可以将希望得到的变量值以错误提示的形式在FLAC3D窗口中得到。常用模型信息传递指针变量返回目录模型源程序分解返回目录静力本构(Mohr-Coulomb)1.屈服函数剪切屈服拉屈服2.塑性势函数非关联关联①②①－剪切屈服②－拉伸屈服MC本构写成线性函数S：剪切屈服修正写成线性函数S：拉伸屈服修正constchar*UserMohrModel::Initialize(unsigneduDim,State*){if((uDim!=2)&&(uDim!=3))return("IllegaldimensioninUserMohrconstitutivemodel");dE1=dBulk+d4d3*dShear;dE2=dBulk-d2d3*dShear;dG2=2.0*dShear;doubledRsin=sin(dFriction*dDegRad);dNPH=(1.0+dRsin)/(1.0-dRsin);dCSN=2.0*dCohesion*sqrt(dNPH);if(dFriction){doubledApex=dCohesion*cos(dFriction*dDegRad)/dRsin;dTension=dTension<dApex?dTension:dApex;}dRsin=sin(dDilation*dDegRad);dRnps=(1.0+dRsin)/(1.0-dRsin);doubledRa=dE1-dRnps*dE2;doubledRb=dE2-dRnps*dE1;doubledRd=dRa-dRb*dNPH;dSC1=dRa/dRd;dSC3=dRb/dRd;dSC2=dE2*(1.0-dRnps)/dRd;dBISC=sqrt(1.0+dNPH*dNPH)+dNPH;dE21=dE2/dE1;return(0);}初始化各式中的常数constchar*UserMohrModel::Run(unsigneduDim,State*ps){if((uDim!=3)&&(uDim!=2))return("IllegaldimensioninMohrconstitutivemodel");if(ps->dHystDampMult>0.0)HDampInit(ps->dHystDampMult);/*---plasticityindicator:*//*store'now'info.as'past'andturn'now'infooff---*/if(ps->mState&mShearNow)ps->mState=(unsignedlong)(ps->mState|mShearPast);ps->mState=(unsignedlong)(ps->mState&~mShearNow);if(ps->mState&mTensionNow)ps->mState=(unsignedlong)(ps->mState|mTensionPast);ps->mState=(unsignedlong)(ps->mState&~mTensionNow);intiPlas=0;doubledTeTens=dTension;/*---trialelasticstresses---*/doubledE11=ps->stnE.d11;doubledE22=ps->stnE.d22;doubledE33=ps->stnE.d33;ps->stnS.d11+=dE11*dE1+(dE22+dE33)*dE2;ps->stnS.d22+=(dE11+dE33)*dE2+dE22*dE1;ps->stnS.d33+=(dE11+dE22)*dE2+dE33*dE1;ps->stnS.d12+=ps->stnE.d12*dG2;ps->stnS.d13+=ps->stnE.d13*dG2;ps->stnS.d23+=ps->stnE.d23*dG2;/*---calculateandsortps->incips->lstressesandps->incips->ldirections---*/AxesaDir;doubledPrinMin,dPrinMid,dPrinMax,sdif=0.0,psdif=0.0;inticase=0;boolbFast=ps->stnS.Resoltopris(&dPrinMin,&dPrinMid,&dPrinMax,&aDir,uDim,&icase,&sdif,&psdif);doubledPrinMinCopy=dPrinMin;doubledPrinMidCopy=dPrinMid;doubledPrinMaxCopy=dPrinMax;/*---Mohr-Coulombfailurecriterion---*/doubledFsurf=dPrinMin-dNPH*dPrinMax+dCSN;/*---Tensilefailurecriteria---*/doubledTsurf=dTension-dPrinMax;doubledPdiv=-dTsurf+(dPrinMin-dNPH*dTension+dCSN)*dBISC;/*---testsforfailure*/if(dFsurf<0.0&&dPdiv<0.0){iPlas=1;/*---shearfailure:correctiontops->incips->lstresses---*/ps->mState=(unsignedlong)(ps->mState|0x01);dPrinMin-=dFsurf*dSC1;dPrinMid-=dFsurf*dSC2;dPrinMax-=dFsurf*dSC3;}elseif(dTsurf<0.0&&dPdiv>0.0){iPlas=2;/*---tensionfailure:correctiontops->incips->lstresses---*/ps->mState=(unsignedlong)(ps->mState|0x02);doubledTco=dE21*dTsurf;dPrinMin+=dTco;dPrinMid+=dTco;dPrinMax=dTension;}if(iPlas){ps->stnS.Resoltoglob(dPrinMin,dPrinMid,dPrinMax,aDir,dPrinMinCopy,dPrinMidCopy,dPrinMaxCopy,uDim,icase,sdif,psdif,bFast);ps->bViscous=false;//Inhibitstiffness-dampingterms}else{ps->bViscous=true;//Allowstiffness-dampingterms}return(0);}蠕变本构(Viscous)Modelvisc（H-N）一.力-位移定律总位移弹性部分粘性部分uHuNFFH-N体示意图将虎克体和牛顿体的本构关系带入做变换可得Visc本构力-位移定律：Maxwell体－M体ModelVisc（H-N）二.中心差分方案（将应力张量做球应力张量和偏应力张量分解-增量关系）其中:偏应力张量和偏应变张量的本构关式中：球应力张量和球应变张量的本构关系总应力：Modelvisc（H-N）求解系数：dD=如果蠕变指标为真，则返回真值；否则dD=0。三.源程序分解constchar*UserViscousModel::Run(unsigneduDim,State*ps){if((uDim!=3)&&(uDim!=2))return("IllegaldimensioninUserViscousModel");doubledD=dGD2V*(ps->bCreep?ps->dTimeStep:0.0);if(dD>0.5)return("TimesteptoolargeforUserViscousModel");doubledVic1=1.0-dD;doubledVic2=1.0/(1.0+dD);/*volumetricanddeviatoricstrainincrements*/doubledDevol=ps->stnE.d11+ps->stnE.d22+ps->stnE.d33;doubledDevol3=d1d3*dDevol;doubledE11d=ps->stnE.d11-dDevol3;doubledE22d=ps->stnE.d22-dDevol3;doubledE33d=ps->stnE.d33-dDevol3;constchar*UserViscousModel::Initialize(unsigned,State*){dG2=2.0*dShear;if(dViscosity<=0.0)dGD2V=0.0;elsedGD2V=0.5*dShear/dViscosity;return(0);}求解系数：dGD2V=；1.如果粘滞系数<=0，dGD2V=0；2.否则dGD2V为真值。求解系数：C1，C2。球应变增量张量偏应变增量张量Modelvisc（H-N）constchar*UserViscousModel::Run(unsigneduDim,State*ps){/*oldisotropicanddeviatoricstresses*/doubledS0=d1d3*(ps->stnS.d11+ps->stnS.d22+ps->stnS.d33);doubledS11d=ps->stnS.d11-dS0;doubledS22d=ps->stnS.d22-dS0