修正的剑桥模型在flac中加卸载准则是什么

wufan123wang

#include "usercam.h"
#include <math.h>

static const unsigned long mShearNow  = 0x01;
static const unsigned long mShearPast = 0x02;
static const double dC1D3 = 1.0 / 3.0;
static const double dC2D3 = 2.0 / 3.0;
static const double dC4D3 = 4.0 / 3.0;
static UserCamClayModel usercamclaymodel(true);

UserCamClayModel::UserCamClayModel(bool bRegister)
             :ConstitutiveModel(mnUserCamClayModel,bRegister),
              dBulk(0.0),dBulkB(0.0),dShear(0.0),dPoisson(0.0),
              dKappa(0.0),dLambda(0.0),dMM(0.0),dMPC(0.0),dMP1(0.0),
              dMV_L(0.0),dMV(0.0),dEV(0.0),dEV_P(0.0),dMP(0.0),dMQ(0.0) {
}

const char *UserCamClayModel::Keyword(void) const { return("usercam-clay"); }

const char *UserCamClayModel::Name(void) const { return("User-Cam-Clay"); }

const char **UserCamClayModel:roperties(void) const {
  static const char *strKey[] = {
    "bulk", "shear", "bulk_bound", "poisson", "kappa", "lambda",
    "mm", "mpc", "mp1", "mv_l", "cv", "cam_ev", "camev_p", "cam_cp", "cq", 0
  };
  return(strKey);
}

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

double UserCamClayModel::GetProperty(unsigned ul) const {
  switch (ul) {
    case 1:  return(dBulk);
    case 2:  return(dShear);
    case 3:  return(dBulkB);
    case 4:  return(dPoisson);
    case 5:  return(dKappa);
    case 6:  return(dLambda);
    case 7:  return(dMM);
    case 8:  return(dMPC);
    case 9:  return(dMP1);
    case 10: return(dMV_L);
    case 11: return(dMV);
    case 12: return(dEV);
    case 13: return(dEV_P);
    case 14: return(dMP);
    case 15: return(dMQ);
  }
  return(0.0);
}

void UserCamClayModel::SetProperty(unsigned ul,const double &dVal) {
  switch (ul) {
    case 1:  dBulk = dVal;  break;
    case 2:  dShear = dVal;  break;
    case 3:  dBulkB = dVal;  break;
    case 4:  dPoisson = dVal;  break;
    case 5:  dKappa = dVal;  break;
    case 6:  dLambda = dVal;  break;
    case 7:  dMM = dVal;  break;
    case 8:  dMPC = dVal;  break;
    case 9:  dMP1 = dVal;  break;
    case 10: dMV_L = dVal;  break;
    case 11: dMV = dVal;  break;
    case 12: dEV = dVal;  break;
    case 13: dEV_P = dVal;  break;
    case 14: dMP = dVal;  break;
    case 15: dMQ = dVal;  break;
  }
}

const char *UserCamClayModel::Copy(const ConstitutiveModel *m) {
  const char *str = ConstitutiveModel::Copy(m);
  if (str) return(str);
  UserCamClayModel *em = (UserCamClayModel *)m;
  dBulk    = em->dBulk;
  dShear   = em->dShear;
  dBulkB   = em->dBulkB;
  dPoisson = em->dPoisson;
  dKappa   = em->dKappa;
  dLambda  = em->dLambda;
  dMM      = em->dMM;
  dMPC     = em->dMPC;
  dMP1     = em->dMP1;
  dMV_L    = em->dMV_L;
  dMV      = em->dMV;
  dEV      = em->dEV;
  dEV_P    = em->dEV_P;
  dMP      = em->dMP;
  dMQ      = em->dMQ;
  return(0);
}

const char *UserCamClayModel::Initialize(unsigned uDim,State *) {
  if ((uDim!=2)&&(uDim!=3)) return("Illegal dimension in UserCamClayModel");
  // Initialize mean pressure
  if (dMP1 == 0.0) dMP1 = 1.0;
  if (dMP1 < 0.0 || dMPC < 0.0) return("roperty bad");
  if (dMV_L < 0.0 || dKappa == 0.0) return("Property bad");
  if (dLambda == dKappa) return("Property bad");
  if (dBulkB == 0.0) return("Property bad");
  /* --- initialize specific volume --- */
  if (dMV == 0.0) {
    double dP0 = dMP;           // initialized in camclay_ini_p
    if (dP0 <= 0.0) return("Cam-clay: Mean effective pressure is negative");
    else {
      dMV = dMV_L-dLambda*log(dMPC/dMP1)+dKappa*log(dMPC/dP0);
      if (dMV < 0.0) return("Cam-clay:  Specific volume is negative");
    }
  }
  /* --- initialize current bulk modulus --- */
  if (dBulk == 0.0) {
    double dP0 = dMP;
    if (dP0 <= 0.0) return("Cam-clay: Mean effective pressure is negative");
    else {
      dBulk = dMV * dP0 / dKappa;
      if (dPoisson)
        dShear = 1.5*dBulk*(1.0-2.0*dPoisson)/(1.0+dPoisson);
    }
  }
  /* --- initialize shear modulus --- */
  if (!dShear && dPoisson)
    dShear = 1.5*dBulk*(1.0-2.0*dPoisson)/(1.0+dPoisson);
  return(0);
}

