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