flac3d中邓肯张自定义本构模型

// Modified at 2006-01-05 // 1. Tension correction should be modified // 2. Shear correction sould be modified // // Modified at 2006-08-15 // 3. To fit the non-initial stress calculation // 4. To 2D FLAC calculation #include "duncan.h" #include "conmodel.h" #include <math.h> #include <stdio.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; static const double dPa = 101325.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; static const unsigned long mKurNow = 0x16; // add 2007-01-17 static const unsigned long mKurPast = 0x32; static const unsigned long mEtNow = 0x64; static const unsigned long mEtPast = 0x128; // One static instance is neccessary as a part of internal registration process of the model with FLAC/FLAC3D static DuncanChang duncanchang(true); DuncanChang::DuncanChang(bool bRegister) :ConstitutiveModel(mnDuncanChang,bRegister), dCohesion(0.0),dFriction(0.0), dFailRatio(0.0), dFricDel(0.0), dKe(0.0),dNe(0.0),dKb(0.0),dMb(0.0), dKur(0.0), dF(0.0),dFcos(0.0),dFsin(0.0),dShear(0.0),dBulk(0.0), dSDHistroy(0.0), dS3Histroy(0.0), dSLHistroy(0.0) { }// const char **DuncanChang::Properties(void) const { static const char *strKey[] = { //double dCohesion,dFriction,dFailRatio,dK,dNe,dKB,dMb,dKur; "cohesion","friction","fricdel","ratiofail","ke","ne", "kb","mb","kur",0 }; return(strKey); } const char **DuncanChang::States(void) const { static const char *strKey[] = { "shear-n","tension-n","shear-p","tension-p", "Kur-n","Kur-p", "Et-n","Et-p",// add new state 2007-01-17 0 }; return(strKey); } /* * Note: Maintain order of property input/output */ double DuncanChang::GetProperty(unsigned ul) const { switch (ul) { ////double dCohesion,dFriction,dFailRatio,dKe,dNe,dKb,dM,dKur; case 1: return(dCohesion); case 2: return(dFriction); case 3: return(dFricDel); case 4: return(dFailRatio); case 5: return(dKe); case 6: return(dNe); case 7: return(dKb); case 8: return(dMb); case 9: return(dKur); } return(0.0); } void DuncanChang::SetProperty(unsigned ul,const double &dVal) { switch (ul) { ////double dCohesion,dFriction,dFailRatio,dKe,dN,dKB,dMb,dKur; case 1: dCohesion = dVal; break; case 2: dFriction = dVal; break; case 3: dFricDel = dVal; break; case 4: dFailRatio = dVal; break; case 5: dKe = dVal; break; case 6: dNe = dVal; break; case 7: dKb = dVal; break; case 8: dMb = dVal; break; case 9: dKur = dVal; break; } } const char *DuncanChang::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); DuncanChang *mm = (DuncanChang *)cm; dCohesion = mm->dCohesion; dFriction = mm->dFriction; dFricDel = mm->dFricDel; dFailRatio = mm->dFailRatio; dKe = mm->dKe; dNe = mm->dNe; dKb = mm->dKb; dMb = mm->dMb; dKur = mm->dKur; return(0); } const char *DuncanChang::Initialize(unsigned uDim,State *ps) { if ((uDim!=2)&&(uDim!=3)) return("Illegal dimension in Duncan-Chang constitutive model"); //(1)求主应力分量 /* --- calculate and sort ps->incips->l stresses and ps->incips->l directions --- */ //Caution:dPrinMin,dPrinMid,dPrinMax 分别是通常意义的三个主应力σ1，σ2，σ3，但均是负值 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 dSD = dPrinMax - dPrinMin;//主应力差>0 if (dSD > dSDHistroy) dSDHistroy = dSD; if (-1.0*dPrinMax > dS3Histroy) dS3Histroy = -1.0*dPrinMax; //考虑摩擦角随围压的影响 double dlog = (dS3Histroy != 0.0)? log10(dS3Histroy / dPa) : 0.0; dF = dFriction - dFricDel * dlog; dFcos = cos(dF * dDegRad); dFsin = sin(dF * dDegRad); /*/下面语句为调试所用 char temp[200]; sprintf(temp,"%f",dF); return(temp); //调试结束*/ //(2)计算应力水平 double dSF = 2* dCohesion * dFcos + (-2.0) * dPrinMax * dFsin;//破坏摩尔圆直径 double dSL = (dSF != 0)? dSD * (1-dFsin) / dSF : 0.0; //除数(no p) if (dSL > dSLHistroy) dSLHistroy = dSL; if (dSL >= 0.99) dSL = 0.99; //(3)加卸载判断，计算模量 double dSPa = dS3Histroy / dPa; double dEi = dKe *dPa * pow(dSPa,dNe); double dEt; if (dSD < dSDHistroy && dSL < dSLHistroy) { dEt = dKur *dPa * pow(dSPa,dNe); } else dEt = dEi * (1 - dFailRatio * dSL) * (1 - dFailRatio * dSL); //2006-06-11 add Emin double dEtmin = 0.25 * dKe * dPa * pow(0.02, dNe); if (dEt <= dEtmin) dEt = dEtmin; //根据0<v<0.49确定体积模量的范围0.33~17Et dBulk = dKb * dPa * pow(dSPa,dMb); dBulk = (dBulk < 0.33*dEt)? 0.33 * dEt : dBulk; dBulk = (dBulk > 17.0*dEt)? 17.0 * dEt : dBulk; dShear = 3 * dBulk * dEt / (9 * dBulk - dEt); // Divide! return(0); } const char *DuncanChang::Run(unsigned uDim,State *ps) { if ((uDim!=3)&&(uDim!=2)) return("Illegal dimension in Duncan-Chang constitutive model"); //拉裂和剪切破坏指示 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); if (ps->mState & mKurNow) ps->mState = (unsigned long)(ps->mState | mKurPast); ps->mState = (unsigned long)(ps->mState & ~mKurNow); if (ps->mState & mEtNow) ps->mState = (unsigned long)(ps->mState | mEtPast); ps->mState = (unsigned long)(ps->mState & ~mEtNow); 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 dSD = dPrinMax - dPrinMin;//主应力差>0 if (dSD > dSDHistroy) dSDHistroy = dSD; if (-1.0*dPrinMax > dS3Histroy) dS3Histroy = -1.0*dPrinMax; //考虑摩擦角随围压的影响 double dlog = (dS3Histroy != 0.0)? log10(dS3Histroy / dPa) : 0.0; dF = dFriction - dFricDel * dlog; dFcos = cos(dF * dDegRad); dFsin = sin(dF * dDegRad); //(2)计算应力水平 double dSF = 2* dCohesion * dFcos + (-2.0) * dPrinMax * dFsin;//破坏摩尔圆直径 double dSL = (dSF != 0)? dSD * (1-dFsin) / dSF : 0.0; //除数(no p) if (dSL > dSLHistroy) dSLHistroy = dSL; if (dSL >= 0.99) dSL = 0.99; //(3)加卸载判断，计算模量 double dSPa = dS3Histroy / dPa; double dEi = dKe *dPa * pow(dSPa,dNe); double dEt; if (dSD < dSDHistroy && dSL < dSLHistroy) { dEt = dKur *dPa * pow(dSPa,dNe); ps->mState = (unsigned long)(ps->mState | 0x16); } else { dEt = dEi * (1 - dFailRatio * dSL) * (1 - dFailRatio * dSL); ps->mState = (unsigned long)(ps->mState | 0x64); /*/下面语句为调试所用 char temp[200]; sprintf(temp,"%f",dEt); return(temp); //调试结束*/ } //2006-06-11 add Emin double dEtmin = 0.25 * dKe * dPa * pow(0.02, dNe); if (dEt <= dEtmin) dEt = dEtmin; //根据0<v<0.49确定体积模量的范围0.33~17Et dBulk = dKb * dPa * pow(dSPa,dMb); dBulk = (dBulk < 0.33*dEt)? 0.33 * dEt : dBulk; dBulk = (dBulk > 17.0*dEt)? 17.0 * dEt : dBulk; dShear = 3 * dBulk * dEt / (9 * dBulk - d