FLAC3D中的自定义本构模型问题

张胤

请教大家，在FLAC3D自定义本构模型中，我用一个很简单的算例验证自己的模型，但在调试时，发现d11,d22,d33,d12,d23,d13,开始时全是零，导致应变也全是零，算例中我是加了三轴应力的，很奇怪！希望高手能给予帮助！谢谢！

aizhongluo

能看看你的模型编译吗

张胤

aizhongluo 发表于 2011-9-26 10:44
能看看你的模型编译吗

#include "HTC.h"
#include "math.h"
#include "contable.h"

//variables used by all model objects. Hence only one copy is maintained for all objects

static const double d2d3 = 2.0 / 3.0;
static const double dPi  = 3.141592653589793238462643383279502884197169399;
static const double dDegRad = dPi / 180.0;
// Plasticity Indicators
static const unsigned long mPlastNow    = 0x01;  /* state logic */
//static const unsigned long mCompressNow  = 0x02;
static const unsigned long mPlastPast   = 0x04;
//static const unsigned long mCompressPast = 0x08;

// One static instance is neccessary as a part of internal registration process of the model with FLAC/FLAC3D
static HTC HTCmodel(true);


HTC::HTC(bool bRegister)
          :ConstitutiveModel(mnHTC,bRegister), dBulk(0.0),
           dShear(0.0), dFtfc(0.0), dFc(0.0), dYoung(0.0), dPoisson(0.0),
                   dFsin(0.0), dFcos(0.0),                                                                                                                                                                                                                    
                   dE1(0.0),dE2(0.0),dG2(0.0),
           dFsurf(0.0),Lode(0.0),I1(0.0)                                                                                         {
}


const char **HTC:roperties(void) const {
  static const char *strKey[] = {
    "bulk",   "shear","ftfc","fc",
    "young","poisson" , 0
  };
  return(strKey);
}


const char **HTC::States(void) const {
  static const char *strKey[] = {
    "plast-n","plast-p",0
  };
  return(strKey);
}

/*  * Note: Maintain order of property input/output
*/
double HTC::GetProperty(unsigned ul) const {
  switch (ul) {
    case 1:  return(dBulk);
    case 2:  return(dShear);
    case 3:  return(dFtfc);
    case 4:  return(dFc);
    case 5:  return(dYoung);
    case 6:  return(dPoisson);
  }
  return(0.0);
}
void HTC::SetProperty(unsigned ul,const double &dVal) {
  switch (ul) {
    case 1: {
      dBulk = dVal;
      YoungPoissonFromBulkShear(&dYoung,&dPoisson,dBulk,dShear);
      break;
    }
    case 2: {
      dShear = dVal;
      YoungPoissonFromBulkShear(&dYoung,&dPoisson,dBulk,dShear);
      break;
    }
    case 3: dFtfc = dVal;  break;
    case 4: dFc = dVal;    break;
    case 5: {
      dYoung = dVal;
      BulkShearFromYoungPoisson(&dBulk,&dShear,dYoung,dPoisson);
      break;
    }
    case 6: {
      if ((dVal==0.5)||(dVal==-1.0)) return;
      dPoisson = dVal;
      BulkShearFromYoungPoisson(&dBulk,&dShear,dYoung,dPoisson);
      break;
    }
  }
}


const char *HTC::Copy(const ConstitutiveModel *cm) {
  //Detects type mismatch error and returns error string. otherwise returns 0
  const char *str = ConstitutiveModel::Copy(cm);
  if (str) return(str);
  HTC *mm = (HTC *)cm;
  dBulk     = mm->dBulk;
  dShear    = mm->dShear;
  dFtfc     = mm->dFtfc;
  dFc       = mm->dFc;
  dYoung    = mm->dYoung;
  dPoisson  = mm->dPoisson;
  return(0);
}


const char *HTC::Initialize(unsigned uDim,State *) {
  if ((uDim!=2)&&(uDim!=3)) return("Illegal dimension in HTC constitutive model");

      dE1   = dBulk + d4d3 * dShear;
      dE2   = dBulk - d2d3 * dShear;
      dG2   = 2.0 * dShear;
//以下是参数选用（ftfc=1/2/3分别对应ft/fc=0.08、0.10、0.12）

if (dFtfc <1.5){
    dPopa = 2.8367;
    dPopb = 0.5659;
    dPopc = 11.8254;
    dPopd = 0.2723;
  }

else if ((dFtfc <2.5)&&(dFtfc >1.5)){
    dPopa = 2.0107;
    dPopb = 0.9713;
    dPopc = 9.1413;
    dPopd = 0.2310;
  }

else if (dFtfc > 2.5){
    dPopa = 1.4608;
    dPopb = 1.2411;
    dPopc = 7.3549;
    dPopd = 0.2035;
  }
  
  return(0);

}