static const int Pav    = 0;
static const int Qav    = 1;
static const int Evav   = 2;
static const int EvdPav = 3;
static const int Pind = 0;
// The above are now indices into ps->working[] and ps->iworkingp[].
// This was necessary for thread safety.
const char *UserCamClayModel::Run(unsigned uDim,State *ps) {
  if ((uDim!=2)&&(uDim!=3)) return("Illegal dimension in UserCamClayModel");
  /* --- plasticity indicator:                                  */
  /*     store 'now' info. as 'past' and turn 'now' info off ---*/
  if (ps->mState & mShearNow) ps->mState |= mShearPast;
  ps->mState &= ~mShearNow;
  /* --- initialize stacks --- */
  if (!ps->bySubZone) {
    ps->iworking[Pind] = 0;
    ps->working[Pav]   = 0.0;
    ps->working[Qav]   = 0.0;
    ps->working[Evav]  = 0.0;
    ps->working[EvdPav] = 0.0;
  }
  /* --- trial elastic stresses --- */
  double dA1 = dBulk + dC4D3 * dShear;
  double dA2 = dBulk - dC2D3 * dShear;
  double dG2 = 2.0 * dShear;
  ps->stnS.d11 += ps->stnE.d11 * dA1 + (ps->stnE.d22 + ps->stnE.d33) * dA2;
  ps->stnS.d22 += ps->stnE.d22 * dA1 + (ps->stnE.d11 + ps->stnE.d33) * dA2;
  ps->stnS.d33 += ps->stnE.d33 * dA1 + (ps->stnE.d11 + ps->stnE.d22) * dA2;
  ps->stnS.d12 += ps->stnE.d12 * dG2;
  ps->stnS.d13 += ps->stnE.d13 * dG2;
  ps->stnS.d23 += ps->stnE.d23 * dG2;
  /* --- mean pressure --- */
  double dPVal = - (ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33) * dC1D3;
  /* --- deviatoric stresses --- */
  double dDS11 = ps->stnS.d11 + dPVal;
  double dDS22 = ps->stnS.d22 + dPVal;
  double dDS33 = ps->stnS.d33 + dPVal;
  double dDS12 = ps->stnS.d12;
  double dDS13 = ps->stnS.d13;
  double dDS23 = ps->stnS.d23;
  /* --- second deviatoric stress invariant --- */
  double dJ2   = 0.5*(dDS11*dDS11+dDS22*dDS22+dDS33*dDS33)+ps->stnS.d12*ps->stnS.d12+ps->stnS.d13*ps->stnS.d13+ps->stnS.d23*ps->stnS.d23;
  double dQVal = sqrt(3.*dJ2);
  /* --- Cam-clay yield criterion --- */
  double dM2 = dMM*dMM;
  double dFS = dQVal*dQVal + dM2 * dPVal * (dPVal - dMPC);
  double dPNew = 0.0;
  double dQNew = 0.0;
  double dEVVal = 0.0;
  double dEvdPVal = 0.0;
  /* --- detect small pc --- */
  if (dMPC < dMP1*1.e-3) {
    ps->iworking[Pind] = 1;
    ps->mState |= mShearNow;
    ps->stnS.d11 = 0.0;
    ps->stnS.d22 = 0.0;
    ps->stnS.d33 = 0.0;
    ps->stnS.d12 = 0.0;
    ps->stnS.d13 = 0.0;
    ps->stnS.d23 = 0.0;
    dPNew   = dPVal;
    dQNew   = dQVal;
    dEVVal  = ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33;
    dEvdPVal = dEVVal;
  } else {
    if (dFS > 0.0) {
      /* ---                        shear failure --- */
      ps->iworking[Pind] = 1;
      ps->mState |= mShearNow;
      double dSa    = 6.0 * dShear * dQVal;
      double dSc    = dM2 * (2.0 * dPVal - dMPC);
      double dSb    = dBulk * dSc;
      double dBa    = dSa*dSa + dM2 * dSb*dSb;
      double dBb    = dSa * dQVal + 0.5 * dSb * dSc;
      double dBc    = dFS;
      double dAlam  = 0.0;
      double dAlam1 = 0.0;
      if (dBa != 0.0) {
        double dBoa   = dBb / dBa;
        double dVal   = dBoa*dBoa - dBc / dBa;
        if (dVal < 0.0) return("Cam-clay:  Yield envelope cannot be reached");
        double dVal1   = sqrt(dVal);
        dAlam  = dBoa + dVal1;
        dAlam1 = dBoa - dVal1;
      } else {
        if (dBb != 0.0) {
           dAlam  = 0.5 * dBc / dBb;
           dAlam1 = dAlam;
        } else return("Cam-clay:  Yield envelope cannot be reached");
      }
      if (fabs(dAlam1) < fabs(dAlam)) dAlam = dAlam1;
      if (dAlam < 0.0) dAlam = 0.0;
      dQNew  = dQVal - dSa * dAlam;
      dPNew  = dPVal - dSb * dAlam;
      /* ---                        update stresses --- */
      if (dQVal == 0.0) {
         double dAdd = dPNew - dPVal;
         ps->stnS.d11 -= dAdd;
         ps->stnS.d22 -= dAdd;
         ps->stnS.d33 -= dAdd;
      } else {
         ps->stnS.d11 = (dDS11 / dQVal) * dQNew - dPNew;
         ps->stnS.d22 = (dDS22 / dQVal) * dQNew - dPNew;
         ps->stnS.d33 = (dDS33 / dQVal) * dQNew - dPNew;
         ps->stnS.d12 = (dDS12 / dQVal) * dQNew;
         ps->stnS.d23 = (dDS23 / dQVal) * dQNew;
         ps->stnS.d13 = (dDS13 / dQVal) * dQNew;
      }
      dEVVal  = ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33;
      dEvdPVal = dAlam * dSc;
    } else {
      /* ---                        no failure --- */
      dPNew   = dPVal;
      dQNew   = dQVal;
      dEVVal  = ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33;
      dEvdPVal = 0.0;
    }
  }
  /* update zone stresses and strains --- */
  double dVol = ps->dSubZoneVolume;
  ps->working[Pav]   += dPNew * dVol;
  ps->working[Qav]   += dQNew * dVol;
  ps->working[Evav]  += dEVVal * dVol;
  ps->working[EvdPav] += dEvdPVal * dVol;
  /* --- the last zone has been processed, update parameters: --- */
  if (ps->bySubZone==ps->byTotSubZones-1) {
    dVol = 1.0 / ps->dZoneVolume;
    if (ps->byOverlay==2) dVol *= 0.5;
    ps->working[Pav]   = ps->working[Pav] * dVol;
    ps->working[Qav]   = ps->working[Qav] * dVol;
    ps->working[Evav]  = ps->working[Evav] * dVol;
    ps->working[EvdPav] = ps->working[EvdPav] * dVol;
    dMP    = ps->working[Pav];
    dMQ    = ps->working[Qav];
    dEV   += ps->working[Evav];
    dEV_P += ps->working[EvdPav];
    /* ---                  specific dVolume --- */
    double dVold  = dMV;
    dMV    = dVold * (1.0 + ps->working[Evav]);
    if (dMV < 0.0) return("Cam-clay:  Specific volume is negative");
    /* ---                  bulk modulus --- */
    dBulk  = dMV * ps->working[Pav] / dKappa;
    if (dBulk > dBulkB) return("Cam-clay:  Bulk modulus bound must be increased for stability");
    /* ---                  cap pressure --- */
    if (ps->iworking[Pind] == 1) {
      if (ps->working[Pav] > dMPC * 0.4) {
        double dVal = dMV + dKappa * log(ps->working[Pav] / dMP1);
        dMPC = dMP1 * exp((dMV_L - dVal)/(dLambda-dKappa));
      } else {
        dMPC = dMPC * (1.0+ps->working[EvdPav]*dVold/(dLambda-dKappa));
      }
      if (dMPC <= 0.0) {
        dMPC  = 0.0;
        dBulk = 0.0;
      }
    }
    /* ---                  shear modulus --- */
    if (dPoisson != 0.0) {
      dShear = 1.5*dBulk*(1.0-2.0*dPoisson)/(1.0+dPoisson);
    } else {
      double dMaxg = 1.5 * dBulk;
      double dMing = 0.0;
      dShear = dShear < dMaxg ? dShear : dMaxg;
      dShear = dShear > dMing ? dShear : dMing;
    }
  }
  return(0);
}

const char *UserCamClayModel::SaveRestore(ModelSaveObject *mso) {
  const char *str = ConstitutiveModel::SaveRestore(mso);
  if (str) return(str);
  mso->Initialize(15,0);
  mso->Save( 0,dBulk);
  mso->Save( 1,dShear);
  mso->Save( 2,dBulkB);
  mso->Save( 3,dPoisson);
  mso->Save( 4,dKappa);
  mso->Save( 5,dLambda);
  mso->Save( 6,dMM);
  mso->Save( 7,dMPC);
  mso->Save( 8,dMP1);
  mso->Save( 9,dMV_L);
  mso->Save(10,dMV);
  mso->Save(11,dEV);
  mso->Save(12,dEV_P);
  mso->Save(13,dMP);
  mso->Save(14,dMQ);
  return(0);
}
// EOF
请哪位高手能帮在下看一下

wufan123wang

有牛人知道吗?