MD_PS01.cppGo to the documentation of this file.00001 00002 //================================================================================ 00003 //# COPYRIGHT (C): :-)) # 00004 //# PROJECT: Object Oriented Finite Element Program # 00005 //# PURPOSE: Manzari-Dafalias potential surface 01(with Pc) # 00006 //# (Ref. Geotechnique v.47 No.2 255-272, 1997) # 00007 //# CLASS: MDPotentialSurface01 # 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: August 03 '93 # 00016 //# UPDATE HISTORY: August 08 '00 # 00017 //# # 00018 //# # 00019 //# # 00020 //# # 00021 //# # 00022 //================================================================================ 00023 00024 00025 #ifndef MD_PS01_CPP 00026 #define MD_PS01_CPP 00027 00028 #include "MD_PS01.h" 00029 #include <basics.h> 00030 #include <math.h> 00031 00032 //================================================================================ 00033 // Normal constructor 00034 //================================================================================ 00035 00036 MDPotentialSurface01::MDPotentialSurface01(double pc ) 00037 { 00038 Pc = pc; 00039 } 00040 00041 00042 //================================================================================ 00043 //create a colne of itself 00044 //================================================================================ 00045 00046 PotentialSurface * MDPotentialSurface01::newObj() { 00047 00048 PotentialSurface *new_PS = new MDPotentialSurface01(Pc); 00049 return new_PS; 00050 00051 } 00052 00053 //================================================================================ 00054 // tensor dQ/dsigma_ij: the normal to the potential surface 00055 //================================================================================ 00056 00057 tensor MDPotentialSurface01::dQods(const EPState *EPS) const { 00058 00059 tensor dQoverds( 2, def_dim_2, 0.0); 00060 tensor I2("I", 2, def_dim_2); 00061 00062 tensor S = EPS->getStress().deviator(); 00063 double p = EPS->getStress().p_hydrostatic(); 00064 p = p - Pc; 00065 00066 tensor alpha = EPS->getTensorVar( 1 ); //the first tensor hardening var is alpha_ij 00067 00068 //double m = EPS->getScalarVar( 1 ); //the first scalar hardening var is m 00069 double D = EPS->getScalarVar( 2 ); // D is stored in scalar var array's second cell 00070 //printf("Here we go! D= %e\n", D); 00071 00072 tensor r = S *(1.0/p); 00073 tensor temp_f1 = r - alpha; 00074 tensor temp_f2 = temp_f1("ij") * temp_f1("ij"); 00075 double temp_f3 = sqrt( temp_f2.trace() ); 00076 00077 stresstensor n; 00078 if ( temp_f3 >= d_macheps() ){ 00079 n = temp_f1 * (1.0 / temp_f3 ); 00080 } 00081 else { 00082 opserr << "MDPotentialSurface01::dQods |n_ij| = 0, divide by zero! Program exits.\n"; 00083 exit(-1); 00084 //::printf(" \n\n n_ij not defined!!!! Program exits\n"); 00085 //exit(1); 00086 } 00087 00088 dQoverds = n + I2 * (D *(1.0/3.0)); 00089 00090 return dQoverds; 00091 } 00092 00093 00094 //================================================================================ 00095 // tensor d2Q/dsigma_ij2: the second derivatives to the potential surface 00096 //================================================================================ 00097 00098 tensor MDPotentialSurface01::d2Qods2(const EPState *EPS) const 00099 { 00100 00101 tensor d2Qoverds2; 00102 tensor I2("I", 2, def_dim_2); 00103 00104 stresstensor stress = EPS->getStress(); 00105 00106 //double J2D = stress.Jinvariant2(); 00107 //double q = stress.q_deviatoric(); 00108 double theta = stress.theta(); 00109 00110 tensor S = stress.deviator(); 00111 double p = stress.p_hydrostatic(); 00112 p = p - Pc; 00113 opserr << " p = " << p << endln; 00114 00115 tensor alpha = EPS->getTensorVar( 1 ); //the first tensor hardening var is alpha_ij 00116 00117 tensor r = S *(1.0/p); 00118 tensor temp_f1 = r - alpha; 00119 tensor temp_f2 = temp_f1("ij") * temp_f1("ij"); 00120 double temp_f3 = sqrt( temp_f2.trace() ); 00121 00122 stresstensor n; 00123 if ( temp_f3 >= d_macheps() ){ 00124 n = temp_f1 * (1.0 / temp_f3 ); 00125 } 00126 else { 00127 opserr << "MDPotentialSurface01::dQods |n_ij| = 0, divide by zero! Program exits.\n"; 00128 exit(-1); 00129 //::printf(" \n\n n_ij not defined!!!! Program exits\n"); 00130 //exit(1); 00131 } 00132 00133 00134 //tensor d2Qoverds2( 2, def_dim_2, 0.0); // dummy second derivatives. To be redefined. 00135 tensor dnds = dnods( EPS); 00136 double A = 2.64; 00137 double Mc = 1.14; 00138 double kcd = 4.2; 00139 00140 if (p < 0.0) 00141 { 00142 opserr << "MDPotentialSurface01::d2Qods2 p < 0, Program exits.\n"; 00143 exit(-1); 00144 } 00145 double ec = 0.80 - 0.025 * log( p / 160 ); 00146 00147 double e = EPS->getScalarVar(3); 00148 double xi = e - ec; 00149 00150 double c = 1.0; 00151 double cd = 0.0167; 00152 00153 double dgodthetac = dgoverdt(theta, c); 00154 double dgodthetacd = dgoverdt(theta, cd); 00155 tensor dthetaods = dthetaoverds(EPS); 00156 //stresstensor d = EPS->getTensorVar( 3 ); // getting d_ij from EPState 00157 00158 //tensor dtods = dthetaods("pq"); 00159 00160 tensor dDds1 = A*dthetaods*(dgodthetac*Mc+dgodthetacd*kcd*xi)*sqrt(2.0/3.0); 00161 tensor dDds2 = A*apqdnods(EPS); 00162 tensor dDds = dDds1 - dDds2; 00163 dDds.null_indices(); 00164 00165 d2Qoverds2 = dnds + I2("pq")*dDds("mn")*(0.33333333); 00166 d2Qoverds2.null_indices(); 00167 00168 return d2Qoverds2; 00169 } 00170 00171 //================================================================================ 00172 // tensor dn_pq/dsigma_mn 00173 //================================================================================ 00174 tensor MDPotentialSurface01::dnods(const EPState *EPS) const 00175 { 00176 tensor dnods( 2, def_dim_2, 0.0); 00177 tensor I2("I", 2, def_dim_2); 00178 00179 stresstensor S = EPS->getStress().deviator(); 00180 //S.reportshort("S"); 00181 stresstensor alpha = EPS->getTensorVar( 1 ); 00182 00183 double p = EPS->getStress().p_hydrostatic(); 00184 p = p - Pc; 00185 00186 //double m = EPS->getScalarVar( 1 ); 00187 //double fac = 1.0/(sqrt(2.0/3.0)*m); 00188 00189 tensor Ipmnq = I2("pm") * I2("nq"); 00190 tensor Imnpq = I2("mn") * I2("pq"); 00191 tensor Apqmn = alpha("pq") * I2("mn"); 00192 tensor X = Ipmnq - Imnpq * (1.0/3.0) - Apqmn * (1.0/3.0); 00193 00194 //Ipmnq.print("in MD_PS01", ""); 00195 00196 tensor sa = S("ij") * alpha("ij"); 00197 double sad =sa.trace(); 00198 tensor aa = alpha("ij") * alpha("ij"); 00199 double aad =aa.trace(); 00200 tensor Y = S - ( p + 0.333333333*(sad-p*aad) )*I2; 00201 Y.null_indices(); 00202 00203 tensor s_bar = S("ij") - p*alpha("ij"); 00204 s_bar.null_indices(); 00205 tensor norm2 = s_bar("ij") * s_bar("ij"); 00206 double norm = norm2.trace(); 00207 00208 s_bar.null_indices(); 00209 tensor tmp = s_bar("pq")*Y("mn"); 00210 dnods = ( X - 2.0 * tmp)*(1.0/norm); 00211 00212 return dnods; 00213 } 00214 00215 //================================================================================ 00216 // tensor alpha_pq*dnods 00217 //================================================================================ 00218 tensor MDPotentialSurface01::apqdnods(const EPState *EPS) const 00219 { 00220 tensor ddnods( 2, def_dim_2, 0.0); 00221 tensor I2("I", 2, def_dim_2); 00222 00223 stresstensor S = EPS->getStress().deviator(); 00224 //S.reportshort("S"); 00225 00226 double p = EPS->getStress().p_hydrostatic(); 00227 p = p - Pc; 00228 00229 //printf("Here we go! p %f\n", p); 00230 00231 stresstensor alpha = EPS->getTensorVar( 1 ); 00232 //stresstensor d = EPS->getTensorVar( 3 ); // getting d_ij from EPState 00233 00234 tensor akk = alpha("ij") * I2("ij"); 00235 double akkd =akk.trace(); 00236 00237 tensor aa = alpha("ij") * alpha("ij"); 00238 double aad =aa.trace(); 00239 00240 tensor aX = alpha - I2*(akkd +aad)*0.3333333; 00241 00242 //Ymn 00243 tensor sa = S("ij") * alpha("ij"); 00244 double sad =sa.trace(); 00245 tensor Y = S - ( p + 0.333333333*(sad-p*aad) )*I2; 00246 Y.null_indices(); 00247 00248 tensor s_bar = S - p*alpha; 00249 s_bar.null_indices(); 00250 00251 tensor as_bar = alpha("pq") * s_bar("pq"); 00252 double as_bard = as_bar.trace(); 00253 s_bar.null_indices(); 00254 00255 //Norm 00256 tensor norm2 = s_bar("ij") * s_bar("ij"); 00257 double norm = norm2.trace(); 00258 00259 00260 ddnods = (aX - Y*2.0*as_bard)*(1.0/norm); 00261 return ddnods; 00262 } 00263 00264 00265 //================================================================================ 00266 // tensor dg(theta, c)/dsigma_mn 00267 //================================================================================ 00268 double MDPotentialSurface01::dgoverdt(double theta, double c) const 00269 { 00270 //stresstensor stress = EPS->getStress(); 00271 00272 //double J2D = stress.Jinvariant2(); 00273 //double q = stress.q_deviatoric(); 00274 //double theta = stress.theta(); 00275 00276 double temp = (-6*(1 - c)*c*sin(3*theta))/pow(1 + c - (1 - c)*cos(3*theta),2); 00277 return temp; 00278 } 00279 00280 //================================================================================ 00281 // tensor dtheta/dsigma_pq 00282 //================================================================================ 00283 tensor MDPotentialSurface01::dthetaoverds(const EPState *EPS) const 00284 { 00285 tensor ret(2, def_dim_2, 0.0); 00286 stresstensor s( 0.0); 00287 stresstensor t( 0.0); 00288 tensor I2("I", 2, def_dim_2); 00289 00290 //double EPS = pow(d_macheps(),(1./3.)); 00291 stresstensor stress = EPS->getStress(); 00292 00293 double J2D = stress.Jinvariant2(); 00294 double q = stress.q_deviatoric(); 00295 double theta = stress.theta(); 00296 00297 //out while ( theta >= 2.0*PI ) 00298 //out theta = theta - 2.0*PI; // if bigger than full cycle 00299 //out while ( theta >= 4.0*PI/3.0 ) 00300 //out theta = theta - 4.0*PI/3.0; // if bigger than four thirds of half cycle 00301 //out while ( theta >= 2.0*PI/3.0 ) 00302 //out theta = theta - 2.0*PI/3.0; // if bigger than two third of half cycle 00303 //out while ( theta >= PI/3.0 ) 00304 //out theta = 2.0*PI/3.0 - theta; // if bigger than one third of half cycle 00305 //out 00306 //out if ( theta < 0.0001 ) 00307 //out { 00308 //out ::printf("theta = %.10e in Material_Model::dthetaoverds(stress)\n", 00309 //out theta); 00310 //out//>><<>><<>><< return ret; 00311 //out } 00312 00313 double c3t = cos(3.0*theta); 00314 double s3t = sin(3.0*theta); 00315 00316 double tempS = (3.0/2.0)*(c3t/(q*q*s3t)); 00317 double tempT = (9.0/2.0)*(1.0/(q*q*q*s3t)); 00318 00319 s = stress.deviator(); 00320 t = s("qk")*s("kp") - I2*(J2D*(2.0/3.0)); 00321 00322 s.null_indices(); 00323 t.null_indices(); 00324 00325 ret = s*tempS - t*tempT; 00326 ret.null_indices(); 00327 return ret; 00328 } 00329 00330 OPS_Stream& operator<< (OPS_Stream& os, const MDPotentialSurface01 &PS) 00331 { 00332 os << "Manzari-Dafalias Potential Surface 01( with Pc) Parameters: " << endln; 00333 os << "Pc = " << PS.Pc << endln; 00334 return os; 00335 } 00336 00337 00338 00339 #endif 00340
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