Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
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 :
16 : InputParameters
17 23 : ACInterfaceStress::validParams()
18 : {
19 23 : InputParameters params = Kernel::validParams();
20 23 : params.addClassDescription("Interface stress driving force Allen-Cahn Kernel");
21 46 : params.addParam<MaterialPropertyName>("mob_name", "L", "The mobility used with the kernel");
22 46 : params.addParam<std::string>("base_name", "Material property base name");
23 46 : params.addRequiredParam<Real>("stress", "Planar stress");
24 69 : params.addRangeCheckedParam<Real>("op_range",
25 46 : 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 23 : return params;
31 0 : }
32 :
33 12 : ACInterfaceStress::ACInterfaceStress(const InputParameters & parameters)
34 : : Kernel(parameters),
35 12 : _L(getMaterialProperty<Real>("mob_name")),
36 24 : _base_name(isParamValid("base_name") ? getParam<std::string>("base_name") + "_" : ""),
37 24 : _strain(getMaterialPropertyByName<RankTwoTensor>(_base_name + "elastic_strain")),
38 48 : _stress(getParam<Real>("stress") / getParam<Real>("op_range"))
39 : {
40 12 : }
41 :
42 : Real
43 62914560 : ACInterfaceStress::computeQpResidual()
44 : {
45 : // no interface, return zero stress
46 62914560 : const Real grad_norm_sq = _grad_u[_qp].norm_sq();
47 62914560 : if (grad_norm_sq < libMesh::TOLERANCE)
48 : return 0.0;
49 :
50 38491960 : 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 38491960 : const Real s = _stress / grad_norm;
57 38491960 : const Real ds = -_stress / (grad_norm * grad_norm_sq);
58 38491960 : const Real dsx = ds * nx;
59 38491960 : const Real dsy = ds * ny;
60 38491960 : const Real dsz = ds * nz;
61 :
62 : // d/d(grad eta)_x
63 38491960 : _dS(0, 0, 0) = (ny * ny + nz * nz) * dsx; // (ny * ny + nz * nz) * s;
64 38491960 : _dS(0, 1, 0) = _dS(0, 0, 1) = -ny * s - nx * ny * dsx; // -nx * ny * s;
65 38491960 : _dS(0, 1, 1) = 2.0 * nx * s + (nx * nx + nz * nz) * dsx; // (nx * nx + nz * nz) * s;
66 38491960 : _dS(0, 2, 0) = _dS(0, 0, 2) = -nz * s - nx * nz * dsx; // -nx * nz * s;
67 38491960 : _dS(0, 2, 1) = _dS(0, 1, 2) = -ny * nz * dsx; // -ny * nz * s;
68 38491960 : _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 38491960 : _dS(1, 0, 0) = 2.0 * ny * s + (ny * ny + nz * nz) * dsy; // (ny * ny + nz * nz) * s;
72 38491960 : _dS(1, 1, 0) = _dS(1, 0, 1) = -nx * s - nx * ny * dsy; // -nx * ny * s;
73 38491960 : _dS(1, 1, 1) = (nx * nx + nz * nz) * dsy; // (nx * nx + nz * nz) * s;
74 38491960 : _dS(1, 2, 0) = _dS(1, 0, 2) = -nx * nz * dsy; // -nx * nz * s;
75 38491960 : _dS(1, 2, 1) = _dS(1, 1, 2) = -nz * s - ny * nz * dsy; // -ny * nz * s;
76 38491960 : _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 38491960 : _dS(2, 0, 0) = 2.0 * nz * s + (ny * ny + nz * nz) * dsz; // (ny * ny + nz * nz) * s;
80 38491960 : _dS(2, 1, 0) = _dS(2, 0, 1) = -nx * ny * dsz; // -nx * ny * s;
81 38491960 : _dS(2, 1, 1) = 2.0 * nz * s + (nx * nx + nz * nz) * dsz; // (nx * nx + nz * nz) * s;
82 38491960 : _dS(2, 2, 0) = _dS(2, 0, 2) = -nx * s - nx * nz * dsz; // -nx * nz * s;
83 38491960 : _dS(2, 2, 1) = _dS(2, 1, 2) = -ny * s - ny * nz * dsz; // -ny * nz * s;
84 38491960 : _dS(2, 2, 2) = (nx * nx + ny * ny) * dsz; // (nx * nx + ny * ny) * s;
85 :
86 38491960 : return _L[_qp] * 0.5 * _dS.doubleContraction(_strain[_qp]) * _grad_test[_i][_qp];
87 : }
88 :
89 : Real
90 63963136 : ACInterfaceStress::computeQpJacobian()
91 : {
92 : // no interface, return zero stress
93 63963136 : const Real grad_norm_sq = _grad_u[_qp].norm_sq();
94 63963136 : if (grad_norm_sq < libMesh::TOLERANCE)
95 : return 0.0;
96 :
97 39116608 : 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 39116608 : const Real px = _grad_phi[_j][_qp](0);
104 39116608 : const Real py = _grad_phi[_j][_qp](1);
105 39116608 : const Real pz = _grad_phi[_j][_qp](2);
106 :
107 39116608 : const Real s = _stress / grad_norm;
108 39116608 : const Real ds = -_stress / (grad_norm * grad_norm_sq);
109 39116608 : const Real dsx = ds * nx;
110 39116608 : const Real dsy = ds * ny;
111 39116608 : const Real dsz = ds * nz;
112 :
113 39116608 : const Real dus = ds * (nx * px + ny * py + pz * nz);
114 :
115 39116608 : const Real b = -3.0 * nx * px - 3.0 * ny * py - 3.0 * nz * pz;
116 39116608 : const Real dudsx = ds * nx / grad_norm_sq * b + ds * px;
117 39116608 : const Real dudsy = ds * ny / grad_norm_sq * b + ds * py;
118 39116608 : const Real dudsz = ds * nz / grad_norm_sq * b + ds * pz;
119 :
120 : // d/du d/d(grad eta)_x
121 39116608 : _ddS(0, 0, 0) = (2.0 * ny * py + 2.0 * nz * pz) * dsx + (ny * ny + nz * nz) * dudsx;
122 39116608 : _ddS(0, 1, 0) = _ddS(0, 0, 1) =
123 39116608 : -py * s - ny * dus - px * ny * dsx - nx * py * dsx - nx * ny * dudsx;
124 39116608 : _ddS(0, 1, 1) = 2.0 * px * s + 2.0 * nx * dus + (2.0 * nx * px + 2.0 * nz * pz) * dsx +
125 39116608 : (nx * nx + nz * nz) * dudsx;
126 39116608 : _ddS(0, 2, 0) = _ddS(0, 0, 2) =
127 39116608 : -pz * s - nz * dus - px * nz * dsx - nx * pz * dsx - nx * nz * dudsx;
128 39116608 : _ddS(0, 2, 1) = _ddS(0, 1, 2) = -py * nz * dsx - ny * pz * dsx - ny * nz * dudsx;
129 39116608 : _ddS(0, 2, 2) = 2.0 * px * s + 2.0 * nx * dus + (2.0 * nx * px + 2.0 * ny * py) * dsx +
130 39116608 : (nx * nx + ny * ny) * dudsx;
131 :
132 : // d/du d/d(grad eta)_y
133 39116608 : _ddS(1, 0, 0) = 2.0 * py * s + 2.0 * ny * dus + (2.0 * ny * py + 2.0 * nz * pz) * dsy +
134 39116608 : (ny * ny + nz * nz) * dudsy;
135 39116608 : _ddS(1, 1, 0) = _ddS(1, 0, 1) =
136 39116608 : -px * s - nx * dus - px * ny * dsy - nx * py * dsy - nx * ny * dudsy;
137 39116608 : _ddS(1, 1, 1) = (2.0 * nx * px + 2.0 * nz * pz) * dsy + (nx * nx + nz * nz) * dudsy;
138 39116608 : _ddS(1, 2, 0) = _ddS(1, 0, 2) = -px * nz * dsy - nx * pz * dsy - nx * nz * dudsy;
139 39116608 : _ddS(1, 2, 1) = _ddS(1, 1, 2) =
140 39116608 : -pz * s - nz * dus - py * nz * dsy - ny * pz * dsy - ny * nz * dudsy;
141 39116608 : _ddS(1, 2, 2) = 2.0 * py * s + 2.0 * ny * dus + (2.0 * nx * px + 2.0 * ny * py) * dsy +
142 39116608 : (nx * nx + ny * ny) * dudsy;
143 :
144 : // d/du d/d(grad eta)_z
145 39116608 : _ddS(2, 0, 0) = 2.0 * pz * s + 2.0 * nz * dus + (2.0 * ny * py + 2.0 * nz * pz) * dsz +
146 39116608 : (ny * ny + nz * nz) * dudsz;
147 39116608 : _ddS(2, 1, 0) = _ddS(2, 0, 1) = -px * ny * dsz - nx * py * dsz - nx * ny * dudsz;
148 39116608 : _ddS(2, 1, 1) = 2.0 * pz * s + 2.0 * nz * dus + (2.0 * nx * px + 2.0 * nz * pz) * dsz +
149 39116608 : (nx * nx + nz * nz) * dudsz;
150 39116608 : _ddS(2, 2, 0) = _ddS(2, 0, 2) =
151 39116608 : -px * s - nx * dus - px * nz * dsz - nx * pz * dsz - nx * nz * dudsz;
152 39116608 : _ddS(2, 2, 1) = _ddS(2, 1, 2) =
153 39116608 : -py * s - ny * dus - py * nz * dsz - ny * pz * dsz - ny * nz * dudsz;
154 39116608 : _ddS(2, 2, 2) = (2.0 * nx * px + 2.0 * ny * py) * dsz + (nx * nx + ny * ny) * dudsz;
155 :
156 39116608 : return _L[_qp] * 0.5 * _ddS.doubleContraction(_strain[_qp]) * _grad_test[_i][_qp];
157 : }
|
