DP_YS.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 # 00006 //# CLASS: DPYieldSurface # 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 08 '00 # 00015 //# UPDATE HISTORY: 20Aug2004 ZC added kinematic hardening part # 00016 //# # 00017 //# # 00018 //# # 00019 //# # 00020 //# # 00021 //================================================================================ 00022 //*/ 00023 00024 #ifndef DP_YS_CPP 00025 #define DP_YS_CPP 00026 00027 #include "DP_YS.h" 00028 00029 00030 //================================================================================ 00031 // Normal constructor 00032 //================================================================================ 00033 00034 //DPYieldSurface::DPYieldSurface( double a1d = 0.5) : alfa1(a1d) { } 00035 00036 //================================================================================ 00037 // Normal constructor 00038 //================================================================================ 00039 00040 DPYieldSurface::~DPYieldSurface( ) 00041 { 00042 } 00043 00044 //================================================================================ 00045 // Copy constrstructor 00046 //================================================================================ 00047 // 00048 //DPYieldSurface::DPYieldSurface(const DPYieldSurface &DPYS ) { 00049 // 00050 //} 00051 00052 //================================================================================ 00053 //create a colne of itself 00054 //================================================================================ 00055 00056 //alpha machines complains on this 00057 //DPYieldSurface * DPYieldSurface::newObj() { 00058 00059 YieldSurface * DPYieldSurface::newObj() { 00060 00061 YieldSurface *new_YS = new DPYieldSurface(); 00062 return new_YS; 00063 00064 } 00065 00066 //================================================================================ 00067 // Yield criterion evaluation function F(EPState) 00068 // f = a*I1 + 0.5*sqrt(2.0)*||s_ij - p*alpha_ij|| - k =0 00069 //================================================================================ 00070 00071 double DPYieldSurface::f(const EPState *EPS) const { 00072 double temp1 = EPS->getStress().Iinvariant1(); 00073 double temp2 = temp1 * EPS->getScalarVar(1); 00074 00075 stresstensor alpha; 00076 stresstensor s_bar; 00077 stresstensor sigma = EPS->getStress(); 00078 double p = sigma.p_hydrostatic(); 00079 stresstensor sdev = sigma.deviator(); 00080 double halfRt2 = 0.5 * sqrt(2.0); 00081 int nod = EPS->getNTensorVar(); 00082 if ( nod >=1 ) { //May not have kinematic hardening 00083 alpha = EPS->getTensorVar(1); 00084 s_bar = sdev - (alpha*p); 00085 } 00086 else { 00087 s_bar = sdev; 00088 } 00089 stresstensor temp3 = s_bar("ij") * s_bar("ij"); 00090 temp3.null_indices(); 00091 double temp4 = temp3.trace(); 00092 temp4 = sqrt(temp4); 00093 double temp5 = halfRt2 * temp4; 00094 00095 double k = EPS->getScalarVar(2); 00096 00097 double f = temp2 + temp5 - k; 00098 00099 return f; 00100 } 00101 00102 //================================================================================ 00103 // tensor dF/dsigma_ij 00104 //================================================================================ 00105 00106 tensor DPYieldSurface::dFods(const EPState *EPS) const { 00107 stresstensor dFods; 00108 00109 tensor KroneckerI("I", 2, def_dim_2); 00110 double temp1 = EPS->getScalarVar(1); 00111 00112 stresstensor alpha; 00113 stresstensor s_bar; 00114 stresstensor Tnsr0; 00115 double temp6 = 0.0; 00116 stresstensor temp5; 00117 stresstensor sigma = EPS->getStress(); 00118 double p = sigma.p_hydrostatic(); 00119 stresstensor sdev = sigma.deviator(); 00120 double halfRt2 = 0.5 * sqrt(2.0); 00121 int nod = EPS->getNTensorVar(); 00122 if ( nod >=1 ) { //May not have kinematic hardening 00123 alpha = EPS->getTensorVar(1); 00124 s_bar = sdev - (alpha*p); 00125 temp5 = alpha("ij") * s_bar("ij"); 00126 temp5.null_indices(); 00127 temp6 = temp5.trace(); 00128 Tnsr0 = KroneckerI*(temp6/3.0); 00129 } 00130 else { 00131 s_bar = sdev; 00132 } 00133 stresstensor temp3 = s_bar("ij") * s_bar("ij"); 00134 temp3.null_indices(); 00135 double temp4 = temp3.trace(); 00136 temp4 = sqrt(temp4); 00137 Tnsr0 += s_bar; 00138 double eps = pow( d_macheps(), 0.5 ); 00139 if ( fabs(temp4) > eps ) { 00140 Tnsr0 = Tnsr0 * (halfRt2/temp4); 00141 } 00142 00143 dFods = (KroneckerI*temp1) + Tnsr0; 00144 00145 return dFods; 00146 } 00147 00148 //================================================================================ 00149 // double xi_s1 = dF/dS1 = dF/dalfa1 = I1 Derivative in terms of first scalar var 00150 //================================================================================ 00151 00152 double DPYieldSurface::xi_s1( const EPState *EPS ) const { 00153 00154 double temp = EPS->getStress().Iinvariant1(); 00155 return temp; 00156 00157 } 00158 00159 //================================================================================ 00160 // double xi_s2 = dF/dS2 = dF/k = -1.0 Derivative in terms of second scalar var 00161 //================================================================================ 00162 00163 double DPYieldSurface::xi_s2( const EPState *EPS ) const { 00164 00165 return -1.0; 00166 00167 } 00168 00169 00170 //================================================================================ 00171 // double xi_t1 = dF/dt1 = dF/dalpha 00172 //================================================================================ 00173 00174 tensor DPYieldSurface::xi_t1( const EPState *EPS) const { 00175 stresstensor xi; 00176 00177 stresstensor alpha; 00178 stresstensor sigma = EPS->getStress(); 00179 double p = sigma.p_hydrostatic(); 00180 double halfRt2 = 0.5 * sqrt(2.0); 00181 int nod = EPS->getNTensorVar(); 00182 if ( nod >=1 ) { //May not have kinematic hardening 00183 alpha = EPS->getTensorVar(1); 00184 } 00185 stresstensor sigma_bar = sigma - alpha*p; 00186 stresstensor s_bar = sigma_bar.deviator(); 00187 stresstensor temp3 = s_bar("ij") * s_bar("ij"); 00188 temp3.null_indices(); 00189 double temp4 = temp3.trace(); 00190 temp4 = sqrt(temp4); 00191 double eps = pow( d_macheps(), 0.5 ); 00192 if ( fabs(temp4) > eps ) { 00193 xi = s_bar * (-halfRt2*p/temp4); //Note the Negative 1 here 00194 } 00195 00196 return xi; 00197 } 00198 00199 //================================================================================ 00200 OPS_Stream& operator<< (OPS_Stream& os, const DPYieldSurface & YS) 00201 { 00202 os << "Drucker-Prager Yield Surface Parameters: " << endln; 00203 os.precision(4); 00204 //os << "alfa1 = " << YS.getalfa1() << endlnn; 00205 return os; 00206 } 00207 00208 #endif
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