www.mooseframework.org
InclusionProperties.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "InclusionProperties.h"
11 #include "libmesh/utility.h"
12 
13 registerMooseObject("TensorMechanicsApp", InclusionProperties);
14 
15 defineLegacyParams(InclusionProperties);
16 
17 InputParameters
18 InclusionProperties::validParams()
19 {
20  InputParameters params = Material::validParams();
21  params.addRequiredParam<Real>("a", "Ellipse semiaxis");
22  params.addRequiredParam<Real>("b", "Ellipse semiaxis");
23  params.addRequiredParam<Real>("lambda", "Lame's first parameter");
24  params.addRequiredParam<Real>("mu", "Shear modulus (aka Lame's second parameter)");
25  params.addRequiredParam<std::vector<Real>>("misfit_strains",
26  "Vector of misfit strains in order eps_11, eps_22");
27  params.addRequiredParam<MaterialPropertyName>(
28  "stress_name", "Name of the material property where analytical stresses will be stored");
29  params.addRequiredParam<MaterialPropertyName>(
30  "strain_name", "Name of the material property where analytical total strains will be stored");
31  params.addRequiredParam<MaterialPropertyName>(
32  "energy_name",
33  "Name of the material property where analytical elastic energies will be stored");
34 
35  return params;
36 }
37 
38 InclusionProperties::InclusionProperties(const InputParameters & parameters)
39  : Material(parameters),
40  _a(getParam<Real>("a")),
41  _b(getParam<Real>("b")),
42  _lambda(getParam<Real>("lambda")),
43  _mu(getParam<Real>("mu")),
44  _misfit(getParam<std::vector<Real>>("misfit_strains")),
45  _stress(declareProperty<RankTwoTensor>(getParam<MaterialPropertyName>("stress_name"))),
46  _strain(declareProperty<RankTwoTensor>(getParam<MaterialPropertyName>("strain_name"))),
47  _elastic_energy(declareProperty<Real>(getParam<MaterialPropertyName>("energy_name")))
48 {
49  if (_misfit.size() != 2)
50  mooseError("Supply 2 misfit_strains in order eps_11, eps_22 in InclusionProperties.");
51 
52  _nu = _lambda / 2.0 / (_lambda + _mu);
53  _kappa = 3 - 4 * _nu;
54  precomputeInteriorProperties();
55 }
56 
57 void
58 InclusionProperties::computeQpProperties()
59 {
60  Real x = _q_point[_qp](0);
61  Real y = _q_point[_qp](1);
62  if (x * x / _a / _a + y * y / _b / _b < 1)
63  {
64  // Inside the inclusion
65  _stress[_qp] = _stress_int;
66  _strain[_qp] = _total_strain_int;
67  _elastic_energy[_qp] = _elastic_energy_int;
68  }
69  else
70  {
71  // Outside the inclusion
72  Real l = 0.5 * (x * x + y * y - _a * _a - _b * _b // Parameter l called lambda in the paper
73  +
74  std::sqrt(Utility::pow<2>((x * x + y * y - _a * _a + _b * _b)) +
75  4 * (_a * _a - _b * _b) * y * y));
76  Real rho_a = _a / sqrt(_a * _a + l);
77  Real rho_b = _b / sqrt(_b * _b + l);
78  Real m_x = x / (_a * _a + l);
79  Real m_y = y / (_b * _b + l);
80  Real n_x = m_x / std::sqrt(m_x * m_x + m_y * m_y);
81  Real n_y = m_y / std::sqrt(m_x * m_x + m_y * m_y);
82  Real T_6 = rho_a * rho_a + rho_b * rho_b - 4 * rho_a * rho_a * n_x * n_x -
83  4 * rho_b * rho_b * n_y * n_y - 4;
84 
85  Real H11 = rho_a * _b * (_a * rho_b + _b * rho_a + 2 * _a * rho_a * rho_a * rho_b +
86  _b * Utility::pow<3>(rho_a)) /
87  Utility::pow<2>((_a * rho_b + _b * rho_a)) +
88  n_x * n_x * (2 - 6 * rho_a * rho_a + (8 * rho_a * rho_a + T_6) * n_x * n_x);
89 
90  Real H22 = rho_b * _a * (_a * rho_b + _b * rho_a + 2 * _b * rho_a * rho_b * rho_b +
91  _a * Utility::pow<3>(rho_b)) /
92  Utility::pow<2>((_a * rho_b + _b * rho_a)) +
93  n_y * n_y * (2 - 6 * rho_b * rho_b + (8 * rho_b * rho_b + T_6) * n_y * n_y);
94 
95  Real H12 = (_a * _a * rho_a * rho_a * rho_b * rho_b + _b * _b * rho_a * rho_a +
96  _a * _b * rho_a * rho_b) /
97  Utility::pow<2>((_a * rho_b + _b * rho_a)) -
98  rho_b * rho_b * n_x * n_x - rho_a * rho_a * n_y * n_y +
99  (4 * rho_a * rho_a + 4 * rho_b * rho_b + T_6) * n_x * n_x * n_y * n_y;
100 
101  Real H31 = 2 * (_b * rho_a / (_a * rho_b + _b * rho_a) - n_x * n_x);
102  Real H32 = 2 * (_a * rho_b / (_a * rho_b + _b * rho_a) - n_y * n_y);
103 
104  Real H41 = n_x * n_y *
105  (1 - 3 * rho_a * rho_a + (6 * rho_a * rho_a + 2 * rho_b * rho_b + T_6) * n_x * n_x);
106  Real H42 = n_x * n_y *
107  (1 - 3 * rho_b * rho_b + (2 * rho_a * rho_a + 6 * rho_b * rho_b + T_6) * n_y * n_y);
108 
109  _stress[_qp](0, 0) =
110  4 * _mu * rho_a * rho_b / (_kappa + 1) * (H11 * _misfit[0] + H12 * _misfit[1]);
111  _stress[_qp](1, 1) =
112  4 * _mu * rho_a * rho_b / (_kappa + 1) * (H12 * _misfit[0] + H22 * _misfit[1]);
113  _stress[_qp](2, 2) =
114  4 * _mu * rho_a * rho_b / (_kappa + 1) * _nu * (H31 * _misfit[0] + H32 * _misfit[1]);
115  _stress[_qp](0, 1) = _stress[_qp](1, 0) =
116  4 * _mu * rho_a * rho_b / (_kappa + 1) * (H41 * _misfit[0] + H42 * _misfit[1]);
117 
118  Real J1 = rho_a * rho_a * rho_b * _b / (_a * rho_b + _b * rho_a);
119  Real J11 = Utility::pow<4>(rho_a) * rho_b * _b / (3 * _a * _a) * (2 * _a * rho_b + _b * rho_a) /
120  Utility::pow<2>((_a * rho_b + _b * rho_a));
121  Real J12 = Utility::pow<3>(rho_a) * Utility::pow<3>(rho_b) /
122  Utility::pow<2>((_a * rho_b + _b * rho_a));
123  Real J2 = rho_b * rho_b * rho_a * _a / (_a * rho_b + _b * rho_a);
124  Real J22 = Utility::pow<4>(rho_b) * rho_a * _a / (3 * _b * _b) * (2 * _b * rho_a + _a * rho_b) /
125  Utility::pow<2>((_a * rho_b + _b * rho_a));
126 
127  Real G1111 = ((1 - 2 * _nu) * J1 + 3 * _a * _a * J11) / (2 * (1 - _nu)) +
128  rho_a * rho_b * n_x * n_x / (2 * (1 - _nu)) *
129  (2 + 2 * _nu - 6 * rho_a * rho_a + (8 * rho_a * rho_a + T_6) * n_x * n_x);
130  Real G1122 = ((2 * _nu - 1) * J1 + _b * _b * J12) / (2 * (1 - _nu)) +
131  rho_a * rho_b / (2 * (1 - _nu)) *
132  ((1 - rho_a * rho_a) * n_y * n_y + (1 - 2 * _nu - rho_b * rho_b) * n_x * n_x +
133  (4 * rho_a * rho_a + 4 * rho_b * rho_b + T_6) * n_x * n_x * n_y * n_y);
134  Real G2211 = ((2 * _nu - 1) * J2 + _a * _a * J12) / (2 * (1 - _nu)) +
135  rho_a * rho_b / (2 * (1 - _nu)) *
136  ((1 - rho_b * rho_b) * n_x * n_x + (1 - 2 * _nu - rho_a * rho_a) * n_y * n_y +
137  (4 * rho_a * rho_a + 4 * rho_b * rho_b + T_6) * n_x * n_x * n_y * n_y);
138  Real G2222 = ((1 - 2 * _nu) * J2 + 3 * _b * _b * J22) / (2 * (1 - _nu)) +
139  rho_a * rho_b * n_y * n_y / (2 * (1 - _nu)) *
140  (2 + 2 * _nu - 6 * rho_b * rho_b + (8 * rho_a * rho_a + T_6) * n_y * n_y);
141  Real G1211 =
142  rho_a * rho_b * n_x * n_y / (2 * (1 - _nu)) *
143  (1 - 3 * rho_a * rho_a + (6 * rho_a * rho_a + 2 * rho_b * rho_b + T_6) * n_x * n_x);
144  Real G1222 =
145  rho_a * rho_b * n_x * n_y / (2 * (1 - _nu)) *
146  (1 - 3 * rho_b * rho_b + (2 * rho_a * rho_a + 6 * rho_b * rho_b + T_6) * n_y * n_y);
147 
148  // Outside the inclusion, total strain = elastic strain so we only need to
149  // calculate one strain tensor
150  _strain[_qp](0, 0) = G1111 * _misfit[0] + G1122 * _misfit[1];
151  _strain[_qp](1, 1) = G2211 * _misfit[0] + G2222 * _misfit[1];
152  _strain[_qp](0, 1) = _strain[_qp](1, 0) = G1211 * _misfit[0] + G1222 * _misfit[1];
153 
154  _elastic_energy[_qp] = 0.5 * _stress[_qp].doubleContraction(_strain[_qp]);
155  }
156 }
157 
158 void
159 InclusionProperties::precomputeInteriorProperties()
160 {
161  Real t = _b / _a;
162  Real T11 = 1 + 2 / t;
163  Real T12 = 1;
164  Real T21 = 1;
165  Real T22 = 1 + 2 * t;
166  Real T31 = (3 - _kappa) / 2 * (1 + 1 / t);
167  Real T32 = (3 - _kappa) / 2 * (1 + t);
168 
169  _stress_int(0, 0) =
170  -4 * _mu * t / (_kappa + 1) / (1 + t) / (1 + t) * (T11 * _misfit[0] + T12 * _misfit[1]);
171  _stress_int(1, 1) =
172  -4 * _mu * t / (_kappa + 1) / (1 + t) / (1 + t) * (T21 * _misfit[0] + T22 * _misfit[1]);
173  _stress_int(2, 2) =
174  -4 * _mu * t / (_kappa + 1) / (1 + t) / (1 + t) * (T31 * _misfit[0] + T32 * _misfit[1]);
175 
176  Real S11 = t * (t + 3 + _kappa * (1 + t));
177  Real S12 = t * (1 + 3 * t - _kappa * (1 + t));
178  Real S21 = t + 3 - _kappa * (1 + t);
179  Real S22 = 1 + 3 * t + _kappa * (1 + t);
180 
181  _total_strain_int(0, 0) =
182  1 / (_kappa + 1) / (1 + t) / (1 + t) * (S11 * _misfit[0] + S12 * _misfit[1]);
183  _total_strain_int(1, 1) =
184  1 / (_kappa + 1) / (1 + t) / (1 + t) * (S21 * _misfit[0] + S22 * _misfit[1]);
185  // Inside the inclusion, elastic strain = total strain - eigenstrain
186  _elastic_strain_int(0, 0) = _total_strain_int(0, 0) - _misfit[0];
187  _elastic_strain_int(1, 1) = _total_strain_int(1, 1) - _misfit[1];
188 
189  _elastic_energy_int = 0.5 * _stress_int.doubleContraction(_elastic_strain_int);
190 }
InclusionProperties::InclusionProperties
InclusionProperties(const InputParameters &parameters)
Definition: InclusionProperties.C:38
InclusionProperties::_b
const Real _b
Definition: InclusionProperties.h:42
defineLegacyParams
defineLegacyParams(InclusionProperties)
InclusionProperties::_strain
MaterialProperty< RankTwoTensor > & _strain
Definition: InclusionProperties.h:66
InclusionProperties::_nu
Real _nu
Poisson's ratio.
Definition: InclusionProperties.h:52
InclusionProperties::validParams
static InputParameters validParams()
Definition: InclusionProperties.C:18
InclusionProperties::_total_strain_int
RankTwoTensor _total_strain_int
Definition: InclusionProperties.h:61
InclusionProperties::_lambda
const Real _lambda
Elastic constants (isotropic)
Definition: InclusionProperties.h:45
InclusionProperties::_elastic_strain_int
RankTwoTensor _elastic_strain_int
Definition: InclusionProperties.h:62
InclusionProperties::computeQpProperties
virtual void computeQpProperties()
Definition: InclusionProperties.C:58
validParams
InputParameters validParams()
registerMooseObject
registerMooseObject("TensorMechanicsApp", InclusionProperties)
InclusionProperties::_stress
MaterialProperty< RankTwoTensor > & _stress
Definition: InclusionProperties.h:65
InclusionProperties::_a
const Real _a
Semimajor axes of the ellipsoidal inclusion.
Definition: InclusionProperties.h:41
InclusionProperties
This material calculates the stresses, strains, and elastic energies for an ellipsoidal inclusion in ...
Definition: InclusionProperties.h:28
InclusionProperties.h
InclusionProperties::_stress_int
RankTwoTensor _stress_int
Interior stress and strain values are constant so they only need to be calculated once.
Definition: InclusionProperties.h:60
InclusionProperties::_kappa
Real _kappa
Kolosov's first constant.
Definition: InclusionProperties.h:54
RankTwoTensorTempl< Real >
InclusionProperties::precomputeInteriorProperties
virtual void precomputeInteriorProperties()
Definition: InclusionProperties.C:159
InclusionProperties::_elastic_energy
MaterialProperty< Real > & _elastic_energy
Definition: InclusionProperties.h:67
InclusionProperties::_elastic_energy_int
Real _elastic_energy_int
Definition: InclusionProperties.h:63
InclusionProperties::_mu
const Real _mu
Definition: InclusionProperties.h:46
InclusionProperties::_misfit
std::vector< Real > _misfit
Misfit strains.
Definition: InclusionProperties.h:49
Generated on Sat Jan 25 2020 12:04:46 for www.mooseframework.org by   doxygen 1.8.16