CAM_YS.cppGo to the documentation of this file.00001 00002 //================================================================================ 00003 //# COPYRIGHT (C): :-)) # 00004 //# PROJECT: Object Oriented Finite Element Program # 00005 //# PURPOSE: Cam Clay yield criterion # 00006 //# (Ref. Wood p113) # 00007 //# CLASS: CAMYieldSurface(for monotonic loading) # 00008 //# # 00009 //# VERSION: # 00010 //# LANGUAGE: C++.ver >= 2.0 ( Borland C++ ver=3.00, SUN C++ ver=2.1 ) # 00011 //# TARGET OS: DOS || UNIX || . . . # 00012 //# PROGRAMMER(S): Boris Jeremic, ZHaohui Yang # 00013 //# # 00014 //# # 00015 //# DATE: Mar. 28, 2001 # 00016 //# UPDATE HISTORY: # 00017 //# # 00018 //# # 00019 //# # 00020 //# # 00021 //# # 00022 //================================================================================ 00023 00024 00025 #ifndef CAM_YS_CPP 00026 #define CAM_YS_CPP 00027 00028 #include "CAM_YS.h" 00029 #include <basics.h> 00030 00031 00032 //================================================================================ 00033 // Normal constructor 00034 //================================================================================ 00035 00036 CAMYieldSurface::CAMYieldSurface(double Mp ) 00037 { 00038 M = Mp; 00039 } 00040 00041 00042 //================================================================================ 00043 //create a colne of itself 00044 //================================================================================ 00045 00046 YieldSurface * CAMYieldSurface::newObj() { 00047 00048 YieldSurface *new_YS = new CAMYieldSurface(M); 00049 return new_YS; 00050 00051 } 00052 00053 //================================================================================ 00054 // Copy constrstructor 00055 //================================================================================ 00056 // 00057 //CAMYieldSurface::CAMYieldSurface(CAMYieldSurface &CAMYS ) { } 00058 // 00059 00060 00061 //================================================================================ 00062 // Yield criterion evaluation function F(EPState) 00063 //================================================================================ 00064 00065 double CAMYieldSurface::f(const EPState *EPS) const 00066 { 00067 double p =EPS->getStress().p_hydrostatic(); 00068 double q =EPS->getStress().q_deviatoric(); 00069 //cout << "p " << p; 00070 00071 double po = EPS->getScalarVar(1); 00072 double temp3 = q*q - M*M*p*(po - p); 00073 00074 //printf("\n========Inside f = %.4f \n ", temp3); 00075 00076 return temp3; 00077 } 00078 00079 00080 //================================================================================ 00081 // tensor dF/dsigma_ij 00082 //================================================================================ 00083 00084 tensor CAMYieldSurface::dFods(const EPState *EPS) const { 00085 00086 tensor dFoverds( 2, def_dim_2, 0.0); 00087 tensor I2("I", 2, def_dim_2); 00088 00089 double p = EPS->getStress().p_hydrostatic(); 00090 double q = EPS->getStress().q_deviatoric(); 00091 double po = EPS->getScalarVar( 1 ); 00092 tensor DpoDs = EPS->getStress().dpoverds(); 00093 tensor DqoDs = EPS->getStress().dqoverds(); 00094 00095 double dFoverdp = -1.0*M*M*( po - 2.0*p ); 00096 double dFoverdq = 2.0*q; 00097 00098 dFoverds = DpoDs * dFoverdp + 00099 DqoDs * dFoverdq; 00100 00101 return dFoverds; 00102 00103 } 00104 00105 00106 //================================================================================ 00107 // double xi_s1 = dF/dpo = -M*M*p Derivative in terms of first scalar var 00108 //================================================================================ 00109 00110 double CAMYieldSurface::xi_s1( const EPState *EPS ) const { 00111 00112 //double p = EPS->getStress().Iinvariant1()/3.0; 00113 double p = EPS->getStress().p_hydrostatic(); 00114 return -1.0 * M * M * p; 00115 00116 } 00117 00118 00119 // moved to stresstensor Boris Jeremic@ucdavis.edu 21Aug2001 00120 // // I took these functions from mmodel.cpp, programmed by Dr. Jeremic 00121 // //############################################################################# 00122 // tensor CAMYieldSurface::dpoverds( ) const 00123 // { 00124 // tensor ret(2, def_dim_2, 0.0); 00125 // tensor I2("I", 2, def_dim_2); 00126 // ret = I2*(-1.0/3.0); 00127 // ret.null_indices(); 00128 // return ret; 00129 // } 00130 // 00131 // //############################################################################# 00132 // tensor CAMYieldSurface::dqoverds(const EPState *EPS) const 00133 // { 00134 // 00135 // stresstensor stress = EPS->getStress(); 00136 // 00137 // tensor ret(2, def_dim_2, 0.0); 00138 // double q = stress.q_deviatoric(); 00139 // stresstensor s( 0.0); 00140 // //... double J2D = stress.Jinvariant2(); 00141 // //... double temp1 = sqrt(J2D); 00142 // //... double temp2 = sqrt(3.0)/(2.0*temp1); 00143 // double temp2 = (3.0/2.0)*(1/q); 00144 // s = stress.deviator(); 00145 // ret = s*temp2; 00146 // ret.null_indices(); 00147 // return ret; 00148 // } 00149 // 00150 // //############################################################################# 00151 // tensor CAMYieldSurface::dthetaoverds(const EPState *EPS) const 00152 // { 00153 // stresstensor stress = EPS->getStress(); 00154 // 00155 // tensor ret(2, def_dim_2, 0.0); 00156 // stresstensor s( 0.0); 00157 // stresstensor t( 0.0); 00158 // tensor I2("I", 2, def_dim_2); 00159 // 00160 // // double EPS = pow(d_macheps(),(1./3.)); 00161 // double J2D = stress.Jinvariant2(); 00162 // double q = stress.q_deviatoric(); 00163 // double theta = stress.theta(); 00164 // 00165 // //out while ( theta >= 2.0*PI ) 00166 // //out theta = theta - 2.0*PI; // if bigger than full cycle 00167 // //out while ( theta >= 4.0*PI/3.0 ) 00168 // //out theta = theta - 4.0*PI/3.0; // if bigger than four thirds of half cycle 00169 // //out while ( theta >= 2.0*PI/3.0 ) 00170 // //out theta = theta - 2.0*PI/3.0; // if bigger than two third of half cycle 00171 // //out while ( theta >= PI/3.0 ) 00172 // //out theta = 2.0*PI/3.0 - theta; // if bigger than one third of half cycle 00173 // //out 00174 // //out if ( theta < 0.0001 ) 00175 // //out { 00176 // //out ::printf("theta = %.10e in CAMYieldSurface::dthetaoverds(stress)\n", 00177 // //out theta); 00178 // //out//>><<>><<>><< return ret; 00179 // //out } 00180 // 00181 // double c3t = cos(3.0*theta); 00182 // double s3t = sin(3.0*theta); 00183 // 00184 // double tempS = (3.0/2.0)*(c3t/(q*q*s3t)); 00185 // double tempT = (9.0/2.0)*(1.0/(q*q*q*s3t)); 00186 // 00187 // s = stress.deviator(); 00188 // t = s("qk")*s("kp") - I2*(J2D*(2.0/3.0)); 00189 // 00190 // s.null_indices(); 00191 // t.null_indices(); 00192 // 00193 // ret = s*tempS - t*tempT; 00194 // ret.null_indices(); 00195 // return ret; 00196 // } 00197 00198 //================================================================================ 00199 double CAMYieldSurface::getM() const 00200 { 00201 return M; 00202 } 00203 00204 00205 OPS_Stream& operator<< (OPS_Stream& os, const CAMYieldSurface & YS) 00206 { 00207 os << "Cam Clay Yield Surface Parameters: " << endln; 00208 os << "M = " << YS.M << endln; 00209 return os; 00210 } 00211 00212 #endif 00213 |