scaldamtime.cpp

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

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

#define nijac 20
#define limit 1.0e-8


/**
  This constructor inializes attributes to zero values.
  
  TKo
*/
scaldamtime::scaldamtime (void)
{
  indam=0.0;  st=0.0;  uf = 0.0; c = 0.0; kf = 0.0; cde = corr_on;
}

/**
  This destructor is only for the formal purposes.
  
  TKo
*/
scaldamtime::~scaldamtime (void)
{

}

/**
  This function reads material parameters from the opened text file given
  by the parameter in.

  @param in - pointer to the opened input text file
  
  TKo
*/
void scaldamtime::read (FILE *in)
{
  fscanf (in,"%lf %lf %lf %d %d",&indam,&st, &uf, (int *)(&ftype), (int *)(&cde));
  if (ftype == norenergy)
    fscanf(in, "%lf", &c);
  if (ftype == vonmises)
    fscanf(in, "%lf", &kf);

  sra.read (in);
  
}

/**
  This function computes parameters for the damage function
  Different types of the parameter computing are given by the
  attribute ftype.

  @param  ipp - integration point number
  @param  eps - %vector of the strains
  @param  kappa - %vector where the computed parameters are stored
  
  TKo
*/
void scaldamtime::damfuncpar(long ipp, vector &eps, vector &kappa)
{
  double epseq;
  vector poseps;
  vector princsig;
  vector princeps;
  vector posprincsig;
  vector tmp;
  vector sig;
  matrix pvect;
  matrix sigt;
  matrix epst;
  matrix dev;
  matrix d;
  double i1e, j2e, a, b, e, nu, temp; 
  long ncomp, idem = Mm->ip[ipp].gemid();
  long im;

  switch (ftype)
  {
    case norstrain :
      // strain norm
      scprd(eps, eps, epseq);
      kappa[0] = sqrt(epseq);
      break;
    case norenergy :
      // energy norm
      ncomp=Mm->ip[ipp].ncompstr;
      allocm(ncomp, ncomp, d);
      elmatstiff (d,ipp);
      vxm(eps, d, tmp);
      scprd(tmp, eps, epseq);
      kappa[0] = sqrt(epseq);
      break;
    case norposstrain :
      // positive strain norm
      allocv(eps.n,poseps);
      extractposv (eps, poseps);
      scprd(poseps, poseps, epseq);
      kappa[0] = sqrt(epseq);
      break;
    case norposenergy :
      // positive energy norm
      ncomp=Mm->ip[ipp].ncompstr;
      allocm(ncomp, ncomp, d);
      allocv(ncomp, tmp);
      elmatstiff (d,ipp);
      allocv(eps.n,poseps);
      extractposv (eps, poseps);
      vxm(poseps, d, tmp);
      scprd(tmp, poseps, epseq);
      kappa[0] = sqrt(epseq)*c;
      break;
    case norrankine :
    // Rankine norm
    {
      e = Mm->give_actual_ym(ipp);
      ncomp=Mm->ip[ipp].ncompstr;
      allocv(ncomp, sig);
      allocm(ncomp, ncomp, d);
      allocm(3, 3, sigt);
      allocm(3, 3, pvect);
      allocv(3, princsig);
      elmatstiff (d,ipp);
      mxv(d, eps, sig);
      vector_tensor(sig, sigt, Mm->ip[ipp].ssst, stress);
      princ_val (sigt,princsig,pvect,nijac,limit,Mp->zero,3);
      if (Mm->ip[ipp].tm[idem] != elisomat)
      {
        fprintf(stderr, "\n\n Invalid type of elastic material is required");
        fprintf(stderr, "\n  in function scaldamtime::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
      }
      // principal values are sorted -> maximal value of principal stresses is at the last
      // position of vector princsig
      if (princsig[2] > 0.0)
        kappa[0] = princsig[2] / e;
      else
        kappa[0] = 0.0;
      break;
    }
    case norrankinesmooth :
    // smoothed Rankine norm
    {
      e = Mm->give_actual_ym(ipp);
      ncomp=Mm->ip[ipp].ncompstr;
      allocv(ncomp, sig);
      allocm(ncomp, ncomp, d);
      allocm(3, 3, sigt);
      allocm(3, 3, pvect);
      allocv(3, princsig);
      allocv(3, posprincsig);
      elmatstiff (d,ipp);
      mxv(d, eps, sig);
      vector_tensor(sig, sigt, Mm->ip[ipp].ssst, stress);
      princ_val (sigt,princsig,pvect,nijac,limit,Mp->zero,3);
      extractposv (princsig, posprincsig);
      scprd(posprincsig, posprincsig, epseq);
      if (Mm->ip[ipp].tm[idem] != elisomat)
      {
        fprintf(stderr, "\n\n Invalid type of elastic material is required");
        fprintf(stderr, "\n  in function scaldamtime::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
      }
      kappa[0] = sqrt(epseq) / e;
      break;
    }
    case normazar :
    // Mazar's norm
    {
      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, epseq);
      kappa[0] = sqrt(epseq);
      break;
    }
    case vonmises :
      // glasgow modified von Mises norm
      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 scaldamtime::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
      }
      im = Mm->ip[ipp].idm[idem];
      e = Mm->give_actual_ym(ipp);
      nu = Mm->eliso[im].nu;
      a   = (kf-1.0)/(2.0*kf)/(1.0-2.0*nu);
      b   = 3.0/kf/(1+nu)/(1.0+nu);
      temp = a*a*i1e*i1e + b*j2e;
      kappa[0] = a*i1e+sqrt(temp);
/*
      // original von Mises norm
      kappa[0]  = (k-1.0)/(2.0*k)/(1.0-2.0*nu)*i1e;
      tmp = (k-1.0)*(k-1.0)/(1.0-2.0*nu)/(1.0-2.0*nu)*i1e*i1e -12.0*k/(1.0+nu)/(1.0+nu)*j2e;
      kappa[0] += 1.0/(2.0*k)*sqrt(tmp);*/
      break;
    default :
    {
      fprintf(stderr, "\n\n Unknown damage parameter function type");
      fprintf(stderr, "\n  in function scaldamtime::damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
    }
  }
  return;
}

/**
  Function computes derivatives of equivalent strain with respect to
  strain.  Different types of the parameter computing are given by the
  attribute ftype. Function returns resulting values in vector form.

  @param  ipp - integration point number
  @param  eps - %vector of the strains
  @param  kappa - %vector where the computed parameters are stored
  @param  depseq - %vector of resulting derivatives
  
  Returns :
    Resulting values in parameter depseq

  TKo
*/
void scaldamtime::der_damfuncpar(long ipp, vector &eps, vector &kappa, vector &depseq)
{
  vector princeps;
  vector sig;
  vector poseps;
  matrix pvect;
  matrix dev;
  matrix sigt;
  matrix epst;
  matrix d;
  long ncomp, idem = Mm->ip[ipp].gemid();
  matrix depseqt(3,3);
  long i, j, k;
  double i1e, j2e, a, b, e, nu, temp; 
  long im;

  switch (ftype)
  {
    case normazar :
    // Mazar's norm
    {
      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);
      for (i=0; i<3; i++)
      {
        for(j=0; j<3; j++)
	{
          depseqt[i][j] = 0.0; 
          for(k=0; k<3; k++)
	  {
            if (princeps[k] <= 0.0)
              continue;
            temp = pvect[i][k]*pvect[j][k];
            temp *= princeps[k];
            temp /= kappa[0];
            depseqt[i][j] += temp; 
	  }
	}  
      }
      tensor_vector(depseq, depseqt, Mm->ip[ipp].ssst, strain);
      break;
    }
    case vonmises :
      // glasgow modified von Mises norm
      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 scaldamtime::der_damfuncpar, (file %s, line %d)\n", __FILE__, __LINE__);
      }
      im = Mm->ip[ipp].idm[idem];
      e = Mm->give_actual_ym(ipp);
      nu = Mm->eliso[im].nu;
      a   = (kf-1.0)/(2.0*kf)/(1.0-2.0*nu);
      b   = 3.0/kf/(1+nu)/(1.0+nu);
      temp = a*a*i1e*i1e + b*j2e;
      kappa[0] = a*i1e+sqrt(temp);
      copym(dev, depseqt);
      cmulm(b, depseqt);
      depseqt[0][0] = 0.0;
      for (i=0; i<3; i++)
        depseqt[i][i] += i1e*a*a*2.0; 
      temp = a*a*i1e*i1e + b*j2e;
      temp = 1.0/(2.0*sqrt(temp));
      cmulm(temp, depseqt);
      for (i=0; i<3; i++)
        depseqt[i][i] += a; 
      tensor_vector(depseq, depseqt, Mm->ip[ipp].ssst, strain);
      break;
    default :
    {
      fprintf(stderr, "\n\n Unkonwn damage parameter function type");
      fprintf(stderr, "\n  in function scaldamtime::der_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 kappa - %vector of the parameters of damage function
  
  TKo
*/
double scaldamtime::damfunction(long ipp, vector &kappa)
{
  double omega = 0.0, e, ft, err, tmp, dtmp, h, uf_a;
  long idem = Mm->ip[ipp].gemid(), i, ni, eid;
  
  ni=sra.give_ni ();
  err=sra.give_err ();
  // Young modulus
  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 scaldamtime::damfunction (file %s, line %d)\n", __FILE__, __LINE__);
    return 0.0;
  }
  e = Mm->give_actual_ym(ipp);
  // tensile strength
  ft = Mm->give_actual_ft(ipp);
  // initial equivalent strain treshold
  indam = ft/e;
  
  if (kappa[0] < indam)
  // elastic loading
    return 0.0;

  // average size of element
  eid = Mm->elip[ipp];
  switch (Mt->give_dimension(eid))
  {
    case 1 :
      h = Mt->give_length(eid);
      break;
    case 2 :
      h = sqrt(Mt->give_area(eid));
      break;
    case 3 :
      h = pow(Mt->give_volume(eid), 1.0/3.0);
      break;
    default :
      fprintf(stderr, "\n\nError determining dimension of element");
      fprintf(stderr, "\n function damfunction, (file %s, line %d)\n", __FILE__, __LINE__);
  }

  // actual value of initial crack opening
  uf_a = uf/(ft/st);

  // no correction of disipated energy only simple scalar damge
  if (cde == corr_off)
  {
    uf_a /= h;
    uf_a += ft/e;
    omega = 1.0-(ft/e/kappa[0])*exp(-(kappa[0]-ft/e)/(uf_a-ft/e));
    return omega;
  }

  // correction of dispated energy - mesh independent scalar damage
  for (i = 0; i < ni; i++)
  {
    dtmp = -e*kappa[0]+ft/uf_a*h*kappa[0]*exp(-h*omega*kappa[0]/uf_a);
    tmp = (1.0-omega)*e*kappa[0]-ft*exp(-h*omega*kappa[0]/uf_a);
    if (fabs(tmp) < err)
    {
      omega += - tmp/dtmp;
      break;
    }
    omega += - tmp/dtmp;
  }
  if (fabs(tmp) > err)
  {
    fprintf(stderr, "\n\nError - iteration of omega doesn't converge with required error");
    fprintf(stderr, "\n in function damfunction(), (file %s, line %d)\n", __FILE__, __LINE__);
  }
  if ((omega > 1.0) || (omega < 0.0))
  {
    fprintf(stderr, "\n\nError - iteration of omega doesn't converge to range <0,1>");
    fprintf(stderr, "\n in function damfunction(), (file %s, line %d)\n", __FILE__, __LINE__);
    fprintf(stderr, "\n omega = %e",omega);
    fprintf(stderr, "\n h = %e",h);
    fprintf(stderr, "\n uf_a = %e",uf_a);
    fprintf(stderr, "\n kappa[0] = %e",kappa[0]);
    fprintf(stderr, "\n st = %e",st);
    fprintf(stderr, "\n ft = %e",ft);
    fprintf(stderr, "\n e = %e",e);
    fprintf(stderr, "\n tmp = %e",tmp);
    fprintf(stderr, "\n dtmp = %e",dtmp);
  }

  return omega;

}

/**
  This function computes derivative of damage function with respect to equivalent strain.

  @param ipp - integration point number