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