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