usersoft.cpp

硬化软化模型UDM自定义开发源代码

#include "usersoft.h"
#include "contable.h"
#include <math.h>

//variables used by all model objects. Hence only one copy is maintained for all objects
static const double d1d3    = 1.0 / 3.0;
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 mShearNow    = 0x01;  /* state logic */
static const unsigned long mTensionNow  = 0x02;
static const unsigned long mShearPast   = 0x04;
static const unsigned long mTensionPast = 0x08;

// One Static instance isneccessary as a part of internal registration process of the model with FLAC/FLAC3D
static UserStrainSofteningModel mUserStrainSofteningModel(true);

UserStrainSofteningModel::UserStrainSofteningModel(bool bRegister)
             :ConstitutiveModel(mnUserStrainSofteningModel,bRegister),
              dBulk(0.0),dShear(0.0),dYoung(0.0),dPoisson(0.0),
              dCohesion(0.0),dFriction(0.0),dDilation(0.0),
              dTension(0.0),dSHP(0.0),dTHP(0.0),iCTab(0),iFTab(0),
              iDTab(0),iTTab(0),dE1(0.0),dE2(0.0),dG2(0.0),dNPH(0.0),
              dCSN(0.0),dSC1(0.0),dSC2(0.0),dSC3(0.0),dSF1(0.0),
              dSF3(0.0),dBisc(0.0),uiCTab(0),uiFTab(0),uiDTab(0),
              uiTTab(0) { }


const char *UserStrainSofteningModel::Keyword(void) const { return("userssoft"); }


const char *UserStrainSofteningModel::Name(void) const { return("User-StrainSoftening"); }


const char **UserStrainSofteningModel::Properties(void) const {
  static const char *strKey[] = {
    "bulk", "shear", "young", "poisson", "cohesion", "friction",
    "dilation", "tension", "ctable", "ftable", "dtable", "ttable", "es_plastic",
    "et_plastic", 0
  };
  return(strKey);
}


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

/*  * Note: Maintain order of property input/output
*/
double UserStrainSofteningModel::GetProperty(unsigned ul) const {
  switch (ul) {
    case 1:  return(dBulk);
    case 2:  return(dShear);
    case 3:  return(dYoung);
    case 4:  return(dPoisson);
    case 5:  return(dCohesion);
    case 6:  return(dFriction);
    case 7:  return(dDilation);
    case 8:  return(dTension);
    case 9:  return((double)iCTab);
    case 10: return((double)iFTab);
    case 11: return((double)iDTab);
    case 12: return((double)iTTab);
    case 13: return(dSHP);
    case 14: return(dTHP);
  }
  return(0.0);
}

void UserStrainSofteningModel::SetProperty(unsigned ul,const double &dVal) {
  switch (ul) {
    case 1: { //  BULK
      dBulk = dVal;
      YoungPoissonFromBulkShear(&dYoung,&dPoisson,dBulk,dShear);
      break;
    }
    case 2: { // SHEAR
      dShear = dVal;
      YoungPoissonFromBulkShear(&dYoung,&dPoisson,dBulk,dShear);
      break;
    }
    case 3: { // YOUNG
      dYoung = dVal;
      BulkShearFromYoungPoisson(&dBulk,&dShear,dYoung,dPoisson);
      break;
    }
    case 4: { // POISSON
      if ((dVal==0.5)||(dVal==-1.0)) return;
      dPoisson = dVal;
      BulkShearFromYoungPoisson(&dBulk,&dShear,dYoung,dPoisson);
      break;
    }
    case 5:  dCohesion = dVal;  break;
    case 6:  dFriction = dVal;  break;
    case 7:  dDilation = dVal;  break;
    case 8:  dTension  = dVal;  break;
    case 9:  iCTab     = ConTableList::Round(dVal);  break;
    case 10: iFTab     = ConTableList::Round(dVal);  break;
    case 11: iDTab     = ConTableList::Round(dVal);  break;
    case 12: iTTab     = ConTableList::Round(dVal);  break;
    case 13: dSHP      = dVal;  break;
    case 14: dTHP      = dVal;  break;
  }
}


const char *UserStrainSofteningModel::Copy(const ConstitutiveModel *m) {

  //Detects type mismatch error and returns error string. otherwise returns 0
  const char *str = ConstitutiveModel::Copy(m);
  if (str) return(str);

  UserStrainSofteningModel *em = (UserStrainSofteningModel *)m;
  dBulk     = em->dBulk;
  dShear    = em->dShear;
  dYoung    = em->dYoung;
  dPoisson  = em->dPoisson;
  dCohesion = em->dCohesion;
  dFriction = em->dFriction;
  dDilation = em->dDilation;
  dTension  = em->dTension;
  dSHP      = em->dSHP;
  dTHP      = em->dTHP;
  return(0);
}


const char *UserStrainSofteningModel::Initialize(unsigned,State *) {
  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);
  double 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;
  dSF1         = 1.0 / dRd;
  dSF3         = -dRnps / dRd;
  dBisc        = sqrt(1.0 + dNPH * dNPH) + dNPH;
  if (iCTab||iFTab||iDTab||iTTab) {
    const char *str = 0;
    ConTableList *ctl = ConTableList::Instance();
    if (!ctl) return("No table implementation available");
    if (iCTab) str = ctl->IndexFromID(iCTab,&uiCTab);
    if (str) return(str);
    if (iFTab) str = ctl->IndexFromID(iFTab,&uiFTab);
    if (str) return(str);
    if (iDTab) str = ctl->IndexFromID(iDTab,&uiDTab);
    if (str) return(str);
    if (iTTab) str = ctl->IndexFromID(iTTab,&uiTTab);
    if (str) return(str);
  }
  return(0);
}


