camclay.cpp

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

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "camclay.h"
#include "global.h"
#include "vecttens.h"
#include "intpoints.h"
#include "matrix.h"
#include "vector.h"
#include "tensor.h"
#include "elastisomat.h"

/**
  This constructor inializes attributes to zero values.
*/
camclay::camclay (void)
{
  m = 0.0;
  lambda = 0.0;
  kappa = 0.0;
}

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

}

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

  @param in - pointer to the opned text file

*/
void camclay::read (XFILE *in)
{
  xfscanf (in, "%lf %lf %lf", &m, &lambda, &kappa);
  sra.read (in);
}

/**
   This function computes the value of yield functions.

   @param sig - stress tensor
   @param q   - %vector of hardening parameter

   @retval The function returns value of yield function for the given stress tensor

   25.3.2002
*/
double camclay::yieldfunction (matrix &sig, vector &q)
{
  double i1s,j2s,f;

  matrix dev(3,3);
  deviator (sig,dev);

  i1s = first_invar (sig)/3.0;
  j2s = second_invar (dev);
//  f = (3*j2s) + m*m*i1s * (i1s - q[0]);
  f = j2s/(m*m) + i1s * (i1s - q[0]);
  if (f > 0.0) 
    return f;
  
  return f;
}

/**
   This function computes derivatives of as-th yield function
   with respect of vector sigma.

   @param sig - stress tensor
   @param q   - %vector of the hardening parameter
   @param dfds - %matrix where the resulting derivatives are stored


   4.1.2002
*/
void camclay::deryieldfsigma (matrix &sig, vector &q, matrix &dfds)
{
  double volume,i1s;
  matrix dev(3,3);

  i1s = first_invar (sig)/3.0;

  deviator (sig,dfds);
/*  
  pro fci plasticity prenasobenou m*m  
  cmulm(3.0, dfds);
  volume=m*m/3.0*(2.0*i1s - q[0]);
  dfds[0][0]+=volume;
  dfds[1][1]+=volume;
  dfds[2][2]+=volume;*/

  cmulm(1.0/(m*m), dfds);
  volume=1.0/3.0*(2.0*i1s - q[0]);
  dfds[0][0]+=volume;
  dfds[1][1]+=volume;
  dfds[2][2]+=volume;  
}

/**
   function computes the second derivatives of yield function
   with respect of stress tensor sigma, 

   @param ddfds - tensor of the 4-th order where the are derivatives stored

   19.12.2002
*/
void camclay::dderyieldfsigma (matrix &ddfds)
{
  fillm(0.0, ddfds);
  ddfds[0][0] = ddfds[1][1] = ddfds[2][2] = 2.0/(3.0*m*m) + 2.0/9.0;
  ddfds[0][1] = ddfds[0][2] = ddfds[1][0] = ddfds[1][2] = -1.0/(3.0*m*m) + 2.0/9.0;
  ddfds[2][0] = ddfds[2][1] = ddfds[0][1];
  ddfds[3][3] = 1.0/(m*m);
  ddfds[4][4] = 1.0/(m*m);
  ddfds[5][5] = 1.0/(m*m);
/*
  switch (ssst)
  {
    case planestress :
    case planestrain :
      ddfds[0][0] = ddfds[1][1] = 2.0/(3.0*m*m) + 2.0/9.0;
      ddfds[0][1] = ddfds[1][0] = -1.0/(3.0*m*m) + 2.0/9.0;
      ddfds[2][2] = 1.0/(m*m);
      break;
    case axisymm:
      ddfds[0][0] = ddfds[1][1] = ddfds[2][2] = 2.0/(3.0*m*m) + 2.0/9.0;
      ddfds[0][1] = ddfds[0][2] = ddfds[1][0] = ddfds[1][2] = -1.0/(3.0*m*m) + 2.0/9.0;
      ddfds[2][0] = ddfds[2][1] = -1.0/(3.0*m*m) + 2.0/9.0;
      ddfds[3][3] = 1.0/(m*m);
      break;
    case spacestress:
      ddfds[0][0] = ddfds[1][1] = ddfds[2][2] = 2.0/(3.0*m*m) + 2.0/9.0;
      ddfds[0][1] = ddfds[0][2] = ddfds[1][0] = ddfds[1][2] = -1.0/(3.0*m*m) + 2.0/9.0;
      ddfds[2][0] = ddfds[2][1] = ddfds[0][1];
      ddfds[3][3] = ddfds[4][4] = ddfds[5][5] = 1.0/(m*m);
      break;
    default:
      fprintf(stderr, "\n\n Error - unknown type of stress/strain state is required.\n");
      fprintf(stderr, "\n\n Function camclay::dderyieldfsigma, (file %s, line %d)\n", __FILE__, __LINE__);
      break;
  }*/
}

/**
   This function computes derivatives of plastic potential function
   with respect of vector sigma.

   @param sig - stress tensor
   @param q   - %vector of the hardening parameter
   @param dgds - %matrix where the resulting derivatives are stored
*/
void camclay::derpotsigma (matrix &sig, vector &q, matrix &dgds)
{
  deryieldfsigma (sig, q, dgds);
}


/**
   This function computes derivatives of as-th yield function
   with respect of vector of hradening parameters.

   @param sig - stress tensor
   @param dfds - %vector where the resulting derivatives are stored


   4.1.2002
*/
void camclay::deryieldfq(matrix &sig, vector &dfq)
{
  dfq[0] = -first_invar (sig)/3.0;

//  pro fci plasticity prenasobenou m*m  
//    dfq[0] = first_invar (sig)/3.0;
//  dfq[0] *= -m*m;
  return;
}



