shefplast.cpp

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

#include <math.h>
#include <float.h>
#include <stdlib.h>
#include "shefplast.h"
#include "global.h"
#include "vecttens.h"
#include "intpoints.h"
#include "matrix.h"
#include "vector.h"
#include "tensor.h"

#define nijac 20
#define limit 1.0e-8

/**
   This constructor inializes attributes to zero values.
*/
shefplast::shefplast (void)
{
  rc = 0.0;
  rt = 0.0;
  gamma = 0.99;
  alpha = 0.5;
  p = 0.0;
  a = -1.0;
  k0 = 0.1;
  ah = bh = ch = 0.0;
  kh0 = 0.001;
  allocm(3,3,dev);
  khat = rhoc = r = m = xi = ft = 0.0;
  b0 = b1 = c0 = c1 = c2 = c3 = c4 = c5 = theta = d0 = d1 = d2 = 0.0;
}

/**
   This destructor is only for the formal purposes.
*/
shefplast::~shefplast (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 shefplast::read (XFILE *in)
{
  xfscanf (in, "%lf %lf %lf %lf %lf %lf", &rc, &rt, &p, &ah, &bh, &ch);
  sra.read (in);
  
  ft = rt/rc;
  m = 3.0*(1.0-pow(ft,2.0/gamma));
  m /= ft+2.0*pow(ft,1.0/gamma);
}

void shefplast::compute_khat(matrix &sig, vector &q)
{
  double k;

  k = k0 + (1.0 - k0)*sqrt(q[0]*(2.0-q[0]));
  khat = pow(k,2.0*p)*(1.0 - xi*xi*(1.0-k)*(1.0-k)/(a*a));
}

void shefplast::compute_rhoc(matrix &sig)
{
  double tmp;

  tmp = -m+sqrt(m*m-12.0*sqrt(3.0)*m*xi+36.0);
  rhoc = pow(1.0/6.0, gamma)*sqrt(2.0/3.0)*pow(tmp,gamma);
}
 
void shefplast::compute_r(matrix &sig)
{
  double tmp, rhoe;
  double pi=4.0*atan(1.0) ;

  tmp = -m+sqrt(m*m-3.0*sqrt(3.0)*m*xi+9.0);
  rhoe = pow(1.0/3.0, gamma)*sqrt(2.0/3.0)*pow(tmp, gamma);
  b0 = rhoe/rhoc;
  if(b0<=0.5) {
    b0 = 0.5 + 1.0e-6 ;
  }
  else if(b0>1.0) {
    b0 = 1.0 ;
  }
  b1 = sqrt(3.0)*(1.0-alpha)*b0/(1.0+b0)+2.0*alpha*b0/sqrt(3.0);
  c0 = (2.0-sqrt(3.0)*b1)*(2.0*b0-sqrt(3.0)*b1) / ((b1*(1.0+b0)-sqrt(3.0)*b0)*(b1*(1.0+b0)-sqrt(3.0)*b0));
  c1 = 3.0-c0*(1.0+b0)*(1.0+b0);
  c2 = 1.0+3.0*c0*(1.0-b0)*(1.0-b0);
  c3 = 2.0*c0*sqrt(3.0)*(1.0-b0*b0);
  c4 = (1.0+b0)*(1.0-b0*c0);
  c5 = (1.0-b0)*(1.0-3.0*b0*c0);
  if((-sqrt(3.0)/2.0*3.0*j3s/pow(j2s,1.5)) > 1.0) {
    theta = pi/6.0 ;
  }
  else if((-sqrt(3.0)/2.0*3.0*j3s/pow(j2s,1.5)) < -1.0) {
    theta = -pi/6.0 ;
  }
  else
    theta = 1.0/3.0*asin(-sqrt(3.0)/2.0*3.0*j3s/pow(j2s,1.5));
  d0 = c1*cos(theta)*cos(theta) - c2*sin(theta)*sin(theta) + c3*sin(theta)*cos(theta);
  d1 = 2.0*(c4*sqrt(3.0)*cos(theta)-c5*sin(theta));
  d2 = b0*(4.0-3.0*b0*c0);
  r = 2.0*d0/(d1-sqrt(d1*d1-4.0*d0*d2));
}

void shefplast::compute_classvar(matrix &sig, vector &q)
{
  deviator(sig,dev);
  i1s = first_invar(sig);
  j2s = second_invar(dev);
  j3s = third_invar(dev);  
  xi = i1s/(rc*sqrt(3.0));
  compute_khat(sig, q);
  compute_rhoc(sig);
  compute_r(sig);
}

/**
   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

   14.1.2004
*/
double shefplast::yieldfunction (matrix &sig, vector &q)
{
  double rho, f;

  compute_classvar(sig, q);
  rho = pow(2.0*j2s,0.5)/rc;
  f = rho*rho - khat*rhoc*rhoc /(r*r);

  return f;
}

/**
   This function computes derivatives of r parameter
   q = 1/r i.e. q in this case does not mean hardening parameters vector.
   with respect of tensor sigma.

   @param sig - stress tensor
   @param drds - %matrix where the resulting derivatives are stored


   14.1.2004
*/
void shefplast::derqdsigma (matrix &sig, matrix &drds)
{
  long i, j, k;
  double dqdd0, dqdd1, dqdd2, dd0db0, dd1db0, dd2db0;
  double db1db0, dc0db0,dc1db0,dc2db0,dc3db0,dc4db0,dc5db0;
  double dqdb0, dqdth, db0dxi;
  double dd0dth, dd1dth;  
  double dthdj2, dthdj3;
  double tmp = sqrt(d1*d1-4.0*d0*d2);
  double tmp1, tmp2, tmp3;
  double drhocdxi,drhoedxi,rhoe;
  double pi=4.0*atan(1.0) ;
  double epsi=1.0e-6 ;
  matrix dthds(3,3);
  matrix dj3ds(3,3);
  matrix dxids(3,3);
  matrix db0ds(3,3);

  dqdd0 = 4*d2*d0/2.0/tmp - d1+tmp;
  dqdd0 /= (2.0*d0*d0);
  dqdd1 = (1.0-d1/tmp)/2.0/d0;
  dqdd2 = 1.0/tmp;
  db1db0 = sqrt(3.0)*(1.0-alpha)/(1.0+b0)/(1.0+b0)+2.0*alpha/sqrt(3.0);
 
  tmp1 = -sqrt(3.0)*db1db0*(2.0*b0-sqrt(3.0)*b1)+(2.0-sqrt(3.0)*b1)*(2.0-sqrt(3.0)*db1db0);
  tmp2 = b1*(1.0+b0)-sqrt(3.0)*b0;
  tmp3 = 2.0*tmp2*(db1db0*(1.0+b0)+b1-sqrt(3.0))*(2.0-sqrt(3.0)*b1)*(2.0*b0-sqrt(3.0)*b1);
  
  dc0db0 = (tmp1*tmp2*tmp2 - tmp3)/pow(tmp2, 4.0);
  dc1db0 = -dc0db0*(1.0+b0)*(1.0+b0) - c0*2.0*(1.0+b0);
  dc2db0 = 3.0*dc0db0*(1.0-b0)*(1.0-b0)-6.0*c0*(1.0-b0);
  dc3db0 = 2.0*sqrt(3.0)*dc0db0*(1.0-b0*b0)-4.0*sqrt(3.0)*c0*b0;
  dc4db0 = 1.0-b0*c0-(1.0+b0)*(c0+b0*dc0db0);
  dc5db0 = 3.0*b0*c0-1.0-3.0*(1.0-b0)*(c0+b0*dc0db0);

  dd0db0 = dc1db0*cos(theta)*cos(theta)-dc2db0*sin(theta)*sin(theta)+dc3db0*sin(theta)*cos(theta);
  dd1db0 = 2.0*(dc4db0*sqrt(3.0)*cos(theta)-dc5db0*sin(theta));
  dd2db0 = 4.0 - 3.0*b0*c0 - 3.0*b0*(c0+b0*dc0db0);
  
  dqdb0 = dqdd0*dd0db0 + dqdd1*dd1db0 + dqdd2*dd2db0;

  // Compute dB0ds
  tmp1 = sqrt(m*m-12.0*sqrt(3.0)*m*xi+36.0);
  drhocdxi = -gamma*pow((tmp1-m)/6.0, gamma-1.0);  
  drhocdxi *= m*sqrt(2.0)/tmp1;

  tmp1 = sqrt(m*m-3.0*sqrt(3.0)*m*xi+9.0);
  drhoedxi = -gamma*pow((tmp1-m)/3.0, gamma-1.0);  
  drhoedxi *= m/sqrt(2.0)/tmp1;
  
  tmp = -m+sqrt(m*m-3.0*sqrt(3.0)*m*xi+9.0);
  rhoe = pow(1.0/3.0, gamma)*sqrt(2.0/3.0)*pow(tmp, gamma);

  db0dxi = (drhoedxi*rhoc-drhocdxi*rhoe)/rhoc/rhoc;
  for(i=0; i < 3; i++)
    dxids[i][i] = 1.0/rc/sqrt(3.0);
  cmulm(db0dxi,dxids,db0ds);
  // End compute dB0ds
  
  dd0dth = -2.0*(c1+c2)*cos(theta)*sin(theta)+c3*(cos(theta)*cos(theta)-sin(theta)*sin(theta));
  dd1dth = 2.0*(-c4*sqrt(3.0)*sin(theta)-c5*cos(theta));
  dqdth = dqdd0*dd0dth+dqdd1*dd1dth;

  // Compute dthds
  for (i=0; i < 3; i++)
    {
      for (j=0; j < 3; j++)
	{
	  for (k = 0; k < 3; k++)
	    dj3ds[i][j] += dev[i][k]*dev[k][j];
	}
      dj3ds[i][i] -= 2.0/3.0*j2s;
    }
  if(fabs(fabs(theta/(pi/6.0))-1.0)<epsi) {
    for (i = 0; i < 3; i++) {
      for (j = 0; j < 3; j++)
	dthds[i][j] = 0.0 ; 
    }
  }
  else {
    dthdj2 = 3.0 * sqrt(3.0)*j3s/4.0/cos(3.0*theta)/pow(j2s,2.5);
    dthdj3 = -sqrt(3.0)/2.0/cos(3.0*theta)/pow(j2s,1.5); 
    for (i = 0; i < 3; i++) {
      for (j = 0; j < 3; j++)
	dthds[i][j] = dthdj2*dev[i][j] + dthdj3*dj3ds[i][j];
    }
  }
  // End compute dthds
  
  // The final addition
  for (i = 0; i < 3; i++) {
    for (j = 0; j < 3; j++) {
      drds[i][j] = dqdb0*db0ds[i][j]+dqdth*dthds[i][j];
    }
  }
  cmulm((2.0/r)*khat*rhoc*rhoc,drds);
}

/**
   This function computes derivatives of khat parameter
   with respect of tensor sigma.

   @param sig - stress tensor
   @param q   - %vector of the hardening parameter
   @param dkhatds - %matrix where the resulting derivatives are stored
   @param ido - index of internal variables for given material in the ipp other array


   4.1.2002
*/
void shefplast::derkhatds (long ipp,matrix &sig, vector &q, matrix &dkhatds, long ido)
{
  long i;
  double k, dkhatdxi;

  matrix dxids(3,3);

  k = k0 + (1.0-k0)*sqrt(q[0]*(2.0-q[0]));

  dkhatdxi = -2.0*pow(k,2.0*p)*(1.0-k)*(1.0-k)*xi/a/a;

  for(i=0; i < 3; i++)
    dxids[i][i] = 1.0/rc/sqrt(3.0);

  // Compute the full term corresponding to dkhatds in dfds
  cmulm(dkhatdxi, dxids, dkhatds);
  cmulm(rhoc*rhoc/r/r, dkhatds);
}

/**
   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
   @param ido - index of internal variables for given material in the ipp other array


   4.1.2002
*/
void shefplast::deryieldfsigma (long ipp,matrix &sig, vector &q, matrix &dfds,long ido)
{
  long i;
  double tmp, tmp1, rhoc;
  matrix drhods(3,3);
  matrix drhocds(3,3);
  matrix dqds(3,3);
  matrix dkhatds(3,3);
  matrix drhoeds(3,3);

  compute_classvar (sig, q);

  // Compute the full term corresponding to drhods in dfds
  cmulm(2.0/rc/rc,dev,drhods);
  // Compute the full term corresponding to drhocds in dfds
  tmp = sqrt(m*m-12.0*sqrt(3.0)*m*xi+36.0);
  tmp1 = -gamma*pow((tmp-m)/6.0,gamma-1.0);  
  tmp1 *= m*sqrt(2.0)/tmp;
  tmp1 /= (sqrt(3.0)*rc);
  for (i = 0; i < 3; i++)
    drhocds[i][i] = tmp1;
  rhoc = pow((1.0/6.0),gamma)*sqrt(2.0/3.0)*pow((tmp-m),gamma);
  cmulm(2.0*rhoc*khat/r/r,drhocds);
  
  // Compute the full term corresponding to dqds in dfds
  derqdsigma(sig, dqds); /* OK 2 */

  derkhatds(ipp,sig, q, dkhatds, ido);

  addm(dkhatds, drhocds, dfds);
  addm(dfds, dqds, dfds);
  subm(drhods, dfds, dfds);
}

/**
   function computes the second derivatives of yield function
   with respect of vector sigma

   @param sig - stress components
   @param ssst - assumed stress/strain state at the given ip

   19.12.2002
*/
void shefplast::dderyieldfsigma (matrix &ddfds, strastrestate ssst)
{
}

/**
   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
   @param ido - index of internal variables for given material in the ipp other array
*/
void shefplast::derpotsigma (long ipp,matrix &sig, vector &q, matrix &dgds,long ido)
{
  deryieldfsigma (ipp,sig, q, dgds,ido);
}


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

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


   4.1.2002
*/
void shefplast::deryieldfq(matrix &sig, vector &q, vector &dfq)
{
  double k ;
  double dkdkh, dkhatdk ;

  k = k0 + (1.0-k0)*sqrt(q[0]*(2.0-q[0]));
  dkhatdk = 2.0*p*pow(k,2.0*p-1.0)*(1.0-xi*xi*(1.0-k)*(1.0-k)/a/a)
    +2.0*(1.0-k)*(pow(k,p)*xi/a)*(pow(k,p)*xi/a) ;

  dkdkh = (1.0-k0)*(1.0-q[0])/pow(q[0]*(2.0-q[0]),0.5) ;

  dfq[0] = -dkhatdk*dkdkh*rhoc*rhoc/r/r;
  return;
}

/**
   This function computes material stiffnes matrix.

   @param d - allocated matrix structure for material stiffness %matrix
   @param ipp - integration point number

*/
void shefplast::matstiff (matrix &d, long ipp,long ido)
{
  if (Mp->stmat==initial_stiff) {
    //  initial elastic matrix
    Mm->elmatstiff (d,ipp);
  }
  if (Mp->stmat==tangent_stiff) {
    //  tangent stiffness matrix

    matrix ad(d.m,d.n);
    Mm->elmatstiff (ad,ipp);
    tangentstiff (ad,d,ipp,ido);
  }
}

/**
   This function returns the tangent stiffness matrix.
   @param d - allocated matrix structure for material stiffness matrix

*/

void shefplast::tangentstiff (matrix &d,matrix &td,long ipp,long ido)
{
  long ncomp=Mm->ip[ipp].ncompstr;
  long i,j;
  double tmps,gamma,hf,dhdk;
  vector mxi,tmpv,str,dfdk(1),dhds,dfdsdk,dfds,av,q(1);
  matrix td7;
  matrix a,ainv,dinv,sig,dfdstens,dfdsds,am,tmp;

  gamma=Mm->ip[ipp].eqother[ido+ncomp];

  if (gamma<1.0e-10) {
    copym (d,td);
  }
  else {
    allocv (ncomp,str);
    allocv (7,mxi);
    allocv (7,tmpv);
    allocv (6,av);
    allocv (6,dhds);
    allocv (6,dfdsdk);
    allocv (6,dfds);

    allocm (3,3,sig);
    allocm(6,6,dinv);
    allocm(7,7,a);
    allocm(7,7,tmp);
    allocm(7,7,ainv);
    allocm(3,3,dfdstens);
    allocm(6,6,dfdsds);
    allocm(6,6,am);
    allocm(d.m+1,d.n+1,td7);

    invm(d,dinv,Mp->zero);

    q[0] = Mm->ip[ipp].eqother[ido+ncomp+1];

    Mm->givestress (0,ipp,str);
    vector_tensor (str,sig,Mm->ip[ipp].ssst,stress);

    deryieldfsigma (ipp,sig,q,dfdstens,ido);
    tensor_vector (dfds,dfdstens,Mm->ip[ipp].ssst,stress);
    numdiff_dfdsdkc (ipp,sig,q,dfdsdk,ido);
    numdiff_dfdsdsc (ipp,sig,q,dfdsds,ido);

    hf = hardening(ipp,sig,q,ido);

    numdiff_dhdkc (ipp,sig,q,dhdk,ido);
    numdiff_dhdsc (ipp,sig,q,dhds,ido);
    deryieldfq (sig,q,dfdk);

    //  block 1,1 (6*6)
    cmulm(gamma,dfdsds,am);
    addm(am,dinv,am);
    
    for (i=0;i<6;i++) {
      for (j=0;j<6;j++) {
	a[i][j]=am[i][j];
      }
    }
    
    //  block 1,2
    cmulv(gamma,dfdsdk,av);
    
    for (i=0;i<6;i++) {
      a[i][6]=av[i];
      mxi[i] = dfds[i] ;
    }
    mxi[6] = dfdk[0] ;
    
    //  block 2,1
    cmulv(-1.0*gamma,dhds,av);
    
    for (i=0;i<6;i++) {
      a[6][i]=av[i];
    }
    
    //  block 2,2
    a[6][6]=1.0-gamma*dhdk;

    invm(a,ainv,Mp->zero);   
    vxm(mxi,ainv,tmpv);

    // We change the last value of vector mxi to get the vector (m,-h)