- 积分
- 0
- 注册时间
- 2007-5-13
- 仿真币
-
- 最后登录
- 1970-1-1
|
本帖最后由 whutzll 于 2011-2-28 20:49 编辑
如题，希望指教
/*------»ìÄýÍÁËðÉ˱¾¹¹Ä£ÐÍ by ZHANAG-------*/
/*------ÐÞ¸Äʱ¼ä£º 10/08/2009--------------------*/
#include "userdamage.h"
#include "contable.h"
#include <stdlib.h> // for __max definition
#include <math.h>
//variables used by all model objects. Hence only one copy is maintained for all objects
static const double dC4D3 = 4.0 / 3.0;
static const double dC2D3 = 2.0 / 3.0;
static const double dC1D3 = 1.0 / 3.0;
static const double dPi = 3.141592653589793238462643383279502884197169399;
static const double dDegRad = dPi/180.0;
static const double dSq3 = 1.732050807568877293527;
// Plasticity Indicators
static const unsigned long mDShearNow = 0x0400; /* state logic */
static const unsigned long mDTensionNow = 0x0800;
static const unsigned long mDShearPast = 0x1000;
static const unsigned long mDTensionPast = 0x2000;
// One Static instance isneccessary as a part of internal registration process of the model with FLAC/FLAC3D
static UserDamageModel mUserDamageModel(true);
UserDamageModel::UserDamageModel(bool bRegister)
:ConstitutiveModel(mnUserDamageModel,bRegister),
dShear(0.0), dYoung(0.0),dPoisson(0.0),dCohesion(0.0),dFriction(0.0),dDilation(0.0),
dTension(0.0),dCd(0.0),dSHP(0.0),dTHP(0.0),dAlpha(0.0),dBeta(0.0),dGamma(0.0),
uiCohesion(0),uiFriction(0),uiDilation(0),uiCd(0),uiBeta(0) {
}
const char *UserDamageModel::Keyword(void) const { return("DamageSoftening"); }
const char *UserDamageModel::Name(void) const { return("DamageSoftening"); }
const char **UserDamageModel: roperties(void) const {
static const char *strKey[] = {
"bulk", "shear","young","poisson",
"cohesion","friction","dilation" ,"cd","tension","alpha","beta","gamma",
"cohtab","fritab","diltab","tentab","cdtab","betab",
0
};
return(strKey);
}
const char **UserDamageModel::States(void) const {
static const char *strKey[] = {
"damage-shear-n", "damage-tension-n", "damage-shear-p", "damage-tension-p",
0
};
return(strKey);
}
/* * Note: Maintain order of property input/output
*/
double UserDamageModel::GetProperty(unsigned ul) const {
switch (ul) {
case 1: return(dBulk);
case 2: return(dShear);
case 3: return(dYoung);
case 4: return(dPoisson);
case 5: return(dCohesion);
case 6: return(dFriction);
case 7: return(dDilation);
case 8: return(dTension);
case 9: return(dCd);
case 10: return(dAlpha);
case 11: return(dBeta);
case 12: return(dGamma);
case 13: return((double)iCohTab);
case 14: return((double)iFriTab);
case 15: return((double)iDilTab);
case 16: return((double)iTenTab);
case 17: return((double)iCdTab);
case 18: return((double)iBeTab);
case 19: return(dSHP);
case 20: return(dTHP);
}
return(0.0);
}
void UserDamageModel::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: {
dYoung = dVal;
BulkShearFromYoungPoisson(&dBulk,&dShear,dYoung,dPoisson);
break;
}
case 4: {
if ((dVal==0.5)||(dVal==-1.0)) return;
dPoisson = dVal;
BulkShearFromYoungPoisson(&dBulk,&dShear,dYoung,dPoisson);
break;
}
case 5: dCohesion = dVal; break;
case 6: dFriction = dVal; break;
case 7: dDilation = dVal; break;
case 8: dTension = dVal; break;
case 9: dCd = dVal; break;
case 10: dAlpha = dVal; break;
case 11: dBeta = dVal; break;
case 12: dGamma = dVal; break;
case 13: iCohTab = ConTableList::Round(dVal); break;
case 14: iFriTab = ConTableList::Round(dVal); break;
case 15: iDilTab = ConTableList::Round(dVal); break;
case 16: iTenTab = ConTableList::Round(dVal); break;
case 17: iCdTab = ConTableList::Round(dVal); break;
case 18: iBeTab = ConTableList::Round(dVal); break;
case 19: dSHP = dVal; break;
case 20: dTHP = dVal; break;
}
}
const char *UserDamageModel::Copy(const ConstitutiveModel *m) {
//Detects type mismatch error and returns error string. otherwise returns 0
const char *str = ConstitutiveModel::Copy(m);
if (str) return(str);
UserDamageModel *em = (UserDamageModel *)m;
dBulk = em->dBulk;
dShear = em->dShear;
dYoung = em->dYoung;
dPoisson = em->dPoisson;
dCohesion = em->dCohesion;
dFriction = em->dFriction;
dDilation = em->dDilation;
dTension = em->dTension;
dCd = em->dCd;
dAlpha = em->dAlpha ;
dBeta = em->dBeta;
dGamma = em->dGamma;
dSHP = em->dSHP;
dTHP = em->dTHP;
return(0);
}
const char *UserDamageModel::Initialize(unsigned,State *) {
dShearNew = dShear;
if (dFriction) {
if (dFriction < 0.0) dFriction = -dFriction;
double dValue = tan(dFriction * dDegRad);
double dApex = dCohesion / dValue;
dTension = dTension < dApex ? dTension : dApex;
}
if (iCohTab||iFriTab||iDilTab||iCdTab||iTenTab||iBeTab) {
const char *str = 0;
ConTableList *ctl = ConTableList::Instance();
if (!ctl) return("No table implementation available");
if (iCohTab) str = ctl->IndexFromID(iCohTab,&uiCohesion);
if (str) return(str);
if (iFriTab) str = ctl->IndexFromID(iFriTab,&uiFriction);
if (str) return(str);
if (iDilTab) str = ctl->IndexFromID(iDilTab,&uiDilation);
if (str) return(str);
if (iTenTab) str = ctl->IndexFromID(iTenTab,&uiTension);
if (str) return(str);
if (iCdTab) str = ctl->IndexFromID(iCdTab,&uiCd);
if (str) return(str);
if (iBeTab) str = ctl->IndexFromID(iBeTab,&uiBeta);
if (str) return(str);
}
return(0);
}
static const int Dqs = 0;
static const int Dqt = 1;
static const int K = 2;
static const int G = 3;
static const int A1 = 4;
static const int A2 = 5;
static const int G2 = 6;
const char *UserDamageModel::Run(unsigned ulDim,State *ps) {
if ((ulDim!=2)&&(ulDim!=3)) return("Illegal dimension in DamageSoftening model");
if(ps->dHystDampMult > 0.0) HDampInit(ps->dHystDampMult);
/* --- plasticity indicator:store 'now' info. as 'past' and turn 'now' info off ---*/
if (ps->mState & mDShearNow) ps->mState |= mDShearPast;
ps->mState &= ~mDShearNow;
if (ps->mState & mDTensionNow) ps->mState |= mDTensionPast;
ps->mState &= ~mDTensionNow;
int iPlas = 0;
/* --- hardening/softening:initialize stacks to calculate hardening parameters for zone --- */
if (!ps->bySubZone) {
ps->working[Dqs] = 0.0;
ps->working[Dqt] = 0.0;
ps->working[K] = dBulk ;
ps->working[G] = dShearNew;
ps->working[A1] = (dBulk + dC4D3 * dShearNew);
ps->working[A2] = (dBulk - dC2D3 * dShearNew);
ps->working[G2] = 2.0 * dShearNew;
}
/* ---²âÊÔµ¯ÐÔÓ¦Á¦ --- */
ps->stnS.d11 += ps->stnE.d11 * ps->working[A1] + (ps->stnE.d22 + ps->stnE.d33) * ps->working[A2];
ps->stnS.d22 += (ps->stnE.d11 + ps->stnE.d33) * ps->working[A2] + ps->stnE.d22 * ps->working[A1];
ps->stnS.d33 += (ps->stnE.d11 + ps->stnE.d22) * ps->working[A2] + ps->stnE.d33 * ps->working[A1];
ps->stnS.d12 += ps->stnE.d12 * ps->working[G2];
ps->stnS.d13 += ps->stnE.d13 * ps->working[G2];
ps->stnS.d23 += ps->stnE.d23 * ps->working[G2];
/* --- ƽ¾ùÓ¦Á¦ --- */
double dSigi = (ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33) / 3.0;
/* --- Æ«Ó¦Á¦ --- */
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;
/* --- Æ«Ó¦Á¦µÚ¶þ²»±äÁ¿ --- */
double dTaui = sqrt(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);
/* --- ¼ÆËãÖ÷Ó¦Á¦Öµcalculate and sort principal stresses and principal 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,ulDim, &icase, &sdif, &psdif);
double dPrinMinCopy = dPrinMin;
double dPrinMidCopy = dPrinMid;
double dPrinMaxCopy = dPrinMax;
/* --- D-P and tensile failure criteria --- */
double dRaoMax = 0.5*(abs(dPrinMax) + dPrinMax);
double dNRaoMax = 0.5*(abs(-dPrinMax) - dPrinMax);
double dQPhi = tan(dFriction * dDegRad);
/* --- ËðÉËÆÆ»µ×¼Ôò --- */
double dFSurf= -1./(1.-dAlpha) * (3.*dAlpha*dSigi + dSq3*dTaui + dBeta*dRaoMax - dGamma*dNRaoMax) + dCohesion;
/*---ËÜÐÔÐÞÕý--*/
double dTSurf = dTension - dSigi;
double dTc1=0.0, QPhi=0.0, dQPsi=0.0, dPDiv = 0.0;
/*-------------------------------------------------------------------------------------------------------*/
if (dTSurf < 0.0) {
/*------ÐÞ¸Ä by ZHANAG*/
double dBisc = sqrt(1. + dQPhi * dQPhi) - dQPhi;
double dFach = dCohesion - dTension * dQPhi;
double dPDiv = dTaui - dFach + dBisc * dTSurf;
if (dPDiv > 0.) {
iPlas = 1;
dQPsi = tan(dDilation * dDegRad);
dTc1 = (1. - dAlpha) / (3.*ps->working[G] + 3.*ps->working[K] * dQPhi * dAlpha);
double dLams = dFSurf * dTc1;
ps->mState = ps->mState| mDShearNow;
/* --- correct second deviatoric stress invariant --- */
double dTaun = dTaui + dLams * ps->working[G] * dSq3;
/* --- correct volumetric stress --- */
double dSign = dSigi + dLams * ps->working[K] * dQPsi;
/* --- correct deviatoric stresses --- */
double dRat = dTaun / dTaui;
stnD.d11 = stnD.d11 * dRat;
stnD.d22 = stnD.d22 * dRat;
stnD.d33 = stnD.d33 * dRat;
stnD.d12 = stnD.d12 * dRat;
stnD.d13 = stnD.d13 * dRat;
stnD.d23 = stnD.d23 * dRat;
/* --- new stresses µÃµ½ÐµÄÓ¦Á¦Öµ--- */
ps->stnS.d11 = stnD.d11 + dSign * (1-dCd);
ps->stnS.d22 = stnD.d22 + dSign * (1-dCd);
ps->stnS.d33 = stnD.d33 + dSign * (1-dCd);
ps->stnS.d12 = stnD.d12;
ps->stnS.d13 = stnD.d13;
ps->stnS.d23 = stnD.d23;
} else {
/* --- tensile failure À ÉìÆÆ»µ--- */
iPlas = 2;
ps->mState = ps->mState | mDTensionNow;
/*ÐÞ¸Ä by ZHANAG*/
ps->stnS.d11 += dTSurf;//* (1-dCd);
ps->stnS.d22 += dTSurf;// * (1-dCd);
ps->stnS.d33 += dTSurf;// * (1-dCd);
}
} else {
if (dFSurf < 0.0) {
/* --- shear failure ¼ôÇÐÆÆ»µ--- */
iPlas = 1;
dQPsi = tan(dDilation * dDegRad);
dTc1 = (1. - dAlpha) / (3.*ps->working[G] + 3.*ps->working[K] * dQPhi * dAlpha);
double dLams = dFSurf * dTc1;
ps->mState = ps->mState | mDShearNow;
/* --- correct second deviatoric stress invariant --- */
double dTaun = dTaui + dLams * ps->working[G] * dSq3;
/* --- correct volumetric stress --- */
double dSign = dSigi + dLams * ps->working[K] * dQPsi;
/* --- correct deviatoric stresses --- */
double dRat = dTaun / dTaui;
stnD.d11 = stnD.d11 * dRat;
stnD.d22 = stnD.d22 * dRat;
stnD.d33 = stnD.d33 * dRat;
stnD.d12 = stnD.d12 * dRat;
stnD.d13 = stnD.d13 * dRat;
stnD.d23 = stnD.d23 * dRat;
/* --- new stresses µÃµ½ÐµÄÓ¦Á¦Öµ--- */
ps->stnS.d11 = stnD.d11 + dSign * (1-dCd);
ps->stnS.d22 = stnD.d22 + dSign * (1-dCd);
ps->stnS.d33 = stnD.d33 + dSign * (1-dCd);
ps->stnS.d12 = stnD.d12;
ps->stnS.d13 = stnD.d13;
ps->stnS.d23 = stnD.d23;
}
}
/*---¸üÐÂÓ²»¯ÏµÊý---*/
if (iPlas) {
//if (iPlas != 0) {
//ps->stnS.Resoltoglob(dPrinMin,dPrinMid, dPrinMax, aDir, dPrinMinCopy,dPrinMidCopy,dPrinMaxCopy, ulDim, icase, sdif, psdif, bFast);
/* --- hadRdening padRameter accumulation --- */
if (iPlas==1) {
/* --- shear padRameter --- */
double dLams = dFSurf * dTc1;
if (dLams < 0.) dLams = -dLams;
double dValue = sqrt(3.+dC1D3*dQPsi*dQPsi);
ps->working[Dqs] += dLams * ps->dSubZoneVolume;
//ps->working[Dqs] += dLams*dValue * ps->dSubZoneVolume;
}
if (iPlas == 2) {
/* --- tensile padRameter --- */
double dLamt = dTSurf / (ps->working[K]);
if (dLamt < 0.) dLamt = -dLamt;
ps->working[Dqt] += dLamt * ps->dSubZoneVolume;
}
}
/* --- plastic padRameter incrementation and properties update ¸üÐÂÓ²»¯²ÎÊýºÍËðÉËÒò×Ó --- */
if (ps->bySubZone==ps->byTotSubZones-1) {
double dValue = 1.0 / ps->dZoneVolume;
if (ps->byOverlay==2) dValue *= 0.5;
dSHP += ps->working[Dqs] * dValue;
dTHP += ps->working[Dqt] * dValue;
if (ps->working[Dqs] > 0.0) {
const char *str=0;
if (uiCohesion) str = ConTableList::Instance()->YFromX(uiCohesion,dSHP,&dCohesion);
if (str) return(str);
if (uiFriction) str = ConTableList::Instance()->YFromX(uiFriction,dSHP,&dFriction);
if (str) return(str);
if (uiDilation) str = ConTableList::Instance()->YFromX(uiDilation,dSHP,&dDilation);
if (str) return(str);
/*if (uiTension) str = ConTableList::Instance()->YFromX(uiTension,dSHP,&dTension);
if (str) return(str);*/
if (uiCd) str = ConTableList::Instance()->YFromX(uiCd,dSHP,&dCd);
if (str) return(str);
if (uiBeta) str = ConTableList::Instance()->YFromX(uiBeta,dSHP,&dBeta);
if (str) return(str);
}
dValue = sin(dFriction * dDegRad);
if (dFriction) {
double dApex = dCohesion * cos(dFriction * dDegRad) / dValue;
dTension = dTension < dApex ? dTension : dApex;
}
if (ps->working[Dqt] > 0.0) {
if (uiTension) {
double dTemp;
const char *str = ConTableList::Instance()->YFromX(uiTension,dTHP,&dTemp);
if (str) return(str);
dTension = dTension < dTemp ? dTension : dTemp;
}
}
}
if(iPlas == 0) {
ps->bViscous = true; // allow viscous strains
} else {
ps->bViscous = false; // inhibit viscous strains
}
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 *UserDamageModel::SaveRestore(ModelSaveObject *mso) {
// Checks for type mismatch and returns error string. Otherwise 0.
const char *str = ConstitutiveModel::SaveRestore(mso);
if (str) return(str);
// 14 represents 14 properties that are doubles
// and 4 represents 4 properties that are integers
mso->Initialize(14,5);
mso->Save(0,dBulk);
mso->Save(1,dShear);
mso->Save(2,dYoung);
mso->Save(3,dPoisson);
mso->Save(4,dCohesion);
mso->Save(5,dFriction);
mso->Save(6,dDilation);
mso->Save(7,dTension);
mso->Save(8,dCd);
mso->Save(9,dAlpha);
mso->Save(10,dBeta);
mso->Save(11,dGamma);
mso->Save(12,dSHP);
mso->Save(13,dTHP);
//Save integers seperately
mso->Save(0,iCohTab);
mso->Save(1,iFriTab);
mso->Save(2,iDilTab);
mso->Save(3,iTenTab);
mso->Save(3,iCdTab);
mso->Save(4,iBeTab);
return(0);
}
void UserDamageModel::HDampInit(const double dHMult)
{
dShearNew = dShear * dHMult;
}
// EOF |
|
IP卡
狗仔卡
发表于 2011-2-18 17:45:56
提升卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主高手
