DP_YS01.cppGo to the documentation of this file.00001 00002 //================================================================================ 00003 //# COPYRIGHT (C): :-)) # 00004 //# PROJECT: Object Oriented Finite Element Program # 00005 //# PURPOSE: Drucker-Prager yield criterion 01 (with Pc) # 00006 //# (Ref. Geotechnique # 00007 //# V.47 No.2 255-272, 1997) # 00008 //# CLASS: DPYieldSurface01 # 00009 //# # 00010 //# VERSION: # 00011 //# LANGUAGE: C++.ver >= 2.0 ( Borland C++ ver=3.00, SUN C++ ver=2.1 ) # 00012 //# TARGET OS: DOS || UNIX || . . . # 00013 //# PROGRAMMER(S): Boris Jeremic, ZHaohui Yang # 00014 //# # 00015 //# # 00016 //# DATE: August 03 '00 # 00017 //# UPDATE HISTORY: December 13, '00 # 00018 //# # 00019 //# # 00020 //# # 00021 //# # 00022 //# # 00023 //================================================================================ 00024 00025 00026 #ifndef DP_YS01_CPP 00027 #define DP_YS01_CPP 00028 00029 #include "DP_YS01.h" 00030 #include <basics.h> 00031 00032 //================================================================================ 00033 // Normal constructor 00034 //================================================================================ 00035 00036 DPYieldSurface01::DPYieldSurface01(double pc) 00037 { 00038 Pc = pc; 00039 } 00040 00041 00042 //================================================================================ 00043 //create a colne of itself 00044 //================================================================================ 00045 00046 YieldSurface * DPYieldSurface01::newObj() { 00047 00048 YieldSurface *new_YS = new DPYieldSurface01( Pc ); 00049 return new_YS; 00050 00051 } 00052 00053 //================================================================================ 00054 // Copy constrstructor 00055 //================================================================================ 00056 // 00057 //DPYieldSurface01::DPYieldSurface01(DPYieldSurface01 &MDYS ) { } 00058 // 00059 00060 00061 //================================================================================ 00062 // Yield criterion evaluation function F(EPState) 00063 //================================================================================ 00064 00065 double DPYieldSurface01::f(const EPState *EPS) const 00066 { 00067 stresstensor S = EPS->getStress().deviator(); 00068 //cout << "s " << S; 00069 //S.print("sigma", "S"); 00070 00071 double p =EPS->getStress().p_hydrostatic(); 00072 p = p - Pc; 00073 opserr << "p " << p; 00074 00075 stresstensor alpha = EPS->getTensorVar( 1 ); 00076 //alpha.print("alpha", " alpha"); 00077 //opserr << "alpha " << alpha; 00078 00079 double m = EPS->getScalarVar(1); 00080 00081 //stresstensor temp1 = S - p*alpha; 00082 stresstensor temp1 = S - p*alpha; // S +p * alpha 00083 temp1.null_indices(); 00084 //temp1.print("temp1 ", " temp1"); 00085 00086 stresstensor temp2 = temp1("ij") * temp1("ij"); 00087 00088 double temp3 = sqrt( temp2.trace() ); 00089 00090 temp3 = temp3 - sqrt(2.0/3.0) * m * p; 00091 //printf("\n========Inside f temp3 = %.4f x = %.4f\n ", temp3, x); 00092 00093 return temp3; 00094 } 00095 00096 00097 //================================================================================ 00098 // tensor dF/dsigma_ij 00099 //================================================================================ 00100 00101 tensor DPYieldSurface01::dFods(const EPState *EPS) const { 00102 00103 tensor dFoverds( 2, def_dim_2, 0.0); 00104 tensor I2("I", 2, def_dim_2); 00105 00106 stresstensor S = EPS->getStress().deviator(); 00107 //S.reportshort("S"); 00108 00109 double p = EPS->getStress().p_hydrostatic(); 00110 p = p - Pc; 00111 //printf("Here we go! p %f\n", p); 00112 00113 stresstensor alpha = EPS->getTensorVar( 1 ); // getting alpha_ij from EPState 00114 //alpha.reportshort("alpha"); 00115 00116 //stresstensor n = EPS->getTensorVar( 3 ); // getting n_ij from EPState 00117 //n.reportshort("n"); 00118 00119 //------------------------------------------------- 00120 // might be moved to Evolution Law 00121 stresstensor r = S * (1.0 / p); 00122 //r.reportshort("r"); 00123 stresstensor r_bar = r - alpha; 00124 //r_bar.reportshort("r_bar"); 00125 stresstensor norm2 = r_bar("ij") * r_bar("ij"); 00126 double norm = sqrt( norm2.trace() ); 00127 00128 //opserr << "d_macheps " << d_macheps() << endlnn; 00129 00130 stresstensor n; 00131 if ( norm >= d_macheps() ){ 00132 n = r_bar*(1.0 / norm ); 00133 } 00134 else { 00135 opserr << "DPYieldSurface01::dFods |n_ij| = 0, divide by zero! Program exits.\n"; 00136 exit(-1); 00137 } 00138 //EPS->setTensorVar( 3, n); //update n_ij// 00139 //------------------------------------------------- 00140 00141 00142 double m = EPS->getScalarVar( 1 ); 00143 00144 00145 //tensorial multiplication produces 1st-order tensor 00146 //tensor temp = n("ij") * n("ij"); 00147 //double temp1 = temp.trace(); 00148 //printf("==== n_ij*n_ij %e\n", temp1); 00149 00151 tensor temp = n("ij") * alpha("ij"); 00152 double N = temp.trace() + sqrt(2.0/3.0)*m; 00153 //printf(" N = %e\n", N); 00154 00155 dFoverds = n - N *I2 *(1.0/3.0); 00156 00157 return dFoverds; 00158 00159 } 00160 00161 00162 //================================================================================ 00163 // double xi_s1 = dF/dm = -(2/3)^0.5 p Derivative in terms of first scalar var 00164 //================================================================================ 00165 00166 double DPYieldSurface01::xi_s1( const EPState *EPS ) const { 00167 00168 double p = EPS->getStress().p_hydrostatic(); 00169 //double p = EPS->getStress().Iinvariant1()/3.0; 00170 p = p - Pc; 00171 return -1.0*sqrt(2.0/3.0) * p; 00172 00173 } 00174 00175 //================================================================================ 00176 // double xi_t1 = dF/dt1 = dF/dalpha= -p*n Derivative in terms of first tensorial var 00177 //================================================================================ 00178 00179 tensor DPYieldSurface01::xi_t1( const EPState *EPS) const { 00180 tensor dFoverds( 2, def_dim_2, 0.0); 00181 tensor I2("I", 2, def_dim_2); 00182 00183 stresstensor S = EPS->getStress().deviator(); 00184 00185 double p = EPS->getStress().p_hydrostatic(); 00186 p = p - Pc; 00187 00188 stresstensor alpha = EPS->getTensorVar( 1 ); // getting alpha_ij from EPState 00189 00190 stresstensor r = S * (1.0 / p); //for p = sig_kk/3 00191 stresstensor r_bar = r - alpha; 00192 stresstensor norm2 = r_bar("ij") * r_bar("ij"); 00193 double norm = sqrt( norm2.trace() ); 00194 00195 stresstensor n; 00196 if ( norm >= d_macheps() ){ 00197 n = r_bar *(1.0 / norm ); 00198 } 00199 else { 00200 opserr << "DPYieldSurface01::dFods |n_ij| = 0, divide by zero! Program exits.\n"; 00201 exit(-1); 00202 } 00203 00204 return (-1.0) * n * p; 00205 } 00206 00207 00208 OPS_Stream& operator<< (OPS_Stream& os, const DPYieldSurface01 & YS) 00209 { 00210 os << "Drucker-Prager Yield Surface 01 Parameters: " << endln; 00211 os << "Pc = " << YS.Pc << endln; 00212 return os; 00213 } 00214 00215 #endif 00216
Generated on Mon Oct 23 15:05:16 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents