creepb.cpp

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

#include "creepb.h"
#include "global.h"
#include "math.h"
#include "globmat.h"
#include "genfile.h"
#include "adaptivity.h"
#include "intpoints.h"
#include "stdlib.h"
#include "elastisomat.h"
#include "therisomat.h"

creepb::creepb (void)
{
  timeMax=Mp->timecon.endtime ();
  double m = log10(1.2*timeMax), mm; long i;
  if ((timeMax >= 0.0) && (timeMax < 110.0))
  {
    nRetTime=5;  tncomp=40;
  }
  if ((timeMax >= 110.0) && (timeMax < 1100.0))
  {
    nRetTime=7;  tncomp=40;
  }
  if ((timeMax >= 1100.0) && (timeMax < 14000.0))
  {
    nRetTime=8;  tncomp=50;
  }
  allocv (nRetTime,retTime);
  allocv (nRetTime,ert);
  retTime[0]=1.0e-9;
  //retTime[1]=1.0e-3;
  retTime[1]=1.0;
                  //  retTime[3]=14.0;  //retTime[4]=100.0; //retTime[5]=800.0;
  mm=1./(nRetTime-2);
  for (i=2;i<nRetTime-1;i++){
   retTime[i]=pow(10.,(i-1)*mm*m);
  }
  retTime[nRetTime -1]=1.2*timeMax;
  type_h=1;
  type_temp=1;
  h_s = 0.4;   // end humidity
  desht=0.0;
  k_s = 1.0;	   //shape factor
  k_d = 1.0;      // thicknes
  timemat=0.0;
  alfa=0.0;
  allocm(nRetTime, nRetTime, apom); 
  ccTime = -1.0e30;
  
  et=0.0;
}
creepb::~creepb (void)
{
}

void creepb::creepinit (long ipp,double val,nonmechquant nmq)
{
  long nc=Mm->ip[ipp].ncompstr;
  if(nmq == moist )Mm->ip[ipp].eqother[nc*(nRetTime+1)+1]=val;
  if(nmq == temperature )Mm->ip[ipp].eqother[nc*(nRetTime+1)+2]=val;
}

void creepb::read (XFILE *in)
{
  xfscanf (in,"%lf %lf %lf %lf %lf %lf %lf %lf %lf %ld %ld",
	 &tb,&t_w,&fc,&wc,&sc,&gc,&c_s,&a1,&k_d,&type_h,&type_temp);
  
  if (type_h==1){xfscanf (in,"%lf ",&h_s); }
  if (type_temp==1){xfscanf (in,"%lf ",&temp_s); }
  

}

long creepb::numberOfCreepb ()
{
  return nRetTime;
}

/**
   function approximates function defined by nodal values

   @param areacoord - vector containing area coordinates
   @param nodval - nodal values

*/
double creepb::approx (vector &areacoord,vector &nodval)
{
  double f;
  scprd (areacoord,nodval,f);
  return f;
}

/**
   function inversion sym. matrix

   @param a - matrix

*/
void creepb::inv_sym (matrix &a)
{
  long i,j,k;
  double diag,an;
  for (k=0;k<a.n;k++){
    diag = a[k][k];
    if (diag<=0.0){
      fprintf (stderr,"\n Creep model: zero diagonal entry in line %ld",k);
      abort();
    }
    for (i=0;i<a.n;i++){
      an = a[i][k]/diag;
      if (i!=k) { 
	for (j=i;j<a.n;j++){
	  if(j!=k) {
	    a[i][j] = a[i][j] - an*a[k][j];
	    a[j][i] = a[i][j];
	  }
	}
      }
      a[i][k] = an;
      a[k][i] = an;
    }
    a[k][k] = - 1.0/diag;
  }
  
  for (i=0;i<a.n;i++){
    for (j=0;j<a.n;j++){
      a[i][j] = - a[i][j];
    }
  }
}

