G03：Itasca（FLAC、UDEC/3DEC、PFC）-霍克-布朗,流程图,模型,学习-霍克-布朗（Hoek-Brown）本构模型学习-附自创流程图，{Simwe仿真论坛(forum.simwe.com)，CAE/CAD/CAM/，FEA/FEM/有限元分析论坛-(手机验证注册)}

whutzll 发表于 2009-6-9 10:39:21

霍克-布朗（Hoek-Brown）本构模型学习-附自创流程图

霍克-布朗模型应用于岩石工程较为理想，我想将来在掌握其本构关系基础上有所突破
恳请大伙批评指正！
------------------------------------------------------RUN-----------------------------------------------------------------------
static const int e3psum = 0;
const char *UserHoekModel::Run(unsigned uDim,State *ps) {

if(ps->dHystDampMult > 0.0) HDampInit(ps->dHystDampMult);

/* --- plasticity indicator:                               */
/* store 'now' info. as 'past' and turn 'now' info off ---*/
double s1n= 0.0, s2n=0.0, s3n=0.0, F0=0.0, F1=0.0, F2=0.0, e0=0.0, e1=0.0, e2=0.0;
double de3p=0.0, repStress=0.0, smallStress=0.0, gammaAF=0.0, gammaCV=-1.0;
double gamma=0.0;
double mult=1.0;
int i=0;
if (ps->mState & mShearNow)
ps->mState = (unsigned long)(ps->mState | mShearPast);
ps->mState = (unsigned long)(ps->mState & ~mShearNow);
/* --- hardening/softening:
initialize stacks to calculate hardening paramters for zone --- */
if (!ps->bySubZone) {
ps->working = 0.0;
}
/* --------------------------------------------------*/
/* --- 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;
double _sMeanOld = ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33;

if (dSigci <= 0.0) { // (no strength)
double sig0 = (ps->stnS.d11+ps->stnS.d22+ps->stnS.d33) / 3.0;
ps->stnS.d11 = sig0;
ps->stnS.d22 = sig0;
ps->stnS.d33 = sig0;
ps->stnS.d12 = 0.0;
ps->stnS.d13 = 0.0;
ps->stnS.d23 = 0.0;
ps->bViscous = false;// inhibit viscous strains
double _sMeanNew = ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33;
   ps->dMeanPlasticStressChange = (_sMeanOld - _sMeanNew) / 3.0;
return(0);
}
/* --- 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;
/* --- Hoek-Brown failure criterion --- */
double s1 = -dPrinMin; // Use positive compression
double s2 = -dPrinMid;
double s3 = -dPrinMax;
double fac = dMb * s3 / dSigci + dS;
double fac2 = 0.0;
if (fac >= 0.0) fac2 = dSigci * pow(fac,dA);
else          fac2 = -dSigci * pow(fabs(fac),dA); // Ensure continuity
double dFsurf = s1 - s3 - fac2;
dF = dFsurf;
if (dFsurf > 0.0) { //以下程式全是在dFsurf > 0时成立
// yield ...屈服时
//double SigUCS = dSigci * pow(dS,dA);// Unconfined compressive strength单轴抗压强度
//double SigTen = -dS * dSigci / dMb; // Tensile stress at apex顶部抗拉应力
ps->mState = (unsigned long)(ps->mState | mShearNow);
ps->bViscous = false; // inhibit viscous strains 禁止粘性应变

repStress = __max(fabs(s1),fabs(s3));
smallStress = __max(repStress,dSigci) * Ftol; // A small stress取一小应变
if (dFsurf <= smallStress) return(0);
// Stress adjustment ...应力调整
de3p = __max(fabs(dE11),fabs(dE22)); // put starting value into de3p
de3p = __max(de3p,fabs(dE33));
de3p = __max(de3p,fabs(ps->stnE.d12));
de3p = __max(de3p,fabs(ps->stnE.d13));
de3p = __max(de3p,fabs(ps->stnE.d23));// 对de3p初始赋值，取应力最大值
bool converged = false;
double dConst1 = 0.0;
double dConst2 = 0.0;
double dConst3 = 0.0;

if(de3p != 0.0)
{ //计算流动系数
   gammaAF = -1.0/(1.0+dA*dSigci*pow(fabs(dS+dMb*s3/dSigci),dA-1.0)*dMb/dSigci);
   gammaCV = -1.0;
   if ((s3<0.0) && (s1<0.0)) {//当主应力均小于0时
   gamma = s1 / s3;    //径向流动法则
   } else if ((s3<=0.0) && (s1>0.0)) { //当最大主应力<=0且最小主应力>0时
   gamma = gammaAF;          //相关联流动法则
   } else if ((s3>0.0) && (s3 < dS3cv)) { //当最小主应力>0且小于sig3max时    gamma = 1.0/((1.0/gammaAF)+((1.0/gammaCV)-(1.0/gammaAF))*s3/dS3cv);// 复合流动法则
   } else {
   gamma = gammaCV; //恒定体积流动法则
   }
   dConst1 = gamma * dE1 + dE2;//
   dConst3 = gamma * dE2 + dE1;//
   dConst2 = dE2 * (1.0 + gamma);//
   s1n= s1 - de3p * dConst1; //最终应力
   s3n= s3 - de3p * dConst3;
   s2n= s2 - de3p * dConst2;
   fac = dMb * s3n / dSigci + dS;
   if (fac >= 0.0) fac2 = dSigci * pow(fac,dA);
   else          fac2 = -dSigci * pow(fabs(fac),dA);
   F1 = s1n - s3n - fac2; //重新计算F值
   F0 = dFsurf;
   e0 = 0.0;
   e1 = de3p;
   dIts = 0.0;

   for (i=0 ; i<15 ; i++) {//运算15步
   if (F1 == F0) break;
   e2= (F1 * e0 - F0 * e1) / (F1 - F0);//利用牛顿法插值
   s1n = s1 - e2 * dConst1;
   s3n = s3 - e2 * dConst3;
   s2n = s2 - e2 * dConst2;
   fac = dMb * s3n / dSigci + dS;
   if (fac >= 0.0) fac2 = dSigci * pow(fac,dA);
   else          fac2 = -dSigci * pow(fabs(fac),dA);
   F2 = s1n - s3n - fac2;
   dF = F2;
   dIts = dIts + 1.0;
   if (fabs(F2) < smallStress) //判断新F
   {
      converged = true;
      break;
   }
   F0 = F1;
   F1 = F2;
   e0 = e1;
   e1 = e2;
   }
}

// Try a bisection algorithm if solution did not converge.二等分算法
if(!converged)
{
   double e2in = (s1 - s3) / ((1.0-gamma) * (dE2 - dE1));// F <= 0
   double e2out = 0.0; // F > 0

   for(i=0;i<100;i++)
   { //运算100步
   double e2mid = (e2in + e2out) / 2.0;
   s1n = s1 - e2mid * dConst1;
   s3n = s3 - e2mid * dConst3;
   fac = dMb * s3n / dSigci + dS;
   if (fac >= 0.0) fac2 =dSigci * pow(fac,dA);
   else          fac2 = -dSigci * pow(fabs(fac),dA);
   F2 = s1n - s3n - fac2;
   if (fabs(F2) < smallStress)
   {
      dF = F2;
      s1n = s1 - e2mid * dConst1;
      s3n = s3 - e2mid * dConst3;
      s2n = s2 - e2mid * dConst2;
      converged = true;
      dIts = i;
      e2 = e2mid;
      break;
   }
   if(F2 < 0.0)
      e2in = e2mid;
   else if(F2 > 0.0)
      e2out = e2mid;
   }
}

if (!converged)
   return("Hoek-Brown model failed to converge.");

/* ---- sum e2 for all sub-zones 对所有子域e2倍递加-----------------------*/
ps->working += e2 * ps->dSubZoneVolume;
if (s3n > s1n) { // Crossed bisector ... set to apex交叉等分线，在顶点设置
   s1n = -dS * dSigci / dMb;
   s2n = s1n;
   s3n = s1n;
}
dPrinMin = -s1n;
dPrinMid = -s2n;
dPrinMax = -s3n;
//
ps->stnS.Resoltoglob(dPrinMin, dPrinMid, dPrinMax, aDir, dPrinMinCopy, dPrinMidCopy, dPrinMaxCopy, uDim, icase, sdif,psdif, bFast);
ps->bViscous = false; // inhibit viscous strains禁止粘性应变
//若存在硬化增量则更新
if (ps->bySubZone==ps->byTotSubZones-1) {
   double dAux = 1.0 / ps->dZoneVolume;
   if (ps->byOverlay==2) dAux *= 0.5;
   dDEP += ps->working * dAux;
   /* soften properties if plastic strain e2 is non-zero */
   if (ps->working != 0.0) {
   const char *str=0;
   if (uiMULTab) { // update mult更新因子
      str = ConTableList::Instance()->YFromX(uiMULTab,fabs(s3n),&mult);
      if (str) return(str);
   }
   if (uiMTab) { // Update m更新m值
      str = ConTableList::Instance()->YFromX(uiMTab,mult*fabs(dDEP),&dMb);
      if (str) return(str);
   }
   if (uiSTab) { // update s更新s值
      str = ConTableList::Instance()->YFromX(uiSTab,mult*fabs(dDEP),&dS);
      if (str) return(str);
   }
   if (uiATab) { // Update a更新a值
      str = ConTableList::Instance()->YFromX(uiATab,mult*fabs(dDEP),&dA);
      if (str) return(str);
   }
   if (uiCTab) { // update sigci更新sigci值
      str = ConTableList::Instance()->YFromX(uiCTab,mult*fabs(dDEP),&dSigci);
      if (str) return(str);
   }
   }
}
}

if (ps->mState & mShearNow)
{
ps->bViscous = false; // inhibit viscous strains禁止粘性应变
double _sMeanNew = ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33;
   ps->dMeanPlasticStressChange = (_sMeanOld - _sMeanNew) / 3.0;
}
else
{
ps->bViscous = true;// allow viscous strains允许粘性应变
}

return(0);
}

