www.mooseframework.org
NonlinearPlaneStrain.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 "NonlinearPlaneStrain.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  : Nonlinear(solid_model, name, parameters),
22  ScalarCoupleable(&solid_model),
23  _grad_disp_x(coupledGradient("disp_x")),
24  _grad_disp_y(coupledGradient("disp_y")),
25  _have_strain_zz(isCoupled("strain_zz")),
26  _strain_zz(_have_strain_zz ? coupledValue("strain_zz") : _zero),
27  _have_scalar_strain_zz(isCoupledScalar("scalar_strain_zz")),
28  _scalar_strain_zz(_have_scalar_strain_zz ? coupledScalarValue("scalar_strain_zz") : _zero),
29  _grad_disp_x_old(coupledGradientOld("disp_x")),
30  _grad_disp_y_old(coupledGradientOld("disp_y")),
31  _strain_zz_old(_have_strain_zz ? coupledValueOld("strain_zz") : _zero),
32  _scalar_strain_zz_old(_have_scalar_strain_zz ? coupledScalarValueOld("scalar_strain_zz")
33  : _zero),
34  _volumetric_locking_correction(_solid_model.getParam<bool>("volumetric_locking_correction"))
35 {
36 }
37 
39 
41 
43 
44 void
46 {
47  // A = grad(u(k+1) - u(k))
48  // Fbar = 1 + grad(u(k))
49  // Fhat = 1 + A*(Fbar^-1)
50  ColumnMajorMatrix A;
51  ColumnMajorMatrix Fbar;
52  ColumnMajorMatrix Fbar_inverse;
53  ColumnMajorMatrix Fhat_average;
54  Real volume(0);
55 
56  _Fbar.resize(_solid_model.qrule()->n_points());
57 
58  for (unsigned qp = 0; qp < _solid_model.qrule()->n_points(); ++qp)
59  {
60  Real strain_zz, strain_zz_old;
61  if (_have_strain_zz)
62  {
63  strain_zz = _strain_zz[qp];
64  strain_zz_old = _strain_zz_old[qp];
65  }
66  else if (_have_scalar_strain_zz && _scalar_strain_zz.size() > 0)
67  {
68  strain_zz = _scalar_strain_zz[qp];
69  strain_zz_old = _scalar_strain_zz_old[qp];
70  }
71  else
72  {
73  strain_zz = 0.0;
74  strain_zz_old = 0.0;
75  }
76 
77  fillMatrix(qp, _grad_disp_x, _grad_disp_y, strain_zz, A);
78  fillMatrix(qp, _grad_disp_x_old, _grad_disp_y_old, strain_zz_old, Fbar);
79 
80  A -= Fbar;
81 
82  Fbar.addDiag(1);
83 
84  _Fbar[qp] = Fbar;
85 
86  // Get Fbar^(-1)
87  // Computing the inverse is generally a bad idea.
88  // It's better to compute LU factors. For now at least, we'll take
89  // a direct route.
90 
91  invertMatrix(Fbar, Fbar_inverse);
92 
93  Fhat[qp] = A * Fbar_inverse;
94  Fhat[qp].addDiag(1);
95 
97  {
98  // Now include the contribution for the integration of Fhat over the element
99  Fhat_average += Fhat[qp] * _solid_model.JxW(qp);
100 
101  volume += _solid_model.JxW(qp); // Accumulate original configuration volume
102  }
103  }
104 
106  {
107  Fhat_average /= volume;
108  const Real det_Fhat_average(detMatrix(Fhat_average));
109 
110  // Finalize volumetric locking correction
111  for (unsigned qp = 0; qp < _solid_model.qrule()->n_points(); ++qp)
112  {
113  const Real det_Fhat(detMatrix(Fhat[qp]));
114  const Real factor(std::cbrt(det_Fhat_average / det_Fhat));
115 
116  Fhat[qp] *= factor;
117  }
118  }
119  // Moose::out << "Fhat(0,0)" << Fhat[0](0,0) << std::endl;
120 }
121 
123 
124 void
125 NonlinearPlaneStrain::computeDeformationGradient(unsigned int qp, ColumnMajorMatrix & F)
126 {
127  mooseAssert(F.n() == 3 && F.m() == 3, "computeDefGrad requires 3x3 matrix");
128 
129  F(0, 0) = _grad_disp_x[qp](0) + 1;
130  F(0, 1) = _grad_disp_x[qp](1);
131  F(0, 2) = 0;
132  F(1, 0) = _grad_disp_y[qp](0);
133  F(1, 1) = _grad_disp_y[qp](1) + 1;
134  F(1, 2) = 0;
135  F(2, 0) = 0;
136  F(2, 1) = 0;
137  if (_have_strain_zz)
138  F(2, 2) = _strain_zz[qp] + 1;
139  else if (_have_scalar_strain_zz && _scalar_strain_zz.size() > 0)
140  F(2, 2) = _scalar_strain_zz[qp] + 1;
141  else
142  F(2, 2) = 1;
143 }
144 
146 
147 void
149  const VariableGradient & grad_x,
150  const VariableGradient & grad_y,
151  const Real & strain_zz,
152  ColumnMajorMatrix & A) const
153 {
154  mooseAssert(A.n() == 3 && A.m() == 3, "computeDefGrad requires 3x3 matrix");
155 
156  A(0, 0) = grad_x[qp](0);
157  A(0, 1) = grad_x[qp](1);
158  A(0, 2) = 0;
159  A(1, 0) = grad_y[qp](0);
160  A(1, 1) = grad_y[qp](1);
161  A(1, 2) = 0;
162  A(2, 0) = 0;
163  A(2, 1) = 0;
164  A(2, 2) = strain_zz;
165 }
166 
168 
169 Real
171 {
172  Real strain_zz_old;
173  if (_have_strain_zz)
174  strain_zz_old = _strain_zz_old[qp];
175  else if (_have_scalar_strain_zz && _scalar_strain_zz.size() > 0)
176  strain_zz_old = _scalar_strain_zz_old[qp];
177  else
178  strain_zz_old = 0.0;
179 
180  ColumnMajorMatrix Fnm1T;
181  fillMatrix(qp, _grad_disp_x_old, _grad_disp_y_old, strain_zz_old, Fnm1T);
182  Fnm1T.addDiag(1);
183 
184  return detMatrix(Fnm1T);
185 }
186 
188 }
Nonlinear is the base class for all large strain/rotation models.
Definition: Nonlinear.h:22
virtual void computeIncrementalDeformationGradient(std::vector< ColumnMajorMatrix > &Fhat)
std::vector< ColumnMajorMatrix > _Fbar
Definition: Nonlinear.h:46
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
NonlinearPlaneStrain(SolidModel &solid_model, const std::string &name, const InputParameters &parameters)
virtual Real volumeRatioOld(unsigned qp) const
virtual void fillMatrix(unsigned int qp, const VariableGradient &grad_x, const VariableGradient &grad_y, const Real &strain_zz, ColumnMajorMatrix &A) const
static Real detMatrix(const ColumnMajorMatrix &A)
Definition: Element.C:28
Real JxW(unsigned i) const
Definition: SolidModel.h:55
const std::vector< ColumnMajorMatrix > & Fhat() const
Definition: Nonlinear.h:31
SolidModel & _solid_model
Definition: Element.h:72
const VariableGradient & _grad_disp_y_old
static void invertMatrix(const ColumnMajorMatrix &A, ColumnMajorMatrix &Ainv)
Definition: Element.C:49
virtual void computeDeformationGradient(unsigned int qp, ColumnMajorMatrix &F)
const VariableGradient & _grad_disp_x_old