const char *UserStrainSofteningModel::Run(unsigned ulDim,State *ps) {
  if ((ulDim!=2)&&(ulDim!=3)) return("Illegal dimension in StrainSoftening model");
  static double dDqs,dDqt;
  /* --- plasticity indicator:                                  */
  /*     store 'now' info. as 'past' and turn 'now' info off ---*/
  if (ps->mState & mShearNow) ps->mState |= mShearPast;
  ps->mState &= ~mShearNow;
  if (ps->mState & mTensionNow) ps->mState |= mTensionPast;
  ps->mState &= ~mTensionNow;
  int iPlas = 0;
  /* --- hardening/softening:
       initialize stacks to calculate hardening parameters for zone --- */
  if (!ps->bySubZone) {
    dDqs = 0.0;
    dDqt = 0.0;
  }
  /* --- trial elastic stresses --- */
  double stnEd11 = ps->stnE.d11;
  double stnEd22 = ps->stnE.d22;
  double stnEd33 = ps->stnE.d33;

  ps->stnS.d11 += stnEd11 * dE1 + (stnEd22 + stnEd33) * dE2;
  ps->stnS.d22 += (stnEd11 + stnEd33) * dE2 + stnEd22 * dE1;
  ps->stnS.d33 += (stnEd11 + stnEd22) * dE2 + stnEd33 * 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 principal stresses and principal directions --- */
  Axes axes;
  double dPrinMin,dPrinMid,dPrinMax,sdif=0.0,psdif=0.0;
  int icase=0;

  bool bFast=ps->stnS.Resoltopris(&dPrinMin,&dPrinMid,&dPrinMax,&axes,ulDim, &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 principal stresses ---*/
    ps->mState |= mShearNow;
    dPrinMin   -= dFsurf * dSC1;
    dPrinMid   -= dFsurf * dSC2;
    dPrinMax   -= dFsurf * dSC3;
  } else if (dTsurf < 0.0 && dPdiv > 0.0) {
    iPlas = 2;
    /* --- tension failure: correction to principal stresses ---*/
    ps->mState |= mTensionNow;
    double dTco = dE2 * dTsurf / dE1;
    dPrinMin   += dTco;
    dPrinMid   += dTco;
    dPrinMax    = dTension;
  }
  if (iPlas != 0) {
	ps->stnS.Resoltoglob(dPrinMin,dPrinMid, dPrinMax, axes, dPrinMinCopy,dPrinMidCopy,dPrinMaxCopy, ulDim, icase, sdif,  psdif, bFast);
    /* --- hadRdening padRameter accumulation --- */
    if (iPlas==1) {
      /* ---                      shear padRameter --- */
      double dDe1p = dFsurf * dSF1;
      double dDe3p = dFsurf * dSF3;
      double dDepa = d1d3 * (dDe1p + dDe3p);
      dDe1p       -= dDepa;
      dDe3p       -= dDepa;
      dDqs        += sqrt(0.5 * (dDe1p*dDe1p+dDepa*dDepa+dDe3p*dDe3p)) * ps->dSubZoneVolume;
    }
    if (iPlas == 2) {
      /* ---                      tensile padRameter --- */
      double dAux = dTsurf / dE1;
      if (dAux < 0.) dAux = -dAux;
      dDqt  += dAux * ps->dSubZoneVolume;
    }
  }
  /* --- plastic padRameter incrementation and properties update --- */
  if (ps->bySubZone==ps->byTotSubZones-1) {
    double dAux = 1.0 / ps->dZoneVolume;
    if (ps->byOverlay==2) dAux *= 0.5;
    dSHP += dDqs * dAux;
    dTHP += dDqt * dAux;
    if (dDqs > 0.0) {
      const char *str=0;
      if (uiCTab) { // Update Cohesion
        str = ConTableList::Instance()->YFromX(uiCTab,dSHP,&dCohesion);
        if (str) return(str);
      }
      if (uiFTab) { // update friction
        str = ConTableList::Instance()->YFromX(uiFTab,dSHP,&dFriction);
        if (str) return(str);
      }
      if (uiDTab) { // update dilation
        str = ConTableList::Instance()->YFromX(uiDTab,dSHP,&dDilation);
        if (str) return(str);
      }
    }
    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);
    double 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;
    dSF1         = 1.0 / dRd;
    dSF3         = -dRnps / dRd;
    dBisc        = sqrt(1.0 + dNPH * dNPH) + dNPH;
    if (dDqt > 0.0) {
      if (uiTTab) { // update tensile strength
        const char *str = ConTableList::Instance()->YFromX(uiTTab,dTHP,&dAux);
        dTension        = dTension < dAux ? dTension : dAux;
      }
    }
  }
  return(0);
}

/* * Save all properties for the model
   * Note: It is not necessary to save and restore member variables that would be
     initialized. This reduces the size of save file.
*/
const char *UserStrainSofteningModel::SaveRestore(ModelSaveObject *mso) {
  // Checks for type mismatch and returns error string. Otherwise 0.
  const char *str = ConstitutiveModel::SaveRestore(mso);
  if (str) return(str);

  // 10 represents 10 properties that are doubles
  // and 4 represents 4 properties that are integers
  mso->Initialize(10,4);
  mso->Save(0,dBulk);
  mso->Save(1,dShear);
  mso->Save(2,dYoung);
  mso->Save(3,dPoisson);
  mso->Save(4,dCohesion);
  mso->Save(5,dFriction);
  mso->Save(6,dDilation);
  mso->Save(7,dTension);
  mso->Save(8,dSHP);
  mso->Save(9,dTHP);
  //Save integers seperately
  mso->Save(0,iCTab);
  mso->Save(1,iFTab);
  mso->Save(2,iDTab);
  mso->Save(3,iTTab);

  return(0);
}

// EOF