glasgmech.cpp

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

#include <math.h>
#include "glasgmech.h"
#include "vector.h"
#include "matrix.h"
#include "vecttens.h"
#include "tensor.h"
#include "global.h"
#include "intpoints.h"
#include "elastisomat.h"


#define nijac 20
#define limit 1.0e-8

glasgmech::glasgmech (void)
{
  tkappa0=0.0;
  st=0.0;
  k=0.0;
  gf0=0.0;
  lc=0.0;
  a=0.0;
  b=0.0;
  c=0.0;
}

glasgmech::~glasgmech (void)
{

}

/**
   function reads input data from the input text file
*/
void glasgmech::read (XFILE *in)
{
  xfscanf (in,"%lf %lf %lf %lf %d", &tkappa0, &st, &gf0, &lc, (int *)(&ft));
  xfscanf(in, "%lf", &k);
  xfscanf (in,"%lf %lf %lf",&a,&b,&c);
}

/**
   function computes coefficient beta
   
   @param ipp - number of integration point
   @param t - temperature
   
   13.7.2004
*/
double glasgmech::betacoeff (long ipp,double t)
{
  double tstar,tbar,beta;
  
  tstar=(470.0-20.0)/100.0;
  tbar=(t-20.0)/100.0;
  
  if (0.0 > tbar){
    fprintf (stderr,"\n\n Warning: normalized temperature t_bar is less than 0.0 at integration point %ld",ipp);
    return 0.0;
  }
  if (0.0<=tbar && tbar<=tstar){
    beta=(2.0*a*tbar+b)*0.01;
  }
  else{
    beta=(2.0*c*(tbar-tstar)+2.0*a*tstar+b)*0.01;
  }
  
  return beta;
}


/**
   function computes norm of equivalent strains
   
   @param ipp - number of integration point
   @param  eps - %vector of the strains
   @param  kappa - %vector where the computed parameters are stored

*/
void glasgmech::damfuncpar(long ipp, vector &eps, vector &kappa)
{
  vector poseps;
  vector princeps;
  matrix pvect;
  matrix epst;
  matrix dev;
  double i1e, j2e, nu, e, tmp, a, b;
  long ncomp, idem;

  switch (ft)
  {
    case vonmises :
      ncomp=Mm->ip[ipp].ncompstr;
      allocm(3, 3, epst);
      allocm(3, 3, dev);
      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::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
      }
      e  = Mm->eliso[idem].e;
      nu = Mm->eliso[idem].nu;
      a   = (k-1.0)/(2.0*k)/(1.0-2.0*nu);
      b   = 3.0/k/(1+nu);
      tmp = a*a*i1e*i1e + b*j2e;
      kappa[0] = a*i1e+sqrt(tmp);
      break;
    case normazar :
      allocm(3, 3, epst);
      allocm(3, 3, pvect);
      allocv(3, princeps);
      allocv(3,poseps);
      fillv(0.0, princeps);
      vector_tensor(eps, epst, Mm->ip[ipp].ssst, strain);
      princ_val (epst,princeps,pvect,nijac,limit,Mp->zero,3);
      extractposv (princeps, poseps);
      scprd(poseps, poseps, tmp);
      kappa[0] = sqrt(tmp);
      break;
    default :
    {
      fprintf(stderr, "\n\n Unknown damage parameter function type");
      fprintf(stderr, "\n  in function glasgmech::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
    }
  }
  return;
}

