creepb.cpp

Finite element program for mechanical problem. It can solve various problem in solid problem

  }
  
  /*
  if (ipp==0){
    fprintf (Out,"\n\n tau faze2");
    //fprintf (Out," %le %le %le %le %le %le",screep[0][0],screep[1][0],screep[2][0],screep[3][0],screep[4][0],screep[5][0]);
    fprintf (Out," %le %le %le %le",screep[0][0],screep[1][0],screep[2][0],screep[3][0]);
    //  Stress data storage first position, Hide (strain) seconde position
    fprintf (Out,"\n\n *****prir sigma phase 2");
    fprintf (Out,"\n %e %e %e ",sig[0],sig[1],sig[2]);
    //fprintf (Out,"\n %e %e %e %e %e %e",sig[0],sig[1],sig[2],sig[3],sig[4],sig[5]);
  }
  */

  Mm->ip[ipp].eqother[nc*(nRetTime+1)+0]+=desht;
  Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]=h_s;
  Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]=temp_s;
  //       fprintf (Out,"\n\n desht hs, temps faze2");
  //       fprintf (Out," %le  %le %le ",desht, h_s, temp_s);

  //  irreversible strain (from t=0)
    for (i=0; i<nc; i++){
      Mm->ip[ipp].eqother[nc*(nRetTime+1)+3+i]+=depstot[i]-deps[i];
      Mm->ip[ipp].eqother[nc*(nRetTime+2)+3+i]+=sig[i];
    }
}

void creepb::get_h (long ipp)
{
  h_slast=1.0;
  if(type_h !=1 ){h_slast=Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]; }
  
  //if (ipp==0)
  //fprintf (Out,"\n h_slast %e",h_slast);
  
}


void creepb::get_temp (long ipp)
{
  temp_slast=283.0;   //10 C
  if(type_temp != 1 ){temp_slast=Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]; }
  
  //if (ipp==0)
  //fprintf (Out,"\n temp_slast %e",temp_slast);
  
}



/**
   function apdates material parameters
   
   @param E - function of creep

   10.10.2002
*/
void creepb::matstiff (matrix &d, long ipp)
//  function returns elastic stiffness matrix
//  d - elastic stiffness matrix
{
  double qq=1.;
  
  nc=Mm->ip[ipp].ncompstr;
  ss=Mm->ip[ipp].ssst;  //planestress=10,planestrain=11,plate=15,
                        //axisymm=20,shell=25,spacestress=30
  imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
   //  initial elastic stiffness matrix
  Mm->elmatstiff (d,ipp);
  //  initial elastic compliance matrix    Mm->elmatcompl (c,ipp);
  e0=Mm->eliso[imat].e;
  mi=Mm->eliso[imat].nu;    //    d[1][0] = nu*e/(1.0-nu*nu);
  
  if(type_h!=1 && type_temp!=1){
    matstiffchange(qq, ipp);
  }
  else{
    matstiffkons(qq);
    //fprintf (Out,"\n ipp %ld   et %le e0 %le",ipp,et,e0);
  }
  //  New material modulus in Time
  

  cmulm( qq, d);
}
void creepb::matstiffkons (double &qq)
//  function returns elastic stiffness matrix
{
  //  New material modulus in Time
  qq=et/e0;
}

