一個FLAC的本構(gòu)模型

ITASCA(武漢)咨詢有限公司依泰斯卡(武漢)內(nèi)容必要性試驗(yàn)總結(jié)的本構(gòu)模型特定條件下的本構(gòu)模型交叉學(xué)科的本構(gòu)模型二次開發(fā)環(huán)境自定義本構(gòu)模型的功能自定義本構(gòu)模型的基本方法常用模型信息傳遞指針變量模型源程序分解二次開發(fā)環(huán)境FLAC3D采用面向?qū)ο蟮恼Z言標(biāo)準(zhǔn)C++編寫本構(gòu)模型都是以動態(tài)連接庫文件(.DLL文件)的形式提供

VC++6.0(SP4)或更高版本的開發(fā)環(huán)境優(yōu)點(diǎn)自定義的本構(gòu)模型和軟件自帶的本構(gòu)模型的執(zhí)行效率處在同一個水平自定義本構(gòu)模型(.DLL文件)適用于高版本的FLAC(2D)、3DEC、UDEC等其他Itasca軟件中返回目錄自定義本構(gòu)模型的功能

主要功能：對給出的應(yīng)變增量得到新的應(yīng)力輔助功能：模型名稱、版本讀寫操作

模型文件的編寫基類(classConstitutiveModel)的描述成員函數(shù)的描述模型的注冊模型與FLAC3D之間的信息交換模型狀態(tài)指示器的描述

