ITASCA(武汉)咨询有限公司 报告人: 朱永生 2007.04.28 依泰斯卡(武汉)

内 容 必要性 试验总结的本构模型 特定条件下的本构模型 交叉学科的本构模型 二次开发环境 自定义本构模型的功能 自定义本构模型的基本方法 内 容 必要性 试验总结的本构模型 特定条件下的本构模型 交叉学科的本构模型 二次开发环境 自定义本构模型的功能 自定义本构模型的基本方法 常用模型信息传递指针变量 模型源程序分解

二次开发环境 FLAC3D采用面向对象的语言标准C++编写 本构模型都是以动态连接库文件(.DLL文件)的形式提供 VC++6.0(SP4)或更高版本的开发环境 优点 自定义的本构模型和软件自带的本构模型的执行效率处在同一个水平 自定义本构模型(.DLL文件)适用于高版本的FLAC(2D)、3DEC、UDEC等其他Itasca软件中 返回目录

自定义本构模型的功能 主要功能：对给出的应变增量得到新的应力 辅助功能： 模型名称、版本 读写操作 模型文件的编写 基类(class Constitutive Model)的描述 成员函数的描述 模型的注册 模型与FLAC3D之间的信息交换 模型状态指示器的描述 返回目录

自定义本构模型的基本方法 头文件(usermodel.h)中进行新的本构模型派生类的声明 修改模型的ID(>100)、名称和版本 修改派生类的私有成员：私有变量和成员函数 C++文件(usermodel.cpp)中修改模型结构 (UserModel::UserModel(bool bRegister): ConstitutiveModel) const char **UserModel::Properties()函数 模型的参数名称字符串：用于在模型中显示 const char **UserModel::States()函数 计算过程中的状态指示器：单元屈服状态 返回目录

double UserModel::GetProperty()和void UserModel:: SetProperty()函数 const char * UserModel::Initialize()函数 参数和状态指示器的初始化，并对派生类声明中定义的私有变量进行赋值 const char * UserModel::Run(unsigned nDin,States *ps) 函数 由应变增量计算得到应力增量，从而获得新的应力 const char * UserModel::SaveRestore()函数 对计算结果进行保存。 程序的调试 在VC++的工程设置中将FLAC3D软件中的EXE文件路径加入到程序的调试范围中，并将FLAC3D自带的DLL文件加入到附加动态链接库(Additional DLLs)中，然后在Initialize()或Run()函数中设置断点，进行调试； 在程序文件中加入return()语句，这样可以将希望得到的变量值以错误提示的形式在FLAC3D窗口中得到。

常用模型信息传递指针变量 返回目录

模型源程序分解 返回目录

静力本构 (Mohr-Coulomb)

MC本构 ① ② ①－剪切屈服 ②－拉伸屈服 1. 屈服函数 2. 塑性势函数 非关联 剪切屈服 拉屈服 关联

剪切屈服修正 写成线性函数S：

拉伸屈服修正 写成线性函数S：

const char *UserMohrModel::Initialize(unsigned uDim,State *) { if ((uDim!=2)&&(uDim!=3)) return("Illegal dimension in UserMohr constitutive model"); dE1 = dBulk + d4d3 * dShear; dE2 = dBulk - d2d3 * dShear; dG2 = 2.0 * dShear; double dRsin = sin(dFriction * dDegRad); dNPH = (1.0 + dRsin) / (1.0 - dRsin); dCSN = 2.0 * dCohesion * sqrt(dNPH); if (dFriction) { double dApex = dCohesion * cos(dFriction * dDegRad) / dRsin; dTension = dTension < dApex ? dTension : dApex; } dRsin = sin(dDilation * dDegRad); dRnps = (1.0 + dRsin) / (1.0 - dRsin); double dRa = dE1 - dRnps * dE2; double dRb = dE2 - dRnps * dE1; double dRd = 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); 初始化各式中的常数

