chen.cpp

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

#include <math.h>
#include "chen.h"
#include "global.h"
#include "genfile.h"
#include "intpoints.h"
#include "vecttens.h"
#include "alias.h"
#include "stochdriver.h"
#include "matrix.h"
#include "hardsoft.h"

chen::chen (void)
{
  //  compression yield stress
  fyc=0.0;
  //  tension yield stress
  fyt=0.0;
  //  bi-axial compression yield stress
  fybc=0.0;
  //  compression ultimate stress
  fc=0.0;
  //  tension ultimate stress
  ft=0.0;
  //  bi-axial compression ultimate stress
  fbc=0.0;
  
  //  auxiliary constants
  ay=0.0;  au=0.0;
  ky=0.0;  ku=0.0;
  
  
  //  general strain/stress state
  state=1;
  //  compression
  //state=2;
  //  tension
  //state=3;
}

chen::~chen (void)
{

}

/**
   function reads material parameters of the model
   
   JK, 21.2.2005
*/
void chen::read (XFILE *in)
{
  //  compression yield stress
  //  tension yield stress
  //  bi-axial compression yield stress
  //  compression ultimate stress
  //  tension ultimate stress
  //  bi-axial compression ultimate stress
  xfscanf (in,"%lf %lf %lf  %lf %lf %lf",&fyc,&fyt,&fybc,&fc,&ft,&fbc);
  //  type of stress return algorithm
  sra.read (in);
  //  type of hardening/softening
  hs.read (in);
  
}


/**
   function returns elastic stiffness %matrix
   
   @param d - elastic stiffness %matrix
   @param ipp - integration point id
   @param ido - first index in the array other
   
   4.8.2001
*/
void chen::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);
  }
}

/**
   function assembles tangent stiffness %matrix
   
   @param d - elastic stiffness %matrix of the material
   @param td - tangent stiffness %matrix
   @param ipp - integration point id
   @param ido - first index in the array other
   
*/
void chen::tangentstiff (matrix &d,matrix &td,long ipp,long ido)
{
  long ncomp=Mm->ip[ipp].ncompstr;
  double denom1,denom2,denom,gamma;
  vector str,av(d.m),q(1),dfdq(1);
  matrix sig(3,3),am(d.m,d.n),dfds(3,3);
  
  gamma=Mm->ip[ipp].eqother[ido+ncomp];
  if (gamma<1.0e-10){
    copym (d,td);
  }
  else{
    allocv (ncomp,str);
    
    //  actual stress
    Mm->givestress (0,ipp,str);
    vector_tensor (str,sig,Mm->ip[ipp].ssst,stress);
    //  actual hardening parameter
    q[0]=Mm->ip[ipp].eqother[ido+ncomp+1];
    
    deryieldfdsigma (sig,q,dfds);
    tensor_vector (str,dfds,Mm->ip[ipp].ssst,stress);
    
    if (Mm->ip[ipp].ssst==planestress){
      vector auxstr(3);
      auxstr[0]=str[0];auxstr[1]=str[1];auxstr[2]=str[2];
      destrv (str);
      allocv (d.m,str);
      str[0]=auxstr[0];str[1]=auxstr[1];str[2]=auxstr[2];
    }
    
    mxv (d,str,av);
    scprd (av,str,denom1);
    
    scprd (str,str,denom2);
    
    deryieldfdq (sig,q,dfdq);
    
    denom2 = sqrt(denom2)*dfdq[0];
    
    denom=denom1-denom2;
    
    //if (fabs(denom)<1.0e-1){
    if (fabs(denom)<1.0e-1 || denom>1.0e5){
      copym (d,td);
    }
    else{
      vxv (str,str,am);
      mxm (d,am,td);
      mxm (td,d,am);
      
      cmulm (1.0/denom,am);
      
      subm (d,am,td);
    }
    
    destrv (str);
  }

}


