Main Page   Class Hierarchy   Alphabetical List   Compound List   File List   Compound Members   File Members  

J2PlaneStress.cpp

Go to the documentation of this file.
00001 /* ****************************************************************** **
00002 **    OpenSees - Open System for Earthquake Engineering Simulation    **
00003 **          Pacific Earthquake Engineering Research Center            **
00004 **                                                                    **
00005 **                                                                    **
00006 ** (C) Copyright 1999, The Regents of the University of California    **
00007 ** All Rights Reserved.                                               **
00008 **                                                                    **
00009 ** Commercial use of this program without express permission of the   **
00010 ** University of California, Berkeley, is strictly prohibited.  See   **
00011 ** file 'COPYRIGHT'  in main directory for information on usage and   **
00012 ** redistribution,  and for a DISCLAIMER OF ALL WARRANTIES.           **
00013 **                                                                    **
00014 ** ****************************************************************** */
00015                                                                         
00016 // $Revision: 1.2 $
00017 // $Date: 2001/01/23 08:46:28 $
00018 // $Source: /usr/local/cvs/OpenSees/SRC/material/nD/J2PlaneStress.cpp,v $
00019 
00020 // Written: Ed "C++" Love
00021 //
00022 // J2PlaneStress isotropic hardening material class
00023 // 
00024 //  Elastic Model
00025 //  sigma = K*trace(epsilion_elastic) + (2*G)*dev(epsilon_elastic)
00026 //
00027 //  Yield Function
00028 //  phi(sigma,q) = || dev(sigma) ||  - sqrt(2/3)*q(xi) 
00029 //
00030 //  Saturation Isotropic Hardening with linear term
00031 //  q(xi) = simga_infty + (sigma_0 - sigma_infty)*exp(-delta*xi) + H*xi 
00032 //
00033 //  Flow Rules
00034 //  \dot{epsilon_p} =  gamma * d_phi/d_sigma
00035 //  \dot{xi}        = -gamma * d_phi/d_q 
00036 //
00037 //  Linear Viscosity 
00038 //  gamma = phi / eta  ( if phi > 0 ) 
00039 //
00040 //  Backward Euler Integration Routine 
00041 //  Yield condition enforced at time n+1 
00042 //
00043 //  Send strains in following format :
00044 // 
00045 //     strain_vec = {   eps_00
00046 //                      eps_11
00047 //                    2 eps_01   }   <--- note the 2
00048 // 
00049 //  set eta := 0 for rate independent case
00050 //
00051 
00052 #include <J2PlaneStress.h>
00053 
00054 Vector J2PlaneStress :: strain_vec(3) ;
00055 Vector J2PlaneStress :: stress_vec(3) ;
00056 Matrix J2PlaneStress :: tangent_matrix(3,3) ;
00057 
00058 //null constructor
00059 J2PlaneStress ::  J2PlaneStress( ) : 
00060 J2Plasticity( ) 
00061 {  }
00062 
00063 
00064 //full constructor
00065 
00066 J2PlaneStress :: J2PlaneStress(   int    tag, 
00067                  double K,
00068                  double G,
00069                  double yield0,
00070                  double yield_infty,
00071                  double d,
00072                  double H,
00073                  double viscosity ) : 
00074 J2Plasticity(tag, ND_TAG_J2PlaneStress, 
00075              K, G, yield0, yield_infty, d, H, viscosity )
00076 { }
00077 
00078 
00079 //elastic constructor
00080 
00081 J2PlaneStress :: J2PlaneStress(   int    tag, 
00082                  double K, 
00083                  double G ) :
00084 J2Plasticity(tag, ND_TAG_J2PlaneStress, K, G )
00085 { }
00086 
00087 
00088 //destructor
00089 J2PlaneStress :: ~J2PlaneStress( ) 
00090 {  } 
00091 
00092 
00093 //make a clone of this material
00094 NDMaterial* J2PlaneStress :: getCopy( ) 
00095 { 
00096   J2PlaneStress  *clone;
00097   clone = new J2PlaneStress( ) ;   //new instance of this class
00098   *clone = *this ;          //asignment to make copy
00099   return clone ;
00100 }
00101 
00102 
00103 //send back type of material
00104 const char* J2PlaneStress :: getType( ) const 
00105 {
00106   return "PlaneStress2D" ;
00107 }
00108 
00109 
00110 //send back order of strain in vector form
00111 int J2PlaneStress :: getOrder( ) const 
00112 { 
00113   return 3 ; 
00114 } 
00115 
00116 //get the strain and integrate plasticity equations
00117 int J2PlaneStress :: setTrialStrain( const Vector &strain_from_element ) 
00118 {
00119   const double tolerance = 1e-12 ;
00120 
00121   const int max_iterations = 25 ;
00122   int iteration_counter  = 0 ;
00123 
00124   int i, j, k, l ;
00125   int ii, jj ;
00126 
00127   double eps22  =  strain(2,2) ;
00128   strain.Zero( ) ;
00129 
00130   strain(0,0) =        strain_from_element(0) ;
00131   strain(1,1) =        strain_from_element(1) ;
00132   strain(0,1) = 0.50 * strain_from_element(2) ;
00133   strain(1,0) =        strain(0,1) ;
00134 
00135   strain(2,2) =        eps22 ; 
00136 
00137   //enforce the plane stress condition sigma_22 = 0 
00138   //solve for epsilon_22 
00139   iteration_counter = 0 ;  
00140   do {
00141 
00142      this->plastic_integrator( ) ;
00143     
00144      strain(2,2) -= stress(2,2) / tangent[2][2][2][2] ;
00145 
00146      iteration_counter++ ;
00147      if ( iteration_counter > max_iterations ) {
00148        cerr << "More than " << max_iterations ;
00149        cerr << " iterations in setTrialStrain of J2PlaneStress \n" ;
00150        break ;
00151      }// end if 
00152 
00153   } while ( fabs(stress(2,2)) > tolerance ) ;
00154 
00155   //modify tangent for plane stress 
00156   for ( ii = 0; ii < 3; ii++ ) {
00157     for ( jj = 0; jj < 3; jj++ )  {
00158 
00159           index_map( ii, i, j ) ;
00160           index_map( jj, k, l ) ;
00161 
00162           tangent[i][j][k][l] -=   tangent[i][j][2][2] 
00163                                  * tangent[2][2][k][l] 
00164                                  / tangent[2][2][2][2] ;
00165 
00166           //minor symmetries 
00167           tangent [j][i][k][l] = tangent[i][j][k][l] ;
00168           tangent [i][j][l][k] = tangent[i][j][k][l] ;
00169           tangent [j][i][l][k] = tangent[i][j][k][l] ;
00170 
00171     } // end for jj
00172   } // end for ii 
00173 
00174   return 0 ;
00175 }
00176 
00177 
00178 //unused trial strain functions
00179 int J2PlaneStress :: setTrialStrain( const Vector &v, const Vector &r )
00180 { 
00181    return this->setTrialStrain( v ) ;
00182 } 
00183 
00184 int J2PlaneStress :: setTrialStrainIncr( const Vector &v ) 
00185 {
00186     return -1 ;
00187 }
00188 
00189 int J2PlaneStress :: setTrialStrainIncr( const Vector &v, const Vector &r ) 
00190 {
00191     return -1 ;
00192 }
00193 
00194 
00195 
00196 //send back the strain
00197 const Vector& J2PlaneStress :: getStrain( ) 
00198 {
00199   strain_vec(0) =       strain(0,0) ;
00200   strain_vec(1) =       strain(1,1) ;
00201   strain_vec(2) = 2.0 * strain(0,1) ;
00202 
00203   return strain_vec ;
00204 } 
00205 
00206 
00207 //send back the stress 
00208 const Vector& J2PlaneStress :: getStress( ) 
00209 {
00210   stress_vec(0) = stress(0,0) ;
00211   stress_vec(1) = stress(1,1) ;
00212   stress_vec(2) = stress(0,1) ;
00213 
00214   return stress_vec ;
00215 }
00216 
00217 //send back the tangent 
00218 const Matrix& J2PlaneStress :: getTangent( ) 
00219 {
00220   // matrix to tensor mapping
00221   //  Matrix      Tensor
00222   // -------     -------
00223   //   0           0 0
00224   //   1           1 1
00225   //   2           0 1  ( or 1 0 ) 
00226   // 
00227        
00228   tangent_matrix(0,0) = tangent [0][0] [0][0] ;
00229   tangent_matrix(1,1) = tangent [1][1] [1][1] ;
00230   tangent_matrix(2,2) = tangent [0][1] [0][1] ;
00231 
00232   tangent_matrix(0,1) = tangent [0][0] [1][1] ;
00233   tangent_matrix(1,0) = tangent [1][1] [0][0] ;
00234 
00235   tangent_matrix(0,2) = tangent [0][0] [0][1] ;
00236   tangent_matrix(2,0) = tangent [0][1] [0][0] ;
00237 
00238   tangent_matrix(1,2) = tangent [1][1] [0][1] ;
00239   tangent_matrix(2,1) = tangent [0][1] [1][1] ;
00240 
00241   return tangent_matrix ;
00242 } 
00243 
00244 //this is mike's problem
00245 int J2PlaneStress :: setTrialStrain(const Tensor &v) 
00246 {
00247   return -1 ;
00248 }
00249 
00250 int J2PlaneStress :: setTrialStrain(const Tensor &v, const Tensor &r)     
00251 {
00252   return -1 ;
00253 }
00254 
00255 int J2PlaneStress :: setTrialStrainIncr(const Tensor &v) 
00256 {
00257   return -1 ;
00258 }
00259 
00260 int J2PlaneStress :: setTrialStrainIncr(const Tensor &v, const Tensor &r) 
00261 {
00262   return -1 ;
00263 }
00264 
00265 const Tensor& J2PlaneStress :: getTangentTensor( ) 
00266 {
00267   return rank4 ;
00268 }
00269 
00270 //jeremic@ucdavis.edu 22jan2001const Tensor& J2PlaneStress :: getStressTensor( ) 
00271 //jeremic@ucdavis.edu 22jan2001{
00272 //jeremic@ucdavis.edu 22jan2001  return rank2 ;
00273 //jeremic@ucdavis.edu 22jan2001}
00274 //jeremic@ucdavis.edu 22jan2001
00275 //jeremic@ucdavis.edu 22jan2001const Tensor& J2PlaneStress :: getStrainTensor( ) 
00276 //jeremic@ucdavis.edu 22jan2001{
00277 //jeremic@ucdavis.edu 22jan2001  return rank2 ;
00278 //jeremic@ucdavis.edu 22jan2001}
00279 
00280 
00281 //this is frank's problem
00282 int J2PlaneStress :: sendSelf(int commitTag, Channel &theChannel)
00283 {
00284   return -1 ;
00285 }
00286 
00287 int J2PlaneStress :: recvSelf(int commitTag, Channel &theChannel, 
00288                        FEM_ObjectBroker &theBroker)
00289 {
00290   return -1 ;
00291 }
00292 
00293 
00294 //matrix_index ---> tensor indices i,j
00295 // plane stress different because of condensation on tangent
00296 // case 3 switched to 1-2 and case 4 to 3-3 
00297 void J2PlaneStress :: index_map( int matrix_index, int &i, int &j )
00298 {
00299   switch ( matrix_index+1 ) { //add 1 for standard tensor indices
00300 
00301     case 1 :
00302       i = 1 ; 
00303       j = 1 ;
00304       break ;
00305  
00306     case 2 :
00307       i = 2 ;
00308       j = 2 ; 
00309       break ;
00310 
00311     case 3 :
00312       i = 1 ;
00313       j = 2 ;
00314       break ;
00315 
00316     case 4 :
00317       i = 3 ;
00318       j = 3 ;
00319       break ;
00320 
00321     case 5 :
00322       i = 2 ;
00323       j = 3 ;
00324       break ;
00325 
00326     case 6 :
00327       i = 3 ;
00328       j = 1 ;
00329       break ;
00330 
00331 
00332     default :
00333       i = 1 ;
00334       j = 1 ;
00335       break ;
00336 
00337   } //end switch
00338 
00339 i-- ; //subtract 1 for C-indexing
00340 j-- ;
00341 
00342 return ; 
00343 }
00344 
Copyright Contact Us