www.pudn.com > Burger.rar > userBurger.cpp, change:2015-05-29,size:7396b


#include "userBurger.h" 
#include <math.h> 
 
static const double dC1d3 = 1.0 / 3.0; 
static UserBurgerModel userburgermodel(true); 
 
 
UserBurgerModel::UserBurgerModel(bool bRegister) 
             :ConstitutiveModel(mnUserBurgerModel,bRegister),  dBulk(0.0), 
              dKshear(0.0), dMshear(0.0), dKviscosity(0.0), dMviscosity(0.0) 
			  { 
				dMekd[0]= 0.0; 
				dMekd[1]= 0.0; 
				dMekd[2]= 0.0; 
				dMekd[3]= 0.0; 
				dMekd[4]= 0.0;  
				dMekd[5]= 0.0;  
			  } 
 
 
const char **UserBurgerModel::Properties(void) const { 
  static const char *strKey[] = { "bulk","kshear","mshear","kviscosity", "mviscosity", 
                                  "k_exx","k_eyy","k_ezz","k_exy","k_exz","k_eyz", 0}; 
  return(strKey); 
} 
 
 
const char **UserBurgerModel::States(void) const { 
  static const char *strKey[] = { 0 }; 
  return(strKey); 
} 
 
 
double UserBurgerModel::GetProperty(unsigned ul) const { 
  switch (ul) { 
    case 1:  return(dBulk); 
    case 2:  return(dKshear); 
	  case 3:  return(dMshear); 
	  case 4:  return(dKviscosity); 
	  case 5:  return(dMviscosity); 
    case 6:  return(dMekd[0]); 
    case 7:  return(dMekd[1]); 
    case 8:  return(dMekd[2]); 
    case 9:  return(dMekd[3]); 
    case 10:  return(dMekd[4]); 
    case 11:  return(dMekd[5]); 
  } 
  return(0.0); 
} 
 
 
void UserBurgerModel::SetProperty(unsigned ul,const double &dVal) { 
  switch (ul) { 
    case 1: { // BULK 
      dBulk = dVal; 
	  break; 
    } 
    case 2: { // KELVIN SHEAR 
      dKshear = dVal; 
	  break; 
    } 
    case 3: { // MAXWELL SHEAR 
      dMshear = dVal; 
	  break; 
    } 
    case 4: {  // KELVIN VISCOSITY 
	    dKviscosity = dVal; 
	    break; 
    } 
    case 5: { 
      dMviscosity = dVal; 
      break; 
    } 
    case 6: {  // KELVIN strain 11 
      dMekd[0] = dVal; 
	    break; 
    } 
    case 7: {  // KELVIN strain 22 
      dMekd[1] = dVal; 
	    break; 
    } 
    case 8: {  // KELVIN strain 33 
      dMekd[2] = dVal; 
	    break; 
    } 
    case 9: {  // KELVIN strain 12 
      dMekd[3] = dVal; 
	    break; 
    } 
    case 10: {  // KELVIN strain 13 
      dMekd[4] = dVal; 
	    break; 
    } 
    case 11: {  // KELVIN strain 23 
      dMekd[5] = dVal; 
	    break; 
    } 
  } 
} 
 
 
const char *UserBurgerModel::Copy(const ConstitutiveModel *cm) { 
  const char *str = ConstitutiveModel::Copy(cm); 
  if (str) return(str); 
  UserBurgerModel *vm = (UserBurgerModel *)cm; 
  dBulk = vm->dBulk; 
  dKshear = vm->dKshear; 
  dMshear = vm->dMshear; 
  dKviscosity = vm->dKviscosity; 
  dMviscosity = vm->dMviscosity; 
  dMekd[0]  = vm->dMekd[0]; 
  dMekd[1]  = vm->dMekd[1]; 
  dMekd[2]  = vm->dMekd[2]; 
  dMekd[3]  = vm->dMekd[3]; 
  dMekd[4]  = vm->dMekd[4];  
  dMekd[5]  = vm->dMekd[5];  
  return(0); 
} 
 
 
const char *UserBurgerModel::Initialize(unsigned,State *) { 
  if(dMshear <= 0.0) dMshear = 1e-20; 
  if(dKshear <= 0.0) dKshear = 0.0; 
  if(dKviscosity <= 0.0) { 
    dKviscosity = 0.0; 
    dKshear     = 0.0; 
  } 
  if(dMviscosity <= 0.0) dMviscosity = 0.0; 
  return(0); 
} 
 
 
const char *UserBurgerModel::Run(unsigned uDim,State *ps) { 
  if ((uDim!=3)&&(uDim!=2)) return("Illegal dimension in UserBurgerModel"); 
 
  // dEkd values now stored in ps->working[] array (necessary for thread safety) 
  double tempk=0, tempm=0; 
  if (!ps->bySubZone) { 
    ps->working[0] = 0.0; 
    ps->working[1] = 0.0; 
    ps->working[2] = 0.0; 
    ps->working[3] = 0.0; 
    ps->working[4] = 0.0; 
    ps->working[5] = 0.0; 
  } 
  // Timestep 
  double dCrtdel =  (ps->bCreep ? ps->dTimeStep : 0.0); 
 
  if (dKviscosity <= 0.0) tempk = 0.0; 
  else tempk = 1.0 / dKviscosity ; 
  if (dMviscosity <= 0.0) tempm = 0.0; 
  else tempm = 1.0 / dMviscosity ; 
  double temp = 0.5 * dKshear * dCrtdel * tempk; 
 
  double dA_con = 1.0 + temp; 
  double dB_con = 1.0 - temp; 
  double dBa    = dB_con/dA_con; 
  double dBal   = dBa - 1.0; 
  temp  = (tempm + tempk / dA_con) * dCrtdel * 0.25 ; 
  double temp1  = 1.0 / (2.0 *dMshear); 
  double dX_con = temp1 + temp; 
  double dY_con = temp1 - temp; 
  double dZ_con = dCrtdel * tempk/(4.0 * dA_con); 
  double c1dxc  = 1.0 / dX_con; 
   // Partition Strains 
  double dDev = ps->stnE.d11 + ps->stnE.d22 + ps->stnE.d33 ; 
  double dDev3 = dC1d3 * dDev; 
  double dE11d = ps->stnE.d11 - dDev3; 
  double dE22d = ps->stnE.d22 - dDev3; 
  double dE33d = ps->stnE.d33 - dDev3; 
  // Partition Stresses 
  double dS0 = dC1d3 * (ps->stnS.d11 + ps->stnS.d22 + ps->stnS.d33); 
  double dS11d = ps->stnS.d11 - dS0; 
  double dS22d = ps->stnS.d22 - dS0; 
  double dS33d = ps->stnS.d33 - dS0; 
  // Remember old stresses 
  double dS11old = dS11d; 
  double dS22old = dS22d; 
  double dS33old = dS33d; 
  double dS12old = ps->stnS.d12; 
  double dS13old = 0.0; 
  double dS23old = 0.0; 
  if (uDim==3) {                     
   dS13old = ps->stnS.d13; 
   dS23old = ps->stnS.d23; 
  } 
  //new deviatoric stresses 
  dS11d = (dS11d * dY_con + dE11d - dMekd[0] * dBal)/dX_con; 
  dS22d = (dS22d * dY_con + dE22d - dMekd[1] * dBal)/dX_con; 
  dS33d = (dS33d * dY_con + dE33d - dMekd[2] * dBal)/dX_con; 
  ps->stnS.d12 = ( ps->stnS.d12 * dY_con + ps->stnE.d12 - dMekd[3] * dBal )*c1dxc; 
  if (uDim==3) {                     
   ps->stnS.d13 = ( ps->stnS.d13 * dY_con + ps->stnE.d13 - dMekd[4] * dBal )*c1dxc; 
   ps->stnS.d23 = ( ps->stnS.d23 * dY_con + ps->stnE.d23 - dMekd[5] * dBal )*c1dxc; 
  } 
  // isotropic stress is elastic 
  dS0 = dS0 + dBulk * dDev; 
  // convert back to x-y components 
  ps->stnS.d11 = dS11d + dS0; 
  ps->stnS.d22 = dS22d + dS0; 
  ps->stnS.d33 = dS33d + dS0; 
  // sub-zone contribution to kelvin-strains 
  ps->working[0] +=   (dMekd[0] * dBa + (dS11d + dS11old) * dZ_con) * ps->dSubZoneVolume; 
  ps->working[1] +=   (dMekd[1] * dBa + (dS22d + dS22old) * dZ_con) * ps->dSubZoneVolume; 
  ps->working[2] +=   (dMekd[2] * dBa + (dS33d + dS33old) * dZ_con) * ps->dSubZoneVolume; 
  ps->working[3] +=   (dMekd[3] * dBa + (ps->stnS.d12 + dS12old) * dZ_con) * ps->dSubZoneVolume; 
  if (uDim==3) {                     
   ps->working[4] +=   (dMekd[4] * dBa + (ps->stnS.d13 + dS13old) * dZ_con) * ps->dSubZoneVolume; 
   ps->working[5] +=   (dMekd[5] * dBa + (ps->stnS.d23 + dS23old) * dZ_con) * ps->dSubZoneVolume; 
  } 
  // update stored kelvin strains 
 if (ps->bySubZone==ps->byTotSubZones-1) 
  { 
	    double Aux = 1./ps->dZoneVolume; 
	    if (ps->byOverlay==2) Aux *= 0.5; 
	    dMekd[0] = ps->working[0]*Aux; 
	    dMekd[1] = ps->working[1]*Aux; 
	    dMekd[2] = ps->working[2]*Aux; 
	    dMekd[3] = ps->working[3]*Aux; 
      if (uDim==3) {                     
	      dMekd[4] = ps->working[4]*Aux; 
	      dMekd[5] = ps->working[5]*Aux; 
      } 
  } 
   
  return(0); 
} 
 
 
const char *UserBurgerModel::SaveRestore(ModelSaveObject *mso) { 
  const char *str = ConstitutiveModel::SaveRestore(mso); 
  if (str) return(str); 
  // Previous versions does not save all properties but saves only the  
  // first 9 properties...fixed 10/3/03 
  if(Version() > 1) 
      mso->Initialize(11,0); 
  else   
      mso->Initialize(9,0); 
 
  mso->Save(0,dBulk); 
  mso->Save(1,dKshear); 
  mso->Save(2,dMshear); 
  mso->Save(3,dKviscosity); 
  mso->Save(4,dMviscosity); 
  mso->Save(5,dMekd[0]); 
  mso->Save(6,dMekd[1]); 
  mso->Save(7,dMekd[2]); 
  mso->Save(8,dMekd[3]); 
  mso->Save(9,dMekd[4]);   
  mso->Save(10,dMekd[5]);  
 
  return(0); 
} 
 
// EOF