www.mooseframework.org
ComputeLayeredCosseratElasticityTensor.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 /****************************************************************/
8 #include "libmesh/utility.h"
9 #include "Function.h"
10 #include "RankTwoTensor.h"
11 
12 template <>
13 InputParameters
15 {
16  InputParameters params = validParams<ComputeElasticityTensorBase>();
17  params.addClassDescription("Computes Cosserat elasticity and flexural bending rigidity tensors "
18  "relevant for simulations with layered materials. The layering "
19  "direction is assumed to be perpendicular to the 'z' direction.");
20  params.addRequiredParam<Real>("young", "The Young's modulus");
21  params.addRequiredParam<Real>("poisson", "The Poisson's ratio");
22  params.addRequiredRangeCheckedParam<Real>(
23  "layer_thickness", "layer_thickness>=0", "The layer thickness");
24  params.addRequiredRangeCheckedParam<Real>(
25  "joint_normal_stiffness", "joint_normal_stiffness>=0", "The joint normal stiffness");
26  params.addRequiredRangeCheckedParam<Real>(
27  "joint_shear_stiffness", "joint_shear_stiffness>=0", "The joint shear stiffness");
28  return params;
29 }
30 
32  const InputParameters & parameters)
33  : ComputeElasticityTensorBase(parameters),
34  _Eijkl(RankFourTensor()),
35  _Bijkl(RankFourTensor()),
36  _Cijkl(RankFourTensor()),
37  _elastic_flexural_rigidity_tensor(
38  declareProperty<RankFourTensor>("elastic_flexural_rigidity_tensor")),
39  _compliance(declareProperty<RankFourTensor>(_base_name + "compliance_tensor"))
40 {
41  if (!isParamValid("elasticity_tensor_prefactor"))
43 
44  const Real E = getParam<Real>("young");
45  const Real nu = getParam<Real>("poisson");
46  const Real b = getParam<Real>("layer_thickness");
47  const Real kn = getParam<Real>("joint_normal_stiffness");
48  const Real ks = getParam<Real>("joint_shear_stiffness");
49 
50  // shear modulus of solid
51  const Real G = 0.5 * E / (1.0 + nu);
52  // shear modulus of jointed material
53  const Real Gprime = G * b * ks / (b * ks + G);
54 
55  const Real a0000 =
56  (b * kn > 0.0)
57  ? E / (1.0 - nu * nu - Utility::pow<2>(nu * (1.0 + nu)) / (1.0 - nu * nu + E / b / kn))
58  : E / (1.0 - nu * nu);
59  const Real a0011 = nu * a0000 / (1.0 - nu);
60  const Real a2222 =
61  (b * kn > 0.0) ? 1.0 / ((1.0 + nu) * (1.0 - 2.0 * nu) / E / (1.0 - nu) + 1.0 / b / kn) : 0.0;
62  const Real a0022 = nu * a2222 / (1.0 - nu);
63  const Real a0101 = G;
64  const Real a66 = Gprime;
65  const Real a77 = 0.5 * (G + Gprime);
66 
67  // Eijkl does not obey the usual symmetries, viz Eijkl != Ejikl, so must fill manually
68  _Eijkl(0, 0, 0, 0) = _Eijkl(1, 1, 1, 1) = a0000;
69  _Eijkl(0, 0, 1, 1) = _Eijkl(1, 1, 0, 0) = a0011;
70  _Eijkl(2, 2, 2, 2) = a2222;
71  _Eijkl(0, 0, 2, 2) = _Eijkl(1, 1, 2, 2) = _Eijkl(2, 2, 0, 0) = _Eijkl(2, 2, 1, 1) = a0022;
72  _Eijkl(0, 1, 0, 1) = _Eijkl(0, 1, 1, 0) = _Eijkl(1, 0, 0, 1) = _Eijkl(1, 0, 1, 0) = a0101;
73  _Eijkl(0, 2, 0, 2) = _Eijkl(0, 2, 2, 0) = _Eijkl(2, 0, 0, 2) = _Eijkl(1, 2, 1, 2) =
74  _Eijkl(1, 2, 2, 1) = _Eijkl(2, 1, 1, 2) = a66;
75  _Eijkl(2, 0, 2, 0) = _Eijkl(2, 1, 2, 1) = a77;
76 
77  // most of Bijkl is zero since the only nonzero moment stresses are m01 and m10.
78  // It also does not have the usual symmetries.
79  const Real D0 = E * Utility::pow<3>(b) / 12.0 / (1.0 - nu * nu); // bending rigidity of a layer
80  const Real b0101 = D0 / b * G / (2.0 * b * ks + G);
81  const Real b0110 = -nu * b0101;
82  _Bijkl(0, 1, 0, 1) = _Bijkl(1, 0, 1, 0) = b0101;
83  _Bijkl(0, 1, 1, 0) = _Bijkl(1, 0, 0, 1) = b0110;
84 
85  // The compliance tensor also does not obey the usual symmetries, and
86  // this is the main reason it is calculated here, since we can't use
87  // _Eijkl.invSymm()
88  const Real pre = (nu - 1.0) / (a0000 * (nu - 1.0) + 2.0 * a2222 * Utility::pow<2>(nu));
89  const Real cp0000 =
90  ((a2222 - a0000) * nu * nu + 2.0 * a0000 * nu - a0000) / (2.0 * a0000 * nu - a0000);
91  const Real cp0011 = -((a0000 + a2222) * nu * nu - a0000 * nu) / (2.0 * a0000 * nu - a0000);
92  _Cijkl(0, 0, 0, 0) = _Cijkl(1, 1, 1, 1) = pre * cp0000;
93  _Cijkl(0, 0, 1, 1) = _Cijkl(1, 1, 0, 0) = pre * cp0011;
94  _Cijkl(2, 2, 2, 2) = pre * a0000 / a2222;
95  _Cijkl(0, 0, 2, 2) = _Cijkl(1, 1, 2, 2) = _Cijkl(2, 2, 0, 0) = _Cijkl(2, 2, 1, 1) = -nu * pre;
96  _Cijkl(0, 1, 0, 1) = _Cijkl(0, 1, 1, 0) = _Cijkl(1, 0, 0, 1) = _Cijkl(1, 0, 1, 0) = 0.25 / a0101;
97  const Real slip = 2.0 / (G - Gprime);
98  _Cijkl(0, 2, 2, 0) = _Cijkl(2, 0, 0, 2) = _Cijkl(1, 2, 2, 1) = _Cijkl(2, 1, 1, 2) = -slip;
99  _Cijkl(0, 2, 0, 2) = _Cijkl(1, 2, 1, 2) = (G + Gprime) * slip / 2.0 / Gprime;
100  _Cijkl(2, 0, 2, 0) = _Cijkl(2, 1, 2, 1) = slip;
101 }
102 
103 void
105 {
106  _elasticity_tensor[_qp] = _Eijkl;
108  _compliance[_qp] = _Cijkl;
109 
111  _compliance[_qp] /= _prefactor_function->value(_t, _q_point[_qp]);
112 }
RankFourTensor _Cijkl
Inverse of elasticity tensor.
ComputeElasticityTensorBase the base class for computing elasticity tensors.
void issueGuarantee(const MaterialPropertyName &prop_name, Guarantee guarantee)
ComputeLayeredCosseratElasticityTensor(const InputParameters &parameters)
MaterialProperty< RankFourTensor > & _elastic_flexural_rigidity_tensor
Flexural rigidity tensor at the qps.
InputParameters validParams< ComputeElasticityTensorBase >()
MaterialProperty< RankFourTensor > & _compliance
Compliance tensor (_Eijkl^-1) at the qps.
Function *const _prefactor_function
prefactor function to multiply the elasticity tensor with
RankFourTensor _Eijkl
Conventional elasticity tensor.
InputParameters validParams< ComputeLayeredCosseratElasticityTensor >()
MaterialProperty< RankFourTensor > & _elasticity_tensor