MechSys: models/subcam.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_SUBCAM_H
00023 #define MECHSYS_SUBCAM_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::Stress_p_q_S_sin3th; // Mean and deviatoric Cambridge stress invariants
00041 using Tensors::AddScaled;          // Z = a*X + b*Y
00042 using Tensors::I;                  // Second order identity tensor
00043 using Tensors::Mult;               // (Matrix) multiplication of two second order tensors
00044 using Tensors::Psd;                // P symmetric deviatoric
00045 using Tensors::IdyI;               // I dyadic I
00046 
00047 class SubCam : public GenericEP<2+2> // 2 => two internal variables + (for debug only) Evp + Evpsy
00048 {
00049 public:
00050 
00051     // Constructor
00052     SubCam(Array<REAL> const & Prms, Array<REAL> const & IniData);
00053 
00054     // Derived Methods
00055     String Name() const { return "SubCam"; }
00056 
00057     // Constants
00058     static REAL   Q_ZERO; // min value of q (used to calc derivatives)
00059     static size_t NPRMS; // number of parameters
00060     static size_t NIDAT; // number of initial data
00061 
00062     // Additional Internal Values output
00063     int  nAdditionalInternalStateValues () const { return 1; } // 1 => L
00064     void AdditionalInternalStateValues  (Array<REAL>   & IntStateVals)  const;
00065     void AdditionalInternalStateNames   (Array<String> & IntStateNames) const;
00066 
00067 private:
00068     // Data
00069     REAL _lam;
00070     REAL _kap;
00071     REAL _phics;
00072     REAL _G;
00073     REAL _c;
00074     REAL _Mcs;
00075     REAL _wcs;
00076     REAL _alp;
00077 
00078     // Methods
00079     void _calc_grads_hards(REAL    const & v  ,        // In:  Specific volume
00080                            Tensor2 const & Sig,        // In:  Stress tensor
00081                            IntVars const & z  ,        // In:  Internal variables
00082                            Tensor2       & r  ,        // Out: Plastic strain direction
00083                            Tensor2       & V  ,        // Out: Derivative of f ("loading" surface) w.r.t Sig (df_dSig)
00084                            T_dfdz        & y  ,        // Out: Derivative of f ("loading" surface) w.r.t z
00085                            HardMod       & HH ,        // Out: Hardening modulae
00086                            REAL          & cp ) const; // Out: Plastic coeficient
00087 
00088     void _calc_De(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & De) const;
00089     void _calc_Ce(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & Ce) const;
00090     void _unloading_update_int_vars();
00091     void _calc_dz(Tensors::Tensor2 const & dSig, Tensors::Tensor2 const & dEps, REAL const & dLam, HardMod const & HH, IntVars & dz) const
00092     { dz = dLam * HH; }
00093 
00094     // Local
00095     REAL _calc_M(REAL const & sin3th) const;
00096 
00097 }; // class SubCam
00098 
00099 REAL   SubCam::Q_ZERO = 1.0e-7; 
00100 size_t SubCam::NPRMS = 5;
00101 size_t SubCam::NIDAT = 5;
00102 
00103 
00105 
00106 
00107 inline SubCam::SubCam(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
00108 {
00109     // ##################################################################### Setup Parameters
00110 
00111     /* example.mat
00112 
00113         #-------------------------------------------------
00114         #            0       1      2    3               4
00115         #          lam     kap  phics    G (kgf/cm²)     c
00116         SubCam  0.0891  0.0196   31.6  100            5000
00117 
00118      */
00119 
00120     // Check number of parameters
00121     const size_t n_prms = 5;
00122     if (Prms.size()!=n_prms)
00123         throw new Fatal(_("SubCam::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     _G     = Prms[3];
00129     _c     = Prms[4];
00130 
00131     // Mohr-Coulomb sin(Phi)
00132     REAL sin_phi = sin(Util::ToRad(_phics));
00133     _Mcs = 6.0*sin_phi/(3.0-sin_phi);
00134     _wcs = pow((3.0-sin_phi)/(3.0+sin_phi),4.0);
00135     _alp = 1.0; //(_Mcs*(_Mcs-9.0)*(_Mcs-3.0)) / (9.0*(6.0-_Mcs)*(1.0-_kap/_lam));
00136     
00137     // ##################################################################### Setup Initial State
00138 
00139     // Check number of initial data
00140     if (IniData.size()!=5) // Layered analysis
00141         throw new Fatal(_("SubCam::Constructor: The number of initial data is not sufficiet (it must be equal to 5). { SigX SigY SigZ vini OCR }"));
00142 
00143     // Read initial data
00144     REAL sig_x = IniData[0]; // X stress component
00145     REAL sig_y = IniData[1]; // Y stress component
00146     REAL sig_z = IniData[2]; // Z stress component
00147             _v = IniData[3]; // Initial specific volume
00148     REAL   OCR = IniData[4]; // Over consolidation rate (not used here)
00149 
00150     // Fill stress and strain tensors
00151     REAL sq2=sqrt(2.0);
00152     _sig = sig_x, sig_y, sig_z, 0.0*sq2, 0.0*sq2, 0.0*sq2; // Mandel representation
00153     _eps = 0.0,0.0,0.0,0.0*sq2,0.0*sq2,0.0*sq2;
00154 
00155     // Calculate internal variables
00156     _unloading_update_int_vars(); // _z(0) => f pass on stress state
00157     _z(1) = _z(0)*OCR;
00158 
00159     // Evp and Evpsy
00160     _z(2) = 0.0; // Evp
00161     _z(3) = 0.0; // Evpsy
00162 
00163 } // }}}
00164 
00165 inline void SubCam::_unloading_update_int_vars() // {{{
00166 {
00167     // Invariants
00168     REAL     p,q,sin3th;
00169     Tensor2  S;
00170     Stress_p_q_S_sin3th(_sig, p,q,S,sin3th);
00171     REAL     M = _calc_M(sin3th);
00172 
00173     // Subloading surface
00174     _z(0) = p + pow(q/M,2.0)/p;
00175 } // }}}
00176 
00177 inline void SubCam::_calc_grads_hards(REAL    const & v  ,       // In:  Specific volume {{{
00178                                       Tensor2 const & Sig,       // In:  Stress tensor
00179                                       IntVars const & z  ,       // In:  Internal variables
00180                                       Tensor2       & r  ,       // Out: Plastic strain direction
00181                                       Tensor2       & V  ,       // Out: Derivative of f ("loading" surface) w.r.t Sig (df_dSig)
00182                                       T_dfdz        & y  ,       // Out: Derivative of f ("loading" surface) w.r.t z
00183                                       HardMod       & HH ,       // Out: Hardening modulae
00184                                       REAL          & cp ) const // Out: Plastic coeficient
00185 {
00186     // Invariants
00187     REAL     p,q,sin3th;
00188     Tensor2  S;
00189     Stress_p_q_S_sin3th(Sig, p,q,S,sin3th);
00190 
00191     // Gradients
00192     REAL  M = _calc_M(sin3th);
00193     REAL MM = M*M;
00194           V = (MM*(2.0*p-z(0))/3.0)*I +  3.0      *S;
00195           r = (MM*(2.0*p-z(0))/3.0)*I + (3.0*_alp)*S;
00196     if (q>Q_ZERO)
00197     {
00198         REAL ss3th = pow(sin3th,2.0);
00199         REAL cos3th;
00200         if (ss3th>1.0) cos3th = 0.0;
00201         else           cos3th = sqrt(1.0-ss3th);
00202         if (cos3th>Q_ZERO)
00203         {
00204             REAL                dM_dth = (0.75*M*(1.0-_wcs)*cos3th) / (1.0+_wcs+(_wcs-1.0)*sin3th);
00205             Tensor2 SS;         Mult(S,S,SS);
00206             Tensor2 dev_SS;     dev_SS = SS - (I * (SS(0)+SS(1)+SS(2))/3.0);
00207             Tensor2 dth_dsig; dth_dsig = (1.5/(q*q*cos3th)) * (dev_SS*(3.0/q) - S*sin3th);
00208                                     V += dth_dsig*(2.0*M*p*(p-z(0))*dM_dth);
00209                                     r += dth_dsig*(2.0*M*p*(p-z(0))*dM_dth);
00210         }
00211     }
00212     y(0) = -MM*p;
00213     y(1) = 0.0;
00214     y(2) = 0.0; // related to Evp   -- NOT used
00215     y(3) = 0.0; // related to Evpsy -- NOT used
00216 
00217     // Hardening
00218     REAL  chi = (_lam - _kap)/v;
00219     REAL    L = _c*pow(chi*log(z(0)/z(1)),2.0)/p;
00220     REAL tr_r = r(0)+r(1)+r(2);
00221         HH(0) = z(0)*(tr_r+L)/chi;
00222         HH(1) = z(1)*tr_r/chi;
00223         HH(2) = tr_r; // related to Evp   = dLam * tr_r
00224         HH(3) = L;    // related to Evpsy = dLam * L
00225 
00226     // Plastic coefficient
00227     cp = -y(0)*HH(0);
00228 } // }}}
00229 
00230 inline void SubCam::_calc_De(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & De) const // {{{
00231 {
00232     // Invariants
00233     REAL p = (Sig(0)+Sig(1)+Sig(2))/3.0;
00234 
00235     // Calculate Bulk modulus
00236     REAL K = p*v/_kap;
00237 
00238     // Calc De
00239     Tensors::AddScaled(2.0*_G,Tensors::Psd, K,Tensors::IdyI, De);
00240 
00241 } // }}}
00242 
00243 inline void SubCam::_calc_Ce(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & Ce) const // {{{
00244 {
00245     // Invariants
00246     REAL p = (Sig(0)+Sig(1)+Sig(2))/3.0;
00247 
00248     // Calculate Bulk modulus
00249     REAL K = p*v/_kap;
00250 
00251     // Calc Ce
00252     AddScaled(1.0/(2.0*_G), Psd, 1.0/(9.0*K), IdyI,  Ce);
00253 
00254 } // }}}
00255 
00256 inline REAL SubCam::_calc_M(REAL const & sin3th) const // {{{
00257 {
00258     return _Mcs*pow( 2.0*_wcs/(1.0+_wcs+(_wcs-1.0)*sin3th) ,0.25);
00259 } // }}}
00260 
00261 inline void SubCam::AdditionalInternalStateValues(Array<REAL> & IntStateVals) const // {{{
00262 {
00263     IntStateVals.resize(nAdditionalInternalStateValues());
00264     REAL   p = (_sig(0)+_sig(1)+_sig(2))/3.0;
00265     REAL chi = (_lam - _kap)/_v;
00266     IntStateVals[0] = _c*pow(chi*log(_z(0)/_z(1)),2.0)/p; // L
00267 } // }}}
00268 
00269 inline void SubCam::AdditionalInternalStateNames(Array<String> & IntStateNames) const // {{{
00270 {
00271     IntStateNames.resize(nAdditionalInternalStateValues());
00272     IntStateNames[0] = "L";
00273 } // }}}
00274 
00275 
00277 
00278 
00279 // Allocate a new SubCam model
00280 EquilibModel * SubCamMaker(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
00281 {
00282     return new SubCam(Prms, IniData);
00283 } // }}}
00284 
00285 // Register an SubCam model into ModelFactory array map
00286 int SubCamRegister() // {{{
00287 {
00288     EquilibModelFactory["SubCam"] = SubCamMaker;
00289     return 0;
00290 } // }}}
00291 
00292 // Execute the autoregistration
00293 int __SubCam_dummy_int = SubCamRegister();
00294 
00295 #endif // MECHSYS_SUBCAM_H
00296 
00297 // vim:fdm=marker