whutzll 发表于 2009-6-9 10:41:14

本帖最后由 whutzll 于 2009-6-19 09:56 编辑

1# whutzll

whutzll 发表于 2009-6-9 10:44:20

以下是流程图中公式

lookcity 发表于 2009-6-9 10:57:23

很有意义的讨论

lookcity 发表于 2009-6-9 11:01:07

sig3cv的取值要依靠实验对位移的分布范围影响比较明显

whutzll 发表于 2009-6-10 20:03:59

5# lookcity
谢谢lookcity的关注！
有关sigci和sig3cv可以参考rocscience.com网站上roclab软件计算。
软件是免费的，较为实用。

lidelin 发表于 2009-7-1 16:44:20

很好，很强大，谢谢楼主分享

lookcity 发表于 2009-7-2 10:31:52

5# lookcity
谢谢lookcity的关注！
有关sigci和sig3cv可以参考rocscience.com网站上roclab软件计算。
软件是免费的，较为实用。
whutzll 发表于 2009-6-10 20:03 http://forum.simwe.com/images/common/back.gif
Rocklab我用过但里面不包括sig3cv的计算。
Rocklab是基于Hoek的理论的，但只包括峰值之前的，FLAC3D增加了流动准则。

luhuan 发表于 2009-7-2 12:31:26

