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 "NodalRankTwoPD.h"
11 : #include "PeridynamicsMesh.h"
12 : #include "RankTwoScalarTools.h"
13 :
14 : registerMooseObject("PeridynamicsApp", NodalRankTwoPD);
15 :
16 : InputParameters
17 234 : NodalRankTwoPD::validParams()
18 : {
19 234 : InputParameters params = AuxKernelBasePD::validParams();
20 234 : params.addClassDescription(
21 : "Class for computing and outputing components and scalar quantities "
22 : "of nodal rank two strain and stress tensors for bond-based and ordinary "
23 : "state-based peridynamic models");
24 :
25 468 : params.addRequiredCoupledVar("displacements", "Nonlinear variable names for the displacements");
26 468 : params.addCoupledVar("temperature", "Nonlinear variable name for the temperature");
27 468 : params.addCoupledVar("scalar_out_of_plane_strain",
28 : "Scalar variable for strain in the out-of-plane direction");
29 468 : params.addParam<bool>(
30 : "plane_stress",
31 468 : false,
32 : "Plane stress problem or not. This option applies to BPD and OSPD models only");
33 468 : params.addParam<Real>("youngs_modulus", "Material constant: Young's modulus");
34 468 : params.addParam<Real>("poissons_ratio", "Material constant: Poisson's ratio");
35 468 : params.addParam<Real>("thermal_expansion_coeff",
36 : "Value of material thermal expansion coefficient");
37 468 : params.addParam<Real>("stress_free_temperature", 0.0, "Stress free temperature");
38 468 : params.addRequiredParam<std::string>(
39 : "rank_two_tensor",
40 : "Parameter to set which rank two tensor: total_strain, mechanical_strain or stress");
41 468 : params.addRequiredParam<std::string>("output_type", "Type of output: component or scalar");
42 468 : params.addParam<MooseEnum>(
43 468 : "scalar_type", RankTwoScalarTools::scalarOptions(), "Type of scalar output");
44 468 : params.addParam<unsigned int>("index_i", "The index i of ij for the tensor to output (0, 1, 2)");
45 468 : params.addParam<unsigned int>("index_j", "The index j of ij for the tensor to output (0, 1, 2)");
46 234 : params.addParam<Point>("point1",
47 234 : Point(0, 0, 0),
48 : "Start point for axis used to calculate some direction dependent material "
49 : "tensor scalar quantities");
50 234 : params.addParam<Point>("point2",
51 234 : Point(1, 0, 0),
52 : "End point for axis used to calculate some direction dependent material "
53 : "tensor scalar quantities");
54 :
55 234 : return params;
56 0 : }
57 :
58 128 : NodalRankTwoPD::NodalRankTwoPD(const InputParameters & parameters)
59 : : AuxKernelBasePD(parameters),
60 128 : _has_temp(isCoupled("temperature")),
61 128 : _temp_var(_has_temp ? getVar("temperature", 0) : nullptr),
62 128 : _bond_status_var(&_subproblem.getStandardVariable(_tid, "bond_status")),
63 128 : _scalar_out_of_plane_strain_coupled(isCoupledScalar("scalar_out_of_plane_strain")),
64 256 : _scalar_out_of_plane_strain(_scalar_out_of_plane_strain_coupled
65 128 : ? coupledScalarValue("scalar_out_of_plane_strain")
66 : : _zero),
67 256 : _plane_stress(getParam<bool>("plane_stress")),
68 256 : _temp_ref(getParam<Real>("stress_free_temperature")),
69 256 : _rank_two_tensor(getParam<std::string>("rank_two_tensor")),
70 256 : _output_type(getParam<std::string>("output_type")),
71 384 : _scalar_type(getParam<MooseEnum>("scalar_type")),
72 128 : _point1(parameters.get<Point>("point1")),
73 256 : _point2(parameters.get<Point>("point2"))
74 : {
75 128 : if (!_var.isNodal())
76 0 : mooseError("NodalRankTwoPD operates on nodal variable!");
77 :
78 128 : if (coupledComponents("displacements") != _dim)
79 0 : mooseError("Number of displacements should be the same as problem dimension!");
80 768 : for (unsigned int i = 0; i < coupledComponents("displacements"); ++i)
81 512 : _disp_var.push_back(getVar("displacements", i));
82 :
83 128 : _Cijkl.zero();
84 128 : if (_rank_two_tensor == "stress")
85 : {
86 296 : if (!isParamValid("youngs_modulus") || !isParamValid("poissons_ratio"))
87 0 : mooseError(
88 : "Both Young's modulus and Poisson's ratio must be provided for stress calculation");
89 : else
90 : {
91 148 : _youngs_modulus = getParam<Real>("youngs_modulus");
92 148 : _poissons_ratio = getParam<Real>("poissons_ratio");
93 : }
94 :
95 74 : std::vector<Real> iso_const(2);
96 74 : iso_const[0] = _youngs_modulus * _poissons_ratio /
97 74 : ((1.0 + _poissons_ratio) * (1.0 - 2.0 * _poissons_ratio));
98 74 : iso_const[1] = _youngs_modulus / (2.0 * (1.0 + _poissons_ratio));
99 :
100 : // fill elasticity tensor
101 74 : _Cijkl.fillFromInputVector(iso_const, RankFourTensor::symmetric_isotropic);
102 74 : }
103 :
104 446 : if (_has_temp && !isParamValid("thermal_expansion_coeff"))
105 0 : mooseError("thermal_expansion_coeff is required when temperature is coupled");
106 128 : else if (_has_temp)
107 212 : _alpha = getParam<Real>("thermal_expansion_coeff");
108 :
109 128 : if (_output_type == "component")
110 : {
111 488 : if (!isParamValid("index_i") || !isParamValid("index_j"))
112 0 : mooseError("The output_type is 'component', both 'index_i' and 'index_j' must be provided!");
113 : else
114 : {
115 244 : _i = getParam<unsigned int>("index_i");
116 244 : _j = getParam<unsigned int>("index_j");
117 : }
118 : }
119 :
120 146 : if (_output_type == "scalar" && !isParamValid("scalar_type"))
121 0 : mooseError("The output_type is 'scalar', but no 'scalar_type' is provided!");
122 128 : }
123 :
124 : Real
125 131120 : NodalRankTwoPD::computeValue()
126 : {
127 :
128 : Real val = 0.0;
129 131120 : RankTwoTensor tensor;
130 : Point p = Point(0, 0, 0);
131 :
132 131120 : computeRankTwoTensors();
133 :
134 131120 : if (_rank_two_tensor == "total_strain")
135 32176 : tensor = _total_strain;
136 98944 : else if (_rank_two_tensor == "mechanical_strain")
137 29184 : tensor = _mechanical_strain;
138 69760 : else if (_rank_two_tensor == "stress")
139 69760 : tensor = _stress;
140 : else
141 0 : mooseError("NodalRankTwoPD Error: Pass valid rank two tensor: total_strain, "
142 : "mechanical_strain or stress");
143 :
144 131120 : if (_output_type == "component")
145 123824 : val = RankTwoScalarTools::component(tensor, _i, _j);
146 7296 : else if (_output_type == "scalar")
147 7296 : val = RankTwoScalarTools::getQuantity(tensor, _scalar_type, _point1, _point2, _q_point[_qp], p);
148 : else
149 0 : mooseError("NodalRankTwoPD Error: Pass valid output type - component or scalar");
150 :
151 131120 : return val;
152 : }
153 :
154 : void
155 131120 : NodalRankTwoPD::computeRankTwoTensors()
156 : {
157 : _total_strain.zero();
158 : _mechanical_strain.zero();
159 : _stress.zero();
160 :
161 131120 : std::vector<dof_id_type> neighbors = _pdmesh.getNeighbors(_current_node->id());
162 131120 : std::vector<dof_id_type> bonds = _pdmesh.getBonds(_current_node->id());
163 131120 : Real horizon = _pdmesh.getHorizon(_current_node->id());
164 :
165 : // calculate the shape tensor and prepare the deformation gradient tensor
166 131120 : RankTwoTensor shape, dgrad, delta(RankTwoTensor::initIdentity), thermal_strain;
167 : shape.zero();
168 : dgrad.zero();
169 131120 : if (_dim == 2)
170 131120 : shape(2, 2) = dgrad(2, 2) = 1.0;
171 :
172 : thermal_strain.zero();
173 :
174 : Real vol_nb, weight;
175 : RealGradient origin_vec, current_vec;
176 :
177 2400388 : for (unsigned int nb = 0; nb < neighbors.size(); ++nb)
178 2269268 : if (_bond_status_var->getElementalValue(_pdmesh.elemPtr(bonds[nb])) > 0.5)
179 : {
180 2264852 : Node * neighbor_nb = _pdmesh.nodePtr(neighbors[nb]);
181 2264852 : vol_nb = _pdmesh.getNodeVolume(neighbors[nb]);
182 2264852 : origin_vec =
183 2264852 : _pdmesh.getNodeCoord(neighbor_nb->id()) - _pdmesh.getNodeCoord(_current_node->id());
184 :
185 6794556 : for (unsigned int k = 0; k < _dim; ++k)
186 4529704 : current_vec(k) = origin_vec(k) + _disp_var[k]->getNodalValue(*neighbor_nb) -
187 4529704 : _disp_var[k]->getNodalValue(*_current_node);
188 :
189 2264852 : weight = horizon / origin_vec.norm();
190 :
191 6794556 : for (unsigned int k = 0; k < _dim; ++k)
192 13589112 : for (unsigned int l = 0; l < _dim; ++l)
193 : {
194 9059408 : shape(k, l) += weight * origin_vec(k) * origin_vec(l) * vol_nb;
195 9059408 : dgrad(k, l) += weight * current_vec(k) * origin_vec(l) * vol_nb;
196 : }
197 : }
198 :
199 : // finalize the deformation gradient tensor
200 131120 : if (shape.det() != 0.)
201 : {
202 130936 : dgrad *= shape.inverse();
203 :
204 : // the green-lagrange strain tensor
205 130936 : _total_strain = 0.5 * (dgrad.transpose() * dgrad - delta);
206 :
207 130936 : if (_scalar_out_of_plane_strain_coupled)
208 36088 : _total_strain(2, 2) = _scalar_out_of_plane_strain[0];
209 94848 : else if (_plane_stress)
210 : {
211 0 : if (_has_temp)
212 : {
213 : Real mstrain00 =
214 0 : (_total_strain(0, 0) - _alpha * (_temp_var->getNodalValue(*_current_node) - _temp_ref));
215 : Real mstrain11 =
216 0 : (_total_strain(1, 1) - _alpha * (_temp_var->getNodalValue(*_current_node) - _temp_ref));
217 : Real mstrain22 =
218 0 : -(_Cijkl(2, 2, 0, 0) * mstrain00 + _Cijkl(2, 2, 1, 1) * mstrain11) / _Cijkl(2, 2, 2, 2);
219 :
220 0 : _total_strain(2, 2) =
221 0 : mstrain22 + _alpha * (_temp_var->getNodalValue(*_current_node) - _temp_ref);
222 : }
223 : else
224 0 : _total_strain(2, 2) =
225 0 : -(_Cijkl(2, 2, 0, 0) * _total_strain(0, 0) + _Cijkl(2, 2, 1, 1) * _total_strain(1, 1)) /
226 0 : _Cijkl(2, 2, 2, 2);
227 : }
228 :
229 130936 : if (_has_temp)
230 122820 : thermal_strain = _alpha * (_temp_var->getNodalValue(*_current_node) - _temp_ref) * delta;
231 :
232 130936 : _mechanical_strain = _total_strain - thermal_strain;
233 :
234 130936 : _stress = _Cijkl * _mechanical_strain;
235 : }
236 131120 : }
|
