- 积分
- 0
- 注册时间
- 2007-5-13
- 仿真币
-
- 最后登录
- 1970-1-1
|
霍克-布朗模型应用于岩石工程较为理想，我想将来在掌握其本构关系基础上有所突破
恳请大伙批评指正！
------------------------------------------------------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[e3psum] = 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[e3psum] += 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[e3psum] * dAux;
/* soften properties if plastic strain e2 is non-zero */
if (ps->working[e3psum] != 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);
} |
|
IP卡
狗仔卡
发表于 2009-6-9 10:39:21
提升卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
