Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://www.mooseframework.org
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "CappedWeakPlaneCosseratStressUpdate.h"
11 :
12 : #include "libmesh/utility.h"
13 :
14 : registerMooseObject("TensorMechanicsApp", CappedWeakPlaneCosseratStressUpdate);
15 :
16 : InputParameters
17 168 : CappedWeakPlaneCosseratStressUpdate::validParams()
18 : {
19 168 : InputParameters params = CappedWeakPlaneStressUpdate::validParams();
20 168 : params.addClassDescription("Capped weak-plane plasticity Cosserat stress calculator");
21 168 : return params;
22 0 : }
23 :
24 126 : CappedWeakPlaneCosseratStressUpdate::CappedWeakPlaneCosseratStressUpdate(
25 126 : const InputParameters & parameters)
26 126 : : CappedWeakPlaneStressUpdate(parameters)
27 : {
28 126 : }
29 :
30 : void
31 3248 : CappedWeakPlaneCosseratStressUpdate::setStressAfterReturn(const RankTwoTensor & stress_trial,
32 : Real p_ok,
33 : Real q_ok,
34 : Real gaE,
35 : const std::vector<Real> & /*intnl*/,
36 : const yieldAndFlow & smoothed_q,
37 : const RankFourTensor & Eijkl,
38 : RankTwoTensor & stress) const
39 : {
40 3248 : stress = stress_trial;
41 3248 : stress(2, 2) = p_ok;
42 : // stress_xx and stress_yy are sitting at their trial-stress values
43 : // so need to bring them back via Poisson's ratio
44 3248 : stress(0, 0) -= Eijkl(2, 2, 0, 0) * gaE / _Epp * smoothed_q.dg[0];
45 3248 : stress(1, 1) -= Eijkl(2, 2, 1, 1) * gaE / _Epp * smoothed_q.dg[0];
46 3248 : if (_in_q_trial == 0.0)
47 672 : stress(0, 2) = stress(1, 2) = 0.0;
48 : else
49 : {
50 2576 : stress(0, 2) = _in_trial02 * q_ok / _in_q_trial;
51 2576 : stress(1, 2) = _in_trial12 * q_ok / _in_q_trial;
52 : }
53 3248 : }
54 :
55 : void
56 80 : CappedWeakPlaneCosseratStressUpdate::consistentTangentOperator(
57 : const RankTwoTensor & /*stress_trial*/,
58 : Real /*p_trial*/,
59 : Real q_trial,
60 : const RankTwoTensor & /*stress*/,
61 : Real /*p*/,
62 : Real q,
63 : Real gaE,
64 : const yieldAndFlow & smoothed_q,
65 : const RankFourTensor & Eijkl,
66 : bool compute_full_tangent_operator,
67 : RankFourTensor & cto) const
68 : {
69 80 : cto = Eijkl;
70 80 : if (!compute_full_tangent_operator)
71 : return;
72 :
73 80 : const Real Ezzzz = Eijkl(2, 2, 2, 2);
74 80 : const Real Exzxz = Eijkl(0, 2, 0, 2);
75 80 : const Real tanpsi = _tan_psi.value(_intnl[_qp][0]);
76 80 : const Real dintnl0_dq = -1.0 / Exzxz;
77 : const Real dintnl0_dqt = 1.0 / Exzxz;
78 80 : const Real dintnl1_dp = -1.0 / Ezzzz;
79 : const Real dintnl1_dpt = 1.0 / Ezzzz;
80 : const Real dintnl1_dq =
81 80 : tanpsi / Exzxz - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dq / Exzxz;
82 : const Real dintnl1_dqt =
83 80 : -tanpsi / Exzxz - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dqt / Exzxz;
84 :
85 320 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
86 : {
87 240 : const Real dpt_depii = Eijkl(2, 2, i, i);
88 240 : cto(2, 2, i, i) = _dp_dpt * dpt_depii;
89 : const Real poisson_effect =
90 240 : Eijkl(2, 2, 0, 0) / Ezzzz *
91 240 : (_dgaE_dpt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dpt +
92 240 : gaE * smoothed_q.d2g[0][1] * _dq_dpt +
93 240 : gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dq * _dq_dpt) +
94 240 : gaE * smoothed_q.d2g_di[0][1] *
95 240 : (dintnl1_dpt + dintnl1_dp * _dp_dpt + dintnl1_dq * _dq_dpt)) *
96 240 : dpt_depii;
97 240 : cto(0, 0, i, i) -= poisson_effect;
98 240 : cto(1, 1, i, i) -= poisson_effect;
99 240 : if (q_trial > 0.0)
100 : {
101 144 : cto(0, 2, i, i) = _in_trial02 / q_trial * _dq_dpt * dpt_depii;
102 144 : cto(1, 2, i, i) = _in_trial12 / q_trial * _dq_dpt * dpt_depii;
103 : }
104 : }
105 :
106 : const Real poisson_effect =
107 80 : -Eijkl(2, 2, 0, 0) / Ezzzz *
108 80 : (_dgaE_dqt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dqt +
109 80 : gaE * smoothed_q.d2g[0][1] * _dq_dqt +
110 80 : gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dqt + dintnl0_dq * _dq_dqt) +
111 80 : gaE * smoothed_q.d2g_di[0][1] * (dintnl1_dqt + dintnl1_dp * _dp_dqt + dintnl1_dq * _dq_dqt));
112 :
113 80 : const Real dqt_dep02 = (q_trial == 0.0 ? 1.0 : _in_trial02 / q_trial) * Eijkl(0, 2, 0, 2);
114 80 : cto(2, 2, 0, 2) = _dp_dqt * dqt_dep02;
115 80 : cto(0, 0, 0, 2) = cto(1, 1, 0, 2) = poisson_effect * dqt_dep02;
116 80 : if (q_trial > 0.0)
117 : {
118 : // for q_trial=0, Jacobian_mult is just given by the elastic case
119 48 : cto(0, 2, 0, 2) = Eijkl(0, 2, 0, 2) * q / q_trial +
120 48 : _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep02;
121 48 : cto(1, 2, 0, 2) = _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep02;
122 : }
123 :
124 80 : const Real dqt_dep20 = (q_trial == 0.0 ? 1.0 : _in_trial02 / q_trial) * Eijkl(0, 2, 2, 0);
125 80 : cto(2, 2, 2, 0) = _dp_dqt * dqt_dep20;
126 80 : cto(0, 0, 2, 0) = cto(1, 1, 2, 0) = poisson_effect * dqt_dep20;
127 80 : if (q_trial > 0.0)
128 : {
129 : // for q_trial=0, Jacobian_mult is just given by the elastic case
130 48 : cto(0, 2, 2, 0) = Eijkl(0, 2, 2, 0) * q / q_trial +
131 48 : _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
132 48 : cto(1, 2, 2, 0) = _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
133 : }
134 :
135 80 : const Real dqt_dep12 = (q_trial == 0.0 ? 1.0 : _in_trial12 / q_trial) * Eijkl(1, 2, 1, 2);
136 80 : cto(2, 2, 1, 2) = _dp_dqt * dqt_dep12;
137 80 : cto(0, 0, 1, 2) = cto(1, 1, 1, 2) = poisson_effect * dqt_dep12;
138 80 : if (q_trial > 0.0)
139 : {
140 : // for q_trial=0, Jacobian_mult is just given by the elastic case
141 48 : cto(0, 2, 1, 2) = _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep12;
142 48 : cto(1, 2, 1, 2) = Eijkl(1, 2, 1, 2) * q / q_trial +
143 48 : _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep12;
144 : }
145 :
146 80 : const Real dqt_dep21 = (q_trial == 0.0 ? 1.0 : _in_trial12 / q_trial) * Eijkl(1, 2, 2, 1);
147 80 : cto(2, 2, 2, 1) = _dp_dqt * dqt_dep21;
148 80 : cto(0, 0, 2, 1) = cto(1, 1, 2, 1) = poisson_effect * dqt_dep21;
149 80 : if (q_trial > 0.0)
150 : {
151 : // for q_trial=0, Jacobian_mult is just given by the elastic case
152 48 : cto(0, 2, 2, 1) = _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
153 48 : cto(1, 2, 2, 1) = Eijkl(1, 2, 2, 1) * q / q_trial +
154 48 : _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
155 : }
156 : }
157 :
158 : RankTwoTensor
159 3248 : CappedWeakPlaneCosseratStressUpdate::dqdstress(const RankTwoTensor & stress) const
160 : {
161 3248 : RankTwoTensor deriv = RankTwoTensor();
162 3248 : const Real q = std::sqrt(Utility::pow<2>(stress(0, 2)) + Utility::pow<2>(stress(1, 2)));
163 3248 : if (q > 0.0)
164 : {
165 2576 : deriv(0, 2) = stress(0, 2) / q;
166 2576 : deriv(1, 2) = stress(1, 2) / q;
167 : }
168 : else
169 : {
170 : // derivative is not defined here. For now i'll set:
171 672 : deriv(0, 2) = 1.0;
172 672 : deriv(1, 2) = 1.0;
173 : }
174 3248 : return deriv;
175 : }
176 :
177 : RankFourTensor
178 0 : CappedWeakPlaneCosseratStressUpdate::d2qdstress2(const RankTwoTensor & stress) const
179 : {
180 0 : RankFourTensor d2 = RankFourTensor();
181 :
182 0 : const Real q = std::sqrt(Utility::pow<2>(stress(0, 2)) + Utility::pow<2>(stress(1, 2)));
183 0 : if (q == 0.0)
184 : return d2;
185 :
186 0 : const Real dq02 = stress(0, 2) / q;
187 0 : const Real dq12 = stress(1, 2) / q;
188 :
189 0 : d2(0, 2, 0, 2) = 1.0 / q - dq02 * dq02 / q;
190 0 : d2(0, 2, 1, 2) = -dq02 * dq12 / q;
191 0 : d2(1, 2, 0, 2) = -dq12 * dq02 / q;
192 0 : d2(1, 2, 1, 2) = 1.0 / q - dq12 * dq12 / q;
193 :
194 0 : return d2;
195 : }
|
