glasgmech.cpp

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

}

/**
   function computes derivative of thermal damage function with respect to temperature
   kappa_td is notation in the notes
   
   @param ipp - integration point number
   @param tempr - actual temperature
   @param kappa - %vector of the parameters of thermal damage function
                  it contains the maximum of either the largest value attained by temperature
		  or the reference temperature

*/
double glasgmech::dtdfdt (long ipp,double tempr,vector &kappa)
{
  double htd;
  
  if (tempr > kappa[0])
    kappa[0] = tempr;
  
  if (kappa[0] > tkappa0)
    htd = 20.0*(1.0-5.0*(kappa[0]-tkappa0)) - 5.0*20.0*(kappa[0] - tkappa0);
  else
    htd = 0.0;
  
  return htd;
}

/**
   function computes increments of load induced thermal strains
   @param ipp - integration point number
   @param epstm - result vector of increments of load induced thermal strains
   @param sigt - stress tensor
   @param told - previous temperature
   @param tnew - actual temperature
   
*/
void glasgmech::compute_lits (long ipp,vector &epstm,matrix &sigt,double told,double tnew)
{
  long idem;
  double dt,i1,beta,nu,fc;
  vector princsig(3);
  matrix pvect(3,3);
  strastrestate ssst;
  
  idem=Mm->ip[ipp].gemid();
  if (Mm->ip[ipp].tm[idem] != elisomat)
    {
      fprintf(stderr, "\n\n Invalid type of elastic material is required");
      fprintf(stderr, "\n  in function glasgmech::compute_lits, (file %s, line %d)\n", __FILE__, __LINE__);
    }
  nu = Mm->eliso[idem].nu;
  
  fillv (0.0,epstm);
  
  //  increment of temperature
  dt=tnew-told;
  if (dt<=0.0){
    return;
  }
  else{
    princ_val (sigt,princsig,pvect,nijac,limit,Mp->zero,3);
    extractnegv (princsig,princsig);
    fillm (0.0,sigt);
    sigt[0][0]=princsig[0];
    sigt[1][1]=princsig[1];
    sigt[2][2]=princsig[2];
    lgmatrixtransf(sigt, pvect);
    
    i1=first_invar (sigt);

    beta = betacoeff (ipp,tnew);
    
    fc = st*k;
    
    cmulm ((1.0+nu),sigt);
    sigt[0][0]-=nu*i1;
    sigt[1][1]-=nu*i1;
    sigt[2][2]-=nu*i1;
    cmulm (beta/fc*dt,sigt);
    
    ssst=Mm->ip[ipp].ssst;
    tensor_vector (epstm,sigt,ssst,strain);
  }
}

/**
   Function computes derivatives of equivalent strain with respect of elastic strains.
   @param ipp - integration point number
   @param epstm - result vector of increments of load induced thermal strains
   @param sigt - stress tensor
   @param told - previous temperature
   @param tnew - actual temperature
   
*/
void glasgmech::depseqdepsel(long ipp, vector &eps, vector &deeqdeel)
{
  double q1, q2, q3, e, nu, i1e, j2e, tmp;
  double r, af, bf;
  long i;
  long idem = Mm->ip[ipp].gemid();
  long ncomp = Mm->ip[ipp].ncompstr;
  vector pv(ncomp), tv(ncomp);
  matrix pm(ncomp,ncomp);
  matrix epst(3,3), dev(3,3);
  
  fillv(0.0, pv); 
  fillm(0.0, pm);
  fillv(0.0,deeqdeel);

  vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
  i1e = first_invar(epst);
  deviator (epst,dev);
  j2e = second_invar (dev);
  idem = Mm->ip[ipp].gemid();
  if (Mm->ip[ipp].tm[idem] != elisomat)
  {
    fprintf(stderr, "\n\n Invalid type of elastic material is required");
    fprintf(stderr, "\n  in function glasgmech::depseqds, (file %s, line %d)\n", __FILE__, __LINE__);
  }
  e  = Mm->eliso[idem].e;
  nu = Mm->eliso[idem].nu;
  af  = (k-1.0)/(2.0*k)/(1.0-2.0*nu);
  bf  = 3.0/k/(1+nu);
  r   = af*af*i1e*i1e + bf*j2e;
  r   = sqrt(r);
  if (r > 0.0)
    r = 0.5 / r;
  else     
    r = 0.0;
  q3 = nu/(1.0-nu);
  q1 = 1.0 + q3 + q3*q3;
  q2 = 1.0 - 2.0*q3 - 2.0*q3*q3;
  switch (Mm->ip[ipp].ssst)
  {
    case planestrain:
    case axisymm:
      pv[0] = 1.0-q3;
      pv[1] = 1.0-q3;
      pv[2] = 0.0;
      pv[3] = 0.0;
      pm[0][0] = 2.0/3.0*q1;
      pm[1][1] = pm[0][0];
      pm[2][2] = 0.0;
      pm[3][3] = 2.0;
      pm[0][1] = -q2/3.0;
      pm[1][0] = pm[0][1];
      break;
    case spacestress:
      for (i = 0; i < 3; i++)
      {
        pv[i] = 1.0;
        pm[i][i] = 2.0/3.0;
        pm[i+3][i+3] = 2.0;
      }
      pm[0][1] = -q2/3.0; 
      pm[0][2] = pm[0][1]; 
      pm[1][2] = pm[0][1]; 
      pm[2][0] = pm[0][1]; 
      pm[2][1] = pm[0][1]; 
      break;
    default:
      fprintf(stderr, "\n\n Invalid stress/strain state is required in function glasgmech::depseqds\n");
      fprintf(stderr, " in integeration point %ld, (file %s, line %d)\n", ipp, __FILE__, __LINE__);           
  }
  for(i = ncomp / 2; i < ncomp; i++)
    eps[i] *= 0.5;
  tmp = af+2.0*af*af*i1e*r;
  cmulv(tmp, pv);
  mxv(pm, pv, tv);
  cmulv(bf*r, tv);
  addv(tv, pv, deeqdeel);
  for(i = ncomp / 2; i < ncomp; i++)
  {
    deeqdeel[i] *= 0.5;
    eps[i] *= 2.0;
  }
}                