/**
  This function computes damage parameter omega which is the result of the
  damage function.

  @param ipp - integration point number
  @param tempr - actual temperature
  @param chi - thermal damage parameter
  @param kappa - %vector of the parameters of damage function
  
*/
double glasgmech::damfunction(long ipp,double tempr,double chi,vector &kappa)
{
  double omega = 0.0, e0, kappa0, tn, gamma, gf,th;
  long idem = Mm->ip[ipp].gemid();
  
  if (Mm->ip[ipp].tm[idem] != elisomat)
  {
    fprintf(stderr, "\n\nError - ip %ld has wrong type of elastic material combined with", ipp);
    fprintf(stderr, "\n scalar damge. The elastic isotropic material is requried");
    fprintf(stderr, "\n Function glasgmech::damfunction (file %s, line %d)\n", __FILE__, __LINE__);
    return 0.0;
  }

  
  // Young modulus
  e0 = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
  
  //  normalized temperature
  tn = (tempr-20.0)/100.0;
  //  maximum reached normalized temperature
  th = (kappa[1]-20.0)/100.0;
  if (tn>th){
    th=tn;
  }
  //  check of maximum reached normalized temperature
  if ((th < 0.0) || (th > 7.9)){
    fprintf(stderr, "\n\nWarning - largest value of normalised temperature\n");
    fprintf(stderr, "          is out of suggested range 0.0-7.9 (ipp %ld)\n",ipp);
  }
  
  //  threshold of the mechanical damage (depends on temperature)
  kappa0 = st/e0*(1.0 - 0.016*th*th)/(1.0 - 0.1*th)/(1.0 - 0.1*th);
  
  //  fracture energy
  gf = gf0*(1.0+0.39*th-0.07*th*th);
  
  //  softening parameter
  gamma = (1.0-chi)*st*lc/gf;
  
  if (kappa[0] < kappa0)
    // elastic loading
    return 0.0;
  
  //  mechanical damage parameter
  omega = 1.0-(kappa0/kappa[0])*exp(-gamma*(kappa[0]-kappa0));

  return omega;
}

/**
  Function computes derivative of mechanical damage parameter with respect to equivalent strain

  @param ipp - integration point number
  @param tempr - actual temperature
  @param chi - thermal damage parameter
  @param kappa - %vector of the parameters of damage function

*/
double glasgmech::domegadkmd (long ipp,double tempr,double chi, vector &kappa)
{
  double domega = 0.0, e0, kappa0, tn, gamma, gf,th;
  long idem = Mm->ip[ipp].gemid();
  
  if (Mm->ip[ipp].tm[idem] != elisomat)
  {
    fprintf(stderr, "\n\nError - ip %ld has wrong type of elastic material combined with", ipp);
    fprintf(stderr, "\n scalar damge. The elastic isotropic material is requried");
    fprintf(stderr, "\n Function glasgmech::damfunction (file %s, line %d)\n", __FILE__, __LINE__);
    return 0.0;
  }

  
  // Young modulus
  e0 = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
  //  normalized temperature
  tn = (tempr-20.0)/100.0;
  //  maximum reached normalized temperature
  th = (kappa[1]-20.0)/100.0;
  if (tn>th){
    th=tn;
  }
  //  check of maximum reached normalized temperature
  if ((th < 0.0) || (th > 7.9)){
    fprintf(stderr, "\n\nWarning - largest value of normalised temperature\n");
    fprintf(stderr, "          is out of suggested range 0.0-7.9 (ipp %ld)\n",ipp);
  }
  
  //  threshold of the mechanical damage (depends on temperature)
  kappa0 = st/e0*(1.0 - 0.016*th*th)/(1.0 - 0.1*th)/(1.0 - 0.1*th);
  
  //  fracture energy
  gf = gf0*(1.0+0.39*th-0.07*th*th);
  
  //  softening parameter
  gamma = (1.0-chi)*st*lc/gf;

  if (kappa[0] < kappa0)
    // elastic loading
    return 0.0;
  
  domega  = -kappa0*exp(-gamma*(kappa[0]-kappa0))/kappa[0];
  domega += kappa0*kappa[0]*gamma*exp(-gamma*(kappa[0]-kappa0));
  domega /= kappa[0]*kappa[0];
  
  return domega;
}


