creepb-f.cpp

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

//  d - elastic stiffness matrix
{
  long i,j,k;
  double jt,et,pomt,pomt1,qq,delYi,delYj,m0,m;
  vector bpom(nRetTime);
//  matrix apom(nRetTime,nRetTime);
  
  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);
  // therm coeff.
  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

  //ddTime=Mp->timefun.getval(Mp->time);
  ddTime=Mp->timecon.actualforwtimeincr ();
  //ccTime=Mp->time;
  ccTime=Mp->timecon.actualtime ();
  pomt=ccTime+ddTime/2.;
  qq=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 od casu pomt do casu 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]),qq)-pow((pomt1/retTime[i]),qq);
	bpom[i]=bpom[i]+(1.-exp(delYi))*jt;
	for (j=0;j<nRetTime;j++){
//	  delYj= pow((pomt/retTime[j]),qq)-pow((pomt1/retTime[j]),qq);
//	  apom[i][j]=apom[i][j]+(1.-exp(delYi))*(1.-exp(delYj));
	}
      }
      k=k+1;
	  pomt1=pow(10.,m0+(m-m0)/40*k);
//      pomt1=pomt1+20*(k);
    }
    
//    inv_sym(apom);
    for (i=0;i<nRetTime;i++){
      //	      ert[i][imat]=0.;
      ert[i]=0.;
      for (j=0;j<nRetTime;j++){
	//	       ert[i][imat]=ert[i][imat]+apom[i][j]*bpom[j];
	ert[i]=ert[i]+apom[i][j]*bpom[j];
      }
      //	      ert[i][imat]=1./ert[i][imat];
      ert[i]=1./ert[i];
      if(ccTime<0.0001){
	delYi=pow((0.0001+ddTime)/retTime[i],qq)-pow(0.0001/retTime[i],qq);}
      else {
	delYi=pow((ccTime+ddTime)/retTime[i],qq)-pow(ccTime/retTime[i],qq);}
      //          lam[i]=1-(1.-(exp(-delYi)))/delYi;
      //          ert[(nRetTime+1][imat]=ert[nRetTime+1][imat]+(delYi-1.+exp(-delYi))/delYi/ert[i][imat];
      et=et+(delYi-1.+exp(-delYi))/delYi/ert[i];
    }
    //	    ert[nRetTime+1][imat]=1./ert[nRetTime+1][imat];
    et=1./et;

    //fprintf (stderr,"\n et %lf",et);
    //fprintf (stderr,"\n apom %lf",apom[0][0]);
    //fprintf (stderr,"\n bpom %lf %lf %lf %lf %lf %lf %lf",bpom[0],bpom[1],bpom[2],bpom[3],
    //bpom[4],bpom[5],bpom[6]);

    //         call Jkontrola(imat)
    
  }
  
  else {
    //        ert[nRetTime+1][imat]=e0/100000.;
    et=e0/100000.;
    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;
  
  cmulm( qq, d);
}

/**
   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 qq,dj,pomt,pp;
  
  qq=2./3.;
  if((ccTime+ddTime)==ddTime) {
  }
  else{
    pomt=ccTime+ddTime;
	if(ss==planestress){
	 for (i=0;i<nRetTime;i++){
		dj=pow(pomt/retTime[i],qq)-pow(ccTime/retTime[i],qq);
		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]*(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],qq)-pow(ccTime/retTime[i],qq);
		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],qq)-pow(ccTime/retTime[i],qq);
		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],qq)-pow(ccTime/retTime[i],qq);
		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 n,c0,cd,  ac,ag,m,x,y,q,q1,q2,q3,q4,q5,esn,kh,tau,ht,ht0,dh_s,dtemp;
  short experiment=0;

/*
  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=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 and drying from computing
     dh_s = h_s-h_slast;   // humidity
     esn = -a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
  //  shrink
     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;
  //  shrink
     des_hn=des_hn+esn*dtemp*1.e-6;
   }
  }

  else {
    des_hn=0.0;
  }
}


/**
   function computes J(t,t')  by Prof. Bazant
   
   @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;
  short experiment=0;
//  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);}
  //      des_hn = esn*kh/1000000.0*(ht-ht0); //increment of eps
  //  drying from huminidy=1 to h_s 
        q5 = 0.757/fc/pow(fabs(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;
//		des_hn=0.0;
	  }
   }
   else {
  //  shrink and drying from computing
     dh_s = h_s-h_slast;   // humidity
     esn = -a1*1.2*(0.019*pow((wc*c_s),2.1)/pow(fc,0.28)+270.);
  //  shrink
     kh = 3.0*h_s*h_s;
  //     des_hn = esn*kh*dh_s*1.e-6;
  //  drying
	 q5 = 0.35/fc;
     cd = -q5*esn*kh*dh_s;
   }
   if ( type_temp ==1 ){
  //  temperature constant
   }
   else {
  //  temperature
     dtemp = temp_s-temp_slast;   
     esn=alfa*1.e6;
  //  shrink
  //   des_hn=des_hn+esn*dtemp*1.e-6;
  //  drying
	 q5= 1.5/fc;
     cd = cd+q5*esn*dtemp;
   }
  }

  else {
    q1 = 0.0;
    c0 = 0.0;
    cd = 0.0;
 //   des_hn=0.0;
  }
  //  measure static modulus - creep modulus  1/Ec=2/(3 Es)
  //  1/E=e-6/psi =e-6/6.895kPa
  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)));
}