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 "HyperElasticPhaseFieldIsoDamage.h"
11 : #include "libmesh/utility.h"
12 :
13 : registerMooseObject("SolidMechanicsApp", HyperElasticPhaseFieldIsoDamage);
14 :
15 : InputParameters
16 0 : HyperElasticPhaseFieldIsoDamage::validParams()
17 : {
18 0 : InputParameters params = FiniteStrainHyperElasticViscoPlastic::validParams();
19 0 : params.addParam<bool>("numerical_stiffness", false, "Flag for numerical stiffness");
20 0 : params.addParam<Real>("damage_stiffness", 1e-8, "Avoid zero after complete damage");
21 0 : params.addParam<Real>("zero_tol", 1e-12, "Tolerance for numerical zero");
22 0 : params.addParam<Real>(
23 0 : "zero_perturb", 1e-8, "Perturbation value when strain value less than numerical zero");
24 0 : params.addParam<Real>("perturbation_scale_factor", 1e-5, "Perturbation scale factor");
25 0 : params.addRequiredCoupledVar("c", "Damage variable");
26 0 : params.addParam<bool>(
27 0 : "use_current_history_variable", false, "Use the current value of the history variable.");
28 0 : params.addParam<MaterialPropertyName>(
29 : "F_name", "E_el", "Name of material property storing the elastic energy");
30 0 : params.addClassDescription(
31 : "Computes damaged stress and energy in the intermediate configuration assuming isotropy");
32 :
33 0 : return params;
34 0 : }
35 :
36 0 : HyperElasticPhaseFieldIsoDamage::HyperElasticPhaseFieldIsoDamage(const InputParameters & parameters)
37 : : FiniteStrainHyperElasticViscoPlastic(parameters),
38 0 : _num_stiffness(getParam<bool>("numerical_stiffness")),
39 0 : _kdamage(getParam<Real>("damage_stiffness")),
40 0 : _use_current_hist(getParam<bool>("use_current_history_variable")),
41 0 : _l(getMaterialProperty<Real>("l")),
42 0 : _gc(getMaterialProperty<Real>("gc_prop")),
43 0 : _zero_tol(getParam<Real>("zero_tol")),
44 0 : _zero_pert(getParam<Real>("zero_perturb")),
45 0 : _pert_val(getParam<Real>("perturbation_scale_factor")),
46 0 : _c(coupledValue("c")),
47 0 : _save_state(false),
48 0 : _dstress_dc(
49 0 : declarePropertyDerivative<RankTwoTensor>(_base_name + "stress", coupledName("c", 0))),
50 0 : _etens(LIBMESH_DIM),
51 0 : _F(declareProperty<Real>(getParam<MaterialPropertyName>("F_name"))),
52 0 : _dFdc(declarePropertyDerivative<Real>(getParam<MaterialPropertyName>("F_name"),
53 0 : coupledName("c", 0))),
54 0 : _d2Fdc2(declarePropertyDerivative<Real>(
55 0 : getParam<MaterialPropertyName>("F_name"), coupledName("c", 0), coupledName("c", 0))),
56 0 : _d2Fdcdstrain(declareProperty<RankTwoTensor>("d2Fdcdstrain")),
57 0 : _hist(declareProperty<Real>("hist")),
58 0 : _hist_old(getMaterialPropertyOld<Real>("hist"))
59 : {
60 0 : }
61 :
62 : void
63 0 : HyperElasticPhaseFieldIsoDamage::computePK2StressAndDerivative()
64 : {
65 0 : computeElasticStrain();
66 :
67 0 : _save_state = true;
68 0 : computeDamageStress();
69 0 : _pk2[_qp] = _pk2_tmp;
70 :
71 0 : _save_state = false;
72 0 : if (_num_stiffness)
73 0 : computeNumStiffness();
74 :
75 0 : if (_num_stiffness)
76 0 : _dpk2_dce = _dpk2_dee * _dee_dce;
77 :
78 0 : _dce_dfe.zero();
79 0 : for (const auto i : make_range(Moose::dim))
80 0 : for (const auto j : make_range(Moose::dim))
81 0 : for (const auto k : make_range(Moose::dim))
82 : {
83 0 : _dce_dfe(i, j, k, i) = _dce_dfe(i, j, k, i) + _fe(k, j);
84 0 : _dce_dfe(i, j, k, j) = _dce_dfe(i, j, k, j) + _fe(k, i);
85 : }
86 :
87 0 : _dpk2_dfe = _dpk2_dce * _dce_dfe;
88 0 : }
89 :
90 : void
91 0 : HyperElasticPhaseFieldIsoDamage::computeDamageStress()
92 : {
93 0 : Real lambda = _elasticity_tensor[_qp](0, 0, 1, 1);
94 0 : Real mu = _elasticity_tensor[_qp](0, 1, 0, 1);
95 :
96 0 : Real c = _c[_qp];
97 0 : Real xfac = Utility::pow<2>(1.0 - c) + _kdamage;
98 :
99 : std::vector<Real> eigval;
100 : RankTwoTensor evec;
101 0 : _ee.symmetricEigenvaluesEigenvectors(eigval, evec);
102 :
103 0 : for (const auto i : make_range(Moose::dim))
104 0 : _etens[i] = RankTwoTensor::selfOuterProduct(evec.column(i));
105 :
106 : Real etr = 0.0;
107 0 : for (const auto i : make_range(Moose::dim))
108 0 : etr += eigval[i];
109 :
110 0 : Real etrpos = (std::abs(etr) + etr) / 2.0;
111 0 : Real etrneg = (std::abs(etr) - etr) / 2.0;
112 :
113 0 : RankTwoTensor pk2pos, pk2neg;
114 :
115 0 : for (const auto i : make_range(Moose::dim))
116 : {
117 0 : pk2pos += _etens[i] * (lambda * etrpos + 2.0 * mu * (std::abs(eigval[i]) + eigval[i]) / 2.0);
118 0 : pk2neg += _etens[i] * (lambda * etrneg + 2.0 * mu * (std::abs(eigval[i]) - eigval[i]) / 2.0);
119 : }
120 :
121 0 : _pk2_tmp = pk2pos * xfac - pk2neg;
122 :
123 0 : if (_save_state)
124 : {
125 0 : std::vector<Real> epos(LIBMESH_DIM), eneg(LIBMESH_DIM);
126 0 : for (const auto i : make_range(Moose::dim))
127 : {
128 0 : epos[i] = (std::abs(eigval[i]) + eigval[i]) / 2.0;
129 0 : eneg[i] = (std::abs(eigval[i]) - eigval[i]) / 2.0;
130 : }
131 :
132 : // sum squares of epos and eneg
133 : Real pval(0.0), nval(0.0);
134 0 : for (const auto i : make_range(Moose::dim))
135 : {
136 0 : pval += epos[i] * epos[i];
137 0 : nval += eneg[i] * eneg[i];
138 : }
139 :
140 : // Energy with positive principal strains
141 0 : const Real G0_pos = lambda * etrpos * etrpos / 2.0 + mu * pval;
142 0 : const Real G0_neg = lambda * etrneg * etrneg / 2.0 + mu * nval;
143 :
144 : // Assign history variable and derivative
145 0 : if (G0_pos > _hist_old[_qp])
146 0 : _hist[_qp] = G0_pos;
147 : else
148 0 : _hist[_qp] = _hist_old[_qp];
149 :
150 0 : Real hist_variable = _hist_old[_qp];
151 0 : if (_use_current_hist)
152 0 : hist_variable = _hist[_qp];
153 :
154 : // Elastic free energy density
155 0 : _F[_qp] = hist_variable * xfac - G0_neg + _gc[_qp] / (2 * _l[_qp]) * c * c;
156 :
157 : // derivative of elastic free energy density wrt c
158 0 : _dFdc[_qp] = -hist_variable * 2.0 * (1.0 - c) * (1 - _kdamage) + _gc[_qp] / _l[_qp] * c;
159 :
160 : // 2nd derivative of elastic free energy density wrt c
161 0 : _d2Fdc2[_qp] = hist_variable * 2.0 * (1 - _kdamage) + _gc[_qp] / _l[_qp];
162 :
163 0 : _dG0_dee = pk2pos;
164 :
165 0 : _dpk2_dc = -pk2pos * 2.0 * (1.0 - c);
166 : }
167 0 : }
168 :
169 : void
170 0 : HyperElasticPhaseFieldIsoDamage::computeNumStiffness()
171 : {
172 0 : RankTwoTensor ee_tmp;
173 :
174 0 : for (const auto i : make_range(Moose::dim))
175 0 : for (unsigned int j = i; j < LIBMESH_DIM; ++j)
176 : {
177 0 : Real ee_pert = _zero_pert;
178 0 : if (std::abs(_ee(i, j)) > _zero_tol)
179 0 : ee_pert = _pert_val * std::abs(_ee(i, j));
180 :
181 0 : ee_tmp = _ee;
182 0 : _ee(i, j) += ee_pert;
183 :
184 0 : computeDamageStress();
185 :
186 0 : for (const auto k : make_range(Moose::dim))
187 0 : for (const auto l : make_range(Moose::dim))
188 : {
189 0 : _dpk2_dee(k, l, i, j) = (_pk2_tmp(k, l) - _pk2[_qp](k, l)) / ee_pert;
190 0 : _dpk2_dee(k, l, j, i) = (_pk2_tmp(k, l) - _pk2[_qp](k, l)) / ee_pert;
191 : }
192 0 : _ee = ee_tmp;
193 : }
194 0 : }
195 :
196 : void
197 0 : HyperElasticPhaseFieldIsoDamage::computeQpJacobian()
198 : {
199 0 : FiniteStrainHyperElasticViscoPlastic::computeQpJacobian();
200 :
201 0 : RankTwoTensor dG0_dce = _dee_dce.innerProductTranspose(_dG0_dee);
202 0 : RankTwoTensor dG0_dfe = _dce_dfe.innerProductTranspose(dG0_dce);
203 0 : RankTwoTensor dG0_df = _dfe_df.innerProductTranspose(dG0_dfe);
204 :
205 : // 2nd derivative wrt c and strain = 0.0 if we used the previous step's history varible
206 0 : if (_use_current_hist)
207 0 : _d2Fdcdstrain[_qp] =
208 0 : -_df_dstretch_inc.innerProductTranspose(dG0_df) * 2.0 * (1.0 - _c[_qp]) * (1 - _kdamage);
209 :
210 : usingTensorIndices(i_, j_, k_, l_);
211 0 : _dstress_dc[_qp] = _fe.times<i_, k_, j_, l_>(_fe) * _dpk2_dc;
212 0 : }
|
