摩尔-库伦模型本构关系学习

whutzll

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);
}

whutzll

1# whutzll

lookcity

谢谢楼主分享~~~~

syx2006

顶一个    顶一下

whutzll

4# syx2006

谢谢鼓励，忘记说明，流程图是看过RUN部分自己做的，有待大伙们的批评指正，多谢，希望对大家有启发！：）

whutzll

5# whutzll
塑性本构关系本质上是增量关系

wuzhide

ｆｆｆｆｆｆ

zhangmingfei

好啊，正在找。做摩尔库伦的一些思考

feigle

谢谢楼主  难度很大啊   佩服