#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