/**
   function evaluates yield function for given stresses
   
   @param sig - stress components stored in 3x3 %matrix
   @param q - internal parameters (hardening)
   
   JK, 4.8.2001
*/
double chen::yieldfunction (matrix &sig,vector &q)
{
  double f,invar,j2,zone;//kappa;
  matrix dev(3,3);
  
  //  computation of stress deviator
  deviator (sig,dev);
  //  first invariant of the stress tensor
  invar=first_invar (sig);
  //  second invariant of the stress deviator
  j2=second_invar(dev);
  
  zone=sqrt(j2)+invar/sqrt(3.0);
  
  if (state==1){
    //  compression-compression
    if(invar<0.0 && zone<0.0){
      ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
      au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
      ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
      ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
      
      if (q[0]>1.0)  q[0]=1.0;
      //  yield function with hardening
      f = j2 + ay/3.0*invar - ky + q[0]*((au-ay)*invar/3.0 - (ku-ky));
    }
    //  otherwise
    else{
      ay = (fyc-fyt)/2.0;
      au = (fc-ft)/2.0;
      ky = fyc*fyt/6.0;
      ku = fc*ft/6.0;
      
      if (q[0]>1.0)  q[0]=1.0;
      //  yield function with hardening
      f = j2 - (invar*invar)/6.0 + ay/3.0*invar - ky +q[0]*((au-ay)*invar/3.0 - (ku-ky));
    }
  }
  if (state==2){
    //  compression-compression
    ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
    au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
    ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
    ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
    
    if (q[0]>1.0)  q[0]=1.0;
    //  yield function with hardening
    f = j2 + ay/3.0*invar - ky + q[0]*((au-ay)*invar/3.0 - (ku-ky));
  }
  if (state==3){
    //  non compression-compression zone
    ay = (fyc-fyt)/2.0;
    au = (fc-ft)/2.0;
    ky = fyc*fyt/6.0;
    ku = fc*ft/6.0;
    
    if (q[0]>1.0)  q[0]=1.0;
    //  yield function with hardening
    f = j2 - (invar*invar)/6.0 + ay/3.0*invar - ky +q[0]*((au-ay)*invar/3.0 - (ku-ky));
  }
  
  return f;
}


/**
   function evaluates derivatives of yield function
   with respect to stress components
   
   @param sig - stress tensor (stored in 3x3 %matrix)
   @param q - hardening parameters
   @param dfds - derivatives of yield function with respect to stress components (stored in 3x3 %matrix)
   
   JK, 20.4.2005
*/
void chen::deryieldfdsigma (matrix &sig,vector &q,matrix &dfds)
{
  double invar,j2,zone;//kappa;
  matrix dev(3,3);
  
  //  computation of stress deviator
  deviator (sig,dev);
  //  first invariant of the stress tensor
  invar=first_invar (sig);
  //  second invariant of the stress deviator
  j2=second_invar(dev);
  
  fillm (0.0,dfds);

  zone=sqrt(j2)+invar/sqrt(3.0);
  
  if (state==1){
    //  compression-compression
    if(invar<0.0 && zone<0.0){
      ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
      au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
      ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
      ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
      
      if (q[0]>1.0)  q[0]=1.0;
      
      dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + ay/3.0 + q[0]/3.0*(au-ay);
      dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + ay/3.0 + q[0]/3.0*(au-ay);
      dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + ay/3.0 + q[0]/3.0*(au-ay);
      
      dfds[0][1]=2*sig[0][1];
      dfds[1][0]=2*sig[0][1];
      dfds[0][2]=2*sig[0][2];
      dfds[2][0]=2*sig[0][2];
      dfds[1][2]=2*sig[1][2];
      dfds[2][1]=2*sig[1][2];
    }
    //  otherwise
    else{ 
      ay = (fyc-fyt)/2.0;
      au = (fc-ft)/2.0;
      ky = fyc*fyt/6.0;
      ku = fc*ft/6.0;
      
      if (q[0]>1.0)  q[0]=1.0;
      
      dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + ay/3.0 - 1.0/3.0*invar + q[0]/3.0*(au-ay);
      dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + ay/3.0 - 1.0/3.0*invar + q[0]/3.0*(au-ay);
      dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + ay/3.0 - 1.0/3.0*invar + q[0]/3.0*(au-ay);
      
      dfds[0][1]=2*sig[0][1];
      dfds[1][0]=2*sig[0][1];
      dfds[0][2]=2*sig[0][2];
      dfds[2][0]=2*sig[0][2];
      dfds[1][2]=2*sig[1][2];
      dfds[2][1]=2*sig[1][2];
    }
  }
  if (state==2){
    ay = (fybc*fybc-fyc*fyc)/(2.0*fybc-fyc);
    au = (fbc*fbc-fc*fc)/(2.0*fbc-fc);
    ky = fyc*fybc*(2*fyc-fybc)/3.0/(2*fybc-fyc);
    ku = fc*fbc*(2*fc-fbc)/3.0/(2*fbc-fc);
    
    if (q[0]>1.0)  q[0]=1.0;
    
    dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + ay/3.0 + q[0]/3.0*(au-ay);
    dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + ay/3.0 + q[0]/3.0*(au-ay);
    dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + ay/3.0 + q[0]/3.0*(au-ay);
    
    dfds[0][1]=2*sig[0][1];
    dfds[1][0]=2*sig[0][1];
    dfds[0][2]=2*sig[0][2];
    dfds[2][0]=2*sig[0][2];
    dfds[1][2]=2*sig[1][2];
    dfds[2][1]=2*sig[1][2];
  }
  if (state==3){
    ay = (fyc-fyt)/2.0;
    au = (fc-ft)/2.0;
    ky = fyc*fyt/6.0;
    ku = fc*ft/6.0;
    
    if (q[0]>1.0)  q[0]=1.0;
    
    dfds[0][0]=1.0/3.0*(2.0*sig[0][0]-sig[1][1]-sig[2][2]) + ay/3.0 - 1.0/3.0*invar + q[0]/3.0*(au-ay);
    dfds[1][1]=1.0/3.0*(2.0*sig[1][1]-sig[0][0]-sig[2][2]) + ay/3.0 - 1.0/3.0*invar + q[0]/3.0*(au-ay);
    dfds[2][2]=1.0/3.0*(2.0*sig[2][2]-sig[0][0]-sig[1][1]) + ay/3.0 - 1.0/3.0*invar + q[0]/3.0*(au-ay);
    
    dfds[0][1]=2*sig[0][1];
    dfds[1][0]=2*sig[0][1];
    dfds[0][2]=2*sig[0][2];
    dfds[2][0]=2*sig[0][2];
    dfds[1][2]=2*sig[1][2];
    dfds[2][1]=2*sig[1][2];
  }
}

