example_src2.rar_udm flac_摩尔库伦_摩尔库伦c++_摩尔库伦模型

#include "usermohr.h" #include <math.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 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 is neccessary as a part of internal registration process of the model with FLAC/FLAC3D static UserMohrModel usermohrmodel(true); UserMohrModel::UserMohrModel(bool bRegister) :ConstitutiveModel(mnUserMohrModel,bRegister), dBulk(0.0), dShear(0.0), dCohesion(0.0), dFriction(0.0), dDilation(0.0), dTension(0.0), dYoung(0.0), dPoisson(0.0), dE1(0.0), dE2(0.0), dG2(0.0), dNPH(0.0), dCSN(0.0), dSC1(0.0), dSC3(0.0), dBISC(0.0), dE21(0.0) { } const char **UserMohrModel::Properties(void) const { static const char *strKey[] = { "bulk", "shear","cohesion","friction","dilation", "tension","young","poisson" , 0 }; return(strKey); } const char **UserMohrModel::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 UserMohrModel::GetProperty(unsigned ul) const { switch (ul) { case 1: return(dBulk); case 2: return(dShear); case 3: return(dCohesion); case 4: return(dFriction); case 5: return(dDilation); case 6: return(dTension); case 7: return(dYoung); case 8: return(dPoisson); } return(0.0); } void UserMohrModel::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: dCohesion = dVal; break; case 4: dFriction = dVal; break; case 5: dDilation = dVal; break; case 6: dTension = dVal; break; case 7: { dYoung = dVal; BulkShearFromYoungPoisson(&dBulk,&dShear,dYoung,dPoisson); break; } case 8: { if ((dVal==0.5)||(dVal==-1.0)) return; dPoisson = dVal; BulkShearFromYoungPoisson(&dBulk,&dShear,dYoung,dPoisson); break; } } } const char *UserMohrModel::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); UserMohrModel *mm = (UserMohrModel *)cm; dBulk = mm->dBulk; dShear = mm->dShear; dCohesion = mm->dCohesion; dFriction = mm->dFriction; dDilation = mm->dDilation; dTension = mm->dTension; dYoung = mm->dYoung; dPoisson = mm->dPoisson; return(0); } const char *UserMohrModel::Initialize(unsigned uDim,State *) { if ((uDim!=2)&&(uDim!=3)) return("Illegal dimension in UserMohr constitutive model"); 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); 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; dBISC = sqrt(1.0 + dNPH * dNPH) + dNPH; dE21 = dE2 / dE1; return(0); } const char *UserMohrModel::Run(unsigned uDim,State *ps){ if ((uDim!=3)&&(uDim!=2)) return("Illegal dimension in Mohr constitutive model"); if(ps->dHystDampMult > 0.0) HDampInit(ps->dHystDampMult); /* --- plasticity indicator: */ /* store 'now' info. as 'past' and turn 'now' info off ---*/ 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 --- */ 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; 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 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) { 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) { 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); } /* * 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 *UserMohrModel::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(8,0); mso->Save(0,dBulk); mso->Save(1,dShear); mso->Save(2,dCohesion); mso->Save(3,dFriction); mso->Save(4,dDilation); mso->Save(5,dTension); mso->Save(6,dYoung); mso->Save(7,dPoisson); return(0); } void UserMohrModel::HDampInit(const double dHMult) { double dShearNew = dShear * dHMult; dE1 = dBulk + d4d3 * dShearNew; dE2 = dBulk - d2d3 * dShearNew; dG2 = 2.0 * dShearNew; 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; dE21 = dE2 / dE1; } // EOF