Tri_a_fail_crit_YS.cppGo to the documentation of this file.00001 00002 //################################################################################ 00003 //# COPYRIGHT (C): :-)) # 00004 //# PROJECT: OpenSees # 00005 //# PURPOSE: Triaxial Failure Criterion for Concrete - yield criterion # 00006 //# CLASS: TriFCYieldSurface # 00007 //# # 00008 //# VERSION: 1.0 # 00009 //# LANGUAGE: C++ (ili tako nesto) # 00010 //# TARGET OS: # 00011 // DESIGNER(S): Boris Jeremic and Zhaohui Yang [jeremic,zhyang]@ucdavis.edu| 00012 // PROGRAMMER(S): Vlado Vukadin | 00013 //# # 00014 //# # 00015 //# DATE: June 01, 2002 # 00016 //# UPDATE HISTORY: bice tako dobr da nece biti potreban update :) # 00017 //# # 00018 //# # 00019 //# # 00020 //# # 00021 //# SHORT EXPLANATION: # 00022 //# # 00023 //# Yield surface is based on article by Menetrey, P. and William, K.J. # 00024 //# published in 1995 in ACI Structural Journal pp 311-318. Purpose of the # 00025 //# Yield surface is to model triaxial strenght of concrete. All the necessary # 00026 //# explanations about functions are there # 00027 //# # 00028 //# # 00029 //################################################################################ 00030 //*/ 00031 00032 00033 #ifndef Tri_a_fail_crit_YS_CPP 00034 #define Tri_a_fail_crit_YS_CPP 00035 00036 #include "Tri_a_fail_crit_YS.h" 00037 00038 00039 //================================================================================ 00040 // Normal constructor 00041 //================================================================================ 00042 00043 TriFCYieldSurface::TriFCYieldSurface (double fc, double ft, double e, double coh ) 00044 { 00045 00046 fcomp = fc; // a way to fix those two factors so they are not 00047 ftens = ft; // in need for evolution law, do not need to be derived and could be 00048 el = e; // suplied as a constansts (how this is done look in definition of Cam Clay model) 00049 c = coh; 00050 } 00051 //================================================================================ 00052 // Destructor 00053 //================================================================================ 00054 00055 TriFCYieldSurface::~TriFCYieldSurface ( ) 00056 { 00057 00058 } 00059 00060 //================================================================================ 00061 //create a clone of itself 00062 //================================================================================ 00063 00064 YieldSurface * TriFCYieldSurface ::newObj() 00065 { 00066 YieldSurface *new_YS = new TriFCYieldSurface (fcomp,ftens, el, c); 00067 return new_YS; 00068 } 00069 00070 //================================================================================ 00071 // Yield criterion evaluation function F(EPState) 00072 //================================================================================ 00073 00074 double TriFCYieldSurface ::f(const EPState *EPS) const 00075 { 00076 double xi = EPS->getStress().xi(); // functions to get Haigh-Westergard 00077 double rho = EPS->getStress().rho(); // stress invariants for 00078 double th = EPS->getStress().theta(); // explanation look Jeremic&Sture,1998 00079 00080 //double el = EPS->getScalarVar(1); // functions to select four parameters needed 00081 //double fcomp = EPS->getScalarVar(3); // to define yield surface. All are called throug 00082 //double ftens = EPS->getScalarVar(4); // EPS state so that scalar evolution laws could 00083 //double c = EPS->getScalarVar(2); // be developed for all of them 00084 00085 double a1 = pow (1.5,(1./2.) ); // constants needed 00086 double a2 = pow (6.,(1./2.) ); // to make calculation of the yield surface 00087 double a3 = pow (3.,(1./2.) ); // more transparent 00088 double a4 = 3.0*((fcomp * fcomp) - (ftens * ftens)); 00089 double a5 = (fcomp * ftens); 00090 double a6 = (el/(1.0+ el)); 00091 double m = a4*a6/a5; //factor m (friction for explanation look in article) 00092 00093 // William Warnke function 00094 double ww1 = ( 4.0 * (1.0 - el * el) * cos(th) * cos(th) + 00095 (2.0 * el - 1.0) * (2.0 * el - 1.0) ); 00096 double ww2 = pow ( ( (-4.0 * el) + (5.0 * el * el) + 00097 (4.0 * (1.0 - el * el) * cos(th) * cos(th)) ), (1./2.) ); 00098 double ww3 = 2.0 * (1.0 - el * el) * cos(th) + 00099 (2.0 * el - 1.0) * ww2; 00100 double ww = ww1/ww3; 00101 00102 00103 // actual yield surface 00104 double yield_func = (pow((a1*rho/fcomp),2.0) + 00105 (m * ( (rho/(a2*fcomp)) * ww + 00106 (xi/(a3*fcomp) ) ) ) - c); 00107 00108 return yield_func; 00109 00110 00111 } 00112 00113 //================================================================================ 00114 // tensor dF/dsigma_ij (this is the place wher a derivativ over yield function is coded) 00115 //================================================================================ 00116 00117 tensor TriFCYieldSurface ::dFods(const EPState *EPS ) const 00118 { 00119 tensor dFoverds( 2, def_dim_2, 0.0); 00120 00121 //double xi = EPS->getStress().xi(); // functions to get Haigh-Westergard 00122 double rho = EPS->getStress().rho(); // stress invariants for 00123 double th = EPS->getStress().theta(); // explanation look Jeremic&Sture,1998 00124 00125 00126 tensor DpoDs = EPS->getStress().dpoverds(); // functions to get already defined 00127 tensor DqoDs = EPS->getStress().dqoverds(); // drivatives of p, q and theta which are defined 00128 tensor DtoDs = EPS->getStress().dthetaoverds(); // in stress.cpp file. With method getStress() 00129 // it is made sure that function is filled with curent stress 00130 00131 // some parametres neccecary to define yield function and make calculation more transparent 00132 00133 double el = EPS->getScalarVar(1); // parameters called through 00134 //double c = EPS->getScalarVar(2); // definitions of 'set EPS' 00135 double b1 = pow (3.,(1./2.) ); 00136 double b2 = pow (2.,(1./2.) ); 00137 //double b3 = pow (1.5,(1./2.) ); 00138 double b4 = pow (6.,(1./2.) ); 00139 00140 double a4 = 3.0*((fcomp * fcomp) - (ftens * ftens)); // variables neccecary to define 00141 double a5 = (fcomp * ftens); // factor m 00142 double a6 = (el/(1.0+ el)); // which is an expresion for 00143 double m = a4*a6/a5; // friction parameter 00144 00145 00146 // William warnke function 00147 00148 double ww1 = 4.0*(1.0-el*el)*cos(th)*cos(th) + (2.0*el-1.0)*(2.0*el-1.0); 00149 double ww2 = 2.0*(1.0-el*el)*cos(th) + 00150 (2.0*el-1)*pow((4.0*(1.0-el*el)*cos(th)*cos(th) + 00151 5.0*el*el -4.0*el),(1./2.)); 00152 double ww = ww1/ww2; 00153 00154 //derivation of the William Warnke fuction over theta 00155 00156 // some factors enabling me to check the function. It might not be 00157 // the example of clearness but it was clear to me at the time I 00158 // wrote this and to be honest function is really ugly 00159 00160 00161 double dWW1 = (-1.0 + 2.0 * el) ; 00162 double dWW2 = (1.0 - el * el) ; 00163 double dWW3 = cos(th); 00164 double dWW4 = cos(th) * cos(th); 00165 double dWW5 = sin(th); 00166 double dWW6 = pow ( ( -4.0 * el + 5.0 * el * el + 4.0 * dWW2 * dWW4 ) , (1.0/2.0 ) ); 00167 double dWW7 = -8.0 * dWW2 * dWW3 * dWW5; 00168 double dWW8 = 2.0 * dWW2 * dWW3; 00169 double dWW9 = dWW1 * dWW6; 00170 double dWW10 = dWW1 * dWW1 + 4.0 * dWW2 * dWW4; 00171 double dWW11 = -2.0 * dWW2 * dWW5; 00172 double dWW12 = (4.0 * dWW1 * dWW2 * dWW3 * dWW5) / dWW6; 00173 double dWW13 = dWW10 * (dWW11 - dWW12); 00174 double dWW14 = (2.0 * dWW2 * dWW3) + (dWW1 * dWW6); 00175 00176 //William Warnke fuction derived over theta 00177 00178 double dWWodth = (dWW7 / (dWW8 + dWW9)) - (dWW13 / (dWW14 * dWW14)); 00179 00180 00181 double dFoverdp = ( (3.0/b1) * m/(b1*fcomp) ); // derivation of yield function over xi 00182 double dFoverdq = ((( 3.0*rho/(fcomp*fcomp)) + ( m*ww/(b4*fcomp)))* b2/b1 ); // derivation of yield function over rho 00183 double dFoverdt = ( (m*rho/(b4*fcomp) )*dWWodth ); // derivation of yield function over theta 00184 00185 dFoverds = DpoDs * dFoverdp + // factors are added to this function to take into 00186 DqoDs * dFoverdq + // account small differences between Haigh - Westergard 00187 DtoDs * dFoverdt; // coordinats and defined p an q functions. Vnesel sem jih v vrste 174-175 00188 // ker definica BJ tenzorja ne dopusca deljenje ce ni v tenzoski obliki 00189 return dFoverds; 00190 00191 00192 00193 } 00194 00195 //================================================================================ 00196 // double xi_s1 = dF/dS1 Derivative in terms of first scalar var 00197 //================================================================================ 00198 00199 //double TriFCYieldSurface ::xi_s1(const EPState *EPS ) const 00200 //{ 00201 // 00202 // double rho = EPS->getStress().rho(); // stress invariants for 00203 // double th = EPS->getStress().theta(); // explanation look Jeremic&Sture,1998 00204 // 00205 // double el = EPS->getScalarVar(1); 00206 // //double c = EPS->getScalarVar(2); 00207 // 00208 // //double a1 = pow (1.5,(1./2.) ); // constants needed 00209 // double a2 = pow (6.,(1./2.) ); // to make calculation of the yield surface 00210 // //double a3 = pow (3.,(1./2.) ); // more transparent 00211 // double a4 = 3*((fcomp * fcomp) - (ftens * ftens)); 00212 // double a5 = (fcomp * ftens); 00213 // double a6 = (el/(1+ el)); 00214 // double m = a4*a6/a5; 00215 // 00216 // //derivation of the William Warnke function over el 00217 // // 00218 // double dWWoel1 = (-1.0 + 2.0 * el) ; 00219 // double dWWoel2 = (1.0 - el * el) ; 00220 // double dWWoel3 = cos(th); 00221 // double dWWoel4 = cos(th) * cos(th); 00222 // double dWWoel5 = (-4.0 + 10.0 * el - 8.0 * el * dWWoel4); 00223 // double dWWoel6 = pow ( ( -4.0 * el + 5.0 * el * el + 4.0 * dWWoel2 * dWWoel4 ) , (1.0/2.0 ) ); 00224 // double dWWoel7 = dWWoel1 + 4.0 * dWWoel2 * dWWoel4; 00225 // double dWWoel8 = 4.0 * el * dWWoel3; 00226 // double dWWoel9 = dWWoel1 * dWWoel5/(2.0 * dWWoel6); 00227 // double dWWoel10 = 2.0 * dWWoel2 * dWWoel3 + dWWoel1 * dWWoel6; 00228 // double dWWoel11 = 4.0 * dWWoel1 - 8.0 * el * dWWoel4; 00229 // 00230 // double dWWoel = -((dWWoel7 * (dWWoel8 + dWWoel9 + 2.0 * dWWoel6))/(dWWoel10 * dWWoel10)) + (dWWoel11/dWWoel10); 00231 // 00232 // double dof1scalar = (m *(rho/(a2*fcomp)) * dWWoel ); 00233 // 00234 // return dof1scalar; 00235 // 00236 // 00237 // } 00238 00239 //================================================================================ 00240 // double xi_s2 = dF/dS2 Derivative in terms of second scalar var 00241 //================================================================================ 00242 00243 // double TriFCYieldSurface ::xi_s2( const EPState *EPS ) const 00244 // { 00245 // 00246 // return -1.0; 00247 // } 00248 00249 00250 //================================================================================ 00251 // double xi_t1 = dF/dt1 Derivative in terms of first tensorial var 00252 //================================================================================ 00253 00254 //tensor TriFCYieldSurface ::xi_t1( ) const 00255 //{ 00256 00257 00258 //} 00259 00260 //================================================================================ 00261 double TriFCYieldSurface::getfcomp() const 00262 { 00263 return fcomp; 00264 } 00265 00266 double TriFCYieldSurface::getftens() const 00267 { 00268 return ftens; 00269 } 00270 00271 double TriFCYieldSurface::getel() const 00272 { 00273 return el; 00274 } 00275 00276 double TriFCYieldSurface::get_c() const 00277 { 00278 return c; 00279 } 00280 00281 00282 OPS_Stream& operator<< (OPS_Stream& os, const TriFCYieldSurface & YS) 00283 { 00284 00285 os << "Triaxial-Failure-Criterion-for-Concrete Yield Surface Parameters: " << endln; 00286 os << "fcomp = " << YS.fcomp << endln; 00287 os << "ftens = " << YS.ftens << endln; 00288 os << "el= " << YS.el << endln; 00289 os << "c= " << YS.c << endln; 00290 os.precision(4); 00291 return os; 00292 00293 00294 } 00295 00296 #endif 00297 |