mohrc.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 "mohrc.h"
#include "global.h"
#include "matrix.h"
#include "vector.h"
#include "intpoints.h"
#include "tensor.h"
#include "vecttens.h"
#include "elastisomat.h"



#define nijac 20
#define limit 1.0e-8

mohrcoulomb::mohrcoulomb (void)
{
  phi=0.0;  c=0.0;  psi=0.0;
}



mohrcoulomb::~mohrcoulomb (void)
{

}



void mohrcoulomb::read (XFILE *in)
{
  xfscanf (in,"%lf %lf %lf",&phi,&c,&psi);
  sra.read (in);
}



/**
   Function computes the value of yield functions

   @param sig - stress components

   10.11.2001
*/
double mohrcoulomb::yieldfunction (vector &psig)
{
  double f,k1,k2, s1, s2, sigma, tau, maxsigt;

  s1 = psig[0];
  s2 = psig[2];
  sigma   =  (s1 + s2) / 2.0;
  tau     = -(s1 - s2) / (2.0 * cos(phi));
  maxsigt =  c / tan(phi);

  k1 = -1.0 + sin(phi);  k2 = 1.0 + sin(phi);
  f = k1*s1 + k2*s2 - 2.0*c*cos(phi);
  return f;
}



/**
   This function computes derivatives of yield function
   with respect of the principal stresses %vector
   
   @param dfds - vector for result derivatives

   10.11.2003
*/
void mohrcoulomb::deryieldfsigma (vector &dfds)
{
  dfds[0] = -1.0 + sin(phi);
  dfds[1] =  0.0;
  dfds[2] =  1.0 + sin(phi);
  return;
}



/**
   This function computes derivatives of plastic potential function
   with respect of the principal stresses %vector.
   
   @param dgds - vector for result derivatives

   10.11.2003
*/
void mohrcoulomb::derplaspotsigma (vector &dgds)
{
  dgds[0] = -1.0 + sin(psi);
  dgds[1] =  0.0;
  dgds[2] =  1.0 + sin(psi);
  return;
}



/**
  This function returns material stiffness %matrix for given integration point.
  Result is stored in the parameter d.

  @param d   - material stiffness matrix
  @param ipp - integration point pointer
  @param ido - index of internal variables for given material in the ipp other array
*/
void mohrcoulomb::matstiff (matrix &d, long ipp,long ido)
{
  matrix ad(d.m,d.n);
  switch (Mp->stmat)
  {
    case initial_stiff:
      //  initial elastic matrix
      Mm->elmatstiff (d,ipp);
      break;
    case tangent_stiff:
      Mm->elmatstiff (ad,ipp);
      tangentstiff (ad,d,ipp,ido);
      break;
    default :
      fprintf(stderr, "\n\nUnknown type of the stiffness matrix is require in the function\n");
      fprintf(stderr, " mohrcoulomb::matstiff(), (file %s, line %d)\n", __FILE__, __LINE__);
      break;
  }
}



void mohrcoulomb::tangentstiff (matrix &d,matrix &td,long ipp,long ido)
{
  long ncomp=Mm->ip[ipp].ncompstr;
  double denom,gamma;
  vector str,av(d.m),q(1),psig(3), b(3);
  matrix sig(3,3),am(d.m,d.n),pvect(3,3);
  
  gamma=Mm->ip[ipp].eqother[ido+ncomp];
  if (gamma<1.0e-10){
    copym (d,td);
  }
  else{
    
    allocv (ncomp,str);
    
    Mm->givestress (0,ipp,str);
    vector_tensor (str,sig,Mm->ip[ipp].ssst,stress);
    princ_val (sig, psig, pvect, nijac, limit, Mp->zero, 3);
    deryieldfsigma (psig);
    mtxv(pvect, psig, b);
    if ((Mm->ip[ipp].ssst==planestress) || (Mm->ip[ipp].ssst==planestrain)){
      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,b,av);
    scprd (av,b,denom);
    
    
    //    q[0] = Mm->ip[ipp].eqother[ido+ncomp+1];
    //    denom+= plasmodscalar(q);
    
    if (fabs(denom)<1.0e-10){
      copym (d,td);
    }
    else{
      vxv (b,b,am);
      mxm (d,am,td);
      mxm (td,d,am);
      
      cmulm (1.0/denom,am);
      
      subm (d,am,td);
    }
  }
  
}



