DM04_alpha_Eij.cppGo to the documentation of this file.00001 00002 // COPYLEFT (C): Woody's viral GPL-like license (by BJ): 00003 // ``This source code is Copyrighted in 00004 // U.S., for an indefinite period, and anybody 00005 // caught using it without our permission, will be 00006 // mighty good friends of ourn, cause we don't give 00007 // a darn. Hack it. Compile it. Debug it. Run it. 00008 // Yodel it. Enjoy it. We wrote it, that's all we 00009 // wanted to do.'' 00010 // 00011 // 00012 // COPYRIGHT (C): :-)) 00013 // PROJECT: Object Oriented Finite Element Program 00014 // FILE: 00015 // CLASS: 00016 // MEMBER FUNCTIONS: 00017 // 00018 // MEMBER VARIABLES 00019 // 00020 // PURPOSE: 00021 // 00022 // RETURN: 00023 // VERSION: 00024 // LANGUAGE: C++ 00025 // TARGET OS: 00026 // DESIGNER: Zhao Cheng, Boris Jeremic 00027 // PROGRAMMER: Zhao Cheng, Mahdi Taiebat 00028 // DATE: Fall 2005 00029 // UPDATE HISTORY: 00030 // 00032 // 00033 00034 // Ref: Dafalias and Manzari 2004: J. Eng. Mech. 130(6), pp 622-634 00035 // Parameters: 00036 // 1- e0: initial void ratio; 00037 // 2- e_r: reference void for critical state line, ec = e_r - lambda_c*(pc/Pat)^xi; 00038 // 3- lambda_c: parameter for critical state line, ec = e_r - lambda_c*(pc/Pat)^xi; 00039 // 4- xi: parameter for critical state line, ec = e_r - lambda_c*(pc/Pat)^xi; 00040 // 5- Pat: atmospherics pressure for critical state line, ec = e0 - lambda_c*(pc/Pat)^xi; 00041 // 6- m: parameter in the yield function; 00042 // 7- M: critical state stress ration; 00043 // 8- cc: tension-compression strength ratio; 00044 // 9- nb: parameter; 00045 // 10- h0: parameter; 00046 // 11- ch: parameter; 00047 // 12- G0: parameter in the elastic part 00048 // 13- alpha: "back-stress ratio" tensor in yield function; (the 1st tensorial internal variable); 00049 // 14- z: fabric dilation internal tensor (the 2nd tensorial internal variable); 00050 00051 #ifndef DM04_alpha_Eij_CPP 00052 #define DM04_alpha_Eij_CPP 00053 00054 #include "DM04_alpha_Eij.h" 00055 00056 stresstensor DM04_alpha_Eij::DM04_alpha_t; 00057 00058 DM04_alpha_Eij::DM04_alpha_Eij(int e0_index_in, 00059 int e_r_index_in, 00060 int lambda_c_index_in, 00061 int xi_index_in, 00062 int Pat_index_in, 00063 int m_index_in, 00064 int M_cal_index_in, 00065 int cc_index_in, 00066 int nb_index_in, 00067 int h0_index_in, 00068 int ch_index_in, 00069 int G0_index_in, 00070 int alpha_index_in, 00071 int z_index_in) 00072 : e0_index(e0_index_in), 00073 e_r_index(e_r_index_in), 00074 lambda_c_index(lambda_c_index_in), 00075 xi_index(xi_index_in), 00076 Pat_index(Pat_index_in), 00077 m_index(m_index_in), 00078 M_cal_index(M_cal_index_in), 00079 cc_index(cc_index_in), 00080 nb_index(nb_index_in), 00081 h0_index(h0_index_in), 00082 ch_index(ch_index_in), 00083 G0_index(G0_index_in), 00084 alpha_index(alpha_index_in), 00085 z_index(z_index_in) 00086 { 00087 a_index = 0; 00088 stresstensor zT; 00089 alpha_in.Initialize(zT); 00090 } 00091 00092 TensorEvolution* DM04_alpha_Eij::newObj() 00093 { 00094 TensorEvolution* nObj = new DM04_alpha_Eij(this->e0_index, 00095 this->e_r_index, 00096 this->lambda_c_index, 00097 this->xi_index, 00098 this->Pat_index, 00099 this->m_index, 00100 this->M_cal_index, 00101 this->cc_index, 00102 this->nb_index, 00103 this->h0_index, 00104 this->ch_index, 00105 this->G0_index, 00106 this->alpha_index, 00107 this->z_index); 00108 return nObj; 00109 } 00110 00111 const straintensor& DM04_alpha_Eij::Hij(const straintensor& plastic_flow, const stresstensor& Stre, 00112 const straintensor& Stra, const MaterialParameter& material_parameter) 00113 { 00114 stresstensor alpha_alpha_in; 00115 double a_in = 0.0; 00116 double h = 0.0; 00117 00118 double e0 = gete0(material_parameter); 00119 double e_r = gete_r(material_parameter); 00120 double lambda_c = getlambda_c(material_parameter); 00121 double xi = getxi(material_parameter); 00122 double Pat = getPat(material_parameter); 00123 double m = getm(material_parameter); 00124 double M_cal = getM_cal(material_parameter); 00125 double cc = getcc(material_parameter); 00126 double nb = getnb(material_parameter); 00127 double h0 = geth0(material_parameter); 00128 double ch = getch(material_parameter); 00129 double G0 = getG0(material_parameter); 00130 00131 stresstensor alpha = getalpha(material_parameter); 00132 stresstensor z = getz(material_parameter); 00133 00134 double p = Stre.p_hydrostatic(); 00135 stresstensor s = Stre.deviator(); 00136 00137 stresstensor n; 00138 stresstensor alpha_b; 00139 stresstensor alpha_b_alpha; 00140 double b0 = 0.0; 00141 double g = 0.0; 00142 double ec = e_r; 00143 double stateParameter = 0.0; 00144 double expnb = 1.0; 00145 double ab = 0.0; 00146 stresstensor s_bar; 00147 double _s_bar_ = 0.0; 00148 00149 double J3D; 00150 double cos3theta = 0.0; 00151 00152 double e = e0 + (1.0 + e0) *Stra.Iinvariant1(); 00153 00154 s_bar = Stre.deviator() - (alpha *p); 00155 _s_bar_ = sqrt( (s_bar("ij")*s_bar("ij")).trace() ); 00156 if (p > 0.0 && _s_bar_ > 0.0) 00157 n = s_bar * (1.0/_s_bar_); 00158 00159 J3D = n.Jinvariant3(); 00160 cos3theta = -3.0*sqrt(6.0) *J3D; 00161 00162 if (p <= 0.0) 00163 cos3theta = 1.0; 00164 00165 if (cos3theta > 1.0) 00166 cos3theta = 1.0; 00167 00168 if (cos3theta < -1.0) 00169 cos3theta = -1.0; 00170 00171 g = getg(cc, cos3theta); 00172 00173 if ( (p/Pat) >= 0.0 ) 00174 ec = getec(e_r, lambda_c, xi, Pat, p); 00175 00176 stateParameter = e - ec; 00177 00178 expnb = exp( -nb *stateParameter ); 00179 00180 // Using the other equation 00181 //ab = sqrt(2.0/3.0) * ( g *M_cal *expnb - m); 00182 //alpha_b = n *ab; 00183 //alpha_b_alpha = alpha_b - alpha; 00184 00185 // Replacing the above equation 00186 ab = sqrt(2.0/3.0) * g *M_cal *expnb; 00187 alpha_b_alpha = n *ab - s *(1.0/p); 00188 00189 00190 if ( (p/Pat) >= 0.01 ) 00191 b0 = G0 *h0 *(1.0-ch*e) *pow(p/Pat, -0.5); 00192 else 00193 b0 = G0 *h0 *(1.0-ch*e) *10.0; 00194 00195 if (a_index == 0) { 00196 alpha_in.Initialize(alpha); 00197 a_in = 0.0; 00198 a_index = 1; 00199 h = 1.0e10 * b0; 00200 } 00201 else { 00202 alpha_alpha_in = alpha - alpha_in; 00203 a_in = (alpha_alpha_in("ij")*n("ij")).trace(); 00204 if ( a_in < 0.0 ) { 00205 a_index = 0; 00206 alpha_in.Initialize(alpha); 00207 } 00208 if ( a_in < 1.0e-10 ) 00209 a_in = 1.0e-10; 00210 h = b0 /a_in; 00211 } 00212 00213 TensorEvolution::TensorEvolutionHij = alpha_b_alpha *(h*2.0/3.0); 00214 00215 return TensorEvolution::TensorEvolutionHij; 00216 } 00217 00218 // to get e0 00219 //================================================================================ 00220 double DM04_alpha_Eij::gete0(const MaterialParameter& material_parameter) const 00221 { 00222 return getParameters(material_parameter, e0_index); 00223 } 00224 00225 // to get e_r 00226 //================================================================================ 00227 double DM04_alpha_Eij::gete_r(const MaterialParameter& material_parameter) const 00228 { 00229 return getParameters(material_parameter, e_r_index); 00230 } 00231 00232 // to get lambda_c 00233 //================================================================================ 00234 double DM04_alpha_Eij::getlambda_c(const MaterialParameter& material_parameter) const 00235 { 00236 return getParameters(material_parameter, lambda_c_index); 00237 } 00238 00239 // to get xi 00240 //================================================================================ 00241 double DM04_alpha_Eij::getxi(const MaterialParameter& material_parameter) const 00242 { 00243 return getParameters(material_parameter, xi_index); 00244 00245 } 00246 00247 // to get Pat 00248 //================================================================================ 00249 double DM04_alpha_Eij::getPat(const MaterialParameter& material_parameter) const 00250 { 00251 return getParameters(material_parameter, Pat_index); 00252 } 00253 00254 // to get m 00255 //================================================================================ 00256 double DM04_alpha_Eij::getm(const MaterialParameter& material_parameter) const 00257 { 00258 return getParameters(material_parameter, m_index); 00259 } 00260 00261 // to get M 00262 //================================================================================ 00263 double DM04_alpha_Eij::getM_cal(const MaterialParameter& material_parameter) const 00264 { 00265 return getParameters(material_parameter, M_cal_index); 00266 } 00267 00268 // to get c 00269 //================================================================================ 00270 double DM04_alpha_Eij::getcc(const MaterialParameter& material_parameter) const 00271 { 00272 return getParameters(material_parameter, cc_index); 00273 } 00274 00275 // to get n_d 00276 //================================================================================ 00277 double DM04_alpha_Eij::getnb(const MaterialParameter& material_parameter) const 00278 { 00279 return getParameters(material_parameter, nb_index); 00280 } 00281 00282 00283 // to get h0 00284 //================================================================================ 00285 double DM04_alpha_Eij::geth0(const MaterialParameter& material_parameter) const 00286 { 00287 return getParameters(material_parameter, h0_index); 00288 00289 } 00290 00291 00292 // to get ch 00293 //================================================================================ 00294 double DM04_alpha_Eij::getch(const MaterialParameter& material_parameter) const 00295 { 00296 return getParameters(material_parameter, ch_index); 00297 00298 } 00299 00300 // to get G0 00301 //================================================================================ 00302 double DM04_alpha_Eij::getG0(const MaterialParameter& material_parameter) const 00303 { 00304 return getParameters(material_parameter, G0_index); 00305 00306 } 00307 00308 // to get alpha 00309 //================================================================================ 00310 const stresstensor& DM04_alpha_Eij::getalpha(const MaterialParameter& material_parameter) const 00311 { 00312 if ( alpha_index <= material_parameter.getNum_Internal_Tensor() && alpha_index > 0) { 00313 DM04_alpha_Eij::DM04_alpha_t = material_parameter.getInternal_Tensor(alpha_index-1); 00314 return DM04_alpha_Eij::DM04_alpha_t; 00315 } 00316 else { 00317 opserr << "DM04_alpha: Invalid Input (alpha) " << endln; 00318 exit (1); 00319 } 00320 } 00321 00322 // to get z 00323 //================================================================================ 00324 const stresstensor& DM04_alpha_Eij::getz(const MaterialParameter& material_parameter) const 00325 { 00326 if ( z_index <= material_parameter.getNum_Internal_Tensor() && z_index > 0) { 00327 DM04_alpha_Eij::DM04_alpha_t = material_parameter.getInternal_Tensor(z_index-1); 00328 return DM04_alpha_Eij::DM04_alpha_t; 00329 } 00330 else { 00331 opserr << "DM04_alpha: Invalid Input (z) " << endln; 00332 exit (1); 00333 } 00334 } 00335 00336 //================================================================================ 00337 double DM04_alpha_Eij::getParameters(const MaterialParameter& material_parameter, int which) const 00338 { 00339 if ( which <= material_parameter.getNum_Material_Parameter() && which > 0) 00340 return material_parameter.getMaterial_Parameter(which-1); 00341 else { 00342 opserr << "DM04_alpha: Invalid Input - #" << which << endln; 00343 exit (1); 00344 } 00345 } 00346 00347 00348 //================================================================================ 00349 double DM04_alpha_Eij::getec(double e_r, double lambda_c, double xi, double Pat, double p_c) const 00350 { 00351 double ee = e_r; 00352 00353 if ( p_c/Pat >= 0.0 ) 00354 ee = e_r - lambda_c * pow(p_c/Pat, xi); 00355 else 00356 opserr << "Warning: DM04_alpha_Eij - 'p_c/Pat' less than zero! " << endln; 00357 00358 return ee; 00359 } 00360 00361 //================================================================================ 00362 double DM04_alpha_Eij::getg(double c, double cos3theta) const 00363 { 00364 return 2.0 * c / ( (1.0+c) - (1.0-c)*cos3theta ); 00365 } 00366 00367 00368 #endif 00369
Generated on Mon Oct 23 15:05:14 2006 for OpenSees by |
opensees-support @ berkeley.edu
Supported by the National Science Foundation
©2006, UC Regents