void creepb::matstiffchange (double &qq, long ipp)
//  function returns elastic stiffness matrix
//  d - elastic stiffness matrix
{
  long i,j,k;
  double jt,pomt,pomt1,q,delYi,m0,m;
  vector bpom(nRetTime);
  
  if(Mm->ip[ipp].hmt & 1) 
	{i=Mm->ip[ipp].idm[Mm->ip[ipp].nm-1]; alfa=Mm->tidilat[i].alpha;}
 // initial and current moist and tempr
  if(type_h!=1) {h_s=Mm->moist[ipp];}   //actual moist
  if(type_temp!=1) {temp_s=Mm->tempr[ipp];}  //ctual tempr
  get_h (ipp);     //initial  moist
  get_temp (ipp);     //initial  tempr

  pomt=ccTime+ddTime/2.;
  q=2./3.;
  jt=0.0;
  //   from concrete starts   tb 
  if(timemat<=pomt) {
    et=0.;
    for (i=0;i<nRetTime;i++){
      bpom[i]=0.;
/*      for (j=0;j<nRetTime;j++){
	    apom[i][j]=0.;
      }*/
    }

    pomt1=pomt+0.0001;   
    m0 = log10(pomt1);
    m = log10(2.*timeMax)*0.999;
    k=0;
    while (pomt1<=2.*timeMax){
      //   creep function from time pomt to time pomt1
      {b3_law(jt, pomt, pomt1);}
      //		fprintf (Out,"\n\n jt v case t   od t0 prirustek dt");
      //		fprintf (Out," %e\t%e\t%e\t%e",jt, pomt, ccTime, ddTime);
      for (i=0;i<nRetTime;i++){
	   delYi=pow((pomt/retTime[i]),q)-pow((pomt1/retTime[i]),q);
	   bpom[i]=bpom[i]+(1.-exp(delYi))*jt;
	   for (j=0;j<nRetTime;j++){
//	  delYj= pow((pomt/retTime[j]),q)-pow((pomt1/retTime[j]),q);
//	  apom[i][j]=apom[i][j]+(1.-exp(delYi))*(1.-exp(delYj));
	   }
      }
      k=k+1;
	  pomt1=pow(10.,m0+(m-m0)/50*k);
//      pomt1=pomt1+20*(k);
    }
    
//    inv_sym(apom);
    for (i=0;i<nRetTime;i++){
      ert[i]=0.;
      for (j=0;j<nRetTime;j++){
	   ert[i]=ert[i]+apom[i][j]*bpom[j];
      }
      ert[i]=1./ert[i];
      if(ccTime<0.0001){
	   delYi=pow((0.0001+ddTime)/retTime[i],q)-pow(0.0001/retTime[i],q);}
      else {
	   delYi=pow((ccTime+ddTime)/retTime[i],q)-pow(ccTime/retTime[i],q);}
       et=et+(delYi-1.+exp(-delYi))/delYi/ert[i];
    }
    if(et>1.e-50)  et=1./et;
    else{fprintf (Out,"\n\n Stiff is zero"); abort();}
  }
  
  else {
    et=1.;
    for (i=0;i<nRetTime;i++){ ert[i]=0.;}
  }
  
  //		fprintf (Out,"\n\n eh v case t   od t0 prirustek dt");
  //		fprintf (Out," %e\t%e\t%e\t%e",et, pomt, ccTime, ddTime);
  
  //  New material modulus in Time
  qq=et/e0;

}

