MechSys: models/subtij.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_SUBTIJ_H
00023 #define MECHSYS_SUBTIJ_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 <cmath>
00034 
00035 #include "models/genericep.h"
00036 #include "tensors/operators.h"
00037 #include "tensors/functions.h"
00038 #include "tensors/tijtensor.h"
00039 #include "util/string.h"
00040 
00041 using Tensors::Tensor2;      // Second order tensor
00042 using Tensors::Stress_tn_ts; // Prof. Nakai's Mean and deviatoric stress invariants
00043 using Tensors::Stress_q;     // Cambridge's deviatoric stress
00044 using Tensors::AddScaled;    // Z = a*X + b*Y
00045 using Tensors::I;            // Second order identity tensor
00046 using Tensors::Mult;         // (Matrix) multiplication of two second order tensors
00047 using Tensors::TijTensor;    // Second order tensor
00048 
00049 class SubTij : public GenericEP<2> // 2 => two internal variables
00050 {
00051     static REAL TIJ_ZERO;
00052     static REAL MIN_Q   ; // min q value to check for hydrostatic state
00053     static REAL DELTA_SX; // increment to Sx to avoid hydrostatic state
00054     /* // PrmSeek
00055     static const REAL TIJ_ZERO = 1.0e-3;
00056     static const REAL MIN_Q    = 1.0e-2; // min q value to check for hydrostatic state
00057     static const REAL DELTA_SX = 1.0e-1; // increment to Sx to avoid hydrostatic state
00058     */
00059 public:
00060     SubTij(Array<REAL> const & Prms, Array<REAL> const & IniData);
00061     String Name() const { return "SubTij"; }
00062 private:
00063     // Data
00064     REAL _lam;
00065     REAL _kap;
00066     REAL _nu ;
00067     REAL _Rcs;
00068     REAL _c  ;
00069     REAL _bet;
00070     REAL _Xcs;
00071     REAL _Ycs;
00072     REAL _mus;
00073 
00074     // Methods
00075     void _calc_grads_hards(REAL    const & v  ,        // In:  Specific volume
00076                            Tensor2 const & Sig,        // In:  Stress tensor
00077                            IntVars const & z  ,        // In:  Internal variables
00078                            Tensor2       & r  ,        // Out: Plastic strain direction
00079                            Tensor2       & V  ,        // Out: Derivative of f ("loading" surface) w.r.t Sig (df_dSig)
00080                            T_dfdz        & y  ,        // Out: Derivative of f ("loading" surface) w.r.t z
00081                            HardMod       & HH ,        // Out: Hardening modulae
00082                            REAL          & cp ) const; // Out: Plastic coeficient
00083 
00084     void _calc_De(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & De) const;
00085     void _calc_Ce(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & Ce) const;
00086     void _unloading_update_int_vars();
00087     void _calc_dz(Tensors::Tensor2 const & dSig, Tensors::Tensor2 const & dEps, REAL const & dLam, HardMod const & HH, IntVars & dz) const
00088     { dz = dLam * HH; }
00089 
00090 }; // class SubTij
00091 
00092 REAL SubTij::TIJ_ZERO = 1.0e-7;
00093 REAL SubTij::MIN_Q    = 1.0e-6;
00094 REAL SubTij::DELTA_SX = 1.0e-5;
00095 
00096 
00098 
00099 
00100 inline SubTij::SubTij(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
00101 {
00102     // ##################################################################### Setup Parameters
00103 
00104     /* example.mat
00105      * #------------------------------------
00106      * #           0      1    2   3   4   5
00107      * #         lam    kap   nu rcs   c bet
00108      * SubTij 0.0891 0.0196 0.20 3.5 500 1.5
00109      */
00110 
00111     // Check number of parameters
00112     const size_t n_prms = 6;
00113     if (Prms.size()!=n_prms)
00114         throw new Fatal(_("SubTij::Constructor: The number of parameters (%d) is incorrect. It must be equal to %d"), Prms.size(), n_prms);
00115 
00116     _lam = Prms[0];
00117     _kap = Prms[1];
00118     _nu  = Prms[2];
00119     _Rcs = Prms[3];
00120     _c   = Prms[4];
00121     _bet = Prms[5];
00122 
00123     // Derived constants
00124     REAL sq2   = sqrt(2.0);
00125     REAL sqRcs = sqrt(_Rcs);
00126     _Xcs = (sqRcs-1.0/sqRcs)*sq2/3.0;
00127     _Ycs = ((1.0-sqRcs)/(0.5+sqRcs))/sq2;
00128     _mus = pow(pow(_Xcs,_bet)+_Ycs*pow(_Xcs,(_bet-1.0)),1.0/_bet);
00129 
00130     // ##################################################################### Setup Initial State
00131 
00132     // Check number of initial data
00133     if (IniData.size()!=5) // Layered analysis
00134         throw new Fatal(_("SubTij::Constructor: The number of initial data is not sufficiet (it must be equal to 5). { SigX SigY SigZ vini OCR }"));
00135 
00136     // Read initial data
00137     REAL sig_x = IniData[0]; // X stress component
00138     REAL sig_y = IniData[1]; // Y stress component
00139     REAL sig_z = IniData[2]; // Z stress component
00140             _v = IniData[3]; // Initial specific volume
00141     REAL   OCR = IniData[4]; // Over consolidation rate (not used here)
00142 
00143     // Fill stress and strain tensors
00144     _sig = sig_x, sig_y, sig_z, 0.0*sq2, 0.0*sq2, 0.0*sq2; // Mandel representation
00145     _eps = 0.0,0.0,0.0,0.0*sq2,0.0*sq2,0.0*sq2;
00146 
00147     // Avoid hydrostatic state
00148     if (Stress_q(_sig)<=MIN_Q) { _sig[0] += DELTA_SX; }
00149 
00150     // Calculate internal variables
00151     _unloading_update_int_vars(); // _z(0) => f pass on stress state
00152     _z(1) = _z(0)*OCR;
00153 
00154 } // }}}
00155 
00156 inline void SubTij::_unloading_update_int_vars() // {{{
00157 {
00158     // Avoid hydrostatic state
00159     if (Stress_q(_sig)<=MIN_Q) { _sig[0] += DELTA_SX; }
00160 
00161     // Invariants
00162     REAL tn,ts;
00163     if (!Stress_tn_ts(_sig,tn,ts))
00164         throw new Fatal(_("SubTij:_unloading_update_int_vars: Stress_tn_ts(_sig,tn,ts) failed"));
00165 
00166     // Subloading surface
00167     _z(0) = tn*exp(pow(ts/(_mus*tn),_bet)/_bet); // tn1
00168 
00169 } // }}}
00170 
00171 inline void SubTij::_calc_grads_hards(REAL    const & v  ,       // In:  Specific volume {{{
00172                                       Tensor2 const & Sig,       // In:  Stress tensor
00173                                       IntVars const & z  ,       // In:  Internal variables
00174                                       Tensor2       & r  ,       // Out: Plastic strain direction
00175                                       Tensor2       & V  ,       // Out: Derivative of f ("loading" surface) w.r.t Sig (df_dSig)
00176                                       T_dfdz        & y  ,       // Out: Derivative of f ("loading" surface) w.r.t z
00177                                       HardMod       & HH ,       // Out: Hardening modulae
00178                                       REAL          & cp ) const // Out: Plastic coeficient
00179 {
00180     // Avoid hydrostatic state
00181     if (Stress_q(Sig)<=MIN_Q) { const_cast<Tensor2&>(Sig)[0] += DELTA_SX; }
00182 
00183     // tij tensor
00184     REAL      tn;      // Professor Nakai mean stress
00185     REAL      ts;      // Professor Nakai deviatoric stress
00186     REAL      Is[3];   // Characteristics invariants (I1,I2,I3)
00187     REAL      sqrt_rz; // sqrt(I3/I2)
00188     Tensor2   tau;     // Tau = sqrt(sig)
00189     Tensor2   inv_tau; // Tau^(-1)
00190     Tensor2   a;       // Professor Nakai tensor normal to SMP
00191     Tensor2   t;       // Professor Nakai tij tensor
00192     Tensor2   dtn_ds;  // derivative of tn w.r.t sig
00193     Tensor2   dts_ds;  // derivative of ts w.r.t sig
00194     Tensor2   dts_dt;  // derivative of ts w.r.t tij
00195     TijTensor tij(TIJ_ZERO);
00196     tij.Derivs(Sig, tn, ts, Is, sqrt_rz, tau, inv_tau, a, t, dtn_ds, dts_ds, dts_dt);
00197 
00198     // Gradients
00199     REAL df_dtn = (1.0-pow(ts/(_mus*tn),_bet))/tn;
00200     REAL df_dts = pow(ts,_bet-1.0)/pow(_mus*tn,_bet);
00201               r =      a*df_dtn + dts_dt*df_dts; // df_dt
00202               V = dtn_ds*df_dtn + dts_ds*df_dts; // df_ds
00203            y(0) = -1.0/z(0);                     // df_dz0
00204            y(1) = 0.0;                           // df_dz1
00205 
00206     // Hardening
00207     REAL    G = _c*pow( (_lam-_kap)*log(z(1)/z(0)) ,2.0);
00208     REAL tr_r = r(0)+r(1)+r(2);
00209     REAL  chi = (_lam - _kap)/v;
00210         HH(0) = z(0)*(tr_r+G/tn)/chi;
00211         HH(1) = z(1)*tr_r/chi;
00212 
00213     // Plastic coefficient
00214     cp = -y(0)*HH(0);
00215 
00216 } // }}}
00217 
00218 inline void SubTij::_calc_De(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & De) const // {{{
00219 {
00220     // Calculate Young modulus
00221     REAL E = (_sig(0)+_sig(1)+_sig(2)) * (1.0-2.0*_nu) * v / _kap;
00222 
00223     // Calculate elastic tangent tensor
00224     AddScaled(     E/  (1.0+_nu)                 , Tensors::IIsym ,
00225                _nu*E/( (1.0+_nu)*(1.0-2.0*_nu) ) , Tensors::IdyI  , De ); // Z = a*X + b*Y
00226 } // }}}
00227 
00228 inline void SubTij::_calc_Ce(REAL const & v, Tensors::Tensor2 const & Sig, IntVars const & z, bool IsUnloading, Tensors::Tensor4 & Ce) const // {{{
00229 {
00230     // Calculate Young modulus
00231     REAL E = (_sig(0)+_sig(1)+_sig(2)) * (1.0-2.0*_nu) * v / _kap;
00232 
00233     // Calculate elastic tangent tensor
00234     AddScaled( (1.0+_nu)/E , Tensors::IIsym ,
00235                    - _nu/E , Tensors::IdyI  , Ce ); // Z = a*X + b*Y
00236 
00237 } // }}}
00238 
00239 
00241 
00242 
00243 // Allocate a new SubTij model
00244 EquilibModel * SubTijMaker(Array<REAL> const & Prms, Array<REAL> const & IniData) // {{{
00245 {
00246     return new SubTij(Prms, IniData);
00247 } // }}}
00248 
00249 // Register an SubTij model into ModelFactory array map
00250 int SubTijRegister() // {{{
00251 {
00252     EquilibModelFactory["SubTij"] = SubTijMaker;
00253     return 0;
00254 } // }}}
00255 
00256 // Execute the autoregistration
00257 int __SubTij_dummy_int = SubTijRegister();
00258 
00259 #endif // MECHSYS_SUBTIJ_H
00260 
00261 // vim:fdm=marker