/**
   function computes second derivatives of yiled function with
   respect to stresses
   
   @param sig - stress tensor (stored in 3x3 %matrix)
   @param dfdsds - second derivatives of yield function with respect to stress components (stored in 6x6 %matrix)

   JK, 20.4.2005
*/
void chen::deryieldfdsigmadsigma (matrix &sig,matrix &dfdsds)
{
  double invar,j2,zone;
  matrix dev(3,3);

  //  computation of stress deviator
  deviator (sig,dev);
  //  first invariant of the stress tensor
  invar=first_invar (sig);
  //  second invariant of the stress deviator
  j2=second_invar(dev);
  
  fillm (0.0,dfdsds);
  
  zone = sqrt(j2)+invar/sqrt(3.0);

  if (state==1){
    //  compression-compression
    if(invar<0.0 && zone<0.0){
      dfdsds[0][0] =  2.0/3.0;
      dfdsds[0][1] = -1.0/3.0;
      dfdsds[0][2] = -1.0/3.0;
      
      dfdsds[1][0] = -1.0/3.0;
      dfdsds[1][1] =  2.0/3.0;
      dfdsds[1][2] = -1.0/3.0;
      
      dfdsds[2][0] = -1.0/3.0;
      dfdsds[2][1] = -1.0/3.0;
      dfdsds[2][2] =  2.0/3.0;
      
      dfdsds[3][3] = 2.0;
      dfdsds[4][4] = 2.0;
      dfdsds[5][5] = 2.0;
    }
    //  otherwise
    else{ 
      dfdsds[0][0] =  1.0/3.0;
      dfdsds[0][1] = -2.0/3.0;
      dfdsds[0][2] = -2.0/3.0;
      
      dfdsds[1][0] = -2.0/3.0;
      dfdsds[1][1] =  1.0/3.0;
      dfdsds[1][2] = -2.0/3.0;
      
      dfdsds[2][0] = -2.0/3.0;
      dfdsds[2][1] = -2.0/3.0;
      dfdsds[2][2] =  1.0/3.0;
      
      dfdsds[3][3] = 2.0;
      dfdsds[4][4] = 2.0;
      dfdsds[5][5] = 2.0;
    }
  }
  if (state==2){
    dfdsds[0][0] =  2.0/3.0;
    dfdsds[0][1] = -1.0/3.0;
    dfdsds[0][2] = -1.0/3.0;
    
    dfdsds[1][0] = -1.0/3.0;
    dfdsds[1][1] =  2.0/3.0;
    dfdsds[1][2] = -1.0/3.0;
    
    dfdsds[2][0] = -1.0/3.0;
    dfdsds[2][1] = -1.0/3.0;
    dfdsds[2][2] =  2.0/3.0;
    
    dfdsds[3][3] = 2.0;
    dfdsds[4][4] = 2.0;
    dfdsds[5][5] = 2.0;
    
  }
  if (state==3){
    dfdsds[0][0] =  1.0/3.0;
    dfdsds[0][1] = -2.0/3.0;
    dfdsds[0][2] = -2.0/3.0;
    
    dfdsds[1][0] = -2.0/3.0;
    dfdsds[1][1] =  1.0/3.0;
    dfdsds[1][2] = -2.0/3.0;
    
    dfdsds[2][0] = -2.0/3.0;
    dfdsds[2][1] = -2.0/3.0;
    dfdsds[2][2] =  1.0/3.0;
    
    dfdsds[3][3] = 2.0;
    dfdsds[4][4] = 2.0;
    dfdsds[5][5] = 2.0;
  }
}