MechSys: models/subcamx.h Source File

00001 /*************************************************************************************
00002  * MechSys - A C++ library to simulate (Continuum) Mechanical Systems                *
00003  * Copyright (C) 2005 Dorival de Moraes Pedroso <dorival.pedroso at gmail.com>       *
00004  * Copyright (C) 2005 Raul Dario Durand Farfan  <raul.durand at gmail.com>           *
00005  *                                                                                   *
00006  * This file is part of MechSys.                                                     *
00007  *                                                                                   *
00008  * MechSys is free software; you can redistribute it and/or modify it under the      *
00009  * terms of the GNU General Public License as published by the Free Software         *
00010  * Foundation; either version 2 of the License, or (at your option) any later        *
00011  * version.                                                                          *
00012  *                                                                                   *
00013  * MechSys is distributed in the hope that it will be useful, but WITHOUT ANY        *
00014  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A   *
00015  * PARTICULAR PURPOSE. See the GNU General Public License for more details.          *
00016  *                                                                                   *
00017  * You should have received a copy of the GNU General Public License along with      *
00018  * MechSys; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, *
00019  * Fifth Floor, Boston, MA 02110-1301, USA                                           *
00020  *************************************************************************************/
00021 
00022 #ifndef MECHSYS_SUBCAMX_H
00023 #define MECHSYS_SUBCAMX_H
00024 
00025 #ifdef HAVE_CONFIG_H
00026   #include "config.h"
00027 #else
00028   #ifndef REAL
00029     #define REAL double
00030   #endif
00031 #endif
00032 
00033 #include "models/genericep.h"
00034 #include "tensors/operators.h"
00035 #include "tensors/functions.h"
00036 #include "util/string.h"
00037 #include "util/util.h"
00038 
00039 using Tensors::Tensor2;             // Second order tensor
00040 using Tensors::Tensor4;             // Fourth order tensor
00041 using Tensors::Stress_p_q;          // Mean and deviatoric Cambridge stress invariants
00042 using Tensors::Stress_p_q_S_sin3th; // Mean and deviatoric Cambridge stress invariants
00043 using Tensors::AddScaled;           // Z = a*X + b*Y
00044 using Tensors::Dyad;                // Dyadic product
00045 using Tensors::Mult;                // (Matrix) multiplication of two second order tensors
00046 using Tensors::Psd;                 // P symmetric deviatoric
00047 using Tensors::IdyI;                // P symmetric deviatoric
00048 using Tensors::I;                   // Second order identity tensor
00049 using Util::ToRad;
00050 using Util::Sq2;
00051 
00052 class SubCamX : public GenericEP<3> // 3 => two internal variables
00053 {
00054     static REAL Q_ZERO; // min value of q (used to calc derivatives)
00055 public:
00056     SubCamX(Array<REAL> const & Prms, Array<REAL> const & IniData);
00057     String Name() const { return "SubCamX"; }
00058 private:
00059     // Data
00060     REAL _lam;
00061     REAL _kap;
00062     REAL _phics;
00063     REAL _c;
00064     REAL _G0;
00065     REAL _psiv0;
00066     REAL _A;
00067     REAL _B;
00068     REAL _alp;
00069     REAL _bet;
00070     REAL _C;
00071 
00072     // Derived
00073     REAL _Mcs;
00074     REAL _wcs;
00075 
00076     // Methods
00077     void _calc_grads_hards(REAL    const & v  ,        // In:  Specific volume
00078                            Tensor2 const & Sig,        // In:  Stress tensor
00079                            IntVars const & z  ,        // In:  Internal variables
00080                            Tensor2       & r  ,        // Out: Plastic strain direction
00081                            Tensor2       & V  ,        // Out: Derivative of f ("loading" surface) w.r.t Sig (df_dSig)
00082                            T_dfdz        & y  ,        // Out: Derivative of f ("loading" surface) w.r.t z
00083                            HardMod       & HH ,        // Out: Hardening modulae
00084                            REAL          & cp ) const; // Out: Plastic coeficient
00085 
00086     void _calc_De(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & De) const;
00087     void _calc_Ce(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & Ce) const;
00088     void _unloading_update_int_vars();
00089     void _calc_dz(Tensors::Tensor2 const & dSig, Tensors::Tensor2 const & dEps, REAL const & dLam, HardMod const & HH, IntVars & dz) const
00090     {
00091         dz(0) = HH(0) *      dLam;
00092         dz(1) = HH(1) *      dLam;
00093         dz(2) = HH(2) * fabs(dLam);
00094     }
00095 
00096     // Local
00097     REAL _calc_M(REAL const & sin3th) const;
00098 
00099 }; // class SubCamX
00100 
00101 REAL SubCamX::Q_ZERO = 1.0e-7;
00102 
00103 
00105 
00106 
00107 inline SubCamX::SubCamX(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
00108 {
00109     // ##################################################################### Setup Parameters
00110 
00111     /* example.mat
00112                                   (kgf/cm²)     (kgf/cm²)
00113        #-----------------------------------------------------------------
00114        #            0      1      2    3      4      5   6  7  8    9  10
00115        #          lam    kap  phics   G0      c  psiv0   A  B  C  alp bet
00116        SubCamX 0.0891 0.0196   31.6  100  14000    500   1  1  1    1   1
00117       
00118      */
00119 
00120     // Check number of parameters
00121     const size_t n_prms = 11;
00122     if (Prms.size()!=n_prms)
00123         throw new Fatal(_("SubCamX::Constructor: The number of parameters (%d) is incorrect. It must be equal to %d"), Prms.size(), n_prms);
00124 
00125     _lam   = Prms[0];
00126     _kap   = Prms[1];
00127     _phics = Prms[2];
00128     _G0    = Prms[3];
00129     _c     = Prms[4];
00130     _psiv0 = Prms[5];
00131     _A     = Prms[6];  // decay of L
00132     _B     = Prms[7]; 
00133     _C     = Prms[8];  // non-associated flow
00134     _alp   = Prms[9];  // hysteresis
00135     _bet   = Prms[10]; // hysteresis
00136 
00137     // Mohr-Coulomb sin(Phi)
00138     REAL sin_phi = sin(Util::ToRad(_phics));
00139     _Mcs = 6.0*sin_phi/(3.0-sin_phi);
00140     _wcs = pow((3.0-sin_phi)/(3.0+sin_phi),4.0);
00141     
00142     // ##################################################################### Setup Initial State
00143 
00144     // Check number of initial data
00145     if (IniData.size()!=5) // Layered analysis
00146         throw new Fatal(_("SubCamX::Constructor: The number of initial data is not sufficiet (it must be equal to 5). { SigX SigY SigZ vini OCR }"));
00147 
00148     // Read initial data
00149     REAL sig_x = IniData[0]; // X stress component
00150     REAL sig_y = IniData[1]; // Y stress component
00151     REAL sig_z = IniData[2]; // Z stress component
00152             _v = IniData[3]; // Initial specific volume
00153     REAL   OCR = IniData[4]; // Over consolidation rate (not used here)
00154 
00155     // Fill stress and strain tensors
00156     _sig = sig_x, sig_y, sig_z, 0.0*Sq2(), 0.0*Sq2(), 0.0*Sq2(); // Mandel representation
00157     _eps = 0.0,0.0,0.0,0.0*Sq2(),0.0*Sq2(),0.0*Sq2();
00158 
00159     // Calculate internal variables
00160     _unloading_update_int_vars(); // _z(0) => f pass on stress state
00161     _z(1) = _z(0)*OCR;
00162     _z(2) = 0.0; // max deviatoric stress "recently" applied (used on uloading only)
00163 
00164 } // }}}
00165 
00166 inline void SubCamX::_unloading_update_int_vars() // {{{
00167 {
00168     // Invariants
00169     REAL     p,q,sin3th;
00170     Tensor2  S;
00171     Stress_p_q_S_sin3th(_sig, p,q,S,sin3th);
00172     REAL     M = _calc_M(sin3th);
00173 
00174     // Subloading surface
00175     _z(0) = p + pow(q/M,2.0)/p;
00176 } // }}}
00177 
00178 inline void SubCamX::_calc_grads_hards(REAL    const & v  ,       // In:  Specific volume {{{
00179                                        Tensor2 const & Sig,       // In:  Stress tensor
00180                                        IntVars const & z  ,       // In:  Internal variables
00181                                        Tensor2       & r  ,       // Out: Plastic strain direction
00182                                        Tensor2       & V  ,       // Out: Derivative of f ("loading" surface) w.r.t Sig (df_dSig)
00183                                        T_dfdz        & y  ,       // Out: Derivative of f ("loading" surface) w.r.t z
00184                                        HardMod       & HH ,       // Out: Hardening modulae
00185                                        REAL          & cp ) const // Out: Plastic coeficient
00186 {
00187     // Invariants
00188     REAL     p,q,sin3th;
00189     Tensor2  S;
00190     Stress_p_q_S_sin3th(Sig, p,q,S,sin3th);
00191 
00192     // Allow Real values only
00193     // (This is IMPORTANT when using prs program because it may become a lot
00194     // slower if the following "if" statements would not be properly computed)
00195     if (Util::IsNanOrInf(q))
00196         throw new Fatal("SubCamX::_calc_grads_hards: q is NAN or INF");
00197 
00198     // Gradients
00199     REAL  M = _calc_M(sin3th);
00200     REAL MM = M*M;
00201           V = (MM*(2.0*p-z(0))/3.0)*I +  3.0    *S;
00202           r = (MM*(2.0*p-z(0))/3.0)*I + (3.0*_C)*S;
00203     if (q>Q_ZERO) // if q is Nan of Inf, prs may become a lot slower
00204     {
00205         REAL  ss3th = pow(sin3th,2.0);
00206         REAL cos3th = (ss3th>1.0 ? 0.0 : sqrt(1.0-ss3th));
00207         if (cos3th>Q_ZERO) // if q is Nan of Inf, prs may become a lot slower
00208         {
00209             REAL                dM_dth = (0.75*M*(1.0-_wcs)*cos3th)/(1.0+_wcs-(1.0-_wcs)*sin3th);
00210             Tensor2 SS;         Mult(S,S,SS);
00211             Tensor2 dev_SS;     dev_SS = SS - (I * (SS(0)+SS(1)+SS(2))/3.0);
00212             Tensor2 dth_dsig; dth_dsig = (1.5/(q*q*cos3th)) * (dev_SS*(3.0/q) - S*sin3th);
00213                                     V += dth_dsig*(2.0*M*p*(p-z(0))*dM_dth);
00214                                     r += dth_dsig*(2.0*M*p*(p-z(0))*dM_dth);
00215         }
00216     }
00217     y(0) = -MM*p;
00218     y(1) = 0.0;
00219     y(2) = 0.0;
00220 
00221     // Hardening
00222     REAL  chi = (_lam - _kap)/v;
00223     REAL    L = (_c*pow(chi*log(z(1)/z(0)),2.0)/p) * exp(_A*z(2)/p);
00224     REAL tr_r = r(0)+r(1)+r(2);
00225         HH(0) = z(0)*(tr_r+L)/chi;
00226         HH(1) = z(1)*tr_r/chi;
00227     
00228     Tensor2 dev_r; dev_r = r - (tr_r/3.0)*I;
00229     //HH(2) = (1.0+z(2))*Norm(dev_r);
00230     HH(2) = exp(z(2))*Norm(dev_r);
00231 
00232     // Plastic coefficient
00233     cp = -y(0)*HH(0);
00234 
00235     // Hysteresis
00236     REAL alp = _alp*log(z(1)/z(0));
00237     REAL bet = _bet*log(z(1)/z(0));
00238     r = (1.0-alp/(1.0+q/p))*dev_r + ((1.0-bet/p)*tr_r/3.0)*I;
00239 
00240 } // }}}
00241 
00242 inline void SubCamX::_calc_De(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & De) const // {{{
00243 {
00244     /*
00245     // Calculate Young modulus
00246     REAL E = (_sig(0)+_sig(1)+_sig(2)) * (1.0-2.0*_nu) * v / _kap;
00247     //REAL E = 4591.0;
00248 
00249     // Calculate elastic tangent tensor
00250     AddScaled(     E/  (1.0+_nu)                 , Tensors::IIsym ,
00251                _nu*E/( (1.0+_nu)*(1.0-2.0*_nu) ) , Tensors::IdyI  , De ); // Z = a*X + b*Y
00252     */
00253     throw new Fatal(_("SubCamX::_calc_De ERROR"));
00254 } // }}}
00255 
00256 inline void SubCamX::_calc_Ce(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & Ce) const // {{{
00257 {
00258     // Invariants
00259     REAL     p,q,sin3th;
00260     Tensor2  S;
00261     Stress_p_q_S_sin3th(Sig, p,q,S,sin3th);
00262     Tensor4  IdyS;  Dyad(I, S, IdyS);
00263 
00264     // Calculate Bulk modulus
00265     REAL K = p*v/_kap;
00266 
00267     if (IsUnloading)
00268     {
00269         AddScaled(0.5/(_G0*(1.0+q/p)), Psd, 1.0/(9.0*K), IdyI,  Ce);
00270         if (q>Q_ZERO)
00271         {
00272             REAL psiv = _psiv0*log(z(1)/z(0))*z(2)*exp(-_B*z(2));
00273             AddScaled(1.0, Ce, -psiv/(2.0*q), IdyS, Ce);
00274         }
00275     }
00276     else
00277         AddScaled(0.5/_G0, Psd, 1.0/(9.0*K), IdyI,  Ce);
00278 
00279 } // }}}
00280 
00281 inline REAL SubCamX::_calc_M(REAL const & sin3th) const // {{{
00282 {
00283     return _Mcs*pow( 2.0*_wcs/(1.0+_wcs-(1.0-_wcs)*sin3th) ,0.25);
00284 } // }}}
00285 
00287 
00288 
00289 // Allocate a new SubCamX model
00290 EquilibModel * SubCamXMaker(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
00291 {
00292     return new SubCamX(Prms, IniData);
00293 } // }}}
00294 
00295 // Register an SubCamX model into ModelFactory array map
00296 int SubCamXRegister() // {{{
00297 {
00298     EquilibModelFactory["SubCamX"] = SubCamXMaker;
00299     return 0;
00300 } // }}}
00301 
00302 // Execute the autoregistration
00303 int __SubCamX_dummy_int = SubCamXRegister();
00304 
00305 #endif // MECHSYS_SUBCAMX_H
00306 
00307 // vim:fdm=marker