chen0407.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"

chen::chen (void)
{
  fyc=0.0;
  fyt=0.0;
  fybc=0.0;
  fc=0.0;
  ft=0.0;
  fbc=0.0;
  epsuc=0.0;
  epsut=0.0;
  
  ay=0.0;  au=0.0;
  ky=0.0;  ku=0.0;
  
  hp=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 (FILE *in)
{
  fscanf (in,"%lf %lf %lf  %lf %lf %lf  %lf %lf  %lf",&fyc,&fyt,&fybc,&fc,&ft,&fbc,&epsuc,&epsut,&hp);
  sra.read (in);
}


/**
   function returns elastic stiffness matrix
   
   @param d - elastic stiffness matrix
   
   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
    //fprintf (stderr,"\n\n tangent stiffness matrix is not implemented yet\n");
    
    matrix ad(d.m,d.n);
    Mm->elmatstiff (ad,ipp);
    tangentstiff (ad,d,ipp,ido);
  }
}


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-10){
      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);
  }

}

void chen::plasmod (matrix &h)
{
  h[0][0]=-1.0*hp;
}

/**
   function evaluates yield function for given stresses
   
   @param sig - stresses
   @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);

  deviator (sig,dev);
  invar=first_invar (sig);
  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){
    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){
    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
   
   JK, 20.4.2005
*/
void chen::deryieldfdsigma (matrix &sig,vector &q,matrix &dfds)
{
  double invar,j2,zone,kappa;
  matrix dev(3,3);
  
  deviator (sig,dev);
  invar=first_invar (sig);
  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;
      
      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
   
   JK, 20.4.2005
*/
void chen::deryieldfdsigmadsigma (matrix &sig,matrix &dfdsds)
{
  double invar,j2,zone;
  matrix dev(3,3);

  deviator (sig,dev);
  invar=first_invar (sig);
  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;