www.mooseframework.org
MaterialTensorCalculatorTools.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 
11 
13 {
14 
15 Real
16 component(const SymmTensor & symm_tensor, unsigned int index)
17 {
18  return symm_tensor.component(index);
19 }
20 
21 Real
22 component(const SymmTensor & symm_tensor, unsigned int index, RealVectorValue & direction)
23 {
24  direction.zero();
25  if (index < 3)
26  direction(index) = 1.0;
27 
28  else if (index == 3) // xy
29  {
30  direction(0) = std::sqrt(0.5);
31  direction(1) = direction(0);
32  }
33  else if (index == 4) // yz
34  {
35  direction(1) = std::sqrt(0.5);
36  direction(2) = direction(1);
37  }
38  else if (index == 5) // zx
39  {
40  direction(0) = std::sqrt(0.5);
41  direction(2) = direction(0);
42  }
43  return symm_tensor.component(index);
44 }
45 
46 Real
47 vonMisesStress(const SymmTensor & symm_stress)
48 {
49  Real value = std::pow(symm_stress.xx() - symm_stress.yy(), 2) +
50  std::pow(symm_stress.yy() - symm_stress.zz(), 2) +
51  std::pow(symm_stress.zz() - symm_stress.xx(), 2) +
52  6 * (std::pow(symm_stress.xy(), 2) + std::pow(symm_stress.yz(), 2) +
53  std::pow(symm_stress.zx(), 2));
54  return std::sqrt(value / 2.0);
55 }
56 
57 Real
58 effectiveStrain(const SymmTensor & symm_strain)
59 {
60  return std::sqrt(2.0 / 3.0 * symm_strain.doubleContraction(symm_strain));
61 }
62 
63 Real
64 hydrostatic(const SymmTensor & symm_tensor)
65 {
66  return symm_tensor.trace() / 3.0;
67 }
68 
69 Real
70 volumetricStrain(const SymmTensor & symm_strain)
71 {
72  // Since the strains are logarithmic strains, which are by definition log(L/L0),
73  // exp(log_strain) = L/L0
74  // The ratio of the volume change of a strained cube to the original volume
75  // (delta V / V) is thus:
76  // exp(log_strain_11) * exp(log_strain_22) * exp(log_strain_33) - 1
77  //
78  // Since eng_strain = exp(log_strain) - 1, the equivalent calculation using
79  // engineering strains would be:
80  // (1 + eng_strain_11) * (1 + eng_strain_22) + (1 + eng_strain_33) - 1
81  // If strains are small, the resulting terms that involve squared and cubed
82  // strains are negligible, resulting in the following approximate form:
83  // strain_11 + strain_22 + strain_33
84  // There is not currently an option to compute this small-strain form of the
85  // volumetric strain, but at small strains, it is approximately equal to the
86  // finite strain form that is computed.
87 
88  return std::exp(symm_strain.xx()) * std::exp(symm_strain.yy()) * std::exp(symm_strain.zz()) - 1.0;
89 }
90 
91 Real
92 firstInvariant(const SymmTensor & symm_tensor)
93 {
94  return symm_tensor.trace();
95 }
96 
97 Real
98 secondInvariant(const SymmTensor & symm_tensor)
99 {
100  Real value = symm_tensor.xx() * symm_tensor.yy() + symm_tensor.yy() * symm_tensor.zz() +
101  symm_tensor.zz() * symm_tensor.xx() - symm_tensor.xy() * symm_tensor.xy() -
102  symm_tensor.yz() * symm_tensor.yz() - symm_tensor.zx() * symm_tensor.zx();
103 
104  return value;
105 }
106 
107 Real
108 thirdInvariant(const SymmTensor & symm_tensor)
109 {
110  Real val = 0.0;
111  val = symm_tensor.xx() * symm_tensor.yy() * symm_tensor.zz() -
112  symm_tensor.xx() * symm_tensor.yz() * symm_tensor.yz() +
113  symm_tensor.xy() * symm_tensor.yz() * symm_tensor.zx() -
114  symm_tensor.xy() * symm_tensor.xy() * symm_tensor.zz() +
115  symm_tensor.zx() * symm_tensor.xy() * symm_tensor.yz() -
116  symm_tensor.zx() * symm_tensor.yy() * symm_tensor.zx();
117 
118  return val;
119 }
120 
121 Real
122 maxPrincipal(const SymmTensor & symm_tensor, RealVectorValue & direction)
123 {
124  Real val = calcPrincipalValues(symm_tensor, (LIBMESH_DIM - 1), direction);
125  return val;
126 }
127 
128 Real
129 midPrincipal(const SymmTensor & symm_tensor, RealVectorValue & direction)
130 {
131  Real val = calcPrincipalValues(symm_tensor, 1, direction);
132  return val;
133 }
134 
135 Real
136 minPrincipal(const SymmTensor & symm_tensor, RealVectorValue & direction)
137 {
138  Real val = calcPrincipalValues(symm_tensor, 0, direction);
139  return val;
140 }
141 
142 // The functions below for maxPrinciple, midPrinciple and minPrinciple are
143 // deprecated. They will be replaced with the correctly spelled versions.
144 
145 Real
146 maxPrinciple(const SymmTensor & symm_tensor, RealVectorValue & direction)
147 {
148  Real val = calcPrincipalValues(symm_tensor, (LIBMESH_DIM - 1), direction);
149  return val;
150 }
151 
152 Real
153 midPrinciple(const SymmTensor & symm_tensor, RealVectorValue & direction)
154 {
155  Real val = calcPrincipalValues(symm_tensor, 1, direction);
156  return val;
157 }
158 
159 Real
160 minPrinciple(const SymmTensor & symm_tensor, RealVectorValue & direction)
161 {
162  Real val = calcPrincipalValues(symm_tensor, 0, direction);
163  return val;
164 }
165 
166 Real
167 calcPrincipalValues(const SymmTensor & symm_tensor, unsigned int index, RealVectorValue & direction)
168 {
169  ColumnMajorMatrix eval(3, 1);
170  ColumnMajorMatrix evec(3, 3);
171  symm_tensor.columnMajorMatrix().eigen(eval, evec);
172  // Eigen computes low to high. We want high first.
173  direction(0) = evec(0, index);
174  direction(1) = evec(1, index);
175  direction(2) = evec(2, index);
176  return eval(index);
177 }
178 
179 Real
180 axialStress(const SymmTensor & symm_stress,
181  const Point & point1,
182  const Point & point2,
183  RealVectorValue & direction)
184 {
185  Point axis = point2 - point1;
186  axis /= axis.norm();
187 
188  // Calculate the stress in the direction of the axis specifed by the user
189  Real axial_stress = 0.0;
190  for (unsigned int i = 0; i < 3; ++i)
191  {
192  for (unsigned int j = 0; j < 3; ++j)
193  axial_stress += axis(j) * symm_stress(j, i) * axis(i);
194  }
195  direction = axis;
196  return axial_stress;
197 }
198 
199 Real
200 hoopStress(const SymmTensor & symm_stress,
201  const Point & point1,
202  const Point & point2,
203  const Point & curr_point,
204  RealVectorValue & direction)
205 {
206  // Calculate the cross of the normal to the axis of rotation from the current
207  // location and the axis of rotation
208  Point xp;
209  normalPositionVector(point1, point2, curr_point, xp);
210 
211  Point axis_rotation = point2 - point1;
212  Point yp = axis_rotation / axis_rotation.norm();
213  Point zp = xp.cross(yp);
214 
215  // Calculate the scalar value of the hoop stress
216  Real hoop_stress = 0.0;
217  for (unsigned int i = 0; i < 3; ++i)
218  {
219  for (unsigned int j = 0; j < 3; ++j)
220  hoop_stress += zp(j) * symm_stress(j, i) * zp(i);
221  }
222  direction = zp;
223  return hoop_stress;
224 }
225 
226 Real
227 radialStress(const SymmTensor & symm_stress,
228  const Point & point1,
229  const Point & point2,
230  const Point & curr_point,
231  RealVectorValue & direction)
232 {
233  Point radial_norm;
234  normalPositionVector(point1, point2, curr_point, radial_norm);
235 
236  // Compute the scalar stress component in the direciton of the normal vector from the
237  // user-defined axis of rotation.
238  Real radial_stress = 0.0;
239  for (unsigned int i = 0; i < 3; ++i)
240  {
241  for (unsigned int j = 0; j < 3; ++j)
242  radial_stress += radial_norm(j) * symm_stress(j, i) * radial_norm(i);
243  }
244  direction = radial_norm;
245  return radial_stress;
246 }
247 
248 void
249 normalPositionVector(const Point & point1,
250  const Point & point2,
251  const Point & curr_point,
252  Point & normalPosition)
253 {
254  // Find the nearest point on the axis of rotation (defined by point2 - point1)
255  // to the current position, e.g. the normal to the axis of rotation at the
256  // current position
257  Point axis_rotation = point2 - point1;
258  Point positionWRTpoint1 = point1 - curr_point;
259  Real projection = (axis_rotation * positionWRTpoint1) / axis_rotation.norm_sq();
260  Point normal = point1 - projection * axis_rotation;
261 
262  // Calculate the direction normal to the plane formed by the axis of rotation
263  // and the normal to the axis of rotation from the current position.
264  normalPosition = curr_point - normal;
265  normalPosition /= normalPosition.norm();
266 }
267 
268 Real
269 directionValueTensor(const SymmTensor & symm_tensor, const RealVectorValue & input_direction)
270 {
271  Real tensor_value_in_direction = 0.0;
272  for (unsigned int i = 0; i < 3; ++i)
273  {
274  for (unsigned int j = 0; j < 3; ++j)
275  tensor_value_in_direction += input_direction(j) * symm_tensor(j, i) * input_direction(i);
276  }
277  return tensor_value_in_direction;
278 }
279 
280 Real
281 triaxialityStress(const SymmTensor & symm_stress)
282 {
283  Real val = hydrostatic(symm_stress) / vonMisesStress(symm_stress);
284  return val;
285 }
286 }
Real minPrincipal(const SymmTensor &symm_tensor, RealVectorValue &direction)
Real midPrincipal(const SymmTensor &symm_tensor, RealVectorValue &direction)
Real midPrinciple(const SymmTensor &symm_tensor, RealVectorValue &direction)
Real minPrinciple(const SymmTensor &symm_tensor, RealVectorValue &direction)
Real volumetricStrain(const SymmTensor &symm_strain)
Real calcPrincipalValues(const SymmTensor &symm_tensor, unsigned int index, RealVectorValue &direction)
Real hydrostatic(const SymmTensor &symm_tensor)
Real hoopStress(const SymmTensor &symm_stress, const Point &point1, const Point &point2, const Point &curr_point, RealVectorValue &direction)
Real maxPrincipal(const SymmTensor &symm_tensor, RealVectorValue &direction)
Real component(const SymmTensor &symm_tensor, unsigned int index)
Real yz() const
Definition: SymmTensor.h:135
Real zx() const
Definition: SymmTensor.h:136
Real effectiveStrain(const SymmTensor &symm_strain)
Real radialStress(const SymmTensor &symm_stress, const Point &point1, const Point &point2, const Point &curr_point, RealVectorValue &direction)
Real triaxialityStress(const SymmTensor &symm_stress)
Real secondInvariant(const SymmTensor &symm_tensor)
Real component(unsigned int i) const
Definition: SymmTensor.h:99
Real vonMisesStress(const SymmTensor &symm_stress)
Real yy() const
Definition: SymmTensor.h:132
Real xx() const
Definition: SymmTensor.h:131
ExpressionBuilder::EBTerm pow(const ExpressionBuilder::EBTerm &left, T exponent)
ColumnMajorMatrix columnMajorMatrix() const
Definition: SymmTensor.h:425
void normalPositionVector(const Point &point1, const Point &point2, const Point &curr_point, Point &normalPosition)
Real doubleContraction(const SymmTensor &rhs) const
Definition: SymmTensor.h:259
Real trace() const
Definition: SymmTensor.h:97
Real firstInvariant(const SymmTensor &symm_tensor)
Real xy() const
Definition: SymmTensor.h:134
Real zz() const
Definition: SymmTensor.h:133
Real thirdInvariant(const SymmTensor &symm_tensor)
Real maxPrinciple(const SymmTensor &symm_tensor, RealVectorValue &direction)
Real axialStress(const SymmTensor &symm_stress, const Point &point1, const Point &point2, RealVectorValue &direction)
Real directionValueTensor(const SymmTensor &symm_tensor, const RealVectorValue &input_direction)