void mohrcoulomb::pelmatstiff (long ipp, matrix &d, double e, double nu)
{
  double b;

  switch (Mm->ip[ipp].ssst)
  {
    case planestress :
      b = e/(1.0-nu*nu);
      d[0][0] = b;     d[0][1] = b*nu;  d[0][2] = b*nu;
      d[1][0] = b*nu;  d[1][1] = b;     d[1][2] = b*nu;
      d[2][0] = b*nu;  d[2][1] = b*nu;  d[2][2] = b;
      break;
    case planestrain :
    case spacestress :
      b = e /((1.0 + nu)*(1.0 - 2.0*nu));
      d[0][0] = d[1][1] = d[2][2] = 1.0 - nu;
      d[0][1] = d[0][2] = d[1][0] = d[1][2] = d[2][0] = d[2][1] = nu;
      cmulm(b, d);
      break;
    default :
      fprintf(stderr, "\nUnknown type of stress/strain state is required\n");
      fprintf(stderr, "\n in function  mohrcoulomb::pelmatstiff(), (file %s, line %d)\n", __FILE__, __LINE__);
      break;
  }
}

void mohrcoulomb::nlstresses (long ipp, long ido)
  //
{
  long i,ni,n=Mm->ip[ipp].ncompstr,mu;
  double gamma,err;
  vector epsn(n),epsp(n),q(0);

  //  initial values
  for (i=0; i<n; i++)
  {
    epsn[i]=Mm->ip[ipp].strain[i];
    epsp[i]=Mm->ip[ipp].eqother[ido+i];
  }
  gamma=Mm->ip[ipp].eqother[ido+n];

  //  try single vector stress return algorithm
  if (sra.give_tsra () == cp){
    ni=sra.give_ni ();
    err=sra.give_err ();
    
    mu = cutting_plane (ipp,gamma,epsn,epsp,q,ni,err);
  }
  else{
    fprintf (stderr,"\n\n wrong type of stress return algorithm is required in nlstresses (file %s, line %d).\n",__FILE__,__LINE__);
    abort ();
  }
  
  // check singularity
  if (mu > 0)
  {
    // singularity was found, use multisurface stress return algorithm
    gamma=Mm->ip[ipp].eqother[ido+n];
    
    ni=sra.give_ni ();
    err=sra.give_err ();
    mc_msurf_cp(ipp, gamma, epsn, epsp, q, mu,ni,err);
  }

  //  new data storage
  for (i=0;i<n;i++){
    Mm->ip[ipp].other[ido+i]=epsp[i];
  }
  Mm->ip[ipp].other[ido+n]=gamma;
}

/**
  This function computes stresses at given integration point ipp,
  depending on the reached averaged nonlocal strains.
  The cutting plane algorithm is used. The stress and the other attribute of
  given integration point is actualized.

  @param ipp - integration point number in the mechmat ip array.
  @param ido - index of internal variables for given material in the ipp other array
*/
void mohrcoulomb::nonloc_nlstresses (long ipp,long ido)
{
  long i,ni,n=Mm->ip[ipp].ncompstr,mu;
  double gamma,err;
  vector epsn(n),epsp(n),q(0);

  //  initial values
  for (i=0; i<n; i++)
  {
    epsn[i]=Mm->ip[ipp].strain[i];
    epsp[i]=Mm->ip[ipp].nonloc[i];
  }
  gamma=Mm->ip[ipp].eqother[ido+n];

  //  try single vector stress return algorithm
  if (sra.give_tsra () == cp){
    ni=sra.give_ni ();
    err=sra.give_err ();
    
    mu = cutting_plane (ipp,gamma,epsn,epsp,q,ni,err);
  }
  else{
    fprintf (stderr,"\n\n wrong type of stress return algorithm is required in nlstresses (file %s, line %d).\n",__FILE__,__LINE__);
    abort ();
  }
  
  // check singularity
  if (mu > 0)
    // singularity was found, use multisurface stress return algorithm

    ni=sra.give_ni ();
    err=sra.give_err ();
    mc_msurf_cp(ipp, gamma, epsn, epsp, q, mu,ni,err);

  //  new data storage
  for (i=0;i<n;i++){
    Mm->ip[ipp].other[ido+i]=epsp[i];
  }
  Mm->ip[ipp].other[ido+n]=gamma;

}

