Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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("SolidMechanicsApp", CappedWeakPlaneCosseratStressUpdate);
15 :
16 : InputParameters
17 336 : CappedWeakPlaneCosseratStressUpdate::validParams()
18 : {
19 336 : InputParameters params = CappedWeakPlaneStressUpdate::validParams();
20 336 : params.addClassDescription("Capped weak-plane plasticity Cosserat stress calculator");
21 336 : return params;
22 0 : }
23 :
24 252 : CappedWeakPlaneCosseratStressUpdate::CappedWeakPlaneCosseratStressUpdate(
25 252 : const InputParameters & parameters)
26 252 : : CappedWeakPlaneStressUpdate(parameters)
27 : {
28 252 : }
29 :
30 : void
31 6368 : 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 6368 : stress = stress_trial;
41 6368 : 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 6368 : stress(0, 0) -= Eijkl(2, 2, 0, 0) * gaE / _Epp * smoothed_q.dg[0];
45 6368 : stress(1, 1) -= Eijkl(2, 2, 1, 1) * gaE / _Epp * smoothed_q.dg[0];
46 6368 : if (_in_q_trial == 0.0)
47 1280 : stress(0, 2) = stress(1, 2) = 0.0;
48 : else
49 : {
50 5088 : stress(0, 2) = _in_trial02 * q_ok / _in_q_trial;
51 5088 : stress(1, 2) = _in_trial12 * q_ok / _in_q_trial;
52 : }
53 6368 : }
54 :
55 : void
56 160 : 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 160 : cto = Eijkl;
70 160 : if (!compute_full_tangent_operator)
71 : return;
72 :
73 160 : const Real Ezzzz = Eijkl(2, 2, 2, 2);
74 160 : const Real Exzxz = Eijkl(0, 2, 0, 2);
75 160 : const Real tanpsi = _tan_psi.value(_intnl[_qp][0]);
76 160 : const Real dintnl0_dq = -1.0 / Exzxz;
77 : const Real dintnl0_dqt = 1.0 / Exzxz;
78 160 : const Real dintnl1_dp = -1.0 / Ezzzz;
79 : const Real dintnl1_dpt = 1.0 / Ezzzz;
80 : const Real dintnl1_dq =
81 160 : tanpsi / Exzxz - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dq / Exzxz;
82 : const Real dintnl1_dqt =
83 160 : -tanpsi / Exzxz - (q_trial - q) * _tan_psi.derivative(_intnl[_qp][0]) * dintnl0_dqt / Exzxz;
84 :
85 640 : for (unsigned i = 0; i < _tensor_dimensionality; ++i)
86 : {
87 480 : const Real dpt_depii = Eijkl(2, 2, i, i);
88 480 : cto(2, 2, i, i) = _dp_dpt * dpt_depii;
89 : const Real poisson_effect =
90 480 : Eijkl(2, 2, 0, 0) / Ezzzz *
91 480 : (_dgaE_dpt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dpt +
92 480 : gaE * smoothed_q.d2g[0][1] * _dq_dpt +
93 480 : gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dq * _dq_dpt) +
94 480 : gaE * smoothed_q.d2g_di[0][1] *
95 480 : (dintnl1_dpt + dintnl1_dp * _dp_dpt + dintnl1_dq * _dq_dpt)) *
96 480 : dpt_depii;
97 480 : cto(0, 0, i, i) -= poisson_effect;
98 480 : cto(1, 1, i, i) -= poisson_effect;
99 480 : if (q_trial > 0.0)
100 : {
101 288 : cto(0, 2, i, i) = _in_trial02 / q_trial * _dq_dpt * dpt_depii;
102 288 : cto(1, 2, i, i) = _in_trial12 / q_trial * _dq_dpt * dpt_depii;
103 : }
104 : }
105 :
106 : const Real poisson_effect =
107 160 : -Eijkl(2, 2, 0, 0) / Ezzzz *
108 160 : (_dgaE_dqt * smoothed_q.dg[0] + gaE * smoothed_q.d2g[0][0] * _dp_dqt +
109 160 : gaE * smoothed_q.d2g[0][1] * _dq_dqt +
110 160 : gaE * smoothed_q.d2g_di[0][0] * (dintnl0_dqt + dintnl0_dq * _dq_dqt) +
111 160 : gaE * smoothed_q.d2g_di[0][1] * (dintnl1_dqt + dintnl1_dp * _dp_dqt + dintnl1_dq * _dq_dqt));
112 :
113 160 : const Real dqt_dep02 = (q_trial == 0.0 ? 1.0 : _in_trial02 / q_trial) * Eijkl(0, 2, 0, 2);
114 160 : cto(2, 2, 0, 2) = _dp_dqt * dqt_dep02;
115 160 : cto(0, 0, 0, 2) = cto(1, 1, 0, 2) = poisson_effect * dqt_dep02;
116 160 : if (q_trial > 0.0)
117 : {
118 : // for q_trial=0, Jacobian_mult is just given by the elastic case
119 96 : cto(0, 2, 0, 2) = Eijkl(0, 2, 0, 2) * q / q_trial +
120 96 : _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep02;
121 96 : cto(1, 2, 0, 2) = _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep02;
122 : }
123 :
124 160 : const Real dqt_dep20 = (q_trial == 0.0 ? 1.0 : _in_trial02 / q_trial) * Eijkl(0, 2, 2, 0);
125 160 : cto(2, 2, 2, 0) = _dp_dqt * dqt_dep20;
126 160 : cto(0, 0, 2, 0) = cto(1, 1, 2, 0) = poisson_effect * dqt_dep20;
127 160 : if (q_trial > 0.0)
128 : {
129 : // for q_trial=0, Jacobian_mult is just given by the elastic case
130 96 : cto(0, 2, 2, 0) = Eijkl(0, 2, 2, 0) * q / q_trial +
131 96 : _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
132 96 : cto(1, 2, 2, 0) = _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep20;
133 : }
134 :
135 160 : const Real dqt_dep12 = (q_trial == 0.0 ? 1.0 : _in_trial12 / q_trial) * Eijkl(1, 2, 1, 2);
136 160 : cto(2, 2, 1, 2) = _dp_dqt * dqt_dep12;
137 160 : cto(0, 0, 1, 2) = cto(1, 1, 1, 2) = poisson_effect * dqt_dep12;
138 160 : if (q_trial > 0.0)
139 : {
140 : // for q_trial=0, Jacobian_mult is just given by the elastic case
141 96 : cto(0, 2, 1, 2) = _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep12;
142 96 : cto(1, 2, 1, 2) = Eijkl(1, 2, 1, 2) * q / q_trial +
143 96 : _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep12;
144 : }
145 :
146 160 : const Real dqt_dep21 = (q_trial == 0.0 ? 1.0 : _in_trial12 / q_trial) * Eijkl(1, 2, 2, 1);
147 160 : cto(2, 2, 2, 1) = _dp_dqt * dqt_dep21;
148 160 : cto(0, 0, 2, 1) = cto(1, 1, 2, 1) = poisson_effect * dqt_dep21;
149 160 : if (q_trial > 0.0)
150 : {
151 : // for q_trial=0, Jacobian_mult is just given by the elastic case
152 96 : cto(0, 2, 2, 1) = _in_trial02 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
153 96 : cto(1, 2, 2, 1) = Eijkl(1, 2, 2, 1) * q / q_trial +
154 96 : _in_trial12 * (_dq_dqt - q / q_trial) / q_trial * dqt_dep21;
155 : }
156 : }
157 :
158 : RankTwoTensor
159 6368 : CappedWeakPlaneCosseratStressUpdate::dqdstress(const RankTwoTensor & stress) const
160 : {
161 6368 : RankTwoTensor deriv = RankTwoTensor();
162 6368 : const Real q = std::sqrt(Utility::pow<2>(stress(0, 2)) + Utility::pow<2>(stress(1, 2)));
163 6368 : if (q > 0.0)
164 : {
165 5088 : deriv(0, 2) = stress(0, 2) / q;
166 5088 : deriv(1, 2) = stress(1, 2) / q;
167 : }
168 : else
169 : {
170 : // derivative is not defined here. For now i'll set:
171 1280 : deriv(0, 2) = 1.0;
172 1280 : deriv(1, 2) = 1.0;
173 : }
174 6368 : 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 : }
|
