www.mooseframework.org
ACInterfaceStress.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 "ACInterfaceStress.h"
11 #include "RankTwoTensor.h"
12 #include "RankThreeTensor.h"
13 
14 registerMooseObject("PhaseFieldApp", ACInterfaceStress);
15 
18 {
20  params.addClassDescription("Interface stress driving force Allen-Cahn Kernel");
21  params.addParam<MaterialPropertyName>("mob_name", "L", "The mobility used with the kernel");
22  params.addParam<std::string>("base_name", "Material property base name");
23  params.addRequiredParam<Real>("stress", "Planar stress");
24  params.addRangeCheckedParam<Real>("op_range",
25  1.0,
26  "op_range > 0.0",
27  "Range over which order parameters change across an "
28  "interface. By default order parameters are assumed to "
29  "vary from 0 to 1");
30  return params;
31 }
32 
34  : Kernel(parameters),
35  _L(getMaterialProperty<Real>("mob_name")),
36  _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
37  _strain(getMaterialPropertyByName<RankTwoTensor>(_base_name + "elastic_strain")),
38  _stress(getParam<Real>("stress") / getParam<Real>("op_range"))
39 {
40 }
41 
42 Real
44 {
45  // no interface, return zero stress
46  const Real grad_norm_sq = _grad_u[_qp].norm_sq();
47  if (grad_norm_sq < libMesh::TOLERANCE)
48  return 0.0;
49 
50  const Real grad_norm = std::sqrt(grad_norm_sq);
51 
52  const Real nx = _grad_u[_qp](0);
53  const Real ny = _grad_u[_qp](1);
54  const Real nz = _grad_u[_qp](2);
55 
56  const Real s = _stress / grad_norm;
57  const Real ds = -_stress / (grad_norm * grad_norm_sq);
58  const Real dsx = ds * nx;
59  const Real dsy = ds * ny;
60  const Real dsz = ds * nz;
61 
62  // d/d(grad eta)_x
63  _dS(0, 0, 0) = (ny * ny + nz * nz) * dsx; // (ny * ny + nz * nz) * s;
64  _dS(0, 1, 0) = _dS(0, 0, 1) = -ny * s - nx * ny * dsx; // -nx * ny * s;
65  _dS(0, 1, 1) = 2.0 * nx * s + (nx * nx + nz * nz) * dsx; // (nx * nx + nz * nz) * s;
66  _dS(0, 2, 0) = _dS(0, 0, 2) = -nz * s - nx * nz * dsx; // -nx * nz * s;
67  _dS(0, 2, 1) = _dS(0, 1, 2) = -ny * nz * dsx; // -ny * nz * s;
68  _dS(0, 2, 2) = 2.0 * nx * s + (nx * nx + ny * ny) * dsx; // (nx * nx + ny * ny) * s;
69 
70  // d/d(grad eta)_y
71  _dS(1, 0, 0) = 2.0 * ny * s + (ny * ny + nz * nz) * dsy; // (ny * ny + nz * nz) * s;
72  _dS(1, 1, 0) = _dS(1, 0, 1) = -nx * s - nx * ny * dsy; // -nx * ny * s;
73  _dS(1, 1, 1) = (nx * nx + nz * nz) * dsy; // (nx * nx + nz * nz) * s;
74  _dS(1, 2, 0) = _dS(1, 0, 2) = -nx * nz * dsy; // -nx * nz * s;
75  _dS(1, 2, 1) = _dS(1, 1, 2) = -nz * s - ny * nz * dsy; // -ny * nz * s;
76  _dS(1, 2, 2) = 2.0 * ny * s + (nx * nx + ny * ny) * dsy; // (nx * nx + ny * ny) * s;
77 
78  // d/d(grad eta)_z
79  _dS(2, 0, 0) = 2.0 * nz * s + (ny * ny + nz * nz) * dsz; // (ny * ny + nz * nz) * s;
80  _dS(2, 1, 0) = _dS(2, 0, 1) = -nx * ny * dsz; // -nx * ny * s;
81  _dS(2, 1, 1) = 2.0 * nz * s + (nx * nx + nz * nz) * dsz; // (nx * nx + nz * nz) * s;
82  _dS(2, 2, 0) = _dS(2, 0, 2) = -nx * s - nx * nz * dsz; // -nx * nz * s;
83  _dS(2, 2, 1) = _dS(2, 1, 2) = -ny * s - ny * nz * dsz; // -ny * nz * s;
84  _dS(2, 2, 2) = (nx * nx + ny * ny) * dsz; // (nx * nx + ny * ny) * s;
85 
86  return _L[_qp] * 0.5 * _dS.doubleContraction(_strain[_qp]) * _grad_test[_i][_qp];
87 }
88 
89 Real
91 {
92  // no interface, return zero stress
93  const Real grad_norm_sq = _grad_u[_qp].norm_sq();
94  if (grad_norm_sq < libMesh::TOLERANCE)
95  return 0.0;
96 
97  const Real grad_norm = std::sqrt(grad_norm_sq);
98 
99  const Real nx = _grad_u[_qp](0);
100  const Real ny = _grad_u[_qp](1);
101  const Real nz = _grad_u[_qp](2);
102 
103  const Real px = _grad_phi[_j][_qp](0);
104  const Real py = _grad_phi[_j][_qp](1);
105  const Real pz = _grad_phi[_j][_qp](2);
106 
107  const Real s = _stress / grad_norm;
108  const Real ds = -_stress / (grad_norm * grad_norm_sq);
109  const Real dsx = ds * nx;
110  const Real dsy = ds * ny;
111  const Real dsz = ds * nz;
112 
113  const Real dus = ds * (nx * px + ny * py + pz * nz);
114 
115  const Real b = -3.0 * nx * px - 3.0 * ny * py - 3.0 * nz * pz;
116  const Real dudsx = ds * nx / grad_norm_sq * b + ds * px;
117  const Real dudsy = ds * ny / grad_norm_sq * b + ds * py;
118  const Real dudsz = ds * nz / grad_norm_sq * b + ds * pz;
119 
120  // d/du d/d(grad eta)_x
121  _ddS(0, 0, 0) = (2.0 * ny * py + 2.0 * nz * pz) * dsx + (ny * ny + nz * nz) * dudsx;
122  _ddS(0, 1, 0) = _ddS(0, 0, 1) =
123  -py * s - ny * dus - px * ny * dsx - nx * py * dsx - nx * ny * dudsx;
124  _ddS(0, 1, 1) = 2.0 * px * s + 2.0 * nx * dus + (2.0 * nx * px + 2.0 * nz * pz) * dsx +
125  (nx * nx + nz * nz) * dudsx;
126  _ddS(0, 2, 0) = _ddS(0, 0, 2) =
127  -pz * s - nz * dus - px * nz * dsx - nx * pz * dsx - nx * nz * dudsx;
128  _ddS(0, 2, 1) = _ddS(0, 1, 2) = -py * nz * dsx - ny * pz * dsx - ny * nz * dudsx;
129  _ddS(0, 2, 2) = 2.0 * px * s + 2.0 * nx * dus + (2.0 * nx * px + 2.0 * ny * py) * dsx +
130  (nx * nx + ny * ny) * dudsx;
131 
132  // d/du d/d(grad eta)_y
133  _ddS(1, 0, 0) = 2.0 * py * s + 2.0 * ny * dus + (2.0 * ny * py + 2.0 * nz * pz) * dsy +
134  (ny * ny + nz * nz) * dudsy;
135  _ddS(1, 1, 0) = _ddS(1, 0, 1) =
136  -px * s - nx * dus - px * ny * dsy - nx * py * dsy - nx * ny * dudsy;
137  _ddS(1, 1, 1) = (2.0 * nx * px + 2.0 * nz * pz) * dsy + (nx * nx + nz * nz) * dudsy;
138  _ddS(1, 2, 0) = _ddS(1, 0, 2) = -px * nz * dsy - nx * pz * dsy - nx * nz * dudsy;
139  _ddS(1, 2, 1) = _ddS(1, 1, 2) =
140  -pz * s - nz * dus - py * nz * dsy - ny * pz * dsy - ny * nz * dudsy;
141  _ddS(1, 2, 2) = 2.0 * py * s + 2.0 * ny * dus + (2.0 * nx * px + 2.0 * ny * py) * dsy +
142  (nx * nx + ny * ny) * dudsy;
143 
144  // d/du d/d(grad eta)_z
145  _ddS(2, 0, 0) = 2.0 * pz * s + 2.0 * nz * dus + (2.0 * ny * py + 2.0 * nz * pz) * dsz +
146  (ny * ny + nz * nz) * dudsz;
147  _ddS(2, 1, 0) = _ddS(2, 0, 1) = -px * ny * dsz - nx * py * dsz - nx * ny * dudsz;
148  _ddS(2, 1, 1) = 2.0 * pz * s + 2.0 * nz * dus + (2.0 * nx * px + 2.0 * nz * pz) * dsz +
149  (nx * nx + nz * nz) * dudsz;
150  _ddS(2, 2, 0) = _ddS(2, 0, 2) =
151  -px * s - nx * dus - px * nz * dsz - nx * pz * dsz - nx * nz * dudsz;
152  _ddS(2, 2, 1) = _ddS(2, 1, 2) =
153  -py * s - ny * dus - py * nz * dsz - ny * pz * dsz - ny * nz * dudsz;
154  _ddS(2, 2, 2) = (2.0 * nx * px + 2.0 * ny * py) * dsz + (nx * nx + ny * ny) * dudsz;
155 
156  return _L[_qp] * 0.5 * _ddS.doubleContraction(_strain[_qp]) * _grad_test[_i][_qp];
157 }
ACInterfaceStress(const InputParameters &parameters)
const MaterialProperty< Real > & _L
Mobility.
const VariableGradient & _grad_u
static InputParameters validParams()
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
registerMooseObject("PhaseFieldApp", ACInterfaceStress)
static constexpr Real TOLERANCE
const VariablePhiGradient & _grad_phi
virtual Real computeQpJacobian() override
ADRealEigenVector< T, D, asd > sqrt(const ADRealEigenVector< T, D, asd > &)
void addRequiredParam(const std::string &name, const std::string &doc_string)
static InputParameters validParams()
Compute the Allen-Cahn interface stress driving force contribution .
RankThreeTensor _dS
d sigma/d(grad eta), derivative of interface stress tensor with order parameter gradient ...
virtual Real computeQpResidual() override
unsigned int _i
const Real _stress
interface stress
unsigned int _j
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const VariableTestGradient & _grad_test
void addClassDescription(const std::string &doc_string)
VectorValue< T > doubleContraction(const RankTwoTensorTempl< T > &b) const
RankThreeTensor _ddS
derivative of _dS w.r.t. the finite element coefficients for the Jacobian calculation ...
void addRangeCheckedParam(const std::string &name, const T &value, const std::string &parsed_function, const std::string &doc_string)
const MaterialProperty< RankTwoTensor > & _strain
unsigned int _qp