www.mooseframework.org
RankTwoScalarTools.C
Go to the documentation of this file.
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 "RankTwoScalarTools.h"
11 
12 // MOOSE includes
13 #include "MooseEnum.h"
14 #include "RankTwoTensor.h"
15 #include "MooseError.h"
16 
17 #include "libmesh/point.h"
18 
19 namespace RankTwoScalarTools
20 {
21 
22 MooseEnum
24 {
25  return MooseEnum("VonMisesStress EffectiveStrain Hydrostatic L2norm MaxPrincipal "
26  "MidPrincipal MinPrincipal VolumetricStrain FirstInvariant SecondInvariant "
27  "ThirdInvariant AxialStress HoopStress RadialStress TriaxialityStress "
28  "Direction MaxShear StressIntensity");
29 }
30 
31 Real
32 getQuantity(const RankTwoTensor & tensor,
33  const MooseEnum scalar_type,
34  const Point & point1,
35  const Point & point2,
36  const Point & curr_point,
37  Point & direction)
38 {
39  Real val = 0.0;
40 
41  switch (scalar_type)
42  {
43  case 0:
44  val = vonMisesStress(tensor);
45  break;
46  case 1:
47  val = effectiveStrain(tensor);
48  break;
49  case 2:
50  val = hydrostatic(tensor);
51  break;
52  case 3:
53  val = L2norm(tensor);
54  break;
55  case 4:
56  val = maxPrincipal(tensor, direction);
57  break;
58  case 5:
59  val = midPrincipal(tensor, direction);
60  break;
61  case 6:
62  val = minPrincipal(tensor, direction);
63  break;
64  case 7:
65  val = volumetricStrain(tensor);
66  break;
67  case 8:
68  val = firstInvariant(tensor);
69  break;
70  case 9:
71  val = secondInvariant(tensor);
72  break;
73  case 10:
74  val = thirdInvariant(tensor);
75  break;
76  case 11:
77  val = axialStress(tensor, point1, point2, direction);
78  break;
79  case 12:
80  val = hoopStress(tensor, point1, point2, curr_point, direction);
81  break;
82  case 13:
83  val = radialStress(tensor, point1, point2, curr_point, direction);
84  break;
85  case 14:
86  val = triaxialityStress(tensor);
87  break;
88  case 15:
89  val = directionValueTensor(tensor, direction);
90  break;
91  case 16:
92  val = maxShear(tensor);
93  break;
94  case 17:
95  val = stressIntensity(tensor);
96  break;
97  default:
98  mooseError("RankTwoScalarAux Error: Pass valid scalar type - " +
99  scalarOptions().getRawNames());
100  }
101 
102  return val;
103 }
104 
105 Real
106 component(const RankTwoTensor & r2tensor, unsigned int i, unsigned int j)
107 {
108  return r2tensor(i, j);
109 }
110 
111 Real
112 component(const RankTwoTensor & r2tensor, unsigned int i, unsigned int j, Point & direction)
113 {
114  direction.zero();
115  if (i == j)
116  direction(i) = 1.0;
117  else
118  {
119  direction(i) = std::sqrt(0.5);
120  direction(j) = std::sqrt(0.5);
121  }
122 
123  return r2tensor(i, j);
124 }
125 
126 Real
128 {
129  RankTwoTensor dev_stress = stress.deviatoric();
130 
131  return std::sqrt(3.0 / 2.0 * dev_stress.doubleContraction(dev_stress));
132 }
133 
134 Real
136 {
137  return std::sqrt(2.0 / 3.0 * strain.doubleContraction(strain));
138 }
139 
140 Real
141 hydrostatic(const RankTwoTensor & r2tensor)
142 {
143  return r2tensor.trace() / 3.0;
144 }
145 
146 Real
147 L2norm(const RankTwoTensor & r2tensor)
148 {
149  return r2tensor.L2norm();
150 }
151 
152 Real
154 {
155  // Since the strains are logarithmic strains, which are by definition log(L/L0),
156  // exp(log_strain) = L/L0
157  // The ratio of the volume change of a strained cube to the original volume
158  // (delta V / V) is thus:
159  // exp(log_strain_11) * exp(log_strain_22) * exp(log_strain_33) - 1
160  //
161  // Since eng_strain = exp(log_strain) - 1, the equivalent calculation using
162  // engineering strains would be:
163  // (1 + eng_strain_11) * (1 + eng_strain_22) + (1 + eng_strain_33) - 1
164  // If strains are small, the resulting terms that involve squared and cubed
165  // strains are negligible, resulting in the following approximate form:
166  // strain_11 + strain_22 + strain_33
167  // There is not currently an option to compute this small-strain form of the
168  // volumetric strain, but at small strains, it is approximately equal to the
169  // finite strain form that is computed.
170  //
171  return std::exp(strain(0, 0)) * std::exp(strain(1, 1)) * std::exp(strain(2, 2)) - 1.0;
172 }
173 
174 Real
175 firstInvariant(const RankTwoTensor & r2tensor)
176 {
177  return r2tensor.trace();
178 }
179 
180 Real
182 {
183  Real val = 0.0;
184  for (unsigned int i = 0; i < 2; ++i)
185  for (unsigned int j = i + 1; j < 3; ++j)
186  {
187  val += r2tensor(i, i) * r2tensor(j, j);
188  val -= (r2tensor(i, j) * r2tensor(i, j) + r2tensor(j, i) * r2tensor(j, i)) * 0.5;
189  }
190 
191  return val;
192 }
193 
194 Real
195 thirdInvariant(const RankTwoTensor & r2tensor)
196 {
197  Real val = 0.0;
198  val = r2tensor(0, 0) * r2tensor(1, 1) * r2tensor(2, 2) -
199  r2tensor(0, 0) * r2tensor(1, 2) * r2tensor(2, 1) +
200  r2tensor(0, 1) * r2tensor(1, 2) * r2tensor(2, 0) -
201  r2tensor(0, 1) * r2tensor(1, 0) * r2tensor(2, 2) +
202  r2tensor(0, 2) * r2tensor(1, 0) * r2tensor(2, 1) -
203  r2tensor(0, 2) * r2tensor(1, 1) * r2tensor(2, 0);
204 
205  return val;
206 }
207 
208 Real
209 maxPrincipal(const RankTwoTensor & r2tensor, Point & direction)
210 {
211  return calcEigenValuesEigenVectors(r2tensor, (LIBMESH_DIM - 1), direction);
212 }
213 
214 Real
215 midPrincipal(const RankTwoTensor & r2tensor, Point & direction)
216 {
217  return calcEigenValuesEigenVectors(r2tensor, 1, direction);
218 }
219 
220 Real
221 minPrincipal(const RankTwoTensor & r2tensor, Point & direction)
222 {
223  return calcEigenValuesEigenVectors(r2tensor, 0, direction);
224 }
225 
226 Real
227 calcEigenValuesEigenVectors(const RankTwoTensor & r2tensor, unsigned int index, Point & eigenvec)
228 {
229  std::vector<Real> eigenval(LIBMESH_DIM);
230  RankTwoTensor eigvecs;
231  r2tensor.symmetricEigenvaluesEigenvectors(eigenval, eigvecs);
232 
233  Real val = eigenval[index];
234  eigenvec = eigvecs.column(index);
235 
236  return val;
237 }
238 
239 Real
240 axialStress(const RankTwoTensor & stress,
241  const Point & point1,
242  const Point & point2,
243  Point & direction)
244 {
245  Point axis = point2 - point1;
246  axis /= axis.norm();
247 
248  // Calculate the stress in the direction of the axis specifed by the user
249  Real axial_stress = 0.0;
250  for (unsigned int i = 0; i < 3; ++i)
251  for (unsigned int j = 0; j < 3; ++j)
252  axial_stress += axis(j) * stress(j, i) * axis(i);
253 
254  direction = axis;
255 
256  return axial_stress;
257 }
258 
259 Real
260 hoopStress(const RankTwoTensor & stress,
261  const Point & point1,
262  const Point & point2,
263  const Point & curr_point,
264  Point & direction)
265 {
266  // Calculate the cross of the normal to the axis of rotation from the current
267  // location and the axis of rotation
268  Point xp;
269  normalPositionVector(point1, point2, curr_point, xp);
270  Point axis_rotation = point2 - point1;
271  Point yp = axis_rotation / axis_rotation.norm();
272  Point zp = xp.cross(yp);
273 
274  // Calculate the scalar value of the hoop stress
275  Real hoop_stress = 0.0;
276  for (unsigned int i = 0; i < 3; ++i)
277  for (unsigned int j = 0; j < 3; ++j)
278  hoop_stress += zp(j) * stress(j, i) * zp(i);
279 
280  direction = zp;
281 
282  return hoop_stress;
283 }
284 
285 Real
287  const Point & point1,
288  const Point & point2,
289  const Point & curr_point,
290  Point & direction)
291 {
292  Point radial_norm;
293  normalPositionVector(point1, point2, curr_point, radial_norm);
294 
295  // Compute the scalar stress component in the direction of the normal vector from the
296  // user-defined axis of rotation.
297  Real radial_stress = 0.0;
298  for (unsigned int i = 0; i < 3; ++i)
299  for (unsigned int j = 0; j < 3; ++j)
300  radial_stress += radial_norm(j) * stress(j, i) * radial_norm(i);
301 
302  direction = radial_norm;
303 
304  return radial_stress;
305 }
306 
307 void
308 normalPositionVector(const Point & point1,
309  const Point & point2,
310  const Point & curr_point,
311  Point & normalPosition)
312 {
313  // Find the nearest point on the axis of rotation (defined by point2 - point1)
314  // to the current position, e.g. the normal to the axis of rotation at the
315  // current position
316  Point axis_rotation = point2 - point1;
317  Point positionWRTpoint1 = point1 - curr_point;
318  Real projection = (axis_rotation * positionWRTpoint1) / axis_rotation.norm_sq();
319  Point normal = point1 - projection * axis_rotation;
320 
321  // Calculate the direction normal to the plane formed by the axis of rotation
322  // and the normal to the axis of rotation from the current position.
323  normalPosition = curr_point - normal;
324  normalPosition /= normalPosition.norm();
325 }
326 
327 Real
328 directionValueTensor(const RankTwoTensor & r2tensor, Point & direction)
329 {
330  Real tensor_value_in_direction = 0.0;
331  for (unsigned int i = 0; i < 3; ++i)
332  for (unsigned int j = 0; j < 3; ++j)
333  tensor_value_in_direction += direction(j) * r2tensor(j, i) * direction(i);
334 
335  return tensor_value_in_direction;
336 }
337 
338 Real
340 {
341  return hydrostatic(stress) / vonMisesStress(stress);
342 }
343 
344 Real
345 maxShear(const RankTwoTensor & stress)
346 {
347  Point dummy;
348  return (maxPrincipal(stress, dummy) - minPrincipal(stress, dummy)) / 2.;
349 }
350 
351 Real
353 {
354  return 2. * maxShear(stress);
355 }
356 }
Real firstInvariant(const RankTwoTensor &r2tensor)
Real radialStress(const RankTwoTensor &stress, const Point &point1, const Point &point2, const Point &curr_point, Point &direction)
Real triaxialityStress(const RankTwoTensor &stress)
Real effectiveStrain(const RankTwoTensor &strain)
Real maxShear(const RankTwoTensor &stress)
Real maxPrincipal(const RankTwoTensor &r2tensor, Point &direction)
Real directionValueTensor(const RankTwoTensor &r2tensor, Point &direction)
Real volumetricStrain(const RankTwoTensor &strain)
Real thirdInvariant(const RankTwoTensor &r2tensor)
Real hoopStress(const RankTwoTensor &stress, const Point &point1, const Point &point2, const Point &curr_point, Point &direction)
Real vonMisesStress(const RankTwoTensor &tensor)
Real stressIntensity(const RankTwoTensor &stress)
Real minPrincipal(const RankTwoTensor &r2tensor, Point &direction)
Real component(const RankTwoTensor &r2tensor, unsigned int i, unsigned int j)
Real axialStress(const RankTwoTensor &stress, const Point &point1, const Point &point2, Point &direction)
Real calcEigenValuesEigenVectors(const RankTwoTensor &r2tensor, unsigned int index, Point &eigenvec)
Real midPrincipal(const RankTwoTensor &r2tensor, Point &direction)
Real getQuantity(const RankTwoTensor &tensor, const MooseEnum scalar_type, const Point &point1, const Point &point2, const Point &curr_point, Point &direction)
Real L2norm(const RankTwoTensor &r2tensor)
void normalPositionVector(const Point &point1, const Point &point2, const Point &curr_point, Point &normalPosition)
Real secondInvariant(const RankTwoTensor &r2tensor)
Real hydrostatic(const RankTwoTensor &r2tensor)