谢谢分享:victory:

wuzhide 发表于 2009-7-21 12:00:44

好东西！！！！！！

recally 发表于 2009-7-21 14:39:06

h-b本来就只是个经验性的强度准则，不包含什么应力应变关系等。

hbsyffzzjj 发表于 2010-11-22 10:31:02

下了，看看！

li_lianchong 发表于 2011-3-19 13:11:31

这样的帖多发点

flyeagle009 发表于 2011-3-20 09:10:43

支持一下，希望有所创新

李凯发表于 2011-8-22 23:54:24

支持，很有技术含量！

李凯发表于 2011-8-22 23:54:41

支持，很有技术含量！

qianxunzn 发表于 2011-8-24 19:27:38

楼主这个怎么这么复杂啊我还想着用霍克布朗准则验算围岩的破坏呢，，现在怕了

邓小一 发表于 2011-8-25 08:43:04

学习了，这个本构的计算效率如何？

李凯发表于 2011-8-25 09:27:57

回复 1# whutzll
请问楼主，在FLAC3D模拟时，怎样使用霍克—布朗本构模型进行数值分析呢？您发的这个程序能直接应用吗？非常感谢

李凯发表于 2011-8-25 15:47:16

回复 1# whutzll

楼主能否指点一下

页: [1] 2

Simwe仿真论坛(forum.simwe.com)，CAE/CAD/CAM/，FEA/FEM/有限元分析论坛---(手机验证注册)'s Archiver

霍克-布朗（Hoek-Brown）本构模型学习-附自创流程图