Line data    Source code

       1             : //
       2             : // This models an isotropic elastic-perfectly-plastic, non-time-dependent solid model.
       3             : // 
       4             : // The energy and derivatives are:
       5             : //
       6             : // .. math::
       7             : //    :nowrap:
       8             : //
       9             : //    \begin{gather}
      10             : //    W = \frac{1}{2}(\varepsilon - \varepsilon_p):\mathbb{C}(\varepsilon-\varepsilon_p) \\
      11             : //    DW = \mathbb{C}(\varepsilon-\varepsilon_p) \\
      12             : //    DDW = \mathbb{C}
      13             : //    \end{gather}
      14             : //
      15             : // where :math:`\mathbb{C}` is an isotropic :code:`Set::Matrix4` and :math:`\varepsilon_p` is 
      16             : // is stored in the :code:`F0` eigenstrain.
      17             : //
      18             : // The plastic strain is evolved according to the following:
      19             : //
      20             : // #. Calculate the deviatoric stress :math:`\sigma_v=\sigma - \frac{1}{3}tr(\sigma)\mathbf{I}`
      21             : // #. Calculate :math:`J_2=\sqrt{\frac{3}{2}\sigma_v:\sigma_v}`
      22             : // #. If :math:`J_2<\sigma_0` then quit - no plasticity occurs
      23             : // #. Calculate :math:`\Delta\sigma = (1-\frac{\sigma_0}{J_2})`, which projects the stress 
      24             : //    back on the yield surface.
      25             : // #. Convert to change in plastic strain, :math:`\Delta\varepsilon=\mathbb{C}^{-1}\Delta\sigma`
      26             : // #. Update plastic strain: :math:`\varepsilon_p += \Delta\varepsilon`
      27             : //
      28             : // Notes:
      29             : //
      30             : // * This does not implement any kind of hardening model. Rate hardening, isotropic hardening,
      31             : //   and kinematic hardening have yet to be implemneted.
      32             : //
      33             : #ifndef MODEL_SOLID_AFFINE_J2_H_
      34             : #define MODEL_SOLID_AFFINE_J2_H_
      35             : 
      36             : #include "AMReX.H"
      37             : #include "IO/ParmParse.H"
      38             : #include "Model/Solid/Affine/Isotropic.H"
      39             : #include "Set/Set.H"
      40             : #include <cmath>
      41             : 
      42             : namespace Model
      43             : {
      44             : namespace Solid
      45             : {
      46             : namespace Affine
      47             : {
      48             : class J2 : public Isotropic
      49             : {
      50             : public:
      51             : 
      52       33500 :     J2() {};
      53             : 
      54        7394 :     void Advance(Set::Scalar dt, Set::Matrix eps, Set::Matrix /*sig*/, Set::Scalar /*time*/) override
      55             :     {
      56        7394 :         Set::Matrix sigma = DW(eps);
      57        7394 :         Set::Matrix sig_dev = sigma - sigma.trace() * Set::Matrix::Identity() / 3.0;
      58        7394 :         if ((hardening<0.0) || (ratecoeff<0.0)) // rate independent
      59             :         {
      60        7394 :             Set::Scalar J2 = sqrt(1.5 * (sig_dev.transpose() * sig_dev).trace());
      61        7394 :             if (J2 < sigma0) return; // No plasticity happens.
      62          79 :             Set::Matrix dsigp = (1. - sigma0/J2) * sig_dev;
      63         158 :             F0 += ddw.Inverse() * dsigp; 
      64             :         }
      65             :         else // rate-dependent
      66             :         {
      67           0 :             Set::Matrix eps_dev = eps - eps.trace() * Set::Matrix::Identity() / (double)AMREX_SPACEDIM;
      68           0 :             Set::Scalar m_mu = ddw.Mu();
      69           0 :             Set::Matrix zeta_trial = sig_dev + 2.0*m_mu*eps_dev;
      70           0 :             Set::Scalar f_trial = zeta_trial.norm() - sqrt(1.0 - 1.0/(double)AMREX_SPACEDIM)*(sigma0 + hardening*alpha);
      71           0 :             if (f_trial <= 0.0) return;
      72           0 :             Set::Matrix n_new = zeta_trial / zeta_trial.norm();
      73           0 :             Set::Scalar dGamma = f_trial / (2.0*m_mu*(1 + hardening/(3.0*m_mu)));
      74           0 :             dGamma *= ratecoeff * dt;
      75           0 :             alpha = alpha + sqrt(1.0 - 1.0/(double)AMREX_SPACEDIM)*dGamma;
      76           0 :             F0 += dGamma * n_new;
      77             :         }
      78             :     }
      79             : 
      80             :     // Parameters
      81             :     Set::Scalar sigma0 = 1.0;
      82             :     Set::Scalar hardening = NAN;
      83             :     Set::Scalar ratecoeff = NAN;
      84             : 
      85             :     // Internal variables to track state
      86             :     Set::Scalar alpha = 0.0;
      87             : 
      88         290 :     static J2 Zero()
      89             :     {
      90         290 :         J2 ret;
      91         290 :         ret.ddw = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Isotropic>::Zero();
      92         290 :         ret.F0 = Set::Matrix::Zero();
      93         290 :         ret.sigma0 = 0.0;
      94         290 :         ret.hardening = 0.0;
      95         290 :         ret.ratecoeff = 0.0;
      96         290 :         return ret;
      97             :     }
      98             : 
      99             :     // See also inputs to :ref:`Model::Solid::Affine::Isotropic`
     100           1 :     static void Parse(J2 & value, IO::ParmParse & pp)
     101             :     {
     102           1 :         Isotropic::Parse(value,pp);
     103             :         // J2 Yield criterion
     104           1 :         pp_query_default("sigma0",value.sigma0,1.0);
     105             :         // Hardening coefficient (negative value disables rate hardening)
     106           1 :         pp_query_default("hardening",value.hardening,-1.0);
     107             :         // Rate coefficient (negative value disables rate hardening)
     108           1 :         pp_query_default("ratecoeff",value.ratecoeff,-1.0);
     109           1 :     }
     110             : 
     111             :     #define OP_CLASS J2
     112             :     #define OP_VARS  X(ddw) X(F0) X(sigma0) X(hardening) X(ratecoeff)
     113             :     #include "Model/Solid/InClassOperators.H"
     114             : };
     115             : #include "Model/Solid/ExtClassOperators.H"
     116             : }
     117             : }
     118             : }
     119             : 
     120             : 
     121             : template<>
     122             : AMREX_ATTRIBUTE_WEAK
     123          17 : int Set::Field<Model::Solid::Affine::J2>::NComp() const 
     124             : {
     125          17 :     return AMREX_SPACEDIM*AMREX_SPACEDIM;
     126             : }
     127             : 
     128             : template<>
     129             : AMREX_ATTRIBUTE_WEAK
     130           8 : std::string Set::Field<Model::Solid::Affine::J2>::Name(int i) const 
     131             : {
     132             : #if AMREX_SPACEDIM==2
     133           8 :     if (i==0) return name + ".epsp_xx";
     134           6 :     if (i==1) return name + ".epsp_xy";
     135           4 :     if (i==2) return name + ".epsp_yx";
     136           2 :     if (i==3) return name + ".epsp_yy";
     137             : #endif
     138             : #if AMREX_SPACEDIM==3
     139           0 :     if (i==0) return name + ".epsp_xx";
     140           0 :     if (i==1) return name + ".epsp_xy";
     141           0 :     if (i==2) return name + ".epsp_xz";
     142           0 :     if (i==3) return name + ".epsp_yx";
     143           0 :     if (i==4) return name + ".epsp_yy";
     144           0 :     if (i==5) return name + ".epsp_yz";
     145           0 :     if (i==6) return name + ".epsp_zx";
     146           0 :     if (i==7) return name + ".epsp_zy";
     147           0 :     if (i==8) return name + ".epsp_zz";
     148             : #endif
     149           0 :     return name;
     150             : }
     151             : 
     152             : template<>
     153             : AMREX_ATTRIBUTE_WEAK
     154           5 : void Set::Field<Model::Solid::Affine::J2>::Copy(int a_lev, amrex::MultiFab &a_dst, int a_dstcomp, int a_nghost) const
     155             : {
     156          25 :     for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
     157             :     {
     158          20 :         const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
     159          20 :         if (bx.ok())
     160             :         {
     161          20 :             amrex::Array4<const Model::Solid::Affine::J2> const & src = ((*this)[a_lev])->array(mfi);
     162          20 :             amrex::Array4<Set::Scalar> const & dst = a_dst.array(mfi);
     163         100 :             for (int n = 0; n < AMREX_SPACEDIM*AMREX_SPACEDIM; n++)
     164             :             {
     165          80 :                 amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
     166       20256 :                     dst(i,j,k,a_dstcomp + n) = src(i,j,k).F0(n/AMREX_SPACEDIM,n%AMREX_SPACEDIM);
     167       10128 :                 });
     168             :             }
     169             :         }
     170             :     }    
     171           5 : }
     172             : 
     173             : 
     174             : #endif