const char *UserMohrModel::Run(unsigned uDim,State *ps){ if ((uDim!=3)&&(uDim!=2)) return("Illegal dimension in Mohr constitutive model"); if(ps->dHystDampMult > 0.0) HDampInit(ps->dHystDampMult); /* --- plasticity indicator: */ /* store 'now' info. as 'past' and turn 'now' info off ---*/ if (ps->mState & mShearNow) ps->mState = (unsigned long)(ps->mState | mShearPast); ps->mState = (unsigned long)(ps->mState & ~mShearNow); if (ps->mState & mTensionNow) ps->mState = (unsigned long)(ps->mState | mTensionPast); ps->mState = (unsigned long)(ps->mState & ~mTensionNow); int iPlas = 0; double dTeTens = dTension; /* --- trial elastic stresses --- */ double dE11 = ps->stnE.d11; double dE22 = ps->stnE.d22; double dE33 = ps->stnE.d33;

ps->stnS.d11 += dE11 * dE1 + (dE22 + dE33) * dE2; ps->stnS.d12 += ps->stnE.d12 * dG2; ps->stnS.d13 += ps->stnE.d13 * dG2; ps->stnS.d23 += ps->stnE.d23 * dG2; /* --- calculate and sort ps->incips->l stresses and ps->incips->l directions --- */ Axes aDir; double dPrinMin,dPrinMid,dPrinMax,sdif=0.0,psdif=0.0; int icase=0; bool bFast=ps->stnS.Resoltopris(&dPrinMin,&dPrinMid,&dPrinMax,&aDir,uDim, &icase, &sdif, &psdif); double dPrinMinCopy = dPrinMin; double dPrinMidCopy = dPrinMid; double dPrinMaxCopy = dPrinMax; /* --- Mohr-Coulomb failure criterion --- */ double dFsurf = dPrinMin - dNPH * dPrinMax + dCSN; /* --- Tensile failure criteria --- */ double dTsurf = dTension - dPrinMax; double dPdiv = -dTsurf + (dPrinMin - dNPH * dTension + dCSN) * dBISC;

/* --- tests for failure */ if (dFsurf < 0.0 && dPdiv < 0.0) { iPlas = 1; /* --- shear failure: correction to ps->incips->l stresses ---*/ ps->mState = (unsigned long)(ps->mState | 0x01); dPrinMin -= dFsurf * dSC1; dPrinMid -= dFsurf * dSC2; dPrinMax -= dFsurf * dSC3; } else if (dTsurf < 0.0 && dPdiv > 0.0) { iPlas = 2; /* --- tension failure: correction to ps->incips->l stresses ---*/ ps->mState = (unsigned long)(ps->mState | 0x02); double dTco = 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; // Inhibit stiffness-damping terms } else { ps->bViscous = true; // Allow stiffness-damping terms return(0);

蠕变本构 (Viscous)

Model visc（H-N） Maxwell体－M体 F 一. 力-位移定律 uH uN 总位移 弹性部分 粘性部分

Model Visc（H-N） 总应力： 二. 中心差分方案（将应力张量做球应力张量和偏应力张量分解-增量关系） 其中: 偏应力张量 和偏应变张量 的本构关 式中： 球应力张量 和球应变张量 的本构关系 总应力：

Model visc（H-N） 三. 源程序分解 求解系数：dD= 如果蠕变指标为真，则返回真值； 否则dD=0。 三. 源程序分解 const char *UserViscousModel::Run(unsigned uDim,State *ps) { if ((uDim!=3)&&(uDim!=2)) return("Illegal dimension in UserViscousModel"); double dD = dGD2V * (ps->bCreep ? ps->dTimeStep : 0.0); if (dD > 0.5) return("Timestep too large for UserViscousModel"); double dVic1 = 1.0 - dD; double dVic2 = 1.0 / (1.0 + dD); /* volumetric and deviatoric strain increments */ double dDevol = ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33; double dDevol3 = d1d3 * dDevol; double dE11d = ps->stnE.d11 - dDevol3; double dE22d = ps->stnE.d22 - dDevol3; double dE33d = ps->stnE.d33 - dDevol3; const char *UserViscousModel::Initialize(unsigned,State *) { dG2 = 2.0 * dShear; if (dViscosity <= 0.0) dGD2V = 0.0; else dGD2V = 0.5 * dShear / dViscosity; return(0); } 求解系数：dGD2V= ； 1. 如果粘滞系数 <=0， dGD2V=0； 2. 否则dGD2V为真值。 求解系数：C1，C2。 球应变增量张量 偏应变增量张量

Model visc（H-N） const char *UserViscousModel::Run(unsigned uDim,State *ps) { /* old isotropic and deviatoric stresses */ double dS0 = d1d3 * (ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33); double dS11d = ps->stnS.d11 - dS0; double dS22d = ps->stnS.d22 - dS0; double dS33d = ps->stnS.d33 - dS0; /* new deviatoric stresses */ dS11d = (dS11d * dVic1 + dG2 * dE11d) * dVic2; dS22d = (dS22d * dVic1 + dG2 * dE22d) * dVic2; dS33d = (dS33d * dVic1 + dG2 * dE33d) * dVic2; ps->stnS.d12 = (ps->stnS.d12 * dVic1 + dG2 * ps->stnE.d12 ) * dVic2; ps->stnS.d23 = (ps->stnS.d23 * dVic1 + dG2 * ps->stnE.d23 ) * dVic2; ps->stnS.d13 = (ps->stnS.d13 * dVic1 + dG2 * ps->stnE.d13 ) * dVic2; /* new elastic isotropic stress */ dS0 += dBulk * dDevol; /* total stresses */ ps->stnS.d11 = dS11d + dS0; ps->stnS.d22 = dS22d + dS0; ps->stnS.d33 = dS33d + dS0; return(0); } 更新i步总应力 i-1步偏应力张量 i-1步球应力张量 求解i步本构方程(偏应力总量)： 求解i步本构方程(球应力总量)：