MD_YS.cppGo to the documentation of this file.00001 00002 //================================================================================ 00003 //# COPYRIGHT (C): :-)) # 00004 //# PROJECT: Object Oriented Finite Element Program # 00005 //# PURPOSE: Mazari-Dafalias yield criterion # 00006 //# (Ref. Geotechnique # 00007 //# V.47 No.2 255-272, 1997) # 00008 //# CLASS: MDYieldSurface01 # 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 MD_YS_CPP 00027 #define MD_YS_CPP 00028 00029 #include "MD_YS.h" 00030 #include <basics.h> 00031 00032 #define sqrt23rd 0.8165 00033 00034 //================================================================================ 00035 // Normal constructor 00036 //================================================================================ 00037 00038 MDYieldSurface::MDYieldSurface( ) { } 00039 00040 00041 //================================================================================ 00042 //create a colne of itself 00043 //================================================================================ 00044 00045 YieldSurface * MDYieldSurface::newObj() { 00046 00047 YieldSurface *new_YS = new MDYieldSurface(); 00048 return new_YS; 00049 00050 } 00051 00052 //================================================================================ 00053 // Copy constrstructor 00054 //================================================================================ 00055 // 00056 //MDYieldSurface::MDYieldSurface(MDYieldSurface &MDYS ) { } 00057 // 00058 00059 00060 //================================================================================ 00061 // Yield criterion evaluation function F(EPState) 00062 //================================================================================ 00063 00064 double MDYieldSurface::f(const EPState *EPS) const 00065 { 00066 stresstensor S = EPS->getStress().deviator(); 00067 //cout << "s " << S; 00068 00069 double p =EPS->getStress().p_hydrostatic(); 00070 //cout << "p " << p; 00071 00072 // Get a_ij 00073 stresstensor alpha = EPS->getTensorVar( 1 ); 00074 //cout << "alpha " << alpha; 00075 00076 double m = EPS->getScalarVar(1); 00077 00078 stresstensor r = S - alpha*p; 00079 r.null_indices(); 00080 stresstensor r2 = r("ij") * r("ij"); 00081 00082 stresstensor n; 00083 00084 double norm = sqrt( r2.trace() ); 00085 00087 //if ( norm >= d_macheps() ){ 00088 // n = r *(1.0 / norm ); 00089 //} 00090 //else { 00091 // g3ErrorHandler->fatal("MDYieldSurface::f |n_ij| = 0, divide by zero! Program exits."); 00092 // exit(-1); 00093 //} 00094 //EPS->setTensorVar( 2, n); //update n_ij// 00095 00096 00097 double ff = norm - sqrt23rd * m * p; 00098 //printf("\n========Inside f temp3 = %.4f x = %.4f\n ", temp3, x); 00099 00100 return ff; 00101 } 00102 00103 00104 //================================================================================ 00105 // tensor dF/dsigma_ij 00106 //================================================================================ 00107 00108 tensor MDYieldSurface::dFods(const EPState *EPS) const { 00109 00110 tensor dFoverds( 2, def_dim_2, 0.0); 00111 tensor I2("I", 2, def_dim_2); 00112 00113 stresstensor S = EPS->getStress().deviator(); 00114 //S.reportshort("S"); 00115 00116 double m = EPS->getScalarVar( 1 ); 00117 00118 double p = EPS->getStress().p_hydrostatic(); 00119 //printf("Here we go! p %f\n", p); 00120 00121 stresstensor alpha = EPS->getTensorVar( 1 ); // getting alpha_ij from EPState 00122 //alpha.reportshort("alpha"); 00123 00124 //stresstensor n = EPS->getTensorVar( 3 ); // getting n_ij from EPState 00125 //n.reportshort("n"); 00126 00127 // might be moved to Evolution Law 00128 stresstensor n = EPS->getTensorVar( 2 ); 00129 // cout << n; 00130 //------------------------------------------------- 00131 // //stresstensor n; 00132 // stresstensor r = S * (1.0 / p); 00133 // //r.reportshort("r"); 00134 // stresstensor r_bar = r - alpha; 00135 // //r_bar.reportshort("r_bar"); 00136 // 00137 // stresstensor norm2 = r_bar("ij") * r_bar("ij"); 00138 // double norm = sqrt( norm2.trace() ); 00139 // //cout << "d_macheps " << d_macheps() << endln; 00140 // if ( norm >= d_macheps() ){ 00141 // n = r_bar *(1.0 / norm ); 00142 // } 00143 // else { 00144 // g3ErrorHandler->fatal("MDYieldSurface::dFods |n_ij| = 0, divide by zero! Program exits."); 00145 // exit(-1); 00146 // } 00147 //EPS->setTensorVar( 3, n); //update n_ij// 00152 00153 00154 00155 00156 //tensorial multiplication produces 1st-order tensor 00157 //tensor temp = n("ij") * n("ij"); 00158 //double temp1 = temp.trace(); 00159 //printf("==== n_ij*n_ij %e\n", temp1); 00160 00162 tensor temp = alpha("ij") * n("ij"); 00163 double N = temp.trace() + sqrt23rd*m; 00164 //double N = temp.trace() + 2.0/3.0*m; 00165 //printf(" N = %e\n", N); 00166 00167 dFoverds = n - I2 * N*(1.0/3.0); 00168 00169 return dFoverds; 00170 00171 } 00172 00173 00174 //================================================================================ 00175 // double xi_s1 = dF/dm = -(2/3)^0.5 p Derivative in terms of first scalar var 00176 //================================================================================ 00177 00178 double MDYieldSurface::xi_s1( const EPState *EPS ) const { 00179 00180 double p = EPS->getStress().p_hydrostatic(); 00181 return -1.0*sqrt23rd * p; 00182 00183 } 00184 00185 //================================================================================ 00186 // double xi_t1 = dF/dt1 = dF/dalpha= -p*n Derivative in terms of first tensorial var 00187 //================================================================================ 00188 00189 tensor MDYieldSurface::xi_t1( const EPState *EPS) const { 00190 tensor dFoverds( 2, def_dim_2, 0.0); 00191 tensor I2("I", 2, def_dim_2); 00192 00193 stresstensor S = EPS->getStress().deviator(); 00194 00195 double p = EPS->getStress().p_hydrostatic(); 00196 00197 stresstensor n = EPS->getTensorVar( 2 ); 00198 //stresstensor n; 00199 00201 //double m = EPS->getScalarVar( 1 ); 00202 //stresstensor alpha = EPS->getTensorVar( 1 ); // getting alpha_ij from EPState 00203 //stresstensor r = S * (1.0 / p); 00204 //stresstensor r_bar = r - alpha; 00205 //stresstensor norm2 = r_bar("ij") * r_bar("ij"); 00206 //double norm = sqrt( norm2.trace() ); 00207 // 00208 //if ( norm >= d_macheps() ){ 00209 // n = r_bar *(1.0 / norm ); 00210 //} 00211 //else { 00212 // opserr->fatal("MDYieldSurface::dFods |n_ij| = 0, divide by zero! Program exits."); 00213 // exit(-1); 00214 //} 00215 //n = r_bar *(1.0 / sqrt23rd / m ); 00217 00218 return n * (-1.0)* p; 00219 } 00220 00221 00222 OPS_Stream& operator<< (OPS_Stream& os, const MDYieldSurface & YS) 00223 { 00224 os << "Manzari-Dafalias Yield Surface Parameters: " << endln; 00225 return os; 00226 } 00227 00228 #endif 00229
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