FLAC3D中的自定义本构模型问题
请教大家，在FLAC3D自定义本构模型中，我用一个很简单的算例验证自己的模型，但在调试时，发现d11,d22,d33,d12,d23,d13,开始时全是零，导致应变也全是零，算例中我是加了三轴应力的，很奇怪！希望高手能给予帮助！谢谢！ 能看看你的模型编译吗 aizhongluo 发表于 2011-9-26 10:44 static/image/common/back.gif能看看你的模型编译吗
#include "HTC.h"
#include "math.h"
#include "contable.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 mPlastNow = 0x01;/* state logic */
//static const unsigned long mCompressNow= 0x02;
static const unsigned long mPlastPast = 0x04;
//static const unsigned long mCompressPast = 0x08;
// One static instance is neccessary as a part of internal registration process of the model with FLAC/FLAC3D
static HTC HTCmodel(true);
HTC::HTC(bool bRegister)
:ConstitutiveModel(mnHTC,bRegister), dBulk(0.0),
dShear(0.0), dFtfc(0.0), dFc(0.0), dYoung(0.0), dPoisson(0.0),
dFsin(0.0), dFcos(0.0),
dE1(0.0),dE2(0.0),dG2(0.0),
dFsurf(0.0),Lode(0.0),I1(0.0) {
}
const char **HTC::Properties(void) const {
static const char *strKey[] = {
"bulk", "shear","ftfc","fc",
"young","poisson" , 0
};
return(strKey);
}
const char **HTC::States(void) const {
static const char *strKey[] = {
"plast-n","plast-p",0
};
return(strKey);
}
/** Note: Maintain order of property input/output
*/
double HTC::GetProperty(unsigned ul) const {
switch (ul) {
case 1:return(dBulk);
case 2:return(dShear);
case 3:return(dFtfc);
case 4:return(dFc);
case 5:return(dYoung);
case 6:return(dPoisson);
}
return(0.0);
}
void HTC::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: dFtfc = dVal;break;
case 4: dFc = dVal; break;
case 5: {
dYoung = dVal;
BulkShearFromYoungPoisson(&dBulk,&dShear,dYoung,dPoisson);
break;
}
case 6: {
if ((dVal==0.5)||(dVal==-1.0)) return;
dPoisson = dVal;
BulkShearFromYoungPoisson(&dBulk,&dShear,dYoung,dPoisson);
break;
}
}
}
const char *HTC::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);
HTC *mm = (HTC *)cm;
dBulk = mm->dBulk;
dShear = mm->dShear;
dFtfc = mm->dFtfc;
dFc = mm->dFc;
dYoung = mm->dYoung;
dPoisson= mm->dPoisson;
return(0);
}
const char *HTC::Initialize(unsigned uDim,State *) {
if ((uDim!=2)&&(uDim!=3)) return("Illegal dimension in HTC constitutive model");
dE1 = dBulk + d4d3 * dShear;
dE2 = dBulk - d2d3 * dShear;
dG2 = 2.0 * dShear;
//以下是参数选用（ftfc=1/2/3分别对应ft/fc=0.08、0.10、0.12）
if (dFtfc <1.5){
dPopa = 2.8367;
dPopb = 0.5659;
dPopc = 11.8254;
dPopd = 0.2723;
}
else if ((dFtfc <2.5)&&(dFtfc >1.5)){
dPopa = 2.0107;
dPopb = 0.9713;
dPopc = 9.1413;
dPopd = 0.2310;
}
else if (dFtfc > 2.5){
dPopa = 1.4608;
dPopb = 1.2411;
dPopc = 7.3549;
dPopd = 0.2035;
}
return(0);
}
const char *HTC::Run(unsigned uDim,State *ps){
if ((uDim!=3)&&(uDim!=2)) return("Illegal dimension in HTC constitutive model");
if(ps->dHystDampMult > 0.0) HDampInit(ps->dHystDampMult);
/* --- plasticity indicator: */
/* store 'now' info. as 'past' and turn 'now' info off ---*/
int iPlas = 0;
if (ps->mState & mPlastNow)
ps->mState = (unsigned long)(ps->mState | mPlastPast);
ps->mState = (unsigned long)(ps->mState & ~mPlastNow);
/* --- trial elastic stresses --- */
ps->stnS.d11 += ps->stnE.d11 * dE1 + (ps->stnE.d22 + ps->stnE.d33) * dE2;
ps->stnS.d22 += (ps->stnE.d11 + ps->stnE.d33) * dE2 + ps->stnE.d22 * dE1;
ps->stnS.d33 += (ps->stnE.d11 + ps->stnE.d22) * dE2 + ps->stnE.d33 * dE1;
ps->stnS.d12 += ps->stnE.d12 * dG2;
ps->stnS.d13 += ps->stnE.d13 * dG2;
ps->stnS.d23 += ps->stnE.d23 * dG2;
/* --- mean stress --- */
double dSigi = (ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33) / 3.0;
/* --- deviatoric stresses --- */
STensor stnD;
stnD.d11 = ps->stnS.d11 - dSigi;
stnD.d22 = ps->stnS.d22 - dSigi;
stnD.d33 = ps->stnS.d33 - dSigi;
stnD.d12 = ps->stnS.d12;
stnD.d13 = ps->stnS.d13;
stnD.d23 = ps->stnS.d23;
/* --- second deviatoric stress invariant --- */
double dJ2 = 0.5 * (stnD.d11 * stnD.d11 +
stnD.d22 * stnD.d22 +
stnD.d33 * stnD.d33)
+ ps->stnS.d12 * ps->stnS.d12
+ ps->stnS.d13 * ps->stnS.d13
+ ps->stnS.d23 * ps->stnS.d23;
doubledJ3 = stnD.d11 * stnD.d22 * stnD.d33 + 2.0 * ps->stnS.d12 * ps->stnS.d23 * ps->stnS.d13
- stnD.d11*ps->stnS.d23 * ps->stnS.d23
- stnD.d22*ps->stnS.d13 * ps->stnS.d13
- stnD.d33*ps->stnS.d12 * ps->stnS.d12;
doubleshang=(-3.0)*sqrt(3.0)*dJ3;
doublexia=2.0*dJ2*sqrt(dJ2);
doublebizhi=shang/xia;
Lode = asin(bizhi);
I1 = ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33;
dFsin = sin(Lode/3.0+120.0*dDegRad);
dFcos = cos(Lode/3.0+120.0*dDegRad);
//定义对J2的偏导
double dJ21 = stnD.d11;
double dJ22 = stnD.d22;
double dJ23 = stnD.d33;
double dJ24 = 2.0*stnD.d12;
double dJ25 = 2.0*stnD.d23;
double dJ26 = 2.0*stnD.d13;
//定义对J3的偏导
double dJ31 = stnD.d22*stnD.d33-stnD.d23*stnD.d23+dJ2/3.0;
double dJ32 = stnD.d33*stnD.d11-stnD.d13*stnD.d13+dJ2/3.0;
double dJ33 = stnD.d22*stnD.d11-stnD.d12*stnD.d12+dJ2/3.0;
double dJ34 = 2.0*(stnD.d12*stnD.d13-stnD.d11*stnD.d23);
double dJ35 = 2.0*(stnD.d12*stnD.d23-stnD.d22*stnD.d13);
double dJ36 = 2.0*(stnD.d23*stnD.d13-stnD.d33*stnD.d12);
//定义系数C
double CC1 = dPopc/3.0+dPopd ;
double CC2 = dPopb+dPopc*2.0*dFsin/sqrt(3.0)-tan(Lode)*dFcos*dPopc*2.0/sqrt(3.0);
double CC3 = dPopc*(-1.0)*dFcos/(dJ2*cos(Lode)) ;
// if (sata > 29.0){
// if((Lode/dDegRad)> 0.0)
// {
// CC1 = dPopc/3.0+dPopd ;
// CC2 = dPopb+dPopc*2.0/sqrt(3.0);
// CC3 = 0.0 ;
//}
// else {
// CC1 = dPopc/3.0+dPopd ;
// CC2 = dPopb+dPopc/sqrt(3.0);
// CC3 = 0.0 ;
//}
// }
double dE11 = ps->stnE.d11;
double dE22 = ps->stnE.d22;
double dE33 = ps->stnE.d33;
double dE12 = ps->stnE.d12;
double dE23 = ps->stnE.d23;
double dE13 = ps->stnE.d13;
//定义对F的偏导
doubledF1 = CC1+CC2*dJ21/(2.0*sqrt(dJ2))+CC3*dJ31;
doubledF2 = CC1+CC2*dJ22/(2.0*sqrt(dJ2))+CC3*dJ32;
doubledF3 = CC1+CC2*dJ23/(2.0*sqrt(dJ2))+CC3*dJ33;
doubledF4 = CC2*dJ24/(2.0*sqrt(dJ2))+CC3*dJ34;
doubledF5 = CC2*dJ25/(2.0*sqrt(dJ2))+CC3*dJ35;
doubledF6 = CC2*dJ26/(2.0*sqrt(dJ2))+CC3*dJ36;
//定义应力的塑性变量
doubleP1= dE1*dF1+dE2*(dF2+dF3);
doubleP2= dE1*dF2+dE2*(dF1+dF3);
doubleP3= dE1*dF3+dE2*(dF2+dF1);
doubleP4= dG2*dF4;
doubleP5= dG2*dF5;
doubleP6= dG2*dF6;
double P11= P1*P1*dE11+P1*P2*dE22+P1*P3*dE33+P1*P4*dE12+P1*P5*dE23+P1*P6*dE13;
double P22= P2*P1*dE11+P2*P2*dE22+P2*P3*dE33+P2*P4*dE12+P2*P5*dE23+P2*P6*dE13;
double P33= P3*P1*dE11+P3*P2*dE22+P3*P3*dE33+P3*P4*dE12+P3*P5*dE23+P3*P6*dE13;
double P44= P4*P1*dE11+P4*P2*dE22+P4*P3*dE33+P4*P4*dE12+P4*P5*dE23+P4*P6*dE13;
double P55= P5*P1*dE11+P5*P2*dE22+P5*P3*dE33+P5*P4*dE12+P5*P5*dE23+P5*P6*dE13;
double P66= P6*P1*dE11+P6*P2*dE22+P6*P3*dE33+P6*P4*dE12+P6*P5*dE23+P6*P6*dE13;
double dPP = P1*dF1+P2*dF2+P3*dF3+P4*dF4+P5*dF5+P6*dF6;
double Dep1= P11/dPP;
double Dep2= P22/dPP;
double Dep3= P33/dPP;
double Dep4= P44/dPP;
double Dep5= P55/dPP;
double Dep6= P66/dPP;
/* --- HTC failure criterion --- */
dFsurf = dPopa*dJ2/dFc+dPopb*sqrt(dJ2)+dPopc*((2.0/sqrt(3.0))*sqrt(dJ2)*dFsin+I1/3.0)+dPopd*I1-dFc;
/* --- tests for failure */
if ( dFsurf > 0.0) {
iPlas = 1;
/* --- Plast failure: correction to ps->incips->l stresses ---*/
ps->mState = ps->mState | mPlastNow;
ps->stnS.d11 -= Dep1;
ps->stnS.d22 -= Dep2;
ps->stnS.d33 -= Dep3;
ps->stnS.d12 -= Dep4;
ps->stnS.d23 -= Dep5;
ps->stnS.d13 -= Dep6;
}
return(0);
}
const char *HTC::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(6,0);
mso->Save(0,dBulk);
mso->Save(1,dShear);
mso->Save(2,dFtfc);
mso->Save(3,dFc);
mso->Save(4,dYoung);
mso->Save(5,dPoisson);
return(0);
}
// EOF 我也做本构模型开发的，如果你愿意，可以互相学习下，QQ301069165 dE1 = dBulk + d4d3 * dShear;
dE2 = dBulk - d2d3 * dShear;
dG2 = 2.0 * dShear;
这个玩意在run 里面也做一次试试
页:
[1]