www.mooseframework.org
Nonlinear.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 "Nonlinear.h"
9 #include "SolidModel.h"
10 #include "Problem.h"
12 
13 #include "libmesh/quadrature.h"
14 
15 namespace SolidMechanics
16 {
17 
19  const std::string & name,
20  const InputParameters & parameters)
21  : Element(solid_model, name, parameters),
22  _decomp_method(RashidApprox),
23  _incremental_rotation(3, 3),
24  _Uhat(3, 3)
25 {
26 
27  std::string increment_calculation = solid_model.getParam<std::string>("increment_calculation");
28  std::transform(increment_calculation.begin(),
29  increment_calculation.end(),
30  increment_calculation.begin(),
31  ::tolower);
32  if (increment_calculation == "rashidapprox")
33  {
35  }
36  else if (increment_calculation == "eigen")
37  {
39  }
40  else
41  {
42  mooseError("The options for the increment calculation are RashidApprox and Eigen.");
43  }
44 }
45 
47 
49 
51 
52 void
54  SymmTensor & strain_increment)
55 {
57  {
58  computeStrainIncrement(Fhat, strain_increment);
60  }
61 
62  else if (_decomp_method == Eigen)
63  {
65  }
66 
67  else
68  {
69  mooseError("Unknown polar decomposition type!");
70  }
71 }
72 
74 
75 void
76 Nonlinear::computeStrainIncrement(const ColumnMajorMatrix & Fhat, SymmTensor & strain_increment)
77 {
78 
79  //
80  // A calculation of the strain at the mid-interval is probably more
81  // accurate (second vs. first order). This would require the
82  // incremental deformation gradient at the mid-step, which we
83  // currently don't have. We would then have to calculate the
84  // rotation for the whole step.
85  //
86  //
87  // We are looking for:
88  // log( Uhat )
89  // = log( sqrt( Fhat^T*Fhat ) )
90  // = log( sqrt( Chat ) )
91  // A Taylor series expansion gives:
92  // ( Chat - 0.25 * Chat^T*Chat - 0.75 * I )
93  // = ( - 0.25 * Chat^T*Chat + Chat - 0.75 * I )
94  // = ( (0.25*Chat - 0.75*I) * (Chat - I) )
95  // = ( B * A )
96  // B
97  // = 0.25*Chat - 0.75*I
98  // = 0.25*(Chat - I) - 0.5*I
99  // = 0.25*A - 0.5*I
100  //
101 
102  const Real Uxx = Fhat(0, 0);
103  const Real Uxy = Fhat(0, 1);
104  const Real Uxz = Fhat(0, 2);
105  const Real Uyx = Fhat(1, 0);
106  const Real Uyy = Fhat(1, 1);
107  const Real Uyz = Fhat(1, 2);
108  const Real Uzx = Fhat(2, 0);
109  const Real Uzy = Fhat(2, 1);
110  const Real Uzz = Fhat(2, 2);
111 
112  const Real Axx = Uxx * Uxx + Uyx * Uyx + Uzx * Uzx - 1.0;
113  const Real Axy = Uxx * Uxy + Uyx * Uyy + Uzx * Uzy;
114  const Real Axz = Uxx * Uxz + Uyx * Uyz + Uzx * Uzz;
115  const Real Ayy = Uxy * Uxy + Uyy * Uyy + Uzy * Uzy - 1.0;
116  const Real Ayz = Uxy * Uxz + Uyy * Uyz + Uzy * Uzz;
117  const Real Azz = Uxz * Uxz + Uyz * Uyz + Uzz * Uzz - 1.0;
118 
119  const Real Bxx = 0.25 * Axx - 0.5;
120  const Real Bxy = 0.25 * Axy;
121  const Real Bxz = 0.25 * Axz;
122  const Real Byy = 0.25 * Ayy - 0.5;
123  const Real Byz = 0.25 * Ayz;
124  const Real Bzz = 0.25 * Azz - 0.5;
125 
126  strain_increment.xx(-(Bxx * Axx + Bxy * Axy + Bxz * Axz));
127  strain_increment.xy(-(Bxx * Axy + Bxy * Ayy + Bxz * Ayz));
128  strain_increment.zx(-(Bxx * Axz + Bxy * Ayz + Bxz * Azz));
129  strain_increment.yy(-(Bxy * Axy + Byy * Ayy + Byz * Ayz));
130  strain_increment.yz(-(Bxy * Axz + Byy * Ayz + Byz * Azz));
131  strain_increment.zz(-(Bxz * Axz + Byz * Ayz + Bzz * Azz));
132 }
133 
135 
136 void
137 Nonlinear::computePolarDecomposition(const ColumnMajorMatrix & Fhat)
138 {
139 
140  // From Rashid, 1993.
141 
142  const Real Uxx = Fhat(0, 0);
143  const Real Uxy = Fhat(0, 1);
144  const Real Uxz = Fhat(0, 2);
145  const Real Uyx = Fhat(1, 0);
146  const Real Uyy = Fhat(1, 1);
147  const Real Uyz = Fhat(1, 2);
148  const Real Uzx = Fhat(2, 0);
149  const Real Uzy = Fhat(2, 1);
150  const Real Uzz = Fhat(2, 2);
151 
152  const Real Ax = Uyz - Uzy;
153  const Real Ay = Uzx - Uxz;
154  const Real Az = Uxy - Uyx;
155  const Real Q = 0.25 * (Ax * Ax + Ay * Ay + Az * Az);
156  const Real traceF = Uxx + Uyy + Uzz;
157  const Real P = 0.25 * (traceF - 1) * (traceF - 1);
158  const Real Y = 1 / ((Q + P) * (Q + P) * (Q + P));
159 
160  const Real C1 = std::sqrt(P * (1 + (P * (Q + Q + (Q + P))) * (1 - (Q + P)) * Y));
161  const Real C2 = 0.125 + Q * 0.03125 * (P * P - 12 * (P - 1)) / (P * P);
162  const Real C3 = 0.5 * std::sqrt((P * Q * (3 - Q) + P * P * P + Q * Q) * Y);
163 
164  // Since the input to this routine is the incremental deformation gradient
165  // and not the inverse incremental gradient, this result is the transpose
166  // of the one in Rashid's paper.
167  _incremental_rotation(0, 0) = C1 + (C2 * Ax) * Ax;
168  _incremental_rotation(0, 1) = (C2 * Ay) * Ax + (C3 * Az);
169  _incremental_rotation(0, 2) = (C2 * Az) * Ax - (C3 * Ay);
170  _incremental_rotation(1, 0) = (C2 * Ax) * Ay - (C3 * Az);
171  _incremental_rotation(1, 1) = C1 + (C2 * Ay) * Ay;
172  _incremental_rotation(1, 2) = (C2 * Az) * Ay + (C3 * Ax);
173  _incremental_rotation(2, 0) = (C2 * Ax) * Az + (C3 * Ay);
174  _incremental_rotation(2, 1) = (C2 * Ay) * Az - (C3 * Ax);
175  _incremental_rotation(2, 2) = C1 + (C2 * Az) * Az;
176 }
177 
179 
180 void
181 Nonlinear::finalizeStress(std::vector<SymmTensor *> & t)
182 {
183  // Using the incremental rotation, update the stress to the current configuration (R*T*R^T)
184  for (unsigned i = 0; i < t.size(); ++i)
185  {
187  }
188 }
189 
191 
192 void
193 Nonlinear::computeStrain(const unsigned qp,
194  const SymmTensor & total_strain_old,
195  SymmTensor & total_strain_new,
196  SymmTensor & strain_increment)
197 {
198  computeStrainAndRotationIncrement(_Fhat[qp], strain_increment);
199 
200  total_strain_new = strain_increment;
201  total_strain_new += total_strain_old;
202 }
203 
205 
206 void
208 {
209  _Fhat.resize(_solid_model.qrule()->n_points());
210 
212 }
213 
215 }
DecompMethod _decomp_method
Definition: Nonlinear.h:40
Real yy() const
Definition: SymmTensor.h:130
static void polarDecompositionEigen(const ColumnMajorMatrix &Fhat, ColumnMajorMatrix &Rhat, SymmTensor &strain_increment)
Definition: Element.C:182
Element is the base class for all of this module&#39;s solid mechanics element formulations.
Definition: Element.h:23
void computeStrainIncrement(const ColumnMajorMatrix &Fhat, SymmTensor &strain_increment)
Definition: Nonlinear.C:76
Nonlinear(SolidModel &solid_model, const std::string &name, const InputParameters &parameters)
Definition: Nonlinear.C:18
Real xx() const
Definition: SymmTensor.h:129
void computePolarDecomposition(const ColumnMajorMatrix &Fhat)
Definition: Nonlinear.C:137
static void rotateSymmetricTensor(const ColumnMajorMatrix &R, const RealTensorValue &T, RealTensorValue &result)
Definition: Element.C:78
SolidModel is the base class for all this module&#39;s solid mechanics material models.
Definition: SolidModel.h:31
QBase * qrule()
Definition: SolidModel.h:53
Real zz() const
Definition: SymmTensor.h:131
ColumnMajorMatrix _incremental_rotation
Definition: Nonlinear.h:42
virtual void computeStrain(const unsigned qp, const SymmTensor &total_strain_old, SymmTensor &total_strain_new, SymmTensor &strain_increment)
Definition: Nonlinear.C:193
void computeStrainAndRotationIncrement(const ColumnMajorMatrix &Fhat, SymmTensor &strain_increment)
Definition: Nonlinear.C:53
virtual void computeIncrementalDeformationGradient(std::vector< ColumnMajorMatrix > &Fhat)=0
const std::vector< ColumnMajorMatrix > & Fhat() const
Definition: Nonlinear.h:31
Real xy() const
Definition: SymmTensor.h:132
SolidModel & _solid_model
Definition: Element.h:72
Real yz() const
Definition: SymmTensor.h:133
Real zx() const
Definition: SymmTensor.h:134
std::vector< ColumnMajorMatrix > _Fhat
Definition: Nonlinear.h:45
virtual void init()
Definition: Nonlinear.C:207
virtual void finalizeStress(std::vector< SymmTensor * > &t)
Rotate stress to current configuration.
Definition: Nonlinear.C:181