void creepb::updateval()
{
 if(type_h!=1&&type_temp!=1) updatevalchange();
 else updatevalkons();
 
}
void creepb::updatevalchange()
{
  long i,j,k;
  double pomt,pomt1,qq,delYi,delYj,m0,m;
  
  if (ccTime != Mp->timecon.actualtime ())
  {
    //ddTime=Mp->timefun.getval(Mp->time);
    ddTime=Mp->timecon.actualforwtimeincr ();
    //ccTime=Mp->time;
    ccTime=Mp->timecon.actualtime ();
    pomt=ccTime+ddTime/2.;
    qq=2./3.;
    //   from concrete starts   tb 
    fillm(0.0, apom); 
    if (timemat<=pomt) {
      pomt1=pomt+0.0001;   
      m0 = log10(pomt1);
      m = log10(2.0*timeMax)*0.999;
      k=0;
      while (pomt1<=2.0*timeMax)
      {
        for (i=0;i<nRetTime;i++)
        {  
	  delYi=pow((pomt/retTime[i]),qq)-pow((pomt1/retTime[i]),qq);
	  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)/tncomp*k);
      }
      inv_sym(apom);
    }
    //fprintf(stdout, "\n ### inversion\n");
  }
}
void creepb::updatevalkons()
{
  long i,j,k;
  double pomt,pomt1,qq,jt,delYi,delYj,m0,m;
  vector bpom(nRetTime);
  
  if (ccTime != Mp->timecon.actualtime ())
  {
    //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.0*timeMax)*0.999;
    k=0;
    while (pomt1<=2.0*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]),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)/tncomp*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],qq)-pow(0.0001/retTime[i],qq);}
      else {
	    delYi=pow((ccTime+ddTime)/retTime[i],qq)-pow(ccTime/retTime[i],qq);}
        et=et+(delYi-1.+exp(-delYi))/delYi/ert[i];
    }
    if(et>1.e-55)  et=1./et;
    else{fprintf (Out,"\n\n Stiff is zero"); abort();}
    //         call Jkontrola(imat)
    
   }
  
   else {
//    et=e0/100000.;
    et=1.;
    for (i=0;i<nRetTime;i++){ ert[i]=0.;}
   }
	
  }
  
}


/**
   function returns eps from history
   
   @param screep - vector of history
   @param epsscr - vector deformation of history
   
   10.10.2002
*/
void creepb::nlstresses (long ipp, long im, long ido)
{
  //  number of components od strain/stress arrays
  nc=Mm->ip[ipp].ncompstr;
  
  //  type of stress state
  //bar=1,plbeam=2,spacebeam=5,
  //planestress=10,planestrain=11,plate=15,
  //axisymm=20,shell=25,spacestress=30
  ss=Mm->ip[ipp].ssst;  
  
  //timeMax=Mp->end_time;
  //ddTime=Mp->timefun.getval(Mp->time);
  ddTime=Mp->timecon.actualforwtimeincr ();
  //ccTime=Mp->time;
  ccTime=Mp->timecon.actualtime ();
  
  //  id of elastic material
  imat=Mm->ip[ipp].idm[Mm->ip[ipp].gemid()];
  
  //  actual moisture
  if(type_h!=1){
    h_s=Mm->moist[ipp];
  }
  //if (ipp==0)
  //fprintf (Out,"\n nlstresses   h_s %e",h_s);
  
  //  actual temperature
  if(type_temp!=1){
    temp_s=Mm->tempr[ipp];
  }
  //if (ipp==0)
  //fprintf (Out,"\n nlstresses    temp_s %e",temp_s);
  
  get_h (ipp);     //initial  moist
  get_temp (ipp);     //initial  tempr
  
  if (Mp->phase==1){
    phase1(ipp,im,ido);
  }
  
  if (Mp->phase==2){
    phase2(ipp);
  }
  
}

void creepb::phase1 (long ipp,long im, long ido)
{
  //  right hand side from time dependent value computation
  //  creep function for shrink from pomt to cctime
  long i,j,ii;
  double qq,dj,endTime;
  
  vector epscr(nc),sig(nc);
  matrix d(nc,nc),screep(nc,nRetTime);
  
  //    Mm->matstiff(d,ipp);  //compute desht
  //    e0=Mm->eliso[imat].e;
  //    mi=Mm->eliso[imat].nu;    //    d[1][0] = nu*e/(1.0-nu*nu);
  endTime=ccTime+ddTime;
  get_desht (desht, ccTime, endTime);
  for (i=0; i<nc; i++){
    epscr[i]=0.0;
  }
  qq=2.0/3.0;
  for (i=0;i<nRetTime;i++){
    ii=nc*(i+1);
    dj=pow(endTime/retTime[i],qq)-pow(ccTime/retTime[i],qq);
    for (j=0; j<nc; j++){
      screep[j][i]=Mm->ip[ipp].eqother[ii+j];	
      epscr[j]=epscr[j]+(1.-exp(-dj))*screep[j][i];
    }
  }
  //        fprintf (Out,"\n\n v case t  prirustek dt faze 1");
  //        fprintf (Out," %e\t%e", ccTime, ddTime);
  //        fprintf (Out,"\n\n eps faze1");
  //        fprintf (Out," %le %le %le %le %le %le",epscr[0],epscr[1],epscr[2],epscr[3],epscr[4],epscr[5]);
  //        fprintf (Out,"\n\n tau faze1");
  //        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]);
  //  increment of shrink and temperature
  
  //  jen kvuli debuggovani
  //desht=0.0;

  epscr[0]=epscr[0]+desht;
  epscr[1]=epscr[1]+desht;
  if(ss==planestrain) epscr[2]=epscr[2]+desht;
  else if(ss==axisymm) epscr[2]=epscr[2]+desht;
  else if(ss==spacestress) epscr[2]=epscr[2]+desht;

  //	fprintf (Out,"\n\n desht prir faze 1"); fprintf (Out," %le",desht);
  //	fprintf (Out,"\n h_s prir faze 1"); fprintf (Out," %le %le",h_s, h_slast);
  //	fprintf (Out,"\n temp_s prir faze 1"); fprintf (Out," %le %le",temp_s, temp_slast);
  
  //		fprintf (Out,"\n\n phase 1 eps1, eps2 prir");
  //		fprintf (Out," %e %e",epscr[0], epscr[1]);
  
  //fprintf (Out,"\n bod %3ld  epscr %le",ipp,epscr[0],epscr[1],epscr[2]);
  Mm->stiff_deps_creep (ipp, im, epscr, nc,nc);  // save sig
}



