OpenSees » SRC » material » nD » Template3Dep

DP_PS.cpp

Go to the documentation of this file.
00001 
00002 //################################################################################
00003 //# COPYRIGHT (C):     :-))                                                      #
00004 //# PROJECT:           Object Oriented Finite Element Program                    #
00005 //# PURPOSE:           Drucker - Prager  yield criterion                         #
00006 //# CLASS:             DPPotentialSurface                                        #
00007 //#                                                                              #
00008 //# VERSION:                                                                     #
00009 //# LANGUAGE:          C++.ver >= 2.0 ( Borland C++ ver=3.00, SUN C++ ver=2.1 )  #
00010 //# TARGET OS:         DOS || UNIX || . . .                                      #
00011 //# PROGRAMMER(S):     Boris Jeremic, ZHaohui Yang                               #
00012 //#                                                                              #
00013 //#                                                                              #
00014 //# DATE:              August 03 '93                                             #
00015 //# UPDATE HISTORY:    August 08 '00                                             #
00016 //# UPDATE HISTORY:    20Aug2004 ZC added kinematic hardening part               #
00017 //#                                                                              #
00018 //#                                                                              #
00019 //#                                                                              #
00020 //#                                                                              #
00021 //#                                                                              #
00022 //################################################################################
00023 //*/
00024 
00025 #ifndef DP_PS_CPP
00026 #define DP_PS_CPP
00027 
00028 #include "DP_PS.h"
00029 
00030 
00031 
00032 //================================================================================
00033 // Normal constructor
00034 //================================================================================
00035 
00036 //DPPotentialSurface::DPPotentialSurface( double a2d = 0.3) : alfa2(a2d) 
00037 //{
00038 //
00039 //} 
00040 
00041 
00042 //================================================================================
00043 // Copy constrstructor  === not necessary for no pointer member
00044 //================================================================================
00045 
00046 DPPotentialSurface::DPPotentialSurface(const DPPotentialSurface &DPPS ) {
00047 
00048     alfa2 =  DPPS.getalfa2();
00049 
00050 }
00051 
00052 
00053 //================================================================================
00054 //create a colne of itself
00055 //================================================================================
00056 
00057 PotentialSurface * DPPotentialSurface::newObj() {  
00058 
00059      PotentialSurface  *new_PS = new DPPotentialSurface( this->getalfa2() );
00060      return new_PS;
00061 
00062 }
00063 
00064 
00065 
00066 //================================================================================
00067 // tensor dQ/dsigma_ij
00068 //================================================================================
00069 
00070 tensor DPPotentialSurface::dQods(const EPState *EPS) const {
00071     stresstensor dQods;
00072     
00073         tensor KroneckerI("I", 2, def_dim_2);
00074     //double temp1 =  EPS->getScalarVar(1);    
00075     double temp1 = getalfa2();    
00076     
00077     stresstensor alpha;
00078     stresstensor s_bar;
00079     stresstensor Tnsr0;
00080     double temp6 = 0.0;
00081     stresstensor temp5;
00082     stresstensor sigma = EPS->getStress();   
00083     double p = sigma.p_hydrostatic();
00084     stresstensor sdev = sigma.deviator();
00085     double halfRt2 = 0.5 * sqrt(2.0);
00086     int nod = EPS->getNTensorVar();
00087     if ( nod >=1 )  { //May not have kinematic hardening
00088       alpha = EPS->getTensorVar(1);  
00089       s_bar = sdev - (alpha*p);
00090       temp5 = alpha("ij") * s_bar("ij");
00091       temp5.null_indices();
00092       temp6 = temp5.trace();
00093       Tnsr0 = KroneckerI*(temp6/3.0);
00094     }
00095     else {
00096           s_bar = sdev;
00097     }
00098     stresstensor temp3 = s_bar("ij") * s_bar("ij");
00099     temp3.null_indices();
00100     double temp4 = temp3.trace();
00101     temp4 = sqrt(temp4);
00102     Tnsr0 += s_bar;
00103     double eps = pow( d_macheps(), 0.5 );
00104     if ( fabs(temp4) > eps )  {
00105       Tnsr0 = Tnsr0 * (halfRt2/temp4);
00106     }
00107             
00108     dQods = (KroneckerI * temp1) + Tnsr0;; 
00109         
00110     return dQods;
00111 }
00112 
00113 
00114 //================================================================================
00115 // tensor d2Q/dsigma_ij_2
00116 //================================================================================
00117 
00118 tensor DPPotentialSurface::d2Qods2(const EPState *EPS) const {
00119     tensor d2Qods2(4, def_dim_4, 0.0);
00120     
00121         tensor KroneckerI("I", 2, def_dim_2);
00122         tensor T1 = KroneckerI("ij")*KroneckerI("mn");
00123         T1.null_indices();
00124         tensor T2 = (T1.transpose0110()+T1.transpose0111())*0.5;
00125         tensor T3 = T2 - T1*(1.0/3.0);
00126         
00127     //double temp1 =  EPS->getScalarVar(1);    
00128     //double temp1 = getalfa2();
00129     
00130     tensor T4;
00131     stresstensor alpha;
00132     stresstensor s_bar;
00133     tensor temp9(4, def_dim_4, 0.0);
00134     stresstensor sigma = EPS->getStress();
00135     double p = sigma.p_hydrostatic();
00136     stresstensor sdev = sigma.deviator();
00137     double halfRt2 = 0.5 * sqrt(2.0);
00138     int nod = EPS->getNTensorVar();
00139     if ( nod >=1 )  { //May not have kinematic hardening
00140       alpha = EPS->getTensorVar(1);
00141       temp9 = KroneckerI("ij") * alpha("mn");
00142       temp9.null_indices();
00143       T4 = T2 + temp9*(1.0/3.0);
00144       s_bar = sdev - (alpha*p);
00145     }
00146     else {
00147           s_bar = sdev;
00148           T4 = T2;
00149     }
00150     T4 = T2 - temp9;
00151     tensor temp3 = s_bar("ij") * s_bar("ij");
00152     temp3.null_indices();  
00153     double temp4 = temp3.trace();
00154     temp4 = sqrt(temp4);
00155     tensor temp5 = s_bar("ij")*s_bar("mn");
00156     double eps = pow( d_macheps(), 0.5 );
00157     if ( fabs(temp4) > eps )  {
00158       d2Qods2 = T3 * (halfRt2/temp4) - temp5*(halfRt2/(temp4*temp4*temp4));
00159       d2Qods2 = T4("ijkl") * d2Qods2("klmn");
00160       d2Qods2.null_indices();
00161     }   
00162     
00163     return d2Qods2;
00164 }
00165 
00166 // For Consistent Algorithm, Z Cheng, Jan 2004
00167 tensor DPPotentialSurface::d2Qodsds1(const EPState *EPS) const 
00168 {  
00169   tensor I("I", 2, def_dim_2);
00170   tensor d2Qoverdsds1 = I;
00171   return d2Qoverdsds1;
00172 }
00173 
00174 //================================================================================
00175 
00176 double DPPotentialSurface::getalfa2() const {
00177 
00178     return alfa2; 
00179 
00180 }
00181 
00182 //================================================================================
00183 OPS_Stream& operator<< (OPS_Stream& os, const DPPotentialSurface &PS)
00184 {
00185    os << "Drucker-Prager Potential Surface Parameters: " << endln;
00186    os << "alfa2 = " << PS.getalfa2() << endln;
00187    return os;
00188 }
00189 
00190 
00191 #endif

Generated on Mon Oct 23 15:05:16 2006 for OpenSees by doxygen 1.5.0