返回目錄自定義本構(gòu)模型的基本方法頭文件(usermodel.h)中進(jìn)行新的本構(gòu)模型派生類的聲明修改模型的ID(>100)、名稱和版本修改派生類的私有成員：私有變量和成員函數(shù)C++文件(usermodel.cpp)中修改模型結(jié)構(gòu)(UserModel::UserModel(boolbRegister):ConstitutiveModel)constchar**UserModel::Properties()函數(shù)模型的參數(shù)名稱字符串：用于在模型中顯示constchar**UserModel::States()函數(shù)計(jì)算過程中的狀態(tài)指示器：單元屈服狀態(tài)返回目錄doubleUserModel::GetProperty()和voidUserModel::SetProperty()函數(shù)constchar*UserModel::Initialize()函數(shù)參數(shù)和狀態(tài)指示器的初始化，并對派生類聲明中定義的私有變量進(jìn)行賦值constchar*UserModel::Run(unsignednDin,States*ps)函數(shù)由應(yīng)變增量計(jì)算得到應(yīng)力增量，從而獲得新的應(yīng)力constchar*UserModel::SaveRestore()函數(shù)對計(jì)算結(jié)果進(jìn)行保存。程序的調(diào)試在VC++的工程設(shè)置中將FLAC3D軟件中的EXE文件路徑加入到程序的調(diào)試范圍中，并將FLAC3D自帶的DLL文件加入到附加動態(tài)鏈接庫(AdditionalDLLs)中，然后在Initialize()或Run()函數(shù)中設(shè)置斷點(diǎn)，進(jìn)行調(diào)試；在程序文件中加入return()語句，這樣可以將希望得到的變量值以錯誤提示的形式在FLAC3D窗口中得到。常用模型信息傳遞指針變量返回目錄模型源程序分解返回目錄靜力本構(gòu)(Mohr-Coulomb)1.屈服函數(shù)剪切屈服拉屈服2.塑性勢函數(shù)非關(guān)聯(lián)關(guān)聯(lián)①②①－剪切屈服②－拉伸屈服MC本構(gòu)寫成線性函數(shù)S：剪切屈服修正寫成線性函數(shù)S：拉伸屈服修正constchar*UserMohrModel::Initialize(unsigneduDim,State*){if((uDim!=2)&&(uDim!=3))return("IllegaldimensioninUserMohrconstitutivemodel");dE1=dBulk+d4d3*dShear;dE2=dBulk-d2d3*dShear;dG2=2.0*dShear;doubledRsin=sin(dFriction*dDegRad);dNPH=(1.0+dRsin)/(1.0-dRsin);dCSN=2.0*dCohesion*sqrt(dNPH);if(dFriction){doubledApex=dCohesion*cos(dFriction*dDegRad)/dRsin;dTension=dTension<dApex?dTension:dApex;}dRsin=sin(dDilation*dDegRad);dRnps=(1.0+dRsin)/(1.0-dRsin);doubledRa=dE1-dRnps*dE2;doubledRb=dE2-dRnps*dE1;doubledRd=dRa-dRb*dNPH;dSC1=dRa/dRd;dSC3=dRb/dRd;dSC2=dE2*(1.0-dRnps)/dRd;dBISC=sqrt(1.0+dNPH*dNPH)+dNPH;dE21=dE2/dE1;return(0);}初始化各式中的常數(shù)constchar*UserMohrModel::Run(unsigneduDim,State*ps){if((uDim!=3)&&(uDim!=2))return("IllegaldimensioninMohrconstitutivemodel");if(ps->dHystDampMult>0.0)HDampInit(ps->dHystDampMult);/*---plasticityindicator:*//*store'now'info.as'past'andturn'now'infooff---*/if(ps->mState&mShearNow)ps->mState=(unsignedlong)(ps->mState|mShearPast);ps->mState=(unsignedlong)(ps->mState&~mShearNow);if(ps->mState&mTensionNow)ps->mState=(unsignedlong)(ps->mState|mTensionPast);ps->mState=(unsignedlong)(ps->mState&~mTensionNow);intiPlas=0;doubledTeTens=dTension;/*---trialelasticstresses---*/doubledE11=ps->stnE.d11;doubledE22=ps->stnE.d22;doubledE33=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;/*---calculateandsortps->incips->lstressesandps->incips->ldirections---*/AxesaDir;doubledPrinMin,dPrinMid,dPrinMax,sdif=0.0,psdif=0.0;inticase=0;boolbFast=ps->stnS.Resoltopris(&dPrinMin,&dPrinMid,&dPrinMax,&aDir,uDim,&icase,&sdif,&psdif);doubledPrinMinCopy=dPrinMin;doubledPrinMidCopy=dPrinMid;doubledPrinMaxCopy=dPrinMax;/*---Mohr-Coulombfailurecriterion---*/doubledFsurf=dPrinMin-dNPH*dPrinMax+dCSN;/*---Tensilefailurecriteria---*/doubledTsurf=dTension-dPrinMax;doubledPdiv=-dTsurf+(dPrinMin-dNPH*dTension+dCSN)*dBISC;/*---testsforfailure*/if(dFsurf<0.0&&dPdiv<0.0){iPlas=1;/*---shearfailure:correctiontops->incips->lstresses---*/ps->mState=(unsignedlong)(ps->mState|0x01);dPrinMin-=dFsurf*dSC1;dPrinMid-=dFsurf*dSC2;dPrinMax-=dFsurf*dSC3;}elseif(dTsurf<0.0&&dPdiv>0.0){iPlas=2;/*---tensionfailure:correctiontops->incips->lstresses---*/ps->mState=(unsignedlong)(ps->mState|0x02);doubledTco=dE21*dTsurf;dPrinMin+=dTco;dPrinMid+=dTco;dPrinMax=dTension;}if(iPlas){ps->stnS.Resoltoglob(dPrinMin,dPrinMid,dPrinMax,aDir,dPrinMinCopy,dPrinMidCopy,dPrinMaxCopy,uDim,icase,sdif,psdif,bFast);ps->bViscous=false;//Inhibitstiffness-dampingterms}else{ps->bViscous=true;//Allowstiffness-dampingterms}return(0);}蠕變本構(gòu)(Viscous)Modelvisc（H-N）一.力-位移定律總位移彈性部分粘性部分uHuNFFH-N體示意圖將虎克體和牛頓體的本構(gòu)關(guān)系帶入做變換可得Visc本構(gòu)力-位移定律：Maxwell體－M體ModelVisc（H-N）二.中心差分方案（將應(yīng)力張量做球應(yīng)力張量和偏應(yīng)力張量分解-增量關(guān)系）其中:偏應(yīng)力張量和偏應(yīng)變張量的本構(gòu)關(guān)式中：球應(yīng)力張量和球應(yīng)變張量的本構(gòu)關(guān)系總應(yīng)力：Modelvisc（H-N）求解系數(shù)：dD=如果蠕變指標(biāo)為真，則返回真值；否則dD=0。三.源程序分解constchar*UserViscousModel::Run(unsigneduDim,State*ps){if((uDim!=3)&&(uDim!=2))return("IllegaldimensioninUserViscousModel");doubledD=dGD2V*(ps->bCreep?ps->dTimeStep:0.0);if(dD>0.5)return("TimesteptoolargeforUserViscousModel");doubledVic1=1.0-dD;doubledVic2=1.0/(1.0+dD);/*volumetricanddeviatoricstrainincrements*/doubledDevol=ps->stnE.d11+ps->stnE.d22+ps->stnE.d33;doubledDevol3=d1d3*dDevol;doubledE11d=ps->stnE.d11-dDevol3;doubledE22d=ps->stnE.d22-dDevol3;doubledE33d=ps->stnE.d33-dDevol3;constchar*UserViscousModel::Initialize(unsigned,State*){dG2=2.0*dShear;if(dViscosity<=0.0)dGD2V=0.0;elsedGD2V=0.5*dShear/dViscosity;return(0);}求解系數(shù)：dGD2V=；1.如果粘滯系數(shù)<=0，dGD2V=0；2.否則dGD2V為真值。求解系數(shù)：C1，C2。球應(yīng)變增量張量偏應(yīng)變增量張量Modelvisc（H-N）constchar*UserViscousModel::Run(unsigneduDim,State*ps){/*oldisotropicanddeviatoricstresses*/doubledS0=d1d3*(ps->stnS.d11+ps->stnS.d22+ps->stnS.d33);doubledS11d=ps->stnS.d11-dS0;doubledS22d=ps->stnS.d22-dS0