void creepb::phase2 (long ipp)
{
//  new total strain and stress
  long i,j,ii;
  double qq,dj, endTime;
  vector epscr(nc),deps(nc),depstot(nc),sig(nc);
  matrix d(nc,nc),screep(nc,nRetTime);
  
  endTime=ccTime+ddTime;
   // fprintf (Out,"\n\n phase 2 v case t   prirustek dt");
   // fprintf (Out," %e\t%e", ccTime, ddTime);
  // arr. strain = is total strain (including e from screep) from t=0
  for (i=0; i<nc; i++){
    depstot[i]=Mm->ip[ipp].strain[i]-Mm->ip[ipp].eqother[i];
  } //increment deps
  qq=2.0/3.0;
  for (i=0;i<nRetTime;i++){
    ii=nc*(i+1);
    dj=pow(endTime/retTime[i],qq)-pow(ccTime/retTime[i],qq);
    for (j=0; j<nc; j++){
      screep[j][i]=Mm->ip[ipp].eqother[ii+j];
      deps[j]=depstot[j]-(1.-exp(-dj))*screep[j][i];     // deps=deps-e(screep)
    }
  }
  get_desht (desht, ccTime, endTime);
  //  fprintf (Out,"\n\n desht prir faze 1"); fprintf (Out," %le",desht);
  
  //  debuggovani
  //desht=0.0;

  deps[0]=deps[0]-desht;
  deps[1]=deps[1]-desht;
  if(ss==planestrain) deps[2]=deps[2]-desht;
  else if(ss==axisymm) deps[2]=deps[2]-desht;
  else if(ss==spacestress) deps[2]=deps[2]-desht;
  //  increment of stress components
  //        fprintf (Out,"\n\n eps celkove,   eps celkove minule,  delta eps");
  //        fprintf (Out,"\n %le %le %le %le %le %le",epscr[0],epscr[1],epscr[2],epscr[3],epscr[4],epscr[5]);
  //        fprintf (Out,"\n %le %le %le ",Mm->ip[ipp].eqother[0],Mm->ip[ipp].eqother[1],Mm->ip[ipp].eqother[2]);
  //        fprintf (Out,"\n %le %le %le ",Mm->ip[ipp].eqother[3],Mm->ip[ipp].eqother[4],Mm->ip[ipp].eqother[5]);
  //        fprintf (Out,"\n\n delta eps");
  //        fprintf (Out,"\n %e %e %e",deps[0],deps[1],deps[2]);
  //  stiffness matrix of material
  Mm->matstiff(d,ipp);  
  e0=Mm->eliso[imat].e;
  mi=Mm->eliso[imat].nu;    //    d[1][0] = nu*e/(1.0-nu*nu);
  mxv (d,deps,sig);

  // save total strain (from t=0)
  for (i=0; i<nc; i++){Mm->ip[ipp].eqother[i]=Mm->ip[ipp].strain[i];}// it is last strain in arr. strain
  //  time history
  seps_time (screep,sig);
  for (i=0; i<nRetTime; i++){
    ii=nc*(i+1);
    for (j=0; j<nc; j++){
      Mm->ip[ipp].eqother[ii+j]=screep[j][i];
    }