/**
   function returns history
   
   @param screep - vector of history
   @param sig - vector od add stress
   
   10.10.2002
*/
void creepb::seps_time (matrix &screep,vector &sig)
{
  int i;
  double q,dj,pomt,pp;
  
  q=2./3.;
  if((ccTime+ddTime)==ddTime) {
  }
  else{
    pomt=ccTime+ddTime;
	if(ss==planestress){
	 for (i=0;i<nRetTime;i++){
		dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
		screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
		screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0])*(1.-exp(-dj))/dj/ert[i];
		screep[2][i]=exp(-dj)*screep[2][i]+ sig[2]*2.*(1.+mi)*(1.-exp(-dj))/dj/ert[i];
		screep[3][i]=0.0;
	 }
	}
	else if(ss==planestrain){  //sx,sy,sy,txy
	 pp=1.+mi;
	 for (i=0;i<nRetTime;i++){
		dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
		screep[0][i]=exp(-dj)*screep[0][i]+ pp*(sig[0]*(1.-mi)-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
		screep[1][i]=exp(-dj)*screep[1][i]+ pp*(sig[1]*(1.-mi)-mi*sig[0])*(1.-exp(-dj))/dj/ert[i];
		screep[2][i]=exp(-dj)*screep[2][i]+ pp*sig[2]*2.   *(1.-exp(-dj))/dj/ert[i];
		screep[3][i]=0.0;
	 }
	}
	else if(ss==axisymm){
	 for (i=0;i<nRetTime;i++){
		dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
		screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
		screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
		screep[2][i]=exp(-dj)*screep[2][i]+(sig[2]-mi*sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
		screep[3][i]=exp(-dj)*screep[3][i]+ sig[3]*(1.+mi)   *(1.-exp(-dj))/dj/ert[i];
	 }
	}
	else if(ss==spacestress){
	 for (i=0;i<nRetTime;i++){
		dj=pow(pomt/retTime[i],q)-pow(ccTime/retTime[i],q);
		screep[0][i]=exp(-dj)*screep[0][i]+(sig[0]-mi*sig[1]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
		screep[1][i]=exp(-dj)*screep[1][i]+(sig[1]-mi*sig[0]-mi*sig[2])*(1.-exp(-dj))/dj/ert[i];
		screep[2][i]=exp(-dj)*screep[2][i]+(sig[2]-mi*sig[0]-mi*sig[1])*(1.-exp(-dj))/dj/ert[i];
		screep[3][i]=exp(-dj)*screep[3][i]+ sig[3]*(1.+mi)   *(1.-exp(-dj))/dj/ert[i];
		screep[4][i]=exp(-dj)*screep[4][i]+ sig[4]*(1.+mi)   *(1.-exp(-dj))/dj/ert[i];
		screep[5][i]=exp(-dj)*screep[5][i]+ sig[5]*(1.+mi)   *(1.-exp(-dj))/dj/ert[i];
	 }
	}
  }

}



/**
   function computes eps  by Prof. Bazant
   
   @param des_hn - function computing of shrink and temperature strain

   10.9.2004
*/
void creepb::get_desht (double &des_hn, double t0, double t)
{
  double esn,kh,tau,ht,ht0,dh_s,dtemp;

/*
  from t0  to t 
  k_s   shape factor slab=1.0, cylinder=1.15, sguare prism.=1.25, sphere=1.3, cube=1.55
  tb    from concrete starts
  t_w   age when drying begins 
  (fc') is 28 day average cilinder strenght fc' [ksi] ksi=1000psi=6.895 MPa(f.e.6.454=44.5MPa)***6.381
  (k_d) effective cross section thickness       D=2*vs_s
   h_s  actual relative humidity 
   h_slast   relative humidity at time t0
   temp_s    actual temperature in K 
   temp_slast  temperature at time t0
*/
  
  if (tb<=t0){
 
  //  shrink
   if ( type_h ==1){
      if ( (t0>=t_w) ){        //   if ( (t0>t_w) && (t<=(t_w+600.)) ){
  //  shrink from huminidy   from  1 to h_s
        tau=k_s*k_s*k_d*k_d*85000./pow(t,0.08)/pow(fc,0.25);
        esn=-a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
	    ht0 =tanh(sqrt((t0-t_w)/tau));
	    ht  =tanh(sqrt((t -t_w)/tau));
		if(h_s<=0.98){kh=(1.0-pow(h_s,3.0));}
		else {kh=-10.0*(1.-h_s);}
        des_hn = esn*kh/1000000.0*(ht-ht0); //increment of eps
	  }
	  else {
		des_hn=0.0;
	  }
   }
   else {
  //  shrink from changeing humidity
     dh_s = h_s-h_slast;   // humidity
     esn = a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
     kh = 3.0*h_s*h_s;
     des_hn = esn*kh*dh_s*1.e-6;
   }
   if ( type_temp ==1 ){
  //  temperature constant
   }
   else {
  //  temperature
     dtemp = temp_s-temp_slast;   
     esn=alfa*1.e6;
     des_hn=des_hn+esn*dtemp*1.e-6;
   }
  }

  else {
    des_hn=0.0;
  }
}


/**
   function computes J(t,t')  by Prof. Bazant    (Jt = strain)
   
   @param jt - function of creep

   10.10.2002
*/
void creepb::b3_law (double &jt, double t0, double t)
{
  double n,c0,cd,  ac,ag,m,x,y,q,q1,q2,q3,q4,q5,esn,kh,tau,ht,ht0,dh_s,dtemp;
//  call only from matstiff 
/*
  from t0  to t 
  K_s  shape factor slab=1.0, cylinder=1.15, sguare prism.=1.25, sphere=1.3, cube=1.55 
  tb   from concrete starts
  t_w  age when drying begins 
  (fc') is 28 day average cilinder strenght fc' [ksi] ksi=1000psi=6.895 MPa(f.e.6.454=44.5MPa)***6.381
  (w/c) is water-cement ratio of the mix by weight   ***0.43
  (s/c) is send-cement ratio of the mix by weight    ***3.4
  (g/c) is gravel-cement ratio of the mix by weight g/c=a/c-s/c     ***1.98
  (a/c) is aggregate-cement ratio of the mix by weight a/c=g/c+s/c   
  (a1)  is coef. for cements of type I,II a1=1.00, III a1=0.93, IV a1=1.05   ***1.05
  (ro)  is mass of concrete in [lb/ft3] =16.03 kg/m3 ***156
  (tl)  effective cross section thickness       D=2*vs_s
   cs  cement content in m3  .. kg/m3
     E0=(0.09+1/(1.7*(0.5*ro*ro*fc*1e-4)*(0.5*ro*ro*fc*1e-4)))
     Et=E0*sqrt(t/(4+0.85*t))          podle ACI Commite 209/II

*/
  
   if (tb<t0){
   if (t_w<tb) t_w = tb;
   if (t_w<=0.0) t_w = 1.0;
   ac=sc+gc;   ag=ac/gc;   m=0.5; //m=0.28+47541025/fc/fc;
    
   x=(2.1*ac/pow(sc,1.4)+0.1*pow((fc/6895),1.5)*pow(wc,(1.0/3.0))*pow(ag,2.2))*a1-4;
   if (x<=4){n=0.12;}
   else {n=0.12+(0.07*pow(x,6.0)/(5130.0+pow(x,6.0)));}
   q1=600000.0/4734.0/sqrt(fc);
   x=1.7*pow(t0,0.12)+8.0;
   y=(0.086*pow(t0,(2.0/9.0))+1.21*pow(t0,(4.0/9.0)));
   q=1.0/y/pow( ( 1.0+ pow( (pow(t0,m)/log(1.+pow((t-t0),n))/y),x ) ),(1./x));
   q2=185.4*sqrt(c_s)/pow(fc,0.9);
   q3=0.29*pow(wc,4.0)*q2;
   q4=20.3/pow(ac,0.7);
  
   c0=q2*q+q3*log(1.0+pow((t-t0),n))+q4*log(t/t0);
 
  //  shrink and drying
   if ( type_h ==1){
      if ( (t0>=t_w) ){        //   if ( (t0>t_w) && (t<=(t_w+600.)) ){
  //  shrink from huminidy=1 to h_s
        tau=k_s*k_s*k_d*k_d*85000./pow(t,0.08)/pow(fc,0.25);
        esn=a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
	    ht0 =tanh(sqrt((t0-t_w)/tau));
	    ht  =tanh(sqrt((t -t_w)/tau));
		if(h_s<=0.98){kh=(1.0-pow(h_s,3.0));}
		else {kh=-10.0*(1.-h_s);}
  //  drying from huminidy=1 to h_s 
        q5 = 0.757/fc/pow(esn,0.6);
        cd = q5*sqrt( exp(-8+8.*(1.0-h_s)*ht)-exp(-8.+8.*(1.0-h_s)*ht0) );
	  }
	  else {
		cd = 0.0;
	  }
   }
   else {
  //  shrink and drying from changeing humidity 
     dh_s = h_s-h_slast;   // humidity
     esn = a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
     kh = 3.0*h_s*h_s;
	 q5 = 0.35/fc;
     cd = fabs(q5*esn*kh*dh_s);
   }
   if ( type_temp ==1 ){
  //  temperature constant
   }
   else {
  //  temperature
     dtemp = temp_s-temp_slast;   
     esn=alfa*1.e6;
	 q5= 1.5/fc;
     cd = cd+fabs(q5*esn*dtemp);
   }
  }

  else {
    q1 = 0.0;
    c0 = 0.0;
    cd = 0.0;
  }
  //  measure static modulus - creep modulus  1/Ec=2/(3 Es)
  //  1/E=e-6/psi =e-6/6.895kPa
  //Jt = strain
  jt=(q1+c0+cd)*1.e-6;     
//    fprintf (Out,"\n t0, t,  q1,q2,q3,q4");
//    fprintf (Out," %e %e\t%e\t%e\t%e\t%e",t0,t,q1,q2*q,q3*log(1+pow((t-t0),n)),q4*log((t-tb)/(t0-tb)));
}



/**
  Function returns irreversible plastic strains.

  @param ipp   - integration point number in the mechmat ip array.
  @param ido   - index of the first internal variable for given material in the ipp other array
  @param epscr  - %vector of irreversible strains
 
  Returns vector of irreversible strains via parameter epscr
*/
void creepb::giveirrstrains (long ipp, long ido, vector &epscr)
{
  long i;
  long ncompstr = Mm->ip[ipp].ncompstr;
  long id=ncompstr+ncompstr*nRetTime+3;
  for (i=0;i<epscr.n;i++)
    epscr[i] = Mm->ip[ipp].eqother[ido+id+i];
}