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 "ComputeLayeredCosseratElasticityTensor.h"
11 : #include "libmesh/utility.h"
12 : #include "Function.h"
13 : #include "RankTwoTensor.h"
14 :
15 : registerMooseObject("SolidMechanicsApp", ComputeLayeredCosseratElasticityTensor);
16 :
17 : InputParameters
18 912 : ComputeLayeredCosseratElasticityTensor::validParams()
19 : {
20 912 : InputParameters params = ComputeElasticityTensorBase::validParams();
21 912 : 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 1824 : params.addRequiredParam<Real>("young", "The Young's modulus of the layered materials");
25 1824 : params.addRequiredParam<Real>("poisson", "The Poisson's ratio of the layered materials");
26 1824 : params.addRequiredRangeCheckedParam<Real>(
27 : "layer_thickness", "layer_thickness>=0", "The layer thickness");
28 1824 : params.addRequiredRangeCheckedParam<Real>(
29 : "joint_normal_stiffness", "joint_normal_stiffness>=0", "The joint normal stiffness");
30 1824 : params.addRequiredRangeCheckedParam<Real>(
31 : "joint_shear_stiffness", "joint_shear_stiffness>=0", "The joint shear stiffness");
32 912 : return params;
33 0 : }
34 :
35 684 : ComputeLayeredCosseratElasticityTensor::ComputeLayeredCosseratElasticityTensor(
36 684 : const InputParameters & parameters)
37 : : ComputeElasticityTensorBase(parameters),
38 684 : _Eijkl(RankFourTensor()),
39 684 : _Bijkl(RankFourTensor()),
40 684 : _Cijkl(RankFourTensor()),
41 684 : _elastic_flexural_rigidity_tensor(
42 684 : declareProperty<RankFourTensor>("elastic_flexural_rigidity_tensor")),
43 1368 : _compliance(declareProperty<RankFourTensor>(_base_name + "compliance_tensor"))
44 : {
45 1368 : if (!isParamValid("elasticity_tensor_prefactor"))
46 1296 : issueGuarantee(_elasticity_tensor_name, Guarantee::CONSTANT_IN_TIME);
47 :
48 1368 : const Real E = getParam<Real>("young");
49 1368 : const Real nu = getParam<Real>("poisson");
50 1368 : const Real b = getParam<Real>("layer_thickness");
51 1368 : const Real kn = getParam<Real>("joint_normal_stiffness");
52 1368 : const Real ks = getParam<Real>("joint_shear_stiffness");
53 :
54 : // shear modulus of solid
55 684 : const Real G = 0.5 * E / (1.0 + nu);
56 : // shear modulus of jointed material
57 684 : const Real Gprime = G * b * ks / (b * ks + G);
58 :
59 : const Real a0000 =
60 684 : (b * kn > 0.0)
61 684 : ? 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 684 : const Real a0011 = nu * a0000 / (1.0 - nu);
64 : const Real a2222 =
65 684 : (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 684 : const Real a0022 = nu * a2222 / (1.0 - nu);
67 : const Real a0101 = G;
68 : const Real a66 = Gprime;
69 684 : const Real a77 = 0.5 * (G + Gprime);
70 :
71 : // Eijkl does not obey the usual symmetries, viz Eijkl != Ejikl, so must fill manually
72 684 : _Eijkl(0, 0, 0, 0) = _Eijkl(1, 1, 1, 1) = a0000;
73 684 : _Eijkl(0, 0, 1, 1) = _Eijkl(1, 1, 0, 0) = a0011;
74 684 : _Eijkl(2, 2, 2, 2) = a2222;
75 684 : _Eijkl(0, 0, 2, 2) = _Eijkl(1, 1, 2, 2) = _Eijkl(2, 2, 0, 0) = _Eijkl(2, 2, 1, 1) = a0022;
76 684 : _Eijkl(0, 1, 0, 1) = _Eijkl(0, 1, 1, 0) = _Eijkl(1, 0, 0, 1) = _Eijkl(1, 0, 1, 0) = a0101;
77 684 : _Eijkl(0, 2, 0, 2) = _Eijkl(0, 2, 2, 0) = _Eijkl(2, 0, 0, 2) = _Eijkl(1, 2, 1, 2) =
78 684 : _Eijkl(1, 2, 2, 1) = _Eijkl(2, 1, 1, 2) = a66;
79 684 : _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 684 : const Real D0 = E * Utility::pow<3>(b) / 12.0 / (1.0 - nu * nu); // bending rigidity of a layer
84 684 : const Real b0101 = D0 / b * G / (2.0 * b * ks + G);
85 684 : const Real b0110 = -nu * b0101;
86 684 : _Bijkl(0, 1, 0, 1) = _Bijkl(1, 0, 1, 0) = b0101;
87 684 : _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 684 : const Real pre = (nu - 1.0) / (a0000 * (nu - 1.0) + 2.0 * a2222 * Utility::pow<2>(nu));
93 684 : const Real cp0000 =
94 684 : ((a2222 - a0000) * nu * nu + 2.0 * a0000 * nu - a0000) / (2.0 * a0000 * nu - a0000);
95 684 : const Real cp0011 = -((a0000 + a2222) * nu * nu - a0000 * nu) / (2.0 * a0000 * nu - a0000);
96 684 : _Cijkl(0, 0, 0, 0) = _Cijkl(1, 1, 1, 1) = pre * cp0000;
97 684 : _Cijkl(0, 0, 1, 1) = _Cijkl(1, 1, 0, 0) = pre * cp0011;
98 684 : _Cijkl(2, 2, 2, 2) = pre * a0000 / a2222;
99 684 : _Cijkl(0, 0, 2, 2) = _Cijkl(1, 1, 2, 2) = _Cijkl(2, 2, 0, 0) = _Cijkl(2, 2, 1, 1) = -nu * pre;
100 684 : _Cijkl(0, 1, 0, 1) = _Cijkl(0, 1, 1, 0) = _Cijkl(1, 0, 0, 1) = _Cijkl(1, 0, 1, 0) = 0.25 / a0101;
101 684 : const Real slip = 2.0 / (G - Gprime);
102 684 : _Cijkl(0, 2, 2, 0) = _Cijkl(2, 0, 0, 2) = _Cijkl(1, 2, 2, 1) = _Cijkl(2, 1, 1, 2) = -slip;
103 684 : _Cijkl(0, 2, 0, 2) = _Cijkl(1, 2, 1, 2) = (G + Gprime) * slip / 2.0 / Gprime;
104 684 : _Cijkl(2, 0, 2, 0) = _Cijkl(2, 1, 2, 1) = slip;
105 684 : }
106 :
107 : void
108 251280 : ComputeLayeredCosseratElasticityTensor::computeQpElasticityTensor()
109 : {
110 251280 : _elasticity_tensor[_qp] = _Eijkl;
111 251280 : _elastic_flexural_rigidity_tensor[_qp] = _Bijkl;
112 251280 : _compliance[_qp] = _Cijkl;
113 :
114 251280 : if (_prefactor_function)
115 0 : _compliance[_qp] /= _prefactor_function->value(_t, _q_point[_qp]);
116 251280 : }
|
