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