/**
  Function computes derivative of mechanical damage parameter with respect to temperature

  @param ipp - integration point number
  @param tempr - actual temperature
  @param chi - thermal damage parameter
  @param kappa - %vector of the parameters of damage function

*/
double glasgmech::domegadt (long ipp,double tempr,double chi, vector &kappa)
{
  double domega = 0.0, e0, kappa0, tn, gamma, gf,th, dkmd0dt, dgdt;
  long idem = Mm->ip[ipp].gemid();
  
  if (Mm->ip[ipp].tm[idem] != elisomat)
  {
    fprintf(stderr, "\n\nError - ip %ld has wrong type of elastic material combined with", ipp);
    fprintf(stderr, "\n scalar damge. The elastic isotropic material is requried");
    fprintf(stderr, "\n Function glasgmech::damfunction (file %s, line %d)\n", __FILE__, __LINE__);
    return 0.0;
  }

  
  // Young modulus
  e0 = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
  //  normalized temperature
  tn = (tempr-20.0)/100.0;
  //  maximum reached normalized temperature
  th = (kappa[1]-20.0)/100.0;
  if (tn>th){
    th=tn;
  }
  //  check of maximum reached normalized temperature
  if ((th < 0.0) || (th > 7.9)){
    fprintf(stderr, "\n\nWarning - largest value of normalised temperature\n");
    fprintf(stderr, "          is out of suggested range 0.0-7.9 (ipp %ld)\n",ipp);
  }
  
  //  threshold of the mechanical damage (depends on temperature)
  kappa0 = st/e0*(1.0 - 0.016*th*th)/(1.0 - 0.1*th)/(1.0 - 0.1*th);
  
  //  fracture energy
  gf = gf0*(1.0+0.39*th-0.07*th*th);
  
  //  softening parameter
  gamma = (1.0-chi)*st*lc/gf;

  if (kappa[0] < kappa0)
    // elastic loading
    return 0.0;
  
  // derivative of kappa_md with respect of temperature
  dkmd0dt = st/e0/100.0*(-0.032*th*(1.0-0.1*th)*(1.0-0.1*th)+(1.0-0.016*th*th)*2.0*(1.0-0.1*th)*0.1);
  dkmd0dt /= (1.0-0.1*th)*(1.0-0.1*th)*(1.0-0.1*th)*(1.0-0.1*th);
  // derivative of gamma with respect of temperature
  dgdt = -(1-chi)*st*lc/(100.0*gf*gf)*(0.39-0.14*th)*gf0;

  domega  = dkmd0dt*exp(-gamma*(kappa[0]-kappa0))/kappa[0];
  domega += kappa0/kappa[0]*exp(-gamma*(kappa[0]-kappa0))*(-(dgdt*(kappa[0]-kappa0)-gamma*dkmd0dt));
  domega  = -domega;
  
  return domega;
}

/**
  This function computes material stiffnes %matrix.

  @param d - allocated matrix structure for material stiffness %matrix
  @param ipp - integration point number
  @param ido - index of internal variables for given material in the ipp other array

*/
void glasgmech::matstiff (matrix &d,long ipp,long ido)
{
  double dp;

  switch (Mp->stmat)
  {
    case initial_stiff :
      Mm->elmatstiff (d,ipp);
      break;
    case tangent_stiff :
      Mm->elmatstiff (d,ipp);
      dp=Mm->ip[ipp].eqother[ido+1];
      if (dp > 0.999999)
        dp = 0.999999;
      cmulm (1.0-dp,d);
      break;
    default :
      fprintf(stderr, "\n\nError - unknown type of stifness matrix");
      fprintf(stderr, "\n in function glasgmech::matstiff (%s, line %d)\n", __FILE__, __LINE__);
  }
}

/**
   function computes 
   
   @param t - temeprature
   @param dt - incrment of temperature
   @param eps - thermal strain increment
   
*/
void glasgmech::compute_thermdilat (double t,double dt,matrix &eps)
{
  double auxt,alpha;
  
  auxt = (t-20.0)/100.0;
  if (auxt>=0.0 && auxt<=6.0){
    alpha = 6.0e-5/(7.0-auxt);
  }
  else{  alpha=0.0; }
  
  fillm (0.0,eps);
  
  eps[0][0]=alpha*dt;
  eps[1][1]=alpha*dt;
  eps[2][2]=alpha*dt;
  
}

/**
  This function computes thermal damage parameter chi which is the result of the
  thermal damage function.

  @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::thermdamfunction (long ipp,double tempr,vector &kappa)
{
  double chi = 0.0;
  
  if (tempr > kappa[0])
    kappa[0] = tempr;
  if (kappa[0] > tkappa0)
    chi = 20.0*(kappa[0] - tkappa0)*(1.0-5.0*(kappa[0]-tkappa0));
  else
    chi = 0.0;
  
  return chi;