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 244445 : J2() {};
53 :
54 178966 : void Advance(Set::Scalar dt, Set::Matrix eps, Set::Matrix /*sig*/, Set::Scalar /*time*/) override
55 : {
56 178966 : Set::Matrix sigma = DW(eps);
57 178966 : Set::Matrix sig_dev = sigma - sigma.trace() * Set::Matrix::Identity() / 3.0;
58 178966 : if ((hardening<0.0) || (ratecoeff<0.0)) // rate independent
59 : {
60 178966 : Set::Scalar J2 = sqrt(1.5 * (sig_dev.transpose() * sig_dev).trace());
61 178966 : if (J2 < sigma0) return; // No plasticity happens.
62 13552 : Set::Matrix dsigp = (1. - sigma0/J2) * sig_dev;
63 27104 : F0 += ddw.Inverse() * dsigp;
64 13552 : }
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 706 : static J2 Zero()
89 : {
90 706 : J2 ret;
91 706 : ret.ddw = Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Isotropic>::Zero();
92 706 : ret.F0 = Set::Matrix::Zero();
93 706 : ret.sigma0 = 0.0;
94 706 : ret.hardening = 0.0;
95 706 : ret.ratecoeff = 0.0;
96 706 : return ret;
97 0 : }
98 :
99 : // See also inputs to :ref:`Model::Solid::Affine::Isotropic`
100 3 : static void Parse(J2 & value, IO::ParmParse & pp)
101 : {
102 3 : Isotropic::Parse(value,pp);
103 : // J2 Yield criterion
104 18 : pp_query_default("sigma0",value.sigma0,1.0);
105 : // Hardening coefficient (negative value disables rate hardening)
106 18 : pp_query_default("hardening",value.hardening,-1.0);
107 : // Rate coefficient (negative value disables rate hardening)
108 15 : pp_query_default("ratecoeff",value.ratecoeff,-1.0);
109 3 : }
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 48 : int Set::Field<Model::Solid::Affine::J2>::NComp() const
124 : {
125 48 : return AMREX_SPACEDIM*AMREX_SPACEDIM;
126 : }
127 :
128 : template<>
129 : AMREX_ATTRIBUTE_WEAK
130 30 : std::string Set::Field<Model::Solid::Affine::J2>::Name(int i) const
131 : {
132 : #if AMREX_SPACEDIM==2
133 12 : if (i==0) return name + ".epsp_xx";
134 9 : if (i==1) return name + ".epsp_xy";
135 6 : if (i==2) return name + ".epsp_yx";
136 3 : if (i==3) return name + ".epsp_yy";
137 : #endif
138 : #if AMREX_SPACEDIM==3
139 18 : if (i==0) return name + ".epsp_xx";
140 16 : if (i==1) return name + ".epsp_xy";
141 14 : if (i==2) return name + ".epsp_xz";
142 12 : if (i==3) return name + ".epsp_yx";
143 10 : if (i==4) return name + ".epsp_yy";
144 8 : if (i==5) return name + ".epsp_yz";
145 6 : if (i==6) return name + ".epsp_zx";
146 4 : if (i==7) return name + ".epsp_zy";
147 2 : if (i==8) return name + ".epsp_zz";
148 : #endif
149 0 : return name;
150 : }
151 :
152 : template<>
153 : AMREX_ATTRIBUTE_WEAK
154 8 : void Set::Field<Model::Solid::Affine::J2>::Copy(int a_lev, amrex::MultiFab &a_dst, int a_dstcomp, int a_nghost) const
155 : {
156 31 : for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
157 : {
158 23 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
159 23 : if (bx.ok())
160 : {
161 23 : amrex::Array4<const Model::Solid::Affine::J2> const & src = ((*this)[a_lev])->array(mfi);
162 23 : amrex::Array4<Set::Scalar> const & dst = a_dst.array(mfi);
163 125 : for (int n = 0; n < AMREX_SPACEDIM*AMREX_SPACEDIM; n++)
164 : {
165 102 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
166 27068 : dst(i,j,k,a_dstcomp + n) = src(i,j,k).F0(n/AMREX_SPACEDIM,n%AMREX_SPACEDIM);
167 13534 : });
168 : }
169 : }
170 8 : }
171 8 : }
172 :
173 :
174 : #endif
|