- 积分
- 0
- 注册时间
- 2007-5-13
- 仿真币
-
- 最后登录
- 1970-1-1
|
const char *UserCamClayModel::Run(unsigned uDim,State *ps) {
if ((uDim!=2)&&(uDim!=3)) return("Illegal dimension in UserCamClayModel");
static double dPav,dQav,dEvav,dEvdPav;
static int iPind;
/* --- plasticity indicator:
---*/
/*
store 'now' info. as 'past' and turn 'now' info off ---*/
if (ps->mState & mShearNow) ps->mState |= mShearPast;
ps->mState &= ~mShearNow;
/* --- initialize stacks --- */
if (!ps->bySubZone) {
iPind = 0;
dPav
= 0.0;
dQav
= 0.0;
dEvav
= 0.0;
dEvdPav = 0.0;
}
/* --- trial elastic stresses试验弹性应力 --- */
double dA1 = dBulk + dC4D3 * dShear;
//体积模量K + 4.0 / 3.0*剪切模量G， dA1 =K+
double dA2 = dBulk - dC2D3 * dShear;
// dA2 =K-
double dG2 = 2.0 * dShear;
// dG2 =2
/*--由主应变分量增量得主应力增量 --*/
ps->stnS.d11 += ps->stnE.d11 * dA1 + (ps->stnE.d22 + ps->stnE.d33) * dA2
ps->stnS.d22 += ps->stnE.d22 * dA1 + (ps->stnE.d11 + ps->stnE.d33) * dA2;
ps->stnS.d33 += ps->stnE.d33 * dA1 + (ps->stnE.d11 + ps->stnE.d22) * dA2;
/* 因为在 中 */
ps->stnS.d12 += ps->stnE.d12 * dG2;
ps->stnS.d13 += ps->stnE.d13 * dG2;
ps->stnS.d23 += ps->stnE.d23 * dG2;
/* --- mean pressure平均压力dPVal
=
--- */
double dPVal = - (ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33) * dC1D3;
/* --- deviatoric stresses偏应力 --- */
double dDS11 = ps->stnS.d11 + dPVal;
double dDS22 = ps->stnS.d22 + dPVal;
double dDS33 = ps->stnS.d33 + dPVal;
double dDS12 = ps->stnS.d12;//
double dDS13 = ps->stnS.d13;
double dDS23 = ps->stnS.d23;
/*second deviatoric stress invariant第二偏应力不变量 - */
Double dJ2 =0.5*(dDS11*dDS11+dDS22*dDS22+dDS33*dDS33) + ps->stnS.d12* ps->stnS.d12 + ps->stnS.d13 * ps->stnS.d13 + ps->stnS.d23 * ps->stnS.d23;
/*----广义剪应力 ----*/
double dQVal = sqrt(3.*dJ2);
/* --- Cam-clay yield criterion剑桥屈服准则--- */
double dM2 = dMM*dMM;
//dMM是摩擦常数M
dM2 =M2
double dFS = dQVal*dQVal + dM2 * dPVal * (dPVal - dMPC); /* ----- 先期固结压力在这里是硬化参数Ha.即Ha=p c , f = ------*/
double dPNew = 0.0;
double dQNew = 0.0;
double dEVVal = 0.0;
double dEvdPVal = 0.0;
/* --- detect small pc --- */
if (dMPC < dMP1*1.e-3) {
/*先期固结压力<基准压力 *10-3 的情况下：*/
iPind = 1;
ps->mState |= mShearNow;
ps->stnS.d11 = 0.0;
// Sx=0
ps->stnS.d22 = 0.0;
ps->stnS.d33 = 0.0;
ps->stnS.d12 = 0.0;
ps->stnS.d13 = 0.0;
ps->stnS.d23 = 0.0;
dPNew
= dPVal;
dQNew
= dQVal;
dEVVal
= ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33; //体应变dEVVal=
dEvdPVal = dEVVal;
} else {
if (dFS > 0.0) {
//如果屈服函数大于零，表示剪切破坏
/* --- shear failure剪切破坏 --- */
iPind = 1;
ps->mState |= mShearNow;
double dSa
= 6.0 * dShear * dQVal;
//dSa=
double dSc
= dM2 * (2.0 * dPVal - dMPC);
//dSc=
double dSb
= dBulk * dSc;
//dSb =弹性体积模量
K*ca
double dBa
= dSa*dSa + dM2 * dSb*dSb;
//dBa=
double dBb = dSa * dQVal + 0.5 * dSb * dSc; //dBb =6qG*q+0.5* K* M2*(2p- pc) *M2*(2p- pc)=
double dBc
= dFS;
// dBc
=
double dAlam
= 0.0;
double dAlam1 = 0.0;
if (dBa != 0.0) {
double dBoa
= dBb / dBa;
//dBoa=
double dVal =dBoa*dBoa - dBc / dBa;
//dVal= 因为a ,b ,c 满足式 ，
if (dVal < 0.0) return("Cam-clay:
Yield envelope cannot be reached");//未达到屈服轨迹线
double dVal1
= sqrt(dVal);
//dVal1=
dAlam
= dBoa + dVal1;
// dAlam= +
dAlam1 = dBoa - dVal1;
//dAlam1= -
} else {
if (dBb != 0.0) {
dAlam
= 0.5 * dBc / dBb;
// dAlam=0.5* /
dAlam1 = dAlam;
} else return("Cam-clay:
Yield envelope cannot be reached");
}
if (fabs(dAlam1) < fabs(dAlam)) dAlam = dAlam1;
if (dAlam < 0.0) dAlam = 0.0;
dQNew
= dQVal - dSa * dAlam;
//
dPNew
= dPVal - dSb * dAlam;
//
/* --- update stresses 更新应力--- */
if (dQVal == 0.0) {
//如果 为0
double dAdd = dPNew - dPVal;
//
ps->stnS.d11 -= dAdd;
ps->stnS.d22 -= dAdd;
ps->stnS.d33 -= dAdd;
} else { //如果 不为0,新的应力分量： ,其中，偏应力分量：
ps->stnS.d11 = (dDS11 / dQVal) * dQNew - dPNew;
//
ps->stnS.d22 = (dDS22 / dQVal) * dQNew - dPNew;
ps->stnS.d33 = (dDS33 / dQVal) * dQNew - dPNew;
ps->stnS.d12 = (dDS12 / dQVal) * dQNew;
//
ps->stnS.d23 = (dDS23 / dQVal) * dQNew;
ps->stnS.d13 = (dDS13 / dQVal) * dQNew;
}
dEVVal
= ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33; //体应变dEVVal=
dEvdPVal = dAlam * dSc;
} else {
/* --- no failure --- */
dPNew
= dPVal;
dQNew
= dQVal;
dEVVal
= ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33; //体应变dEVVal=
dEvdPVal = 0.0;
}
}
/* update zone stresses and strains更新单元的应力应变--- */
double dVol = ps->dSubZoneVolume;
dPav
+= dPNew * dVol;
//在前面有/dPav,dQav,dEvav,dEvdPav;且设置初始值为0
dQav
+= dQNew * dVol;
dEvav
+= dEVVal * dVol;
dEvdPav += dEvdPVal * dVol;
/* --- the last zone has been processed, update parameters: 更新变量--- */
if (ps->bySubZone==ps->byTotSubZones-1) {
dVol = 1.0 / ps->dZoneVolume;
if (ps->byOverlay==2) dVol *= 0.5;
dPav
= dPav * dVol;
dQav
= dQav * dVol;
dEvav
= dEvav * dVol;
dEvdPav = dEvdPav * dVol;
dMP
= dPav;
//当前的平均有效应力
dMQ
= dQav;
//当前的平均偏应力
dEV
+= dEvav;
//累计总体积应变+=
dEV_P += dEvdPav;
//累计塑性体积应变
/* --- specific dVolume --- */
double dVold
= dMV;
//初始比容v0
dMV
= dVold * (1.0 + dEvav);
if (dMV < 0.0) return("Cam-clay:
Specific volume is negative");
//初始比容v0
/* ---
bulk modulus --- */
dBulk
= dMV * dPav / dKappa; //弹性体积模量K =初始比容v0*弹性膨胀线的斜率
if (dBulk > dBulkB) return("Cam-clay: Bulk modulus bound must be increased for stability");//弹性体积模量K
/* --- cap pressure硬化—软化准则：屈服曲线的形状取决于固结压力pc的值，此压力是塑性体积变化的函数，且随比容而变化。对于新的 , , ,可以通过 平面内固结线和膨胀线的交点得出，其表达式如下：--- */
if (iPind == 1) {
if (dPav > dMPC * 0.4) {
//先期固结压力
double dVal = dMV + dKappa * log(dPav / dMP1); /*初始比容v0+
弹性膨胀线的斜率
/基准压力 , 式中
dMPC = dMP1 * exp((dMV_L - dVal)/(dLambda-dKappa)); /*---先期固结压力=基准压力 ╳e （正常固结线的比容 - 弹性膨胀线的斜率 - ---*/
} else {
dMPC=dMPC*(1.0+dEvdPav*dVold/(dLambda-dKappa));//先期固结压力=弹性膨胀线斜率
}
if (dMPC <= 0.0) { //先期固结压力小于等于
dMPC
= 0.0;
//那设置先期固结压力0，因为土不能抗拉
dBulk = 0.0;
//弹性体积模量K
}
}
/* --- shear modulus --- */
if (dPoisson != 0.0) {
//泊松比不等于0
dShear=1.5*dBulk*(1.0-2.0*dPoisson)/(1.0+dPoisson); //剪切模量, 弹性体积模量,泊松比
/* ,弹性体积模量 随比容和平均应力的函数变化的*/
} else {
double dMaxg = 1.5 * dBulk;
//弹性体积模量K
double dMing = 0.0;
dShear = dShear < dMaxg ? dShear : dMaxg; //剪切模量
dShear = dShear > dMing ? dShear : dMing;
}
}
return(0);
} |
本帖子中包含更多资源
您需要 登录 才可以下载或查看，没有账号？注册
×
|
IP卡
狗仔卡
发表于 2009-6-22 19:40:55
提升卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