const char *HTC::Run(unsigned uDim,State *ps){
  
if ((uDim!=3)&&(uDim!=2)) return("Illegal dimension in HTC constitutive model");

  if(ps->dHystDampMult > 0.0) HDampInit(ps->dHystDampMult);


  /* --- plasticity indicator:                                  */
  /*     store 'now' info. as 'past' and turn 'now' info off ---*/
  int iPlas = 0;
  if (ps->mState & mPlastNow)
      ps->mState = (unsigned long)(ps->mState | mPlastPast);
    ps->mState = (unsigned long)(ps->mState & ~mPlastNow);

  
/* --- trial elastic stresses --- */
  ps->stnS.d11 += ps->stnE.d11 * dE1 + (ps->stnE.d22 + ps->stnE.d33) * dE2;
  ps->stnS.d22 += (ps->stnE.d11 + ps->stnE.d33) * dE2 + ps->stnE.d22 * dE1;
  ps->stnS.d33 += (ps->stnE.d11 + ps->stnE.d22) * dE2 + ps->stnE.d33 * dE1;
  ps->stnS.d12 += ps->stnE.d12 * dG2;
  ps->stnS.d13 += ps->stnE.d13 * dG2;
  ps->stnS.d23 += ps->stnE.d23 * dG2;
  /* --- mean stress --- */
  double dSigi = (ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33) / 3.0;
  /* --- deviatoric stresses --- */
  STensor stnD;
  stnD.d11 = ps->stnS.d11 - dSigi;
  stnD.d22 = ps->stnS.d22 - dSigi;
  stnD.d33 = ps->stnS.d33 - dSigi;
  stnD.d12 = ps->stnS.d12;
  stnD.d13 = ps->stnS.d13;
  stnD.d23 = ps->stnS.d23;
  /* --- second deviatoric stress invariant --- */
  double dJ2 = 0.5 * (stnD.d11 * stnD.d11 +
                             stnD.d22 * stnD.d22 +
                             stnD.d33 * stnD.d33)
                      + ps->stnS.d12 * ps->stnS.d12
                      + ps->stnS.d13 * ps->stnS.d13
                      + ps->stnS.d23 * ps->stnS.d23;

  double  dJ3 = stnD.d11 * stnD.d22 * stnD.d33 + 2.0 * ps->stnS.d12 * ps->stnS.d23 * ps->stnS.d13
                      - stnD.d11*ps->stnS.d23 * ps->stnS.d23
                                          - stnD.d22*ps->stnS.d13 * ps->stnS.d13
                                          - stnD.d33*ps->stnS.d12 * ps->stnS.d12;
double  shang=(-3.0)*sqrt(3.0)*dJ3;
double  xia=2.0*dJ2*sqrt(dJ2);
double  bizhi=shang/xia;
  Lode = asin(bizhi);
   I1   = ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33;
  dFsin = sin(Lode/3.0+120.0*dDegRad);
  dFcos = cos(Lode/3.0+120.0*dDegRad);
   
     
//定义对J2的偏导
double   dJ21 = stnD.d11;
double   dJ22 = stnD.d22;
double   dJ23 = stnD.d33;
double   dJ24 = 2.0*stnD.d12;
double   dJ25 = 2.0*stnD.d23;
double   dJ26 = 2.0*stnD.d13;
//定义对J3的偏导
double   dJ31 = stnD.d22*stnD.d33-stnD.d23*stnD.d23+dJ2/3.0;
double   dJ32 = stnD.d33*stnD.d11-stnD.d13*stnD.d13+dJ2/3.0;
double   dJ33 = stnD.d22*stnD.d11-stnD.d12*stnD.d12+dJ2/3.0;
double   dJ34 = 2.0*(stnD.d12*stnD.d13-stnD.d11*stnD.d23);
double   dJ35 = 2.0*(stnD.d12*stnD.d23-stnD.d22*stnD.d13);
double   dJ36 = 2.0*(stnD.d23*stnD.d13-stnD.d33*stnD.d12);
//定义系数C
double   CC1 = dPopc/3.0+dPopd ;   
double   CC2 = dPopb+dPopc*2.0*dFsin/sqrt(3.0)-tan(Lode)*dFcos*dPopc*2.0/sqrt(3.0);
double   CC3 = dPopc*(-1.0)*dFcos/(dJ2*cos(Lode)) ;

// if (sata > 29.0){
//         if((Lode/dDegRad)> 0.0)
//         {
//         CC1 = dPopc/3.0+dPopd ;   
//         CC2 = dPopb+dPopc*2.0/sqrt(3.0);
//         CC3 = 0.0 ;
//  }
//     else {
//          CC1 = dPopc/3.0+dPopd ;   
//          CC2 = dPopb+dPopc/sqrt(3.0);
//         CC3 = 0.0 ;
//  }
// }
double dE11 = ps->stnE.d11;
double dE22 = ps->stnE.d22;
double dE33 = ps->stnE.d33;
double dE12 = ps->stnE.d12;
double dE23 = ps->stnE.d23;
double dE13 = ps->stnE.d13;


//定义对F的偏导
  
double  dF1 = CC1+CC2*dJ21/(2.0*sqrt(dJ2))+CC3*dJ31;
double  dF2 = CC1+CC2*dJ22/(2.0*sqrt(dJ2))+CC3*dJ32;
double  dF3 = CC1+CC2*dJ23/(2.0*sqrt(dJ2))+CC3*dJ33;
double  dF4 =     CC2*dJ24/(2.0*sqrt(dJ2))+CC3*dJ34;
double  dF5 =     CC2*dJ25/(2.0*sqrt(dJ2))+CC3*dJ35;
double  dF6 =     CC2*dJ26/(2.0*sqrt(dJ2))+CC3*dJ36;

//定义应力的塑性变量
  
double  P1  = dE1*dF1+dE2*(dF2+dF3);
double  P2  = dE1*dF2+dE2*(dF1+dF3);
double  P3  = dE1*dF3+dE2*(dF2+dF1);
double  P4  = dG2*dF4;
double  P5  = dG2*dF5;
double  P6  = dG2*dF6;
  
  double P11  = P1*P1*dE11+P1*P2*dE22+P1*P3*dE33+P1*P4*dE12+P1*P5*dE23+P1*P6*dE13;
  double P22  = P2*P1*dE11+P2*P2*dE22+P2*P3*dE33+P2*P4*dE12+P2*P5*dE23+P2*P6*dE13;
  double P33  = P3*P1*dE11+P3*P2*dE22+P3*P3*dE33+P3*P4*dE12+P3*P5*dE23+P3*P6*dE13;
  double P44  = P4*P1*dE11+P4*P2*dE22+P4*P3*dE33+P4*P4*dE12+P4*P5*dE23+P4*P6*dE13;
  double P55  = P5*P1*dE11+P5*P2*dE22+P5*P3*dE33+P5*P4*dE12+P5*P5*dE23+P5*P6*dE13;
  double P66  = P6*P1*dE11+P6*P2*dE22+P6*P3*dE33+P6*P4*dE12+P6*P5*dE23+P6*P6*dE13;

  double dPP   = P1*dF1+P2*dF2+P3*dF3+P4*dF4+P5*dF5+P6*dF6;

  double Dep1  = P11/dPP;
  double Dep2  = P22/dPP;
  double Dep3  = P33/dPP;
  double Dep4  = P44/dPP;
  double Dep5  = P55/dPP;
  double Dep6  = P66/dPP;

  /* --- HTC failure criterion --- */
  

  dFsurf = dPopa*dJ2/dFc+dPopb*sqrt(dJ2)+dPopc*((2.0/sqrt(3.0))*sqrt(dJ2)*dFsin+I1/3.0)+dPopd*I1-dFc;

  /* --- tests for failure */
  if ( dFsurf > 0.0) {
    iPlas = 1;

    /* --- Plast failure: correction to ps->incips->l stresses ---*/
        ps->mState = ps->mState | mPlastNow;
    ps->stnS.d11 -= Dep1;
    ps->stnS.d22 -= Dep2;
    ps->stnS.d33 -= Dep3;
        ps->stnS.d12 -= Dep4;
    ps->stnS.d23 -= Dep5;
    ps->stnS.d13 -= Dep6;

  }


   
                return(0);
}


const char *HTC::SaveRestore(ModelSaveObject *mso) {
// Checks for type mismatch and returns error string. Otherwise 0.
  const char *str = ConstitutiveModel::SaveRestore(mso);
  if (str) return(str);
  // 8 represents 8 properties that are doubles
  // and 0 represents 0 properties that are integers
  mso->Initialize(6,0);
  mso->Save(0,dBulk);
  mso->Save(1,dShear);
  mso->Save(2,dFtfc);
  mso->Save(3,dFc);
  mso->Save(4,dYoung);
  mso->Save(5,dPoisson);
  return(0);
}

// EOF

yangxin546

我也做本构模型开发的，如果你愿意，可以互相学习下，QQ301069165

kays

     dE1   = dBulk + d4d3 * dShear;
dE2   = dBulk - d2d3 * dShear;
      dG2   = 2.0 * dShear;

这个玩意在run 里面也做一次试试