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