/**
   function computes correct stresses in the integration point and stores
   them into ip stress array.
   
   @param ipp - integration point pointer
   @param im - index of material
   @param ido - index of the viscous material in the array eqother
   
*/
void glasgmech::nlstresses (long ipp,long im,long ido)
{
  long i,ncompstr;
  
  //  number of strain/stress components in the problem
  ncompstr = Mm->ip[ipp].ncompstr;
  
  double newt,oldt,dt,chi,htd,omega,domega, omegao, kappao, depseq;
  
  vector epsn(ncompstr),epso(ncompstr),epse(ncompstr),epsft(ncompstr),epstm(ncompstr);
  vector epsc(ncompstr),epsirr(ncompstr),epstirr(ncompstr), deeqdeel(ncompstr), deps(ncompstr);
  vector sig(ncompstr),sigtrial(ncompstr);
  vector kappa(2),tkappa(1);
  vector deftdt(ncompstr), s(ncompstr); 
//  vector detmdt(ncompstr);
  matrix eps(3,3),d(ncompstr,ncompstr);
  matrix sigt(3,3);
  
  if (Mp->phase==1){
    /* *******************************/
    //  right hand side computation  //
    /* *******************************/
    
    //  total new strains
    for (i=0;i<ncompstr;i++){
      epsn[i] = Mm->ip[ipp].strain[i];
    }
    //  previous total strains
    for (i=0;i<ncompstr;i++){
      epso[i] = Mm->ip[ipp].eqother[ncompstr+i];
    }
    

    /* ************************/
    //  free thermal strains  //
    /* ************************/
    
    //  actual temperature
    newt = Mm->tempr[ipp];
    //  temperature from the previous step
    oldt = Mm->ip[ipp].eqother[4*ncompstr];
    //  increment of temperature
    dt = newt-oldt;
    
    //  increments of free thermal strains
    compute_thermdilat (newt,dt,eps);
    //  conversion of tensor notation into vector one
    tensor_vector (epsft,eps,Mm->ip[ipp].ssst,strain);
    
    /* ********************************/
    //  load induced thermal strains  //
    /* ********************************/
    vector_tensor(sig, sigt, Mm->ip[ipp].ssst, stress);
    compute_lits (ipp, epstm, sigt, oldt, newt);
   //compute_lits (epstm);
    
    /* ******************/
    //  thermal damage  //
    /* ******************/
    
    //  history parameter of thermal damage
    tkappa[0]=Mm->ip[ipp].eqother[4*ncompstr+3];
    
    //  thermal damage parameter
    chi = thermdamfunction (ipp,newt,tkappa);
    
    /* *********************/
    //  mechanical damage  //
    /* *********************/
    
    //  norm of equivalent strain
    damfuncpar(ipp,epsn,kappa);
    
    // previous equivalent strain norm
    kappao = Mm->ip[ipp].eqother[4*ncompstr+1];
    //  previous damage
    omegao = Mm->ip[ipp].eqother[4*ncompstr+2];
    // previous temperature
    kappa[1] = Mm->ip[ipp].eqother[4*ncompstr];
    //  mechanical damage parameter
    omega = damfunction(ipp,newt,chi,kappa);
    if (omega < kappao)
    {
      omega = omegao;
      domega = 0.0;
    }
    else
    {
      // depseq = kappa[0] - kappao;

      // increment of equvalent strain
      depseqdepsel(ipp, epsn, deeqdeel);
      subv (epsn, epso, deps);
      subv (deps, epsft, deps);
      subv (deps, epstm, deps);
      scprd(deps, deeqdeel, depseq);
      // mechanical damage increment
      domega = domegadt(ipp, newt, chi, kappa)*depseq;
      domega += domegadkmd(ipp, newt, chi, kappa)*(newt-oldt);
    }
    
    
    /* *********/
    //  creep  //
    /* *********/
    
    //give_creep_eps (epsc);
    
    /* *************************************/
    //  increment of irreversible strains  //
    /* *************************************/
    for (i=0;i<ncompstr;i++){
      epsirr[i]=epsft[i]+epsc[i]+epstm[i];
    }
    
    /* ******************************/
    //  storage of computed values  //
    /* ******************************/
    for (i=0;i<ncompstr;i++){
      Mm->ip[ipp].eqother[3*ncompstr+i]=epsirr[i];
    }
    Mm->ip[ipp].eqother[4*ncompstr+1]=kappa[0];
    Mm->ip[ipp].eqother[4*ncompstr+2]=omega;
    Mm->ip[ipp].eqother[4*ncompstr+3]=tkappa[0];
    Mm->ip[ipp].eqother[4*ncompstr+4]=chi;
    Mm->ip[ipp].eqother[4*ncompstr+5]=domega;
    
    //  stiffness matrix of the material
    Mm->elmatstiff(d,ipp);
    
    //  stress increments
    mxv (d,epsirr,sig);
    
    Mm->storestress (0,ipp,sig);
  }
  
  
  if (Mp->phase==2){
    
    for (i=0;i<ncompstr;i++){
      //  new total strain
      epsn[i] = Mm->ip[ipp].strain[i];
      //  previous total strain
      epso[i] = Mm->ip[ipp].eqother[ido+ncompstr+i];
      //  previous total irreversible strains
      epstirr[i] = Mm->ip[ipp].eqother[ido+2*ncompstr+i];
      //  increment of irreversible strain
      epsirr[i] = Mm->ip[ipp].eqother[ido+3*ncompstr+i];
      
      //  total elastic strains
      epse[i]=epsn[i]-epsirr[i];
    }
    omega =Mm->ip[ipp].eqother[4*ncompstr+2];
    chi   =Mm->ip[ipp].eqother[4*ncompstr+4];
    domega=Mm->ip[ipp].eqother[4*ncompstr+5];
    //  actual temperature
    //newt = Mm->tempr[ipp];
    newt = 20.0;
    //  temperature from the previous step
    oldt = Mm->ip[ipp].eqother[4*ncompstr];
    //  increment of temperature
    dt = newt-oldt;
    //  history parameter of thermal damage
    tkappa[0]=Mm->ip[ipp].eqother[4*ncompstr+3];
    //  derivative of thermal damage function with respect to temperature
    htd = dtdfdt (ipp,newt,tkappa);
    
    
    //  stiffness matrix of material
    Mm->matstiff(d,ipp);
    //  effective stress
    mxv (d,epse,sigtrial);

    //  total strain increment
    subv (epsn,epso,epso);
    //  elastic strain increment
    subv (epso,epsirr,epse);
    //  secant stiffness matrix
    cmulm ((1.0-omega)*(1.0-chi),d);
    //  stress increment
    mxv (d,epse,sig);

    cmulv (htd*dt*(1.0-omega)+domega*(1.0-chi),sigtrial);
    
    subv (sig,sigtrial,sig);
    

    //  new data storage
    Mm->ip[ipp].eqother[4*ncompstr+0]=newt;
    for (i=0;i<ncompstr;i++){
      //  new stress
      Mm->ip[ipp].eqother[ido+i] += sig[i];
      //  total strains
      Mm->ip[ipp].eqother[ido+ncompstr+i] = epsn[i];
      //  total irreversible strains
      Mm->ip[ipp].eqother[ido+2*ncompstr+i] += epsirr[i];
      //  stress increment
      Mm->ip[ipp].stress[i]=sig[i];
    }

  }
}