Tri_a_fail_crit_YS.cpp

Go 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 

Generated on Mon Oct 23 15:05:17 2006 for OpenSees by doxygen 1.5.0