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 "BiLinearMixedModeTraction.h"
11 : #include "libmesh/utility.h"
12 :
13 : registerMooseObject("SolidMechanicsApp", BiLinearMixedModeTraction);
14 :
15 : InputParameters
16 70 : BiLinearMixedModeTraction::validParams()
17 : {
18 70 : InputParameters params = CZMComputeLocalTractionTotalBase::validParams();
19 70 : params.addClassDescription("Mixed mode bilinear traction separation law.");
20 140 : params.addRequiredParam<Real>("penalty_stiffness", "Penalty stiffness.");
21 140 : params.addRequiredParam<MaterialPropertyName>(
22 : "GI_c", "Critical energy release rate in normal direction.");
23 140 : params.addRequiredParam<MaterialPropertyName>("GII_c",
24 : "Critical energy release rate in shear direction.");
25 140 : params.addRequiredParam<MaterialPropertyName>("normal_strength",
26 : "Tensile strength in normal direction.");
27 140 : params.addRequiredParam<MaterialPropertyName>("shear_strength",
28 : "Tensile strength in shear direction.");
29 140 : params.addRequiredParam<Real>("eta", "The power law parameter.");
30 140 : MooseEnum criterion("POWER_LAW BK", "BK");
31 140 : params.addParam<Real>("viscosity", 0.0, "Viscosity.");
32 140 : params.addParam<MooseEnum>(
33 : "mixed_mode_criterion", criterion, "Option for mixed mode propagation criterion.");
34 :
35 : // Advanced parameters to improve numerical convergence
36 140 : params.addParam<bool>(
37 140 : "lag_mode_mixity", true, "Whether to use old displacement jumps to compute the mode mixity.");
38 140 : params.addParam<bool>(
39 : "lag_displacement_jump",
40 140 : false,
41 : "Whether to use old displacement jumps to compute the effective displacement jump.");
42 140 : params.addParam<Real>("alpha", 1e-10, "Regularization parameter for the Macaulay bracket.");
43 140 : params.addParamNamesToGroup("lag_mode_mixity lag_displacement_jump alpha", "Advanced");
44 :
45 70 : return params;
46 70 : }
47 :
48 35 : BiLinearMixedModeTraction::BiLinearMixedModeTraction(const InputParameters & parameters)
49 : : CZMComputeLocalTractionTotalBase(parameters),
50 35 : _K(getParam<Real>("penalty_stiffness")),
51 35 : _d(declareProperty<Real>(_base_name + "damage")),
52 70 : _d_old(getMaterialPropertyOld<Real>(_base_name + "damage")),
53 35 : _interface_displacement_jump_old(
54 35 : getMaterialPropertyOld<RealVectorValue>(_base_name + "interface_displacement_jump")),
55 35 : _delta_init(
56 35 : declareProperty<Real>(_base_name + "effective_displacement_jump_at_damage_initiation")),
57 35 : _delta_final(
58 35 : declareProperty<Real>(_base_name + "effective_displacement_jump_at_full_degradation")),
59 35 : _delta_m(declareProperty<Real>(_base_name + "effective_displacement_jump")),
60 105 : _GI_c(getMaterialProperty<Real>(_base_name + getParam<MaterialPropertyName>("GI_c"))),
61 105 : _GII_c(getMaterialProperty<Real>(_base_name + getParam<MaterialPropertyName>("GII_c"))),
62 105 : _N(getMaterialProperty<Real>(_base_name + getParam<MaterialPropertyName>("normal_strength"))),
63 105 : _S(getMaterialProperty<Real>(_base_name + getParam<MaterialPropertyName>("shear_strength"))),
64 70 : _eta(getParam<Real>("eta")),
65 35 : _beta(declareProperty<Real>(_base_name + "mode_mixity_ratio")),
66 70 : _viscosity(getParam<Real>("viscosity")),
67 70 : _criterion(getParam<MooseEnum>("mixed_mode_criterion").getEnum<MixedModeCriterion>()),
68 70 : _lag_mode_mixity(getParam<bool>("lag_mode_mixity")),
69 70 : _lag_disp_jump(getParam<bool>("lag_displacement_jump")),
70 105 : _alpha(getParam<Real>("alpha"))
71 : {
72 35 : }
73 :
74 : void
75 180 : BiLinearMixedModeTraction::initQpStatefulProperties()
76 : {
77 180 : CZMComputeLocalTractionTotalBase::initQpStatefulProperties();
78 :
79 180 : _d[_qp] = 0.0;
80 180 : }
81 :
82 : void
83 6312 : BiLinearMixedModeTraction::computeInterfaceTractionAndDerivatives()
84 : {
85 6312 : _interface_traction[_qp] = computeTraction();
86 6312 : _dinterface_traction_djump[_qp] = computeTractionDerivatives();
87 : }
88 :
89 : RealVectorValue
90 6312 : BiLinearMixedModeTraction::computeTraction()
91 : {
92 6312 : computeModeMixity();
93 6312 : computeCriticalDisplacementJump();
94 6312 : computeFinalDisplacementJump();
95 6312 : computeEffectiveDisplacementJump();
96 6312 : computeDamage();
97 :
98 : // Split displacement jump into active and inactive parts
99 6312 : const RealVectorValue delta = _interface_displacement_jump[_qp];
100 6312 : const RealVectorValue delta_active(std::max(delta(0), 0.), delta(1), delta(2));
101 6312 : const RealVectorValue delta_inactive(std::min(delta(0), 0.), 0, 0);
102 :
103 6312 : return (1 - _d[_qp]) * _K * delta_active + _K * delta_inactive;
104 : }
105 :
106 : RankTwoTensor
107 6312 : BiLinearMixedModeTraction::computeTractionDerivatives()
108 : {
109 : // The displacement jump split depends on the displacement jump, obviously
110 6312 : const RealVectorValue delta = _interface_displacement_jump[_qp];
111 6312 : RankTwoTensor ddelta_active_ddelta, ddelta_inactive_ddelta;
112 7914 : ddelta_active_ddelta.fillFromInputVector({delta(0) > 0 ? 1. : 0., 1., 1.});
113 11202 : ddelta_inactive_ddelta.fillFromInputVector({delta(0) < 0 ? 1. : 0., 0., 0.});
114 :
115 : RankTwoTensor dtraction_ddelta =
116 6312 : (1 - _d[_qp]) * _K * ddelta_active_ddelta + _K * ddelta_inactive_ddelta;
117 :
118 : // The damage may also depend on the displacement jump
119 6312 : const RealVectorValue delta_active(std::max(delta(0), 0.), delta(1), delta(2));
120 6312 : RankTwoTensor A;
121 6312 : A.vectorOuterProduct(delta_active, _dd_ddelta);
122 6312 : dtraction_ddelta -= _K * A;
123 :
124 6312 : return dtraction_ddelta;
125 : }
126 :
127 : void
128 6312 : BiLinearMixedModeTraction::computeModeMixity()
129 : {
130 : const RealVectorValue delta =
131 6312 : _lag_mode_mixity ? _interface_displacement_jump_old[_qp] : _interface_displacement_jump[_qp];
132 :
133 6312 : if (delta(0) > 0)
134 : {
135 4422 : const Real delta_s = std::sqrt(delta(1) * delta(1) + delta(2) * delta(2));
136 4422 : _beta[_qp] = delta_s / delta(0);
137 :
138 4422 : if (!_lag_mode_mixity)
139 0 : _dbeta_ddelta = RealVectorValue(-delta_s / delta(0) / delta(0),
140 0 : delta(1) / delta_s / delta(0),
141 0 : delta(2) / delta_s / delta(0));
142 : }
143 : else
144 : {
145 1890 : _beta[_qp] = 0;
146 1890 : _dbeta_ddelta = RealVectorValue(0, 0, 0);
147 : }
148 6312 : }
149 :
150 : void
151 6312 : BiLinearMixedModeTraction::computeCriticalDisplacementJump()
152 : {
153 : const RealVectorValue delta =
154 6312 : _lag_mode_mixity ? _interface_displacement_jump_old[_qp] : _interface_displacement_jump[_qp];
155 :
156 6312 : const Real delta_normal0 = _N[_qp] / _K;
157 6312 : const Real delta_shear0 = _S[_qp] / _K;
158 :
159 6312 : _delta_init[_qp] = delta_shear0;
160 6312 : _ddelta_init_ddelta = RealVectorValue(0, 0, 0);
161 6312 : if (delta(0) > 0)
162 : {
163 : const Real delta_mixed =
164 4422 : std::sqrt(delta_shear0 * delta_shear0 + Utility::pow<2>(_beta[_qp] * delta_normal0));
165 4422 : _delta_init[_qp] =
166 4422 : delta_normal0 * delta_shear0 * std::sqrt(1 + _beta[_qp] * _beta[_qp]) / delta_mixed;
167 4422 : if (!_lag_mode_mixity)
168 : {
169 : const Real ddelta_init_dbeta =
170 0 : _delta_init[_qp] * _beta[_qp] *
171 0 : (1 / (1 + _beta[_qp] * _beta[_qp]) - Utility::pow<2>(_delta_init[_qp] / delta_mixed));
172 0 : _ddelta_init_ddelta = ddelta_init_dbeta * _dbeta_ddelta;
173 : }
174 : }
175 6312 : }
176 :
177 : void
178 6312 : BiLinearMixedModeTraction::computeFinalDisplacementJump()
179 : {
180 : const RealVectorValue delta =
181 6312 : _lag_mode_mixity ? _interface_displacement_jump_old[_qp] : _interface_displacement_jump[_qp];
182 :
183 6312 : _delta_final[_qp] = std::sqrt(2) * 2 * _GII_c[_qp] / _S[_qp];
184 6312 : _ddelta_final_ddelta = RealVectorValue(0, 0, 0);
185 6312 : if (delta(0) > 0)
186 : {
187 4422 : if (_criterion == MixedModeCriterion::BK)
188 : {
189 3826 : _delta_final[_qp] =
190 3826 : 2 / _K / _delta_init[_qp] *
191 3826 : (_GI_c[_qp] +
192 3826 : (_GII_c[_qp] - _GI_c[_qp]) *
193 3826 : std::pow(_beta[_qp] * _beta[_qp] / (1 + _beta[_qp] * _beta[_qp]), _eta));
194 3826 : if (!_lag_mode_mixity)
195 : {
196 0 : const Real ddelta_final_ddelta_init = -_delta_final[_qp] / _delta_init[_qp];
197 : const Real ddelta_final_dbeta =
198 0 : 2 / _K / _delta_init[_qp] * (_GII_c[_qp] - _GI_c[_qp]) * _eta *
199 0 : std::pow(_beta[_qp] * _beta[_qp] / (1 + _beta[_qp] * _beta[_qp]), _eta - 1) * 2 *
200 0 : _beta[_qp] * (1 - Utility::pow<2>(_beta[_qp] / (1 + _beta[_qp] * _beta[_qp])));
201 0 : _ddelta_final_ddelta =
202 : ddelta_final_ddelta_init * _ddelta_init_ddelta + ddelta_final_dbeta * _dbeta_ddelta;
203 : }
204 : }
205 596 : else if (_criterion == MixedModeCriterion::POWER_LAW)
206 : {
207 : const Real Gc_mixed =
208 596 : std::pow(1 / _GI_c[_qp], _eta) + std::pow(_beta[_qp] * _beta[_qp] / _GII_c[_qp], _eta);
209 596 : _delta_final[_qp] =
210 596 : (2 + 2 * _beta[_qp] * _beta[_qp]) / _K / _delta_init[_qp] * std::pow(Gc_mixed, -1 / _eta);
211 596 : if (!_lag_mode_mixity)
212 : {
213 0 : const Real ddelta_final_ddelta_init = -_delta_final[_qp] / _delta_init[_qp];
214 : const Real ddelta_final_dbeta =
215 0 : _delta_final[_qp] * 2 * _beta[_qp] / (1 + _beta[_qp] * _beta[_qp]) -
216 0 : (2 + 2 * _beta[_qp] * _beta[_qp]) / _K / _delta_init[_qp] *
217 0 : std::pow(Gc_mixed, -1 / _eta - 1) *
218 0 : (std::pow(1 / _GI_c[_qp], _eta - 1) +
219 0 : std::pow(_beta[_qp] * _beta[_qp] / _GII_c[_qp], _eta - 1) * 2 * _beta[_qp] /
220 : _GII_c[_qp]);
221 0 : _ddelta_final_ddelta =
222 : ddelta_final_ddelta_init * _ddelta_init_ddelta + ddelta_final_dbeta * _dbeta_ddelta;
223 : }
224 : }
225 : }
226 6312 : }
227 :
228 : void
229 6312 : BiLinearMixedModeTraction::computeEffectiveDisplacementJump()
230 : {
231 : const RealVectorValue delta =
232 6312 : _lag_disp_jump ? _interface_displacement_jump_old[_qp] : _interface_displacement_jump[_qp];
233 :
234 6312 : const Real delta_normal_pos = MathUtils::regularizedHeavyside(delta(0), _alpha) * delta(0);
235 6312 : _delta_m[_qp] = std::sqrt(Utility::pow<2>(delta(1)) + Utility::pow<2>(delta(2)) +
236 : Utility::pow<2>(delta_normal_pos));
237 6312 : _ddelta_m_ddelta = RealVectorValue(0, 0, 0);
238 6312 : if (!_lag_disp_jump && !MooseUtils::absoluteFuzzyEqual(_delta_m[_qp], 0))
239 : {
240 : const Real ddelta_normal_pos_ddelta_normal =
241 5234 : MathUtils::regularizedHeavysideDerivative(delta(0), 1e-6) * delta(0) +
242 5234 : MathUtils::regularizedHeavyside(delta(0), _alpha);
243 5234 : _ddelta_m_ddelta =
244 5234 : RealVectorValue(delta_normal_pos * ddelta_normal_pos_ddelta_normal, delta(1), delta(2));
245 : _ddelta_m_ddelta /= _delta_m[_qp];
246 : }
247 6312 : }
248 :
249 : void
250 6312 : BiLinearMixedModeTraction::computeDamage()
251 : {
252 6312 : if (_delta_m[_qp] < _delta_init[_qp])
253 3041 : _d[_qp] = 0;
254 3271 : else if (_delta_m[_qp] > _delta_final[_qp])
255 700 : _d[_qp] = 1;
256 : else
257 2571 : _d[_qp] = _delta_final[_qp] * (_delta_m[_qp] - _delta_init[_qp]) / _delta_m[_qp] /
258 2571 : (_delta_final[_qp] - _delta_init[_qp]);
259 :
260 6312 : _dd_ddelta = RealVectorValue(0, 0, 0);
261 6312 : if (_d[_qp] < _d_old[_qp])
262 : // Irreversibility
263 2062 : _d[_qp] = _d_old[_qp];
264 4250 : else if (_delta_m[_qp] >= _delta_init[_qp] && _delta_m[_qp] <= _delta_final[_qp])
265 1230 : _dd_ddelta = (_ddelta_final_ddelta * (_delta_m[_qp] - _delta_init[_qp]) +
266 : _delta_final[_qp] * (_ddelta_m_ddelta - _ddelta_init_ddelta)) /
267 : _delta_m[_qp] / (_delta_final[_qp] - _delta_init[_qp]) -
268 1230 : _delta_final[_qp] * (_delta_m[_qp] - _delta_init[_qp]) *
269 : (_ddelta_m_ddelta * (_delta_final[_qp] - _delta_init[_qp]) +
270 : _delta_m[_qp] * (_ddelta_final_ddelta - _ddelta_init_ddelta)) /
271 1230 : Utility::pow<2>(_delta_m[_qp] * (_delta_final[_qp] - _delta_init[_qp]));
272 :
273 : // Viscous regularization
274 6312 : _d[_qp] = (_d[_qp] + _viscosity * _d_old[_qp] / _dt) / (_viscosity / _dt + 1);
275 6312 : _dd_ddelta /= (_viscosity / _dt + 1);
276 6312 : }
|