/**
   Function returns stresses on sufrace of plasticity
   cutting plane method in principal stresses is used

   gamma,epsp and q will be replaced by new values

   @param ipp   - integration point pointer
   @param gamma - consistency parameter
   @param epsn  - total strain components
   @param epsp  - plastic strain components
   @param q     - hardening parameters
   @param ni    - maximum number of iterations
   @param err   - required error
   
   @retval Function returns indicator of active yield surfaces in case detection of the singularity
           or zero in case of successfull stress return.
   4.8.2001
*/
long mohrcoulomb::cutting_plane(long ipp, double &gamma, vector &epsn, vector &epsp, vector &q, long ni, double err)
{
  long i,j,ncomp=epsn.n,idem,mu,zps;
  double f,denom,dgamma,e,nu;
  vector epsa(ncomp), sig(ncomp), epsptr(ncomp);
  vector psig(3), peps(3), pdepsp(3), pepsa(3), dfds(3), dgds(3), depsp(ncomp);
  matrix d(ncomp,ncomp), sigt(3,3), epsat(3,3), epspt(3,3);
  matrix pd(3,3), pvect(3,3);

  // Initialization
  // checking for the correct type of the elastic material
  idem = Mm->ip[ipp].gemid();
  if (Mm->ip[ipp].tm[idem] != elisomat)
  {
    fprintf(stderr, "\nError - The Mohr-Coulomb material is combined with\n");
    fprintf(stderr, "        the nonsupported elastic material.\n");
    fprintf(stderr, "        Only elisomat type is supported\n");
    abort ();
  }

  dgamma=0.0;
  e  = Mm->eliso[Mm->ip[ipp].idm[idem]].e;
  nu = Mm->eliso[Mm->ip[ipp].idm[idem]].nu;
  // elastic stiffness matrix computation
  Mm->elmatstiff(d, ipp);
  pelmatstiff (ipp, pd, e, nu);
  // computation complementary values for plain stress/strain
  if (Mm->ip[ipp].ssst == planestrain)
  {
     d[0][3] = d[0][1]; d[1][3] = d[1][0];
     d[3][0] = d[1][0]; d[3][1] = d[1][0]; d[3][3] = d[1][1];
  }
  if (Mm->ip[ipp].ssst == planestress)
    epsn[3] = -nu / (1.0 - nu) * (epsn[0]+epsn[1]);
  copyv(epsp, epsptr);

  // Computing principal directions - they do not change during stress return
  // elastic strain
  subv (epsn, epsp, epsa);
  // trial stress computation
  mxv (d, epsa, sig);
  // principal stresses and their directions
  vector_tensor (sig, sigt, Mm->ip[ipp].ssst, stress);
  princ_val (sigt, psig, pvect, nijac, limit, Mp->zero, 3);
  // transformation elastic strain into the computed principal directions
  vector_tensor (epsa, epsat, Mm->ip[ipp].ssst, strain);
  glmatrixtransf(epsat, pvect);
  peps[0] = epsat[0][0];
  peps[1] = epsat[1][1];
  peps[2] = epsat[2][2];
  if (Mm->ip[ipp].ssst == planestress)
  {
    zps = checkzeropsig(psig);
    fillrow(0.0, zps, pd);
    fillcol(0.0, zps, pd);
  }
  //  Main iteration loop
  for (i=0;i<ni;i++){
    //  elastic strain
    subv (peps, pdepsp, pepsa);
    //  trial stress computation
    mxv (pd, pepsa, psig);
    // checking ordering of the reached principal stresses and vertex singularity
    mu = checkpsig(psig);
    if (mu > 0)
      return mu;
    // checking yield function
    f = yieldfunction (psig);
    if (f<err)  break;
    if (i==ni-1 && f>err)
    {
      fprintf (stderr,"\n yield surface was not reached for ipp = %ld", ipp);
      exit(0);
    }
    // compute necessary derivations of yield function and plastic potential function.
    deryieldfsigma (dfds);
    derplaspotsigma (dgds);

    mxv (pd, dgds, psig);
    scprd (dfds, psig, denom);
    denom += plasmodscalar(ipp, q);
    //  new increment of consistency parameter
    dgamma = f/denom;
    //  new internal variables
    gamma+=dgamma;
    // compute new principal plastic strain increments and transform them into global c.s.
    for (j=0; j<3; j++)
    {
      epspt[j][j]+=dgamma*dgds[j];
      pdepsp[j]+=dgamma*dgds[j];
    }
    lgmatrixtransf(epspt, pvect);
    tensor_vector(depsp, epspt, Mm->ip[ipp].ssst, strain);
    addv(epsp, depsp, epsptr);
    // update hardening parameter
    updateq(ipp, epsptr, q);
  }
  // updating plastic strains
  copyv(epsptr, epsp);
  //  elastic strain
  subv (epsn, epsp, epsa);
  //  trial stress computation
  mxv (d, epsa, sig);
  //  storing resulting stresses
  Mm->storestress (0, ipp, sig);
  return 0;
}


/**
   Function returns plastic modulus in case that hardening/softening rule is adopted

   @param ipp   - integration point pointer
   @param q     - hardening parameters

   @retval It returns value of plastic modulus.

   4.11.2003
*/
double mohrcoulomb::plasmodscalar(long ipp, vector &q)