Line data Source code
1 : //
2 : // A basic and relatively untested implementation of finite-deformation FCC crystal plasticity.
3 : // Use with caution.
4 : //
5 : // Inherits from :ref:`PseudoLinear::Cubic` to provide elastic response.
6 : //
7 : // Plastic flow modeled as:
8 : //
9 : // .. math::
10 : //
11 : // \mathbf{F}^p\mathbf{F}^{p-1} =
12 : // \sum \dot{\gamma}_n \mathbf{a}_n\otimes\mathbf{N}_n
13 : //
14 : // where :math:`\mathbf{F}=\mathbf{F}^e\mathbf{F}^p` and :math:`\gamma_n` are
15 : // slips on system :math:`n`.
16 : //
17 : // Power law rate hardening is used in the integration of slips:
18 : //
19 : // .. math::
20 : //
21 : // \dot{\gamma}_n = \dot{\gamma}_0 \Big(\frac{\tau}{\tau_{crss}}\Big)^m \operatorname{sign}(\tau)
22 : //
23 : //
24 :
25 : #ifndef MODEL_SOLID_FINITE_CRYSTALPLASTIC_H_
26 : #define MODEL_SOLID_FINITE_CRYSTALPLASTIC_H_
27 :
28 :
29 : #include "IO/ParmParse.H"
30 : #include "Model/Solid/Solid.H"
31 : #include "PseudoLinear/Cubic.H"
32 :
33 : namespace Model
34 : {
35 : namespace Solid
36 : {
37 : namespace Finite
38 : {
39 : class CrystalPlastic : public PseudoLinear::Cubic
40 : {
41 : public:
42 0 : CrystalPlastic() {};
43 :
44 0 : virtual ~CrystalPlastic() {};
45 :
46 0 : Set::Scalar W(const Set::Matrix & F) const override
47 : {
48 0 : return PseudoLinear::Cubic::W(F * Set::reduce(Fp).inverse());
49 : }
50 0 : Set::Matrix DW(const Set::Matrix & F) const override
51 : {
52 0 : Set::Matrix Fpinv = Set::reduce(Fp).inverse();
53 0 : return PseudoLinear::Cubic::DW(F * Fpinv) * Fpinv.transpose();
54 : }
55 0 : Set::Matrix4<AMREX_SPACEDIM,Set::Sym::Major> DDW(const Set::Matrix & F) const override
56 : {
57 0 : Set::Matrix Fpinv = Set::reduce(Fp).inverse();
58 0 : auto DDW = PseudoLinear::Cubic::DDW(F * Fpinv);
59 0 : auto ret = Set::Matrix4<AMREX_SPACEDIM, Set::Sym::Major>::Zero();
60 0 : for (int i = 0; i < AMREX_SPACEDIM; i++)
61 0 : for (int J = 0; J < AMREX_SPACEDIM; J++)
62 0 : for (int k = 0; k < AMREX_SPACEDIM; k++)
63 0 : for (int L = 0; L < AMREX_SPACEDIM; L++)
64 : {
65 0 : ret(i,J,k,L) = 0.0;
66 0 : for (int Q = 0; Q < AMREX_SPACEDIM; Q++)
67 0 : for (int R = 0; R < AMREX_SPACEDIM; R++)
68 0 : ret(i,J,k,L) += DDW(i,Q,k,R)*Fpinv(J,Q)*Fpinv(L,R);
69 : }
70 0 : return ret;
71 : }
72 :
73 0 : virtual void Advance(Set::Scalar dt, Set::Matrix /*F*/, Set::Matrix P, Set::Scalar time) override
74 : {
75 0 : if (time < tstart) return;
76 :
77 0 : q.normalize();
78 0 : Set::Matrix3d R = q.toRotationMatrix();
79 0 : std::array<std::pair<Set::Vector3d,Set::Vector3d>,12> ss = slipSystems(R);
80 :
81 0 : Set::Matrix3d L = Set::Matrix3d::Zero();
82 :
83 0 : for (int n = 0; n < 12; n++)
84 : {
85 0 : Set::Vector3d A = ss[n].first;
86 0 : Set::Vector3d N = ss[n].second;
87 0 : Set::Scalar tau = A.dot(Set::expand(P)*N);
88 0 : Set::Scalar sign = tau>0 ? 1.0 : -1.0;
89 :
90 : Set::Scalar gammadot
91 0 : = gammadot0 * std::pow(fabs(tau/tau_crss[n]),m_rate_inv) * sign;
92 :
93 0 : L += gammadot * A * N.transpose();
94 :
95 0 : gamma[n] += gammadot * dt;
96 : }
97 :
98 0 : Fp += L * Fp * dt;
99 : }
100 :
101 :
102 : public:
103 : //Set::Matrix4<AMREX_SPACEDIM,Set::Sym::MajorMinor> CC;
104 : //static constexpr KinematicVariable kinvar = KinematicVariable::F;
105 : //Set::Scalar C11=1.68, C12=1.21, C44=0.75;
106 : //mutable Set::Quaternion q = Set::Quaternion(1.0, 0.0, 0.0, 0.0);
107 :
108 : // Plastic state variables
109 : Set::Matrix3d Fp = Set::Matrix3d::Identity();
110 : Eigen::Matrix<Set::Scalar,12,1> gamma = Eigen::Matrix<Set::Scalar,12,1>::Zero();
111 :
112 : // Parameters
113 : Set::Scalar m_rate_inv = NAN;
114 : //Set::Scalar tau_crss = NAN;
115 : Eigen::Matrix<Set::Scalar,12,1> tau_crss = Eigen::Matrix<Set::Scalar,12,1>::Zero();
116 : Set::Scalar gammadot0 = NAN;
117 :
118 : // Control variables
119 : Set::Scalar tstart = NAN;
120 :
121 : private:
122 0 : static std::array<std::pair<Set::Vector3d,Set::Vector3d>,12> slipSystems(Set::Matrix3d R)
123 : {
124 : const std::array<Set::Vector3d,12> n0 =
125 : {
126 0 : Set::Vector3d::Constant( 1, 1, 1) / sqrt(3.0),
127 0 : Set::Vector3d::Constant( 1, 1, 1) / sqrt(3.0),
128 0 : Set::Vector3d::Constant( 1, 1, 1) / sqrt(3.0),
129 0 : Set::Vector3d::Constant(-1, 1,-1) / sqrt(3.0),
130 0 : Set::Vector3d::Constant(-1, 1,-1) / sqrt(3.0),
131 0 : Set::Vector3d::Constant(-1, 1,-1) / sqrt(3.0),
132 0 : Set::Vector3d::Constant( 1,-1,-1) / sqrt(3.0),
133 0 : Set::Vector3d::Constant( 1,-1,-1) / sqrt(3.0),
134 0 : Set::Vector3d::Constant( 1,-1,-1) / sqrt(3.0),
135 0 : Set::Vector3d::Constant(-1,-1, 1) / sqrt(3.0),
136 0 : Set::Vector3d::Constant(-1,-1, 1) / sqrt(3.0),
137 0 : Set::Vector3d::Constant(-1,-1, 1) / sqrt(3.0)
138 0 : };
139 : const std::array<Set::Vector3d,12> a0 =
140 : {
141 0 : Set::Vector3d::Constant( 1, 0,-1) / sqrt(2.0), // 1 1 1
142 0 : Set::Vector3d::Constant( 0,-1, 1) / sqrt(2.0), // 1 1 1
143 0 : Set::Vector3d::Constant( 1,-1, 0) / sqrt(2.0), // 1 1 1
144 0 : Set::Vector3d::Constant( 1, 0,-1) / sqrt(2.0), // -1 1 -1
145 0 : Set::Vector3d::Constant( 1, 1, 0) / sqrt(2.0), // -1 1 -1
146 0 : Set::Vector3d::Constant( 0, 1, 1) / sqrt(2.0), // -1 1 -1
147 0 : Set::Vector3d::Constant( 1, 1, 0) / sqrt(2.0), // 1 -1 -1
148 0 : Set::Vector3d::Constant( 0,-1, 1) / sqrt(2.0), // 1 -1 -1
149 0 : Set::Vector3d::Constant( 1, 0, 1) / sqrt(2.0), // 1 -1 -1
150 0 : Set::Vector3d::Constant( 0, 1, 1) / sqrt(2.0), // -1 -1 1
151 0 : Set::Vector3d::Constant( 1, 0, 1) / sqrt(2.0), // -1 -1 1
152 0 : Set::Vector3d::Constant( 1,-1, 0) / sqrt(2.0) // -1 -1 1
153 0 : };
154 :
155 :
156 0 : static std::array<std::pair<Set::Vector3d,Set::Vector3d>,12> ret;
157 0 : for (int n = 0 ; n < 12 ; n ++)
158 : {
159 0 : ret[n].first = R*a0[n];
160 0 : ret[n].second = R*n0[n];
161 : }
162 0 : return ret;
163 : }
164 :
165 : public:
166 0 : static CrystalPlastic Zero()
167 : {
168 0 : CrystalPlastic ret;
169 0 : ret.C11 = 0.0;
170 0 : ret.C12 = 0.0;
171 0 : ret.C44 = 0.0;
172 0 : ret.q = Set::Quaternion(0.0,0.0,0.0,0.0);
173 0 : ret.Fp = Set::Matrix3d::Zero();
174 0 : ret.gamma = Eigen::Matrix<Set::Scalar,12,1>::Zero();
175 0 : ret.m_rate_inv = 0.0;
176 0 : ret.tau_crss = Eigen::Matrix<Set::Scalar,12,1>::Zero();
177 0 : ret.gammadot0 = 0.0;
178 0 : return ret;
179 : }
180 : static CrystalPlastic Random()
181 : {
182 : return Random(Util::Random(), Util::Random(), Util::Random());
183 : }
184 : static CrystalPlastic Random(Set::Scalar C11, Set::Scalar C12, Set::Scalar C44)
185 : {
186 : CrystalPlastic ret = CrystalPlastic::Zero();
187 : ret.C11 = C11;
188 : ret.C12 = C12;
189 : ret.C44 = C44;
190 : ret.q = Set::Quaternion::UnitRandom();
191 : return ret;
192 : }
193 0 : static void Parse(CrystalPlastic & value, IO::ParmParse & pp)
194 : {
195 0 : PseudoLinear::Cubic::Parse(value,pp);
196 : // TODO Add inputs for the CP values
197 :
198 0 : std::vector<Set::Scalar> tau_crss;
199 : // Critical resolved shear stress :math:`\tau_{crss}`
200 0 : pp_queryarr("tau_crss",tau_crss);
201 0 : if (tau_crss.size() == 1)
202 0 : for (int n = 0; n < 12; n++) value.tau_crss[n] = tau_crss[0];
203 0 : else if (tau_crss.size() == 12)
204 0 : for (int n = 0; n < 12; n++) value.tau_crss[n] = tau_crss[n];
205 : else
206 0 : Util::Exception(INFO,"Expected either 1 or 12 values for tau_crss but got ",tau_crss.size());
207 :
208 : // Rate hardening coefficient :math:`\dot{\gamma}_0`
209 0 : pp_query_default("gammadot0",value.gammadot0,1.0);
210 :
211 : // Inverse of the hardening exponent :math:`\frac{1}{m}`
212 0 : pp_query_default("m_rate_inv",value.m_rate_inv,0.5);
213 :
214 : // Time to activate plastic slip
215 0 : pp_query_default("tstart",value.tstart,0.0);
216 0 : }
217 :
218 : #define OP_CLASS CrystalPlastic
219 : #define OP_VARS X(C11) X(C12) X(C44) X(q) X(Fp) X(gamma) X(m_rate_inv) X(tau_crss) X(gammadot0) X(tstart)
220 : #include "Model/Solid/InClassOperators.H"
221 : };
222 : #include "Model/Solid/ExtClassOperators.H"
223 :
224 :
225 :
226 :
227 :
228 : }
229 : }
230 : }
231 :
232 :
233 : template<>
234 0 : inline int Set::Field<Model::Solid::Finite::CrystalPlastic>::NComp() const
235 : {
236 0 : return 12;
237 : }
238 :
239 : template<>
240 0 : inline std::string Set::Field<Model::Solid::Finite::CrystalPlastic>::Name(int i) const
241 : {
242 0 : if (i==0) return name + "gamma1";
243 0 : if (i==1) return name + "gamma2";
244 0 : if (i==2) return name + "gamma3";
245 0 : if (i==3) return name + "gamma4";
246 0 : if (i==4) return name + "gamma5";
247 0 : if (i==5) return name + "gamma6";
248 0 : if (i==6) return name + "gamma7";
249 0 : if (i==7) return name + "gamma8";
250 0 : if (i==8) return name + "gamma9";
251 0 : if (i==9) return name + "gamma10";
252 0 : if (i==10) return name + "gamma11";
253 0 : if (i==11) return name + "gamma12";
254 0 : return name;
255 : }
256 :
257 : template<>
258 0 : inline void Set::Field<Model::Solid::Finite::CrystalPlastic>::Copy(int a_lev, amrex::MultiFab &a_dst, int a_dstcomp, int a_nghost) const
259 : {
260 0 : for (amrex::MFIter mfi(a_dst, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi)
261 : {
262 0 : const amrex::Box& bx = mfi.growntilebox(amrex::IntVect(a_nghost));
263 0 : if (bx.ok())
264 : {
265 0 : amrex::Array4<const Model::Solid::Finite::CrystalPlastic> const & src = ((*this)[a_lev])->array(mfi);
266 0 : amrex::Array4<Set::Scalar> const & dst = a_dst.array(mfi);
267 0 : for (int n = 0; n < 12; n++)
268 : {
269 0 : amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) {
270 0 : dst(i,j,k,a_dstcomp + n) = src(i,j,k).gamma(n);
271 0 : });
272 : }
273 : }
274 : }
275 0 : }
276 :
277 :
278 :
279 : #endif
|