/**
   This function computes the second derivatives of yield function
   with respect to hradening parameter.

   @param dfds - %matrix, where the resulting derivatives are stored


   12.4.2005
*/
void camclay::deryieldfdqdq(matrix &dfq)
{  
  dfq[0][0] = 0.0;
  return;
}



/**
   This function computes the second derivatives of yield function
   with respect to stresses.

   @param dfds - tensor, where the resulting derivatives are stored.
                 size of dfds = (6,number_of_hardening_param)


   12.4.2005
*/
void camclay::deryieldfdsdq(matrix &dfdsdqt)
{  
  dfdsdqt[0][0] = dfdsdqt[1][0] = dfdsdqt[2][0] = -1.0/3.0;
  return;
}



/**
   This function computes the derivatives of hardening parameters
   with respect to stresses.

  @param sigt - stress tensor
  @param qt   - %vector of hardening variables
  @param dqds - tensor, where the resulting derivatives are stored.
                 size of dqds = (3,3)


   12.4.2005
*/
void camclay::dqdsigma(matrix &sigt, vector &qt, matrix &dqds)
{  
  double i1s;

  i1s = first_invar (sigt)/3.0;

  fillm(0.0, dqds);
  dqds[0][0] = dqds[1][1] = dqds[2][2] = qt[0]*kappa/i1s*1/3.0;
  return;
}



/**
   This function computes the derivatives of hardening function
   with respect to hardening parameters.

   
   @param ipp - integration point number
   @param ido - index of internal variables for given material in the ipp other array
   @param dgamma - increment of consistency parameter gamma
   @param qt   - %vector of hardening variables
   @param dqdq - tensor, where the resulting derivatives are stored.
                 size of dqdq = (number_of_hardening_param, number_of_hardening_param)


   12.4.2005
*/
void camclay::dhardfdq(long ipp, long ido, double dgamma, vector &qt, matrix &dqdq)
{  
  double v_ini;
  long ncompstr = Mm->ip[ipp].ncompstr;

  fillm(0.0, dqdq);
  v_ini = Mm->ip[ipp].eqother[ido+ncompstr+4];
  dqdq[0][0] = -dgamma * v_ini * qt[0]/(kappa-lambda);
  return;
}



/**
   This function computes derivatives of hardening paramters
   with respect of consistency parameter gamma.

   @param ipp - integration point number
   @param ido - index of internal variables for given material in the ipp other array
   @param sig - stress tensor
   @param qtr  - %vector of hardening variables
   @param dqdg - %matrix where the resulting derivatives are stored


   4.1.2002
*/
void camclay::der_q_gamma(long ipp, long ido, matrix &sig, vector &qtr, vector &dqdg)
{
  double v_ini, i1s, v_pc0, v_lambda1;
  double v_kappa1, p1, pc_0, p_ini;
  long ncompstr = Mm->ip[ipp].ncompstr;
  long i;

  // original specific volume before any loading
  if (qtr[0] == 0.0)
  {
    v_kappa1 = Mm->ic[ipp][0];
    p1       = Mm->ic[ipp][1];
    pc_0     = Mm->ic[ipp][2];
    p_ini    = 0.0;
    for (i=0; i < ncompstr; i++)
      p_ini += Mm->ic[ipp][i+3];
    p_ini /= 3.0;
    Mm->ip[ipp].eqother[ido+ncompstr+2] = v_pc0     = v_kappa1 - kappa*log(pc_0/p1);
    Mm->ip[ipp].eqother[ido+ncompstr+3] = v_lambda1 = v_pc0 + lambda*log(pc_0/p1);
    Mm->ip[ipp].eqother[ido+ncompstr+4] = v_ini     = v_pc0 + kappa*log(pc_0/p_ini);
  }
  else
  {
    v_ini = Mm->ip[ipp].eqother[ido+ncompstr+4];
  }
  i1s = first_invar (sig)/3.0;

  dqdg[0] = qtr[0] * (v_ini)/(kappa-lambda) * (2.0*i1s-qtr[0]);

//  pro fci plasticity prenasobenou m*m  
//  dqdg[0] = qtr[0] * (-v_ini)/(lambda-kappa)*m*m*(2.0*i1s-qtr[0]);
  return;
}


/**
  This function computes plastic modulus.

  @param ipp - integration point number
  @param ido - index of internal variables for given material in the ipp other array
  @param sig - stress tensor
  @param qtr - %vector of hardening parameters

  @retval The function returns value of the plastic modulus.
*/
double camclay::plasmodscalar(long ipp, long ido, matrix &sig, vector &qtr)
{
  double ret;
  vector dfq(1);
  vector dqg(1);

  deryieldfq(sig, dfq);
  der_q_gamma(ipp, ido, sig, qtr, dqg);
  scprd(dfq, dqg, ret);

  return -ret;
}

/**
  This function computes new value of teh hardening parameter q.

  @param ipp - integration point pointer
  @param ido - index of internal variables for given material in the ipp other array
  @param eps - %vector of the reached strains
  @param epsp - %vector of the reached plastic strains
  @param sig - reached stress tensor
  @param ssst - stress/strain state parameter (for the used vector_tensor function)
  @param q    - %vector of the hardening parameters

*/
void camclay::updateq(long ipp, long ido, vector &eps, vector &epsp, matrix &sig, vector &q)
{
  matrix epst(3,3);
  double v_kappa1, v_lambda1, p1, pc_0, v_pc0, vk, epsv, i1s, v_ini, p_ini, pc_new;