摩尔-库伦模型本构关系学习
mohr-coulomb本构模型编译流程图//在FLAC3D单元计算的每一个循环、每一个子单元中都要调用，其主要作用是根据应变增量计算应力张量
const char *UserMohrModel::Run(unsigned uDim,State *ps) {
if ((uDim!=3)&&(uDim!=2)) return("Illegal dimension in Mohr constitutive model");
/* --- plasticity indicator: --- */ //塑性指示
/*---store 'now' info. as 'past' and turn 'now' info off ---*//*--- ,(i=1,2,3)----*/
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.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; //
/*calculate and sort ps->incips->l stresses and ps->incips->l directions*///按照胡克定律计算单元的6个应力张量分量，然后按照塑性理论对这6个分量进行塑性修正
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; //胡克定律计算得到的3个主应力分量
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) { // <0 或者h<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) { //* <0或者h>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) {
//按照完成塑性修正得到的3个主应力进行应力状态的恢复
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);
} 1# whutzll 谢谢楼主分享~~~~ 顶一个 顶一下 本帖最后由 whutzll 于 2009-5-26 09:19 编辑
4# syx2006
谢谢鼓励，忘记说明，流程图是看过RUN部分自己做的，有待大伙们的批评指正，多谢，希望对大家有启发！：） 5# whutzll
塑性本构关系本质上是增量关系 ｆｆｆｆｆｆ 好啊，正在找。做摩尔库伦的一些思考 谢谢楼主难度很大啊 佩服
页:
[1]