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 "INSChorinPredictor.h"
11 : #include "MooseMesh.h"
12 :
13 : registerMooseObject("NavierStokesApp", INSChorinPredictor);
14 :
15 : InputParameters
16 82 : INSChorinPredictor::validParams()
17 : {
18 82 : InputParameters params = Kernel::validParams();
19 :
20 82 : params.addClassDescription("This class computes the 'Chorin' Predictor equation in "
21 : "fully-discrete (both time and space) form.");
22 : // Coupled variables
23 164 : params.addRequiredCoupledVar("u", "x-velocity");
24 164 : params.addCoupledVar("v", "y-velocity"); // only required in 2D and 3D
25 164 : params.addCoupledVar("w", "z-velocity"); // only required in 3D
26 :
27 : // Make star also be required, even though we might not use it?
28 164 : params.addRequiredCoupledVar("u_star", "star x-velocity");
29 164 : params.addCoupledVar("v_star", "star y-velocity"); // only required in 2D and 3D
30 164 : params.addCoupledVar("w_star", "star z-velocity"); // only required in 3D
31 :
32 : // Required parameters
33 164 : params.addRequiredRangeCheckedParam<unsigned>(
34 : "component",
35 : "component>=0 & component<=2",
36 : "0,1,2 depending on if we are solving the x,y,z component of the Predictor equation");
37 164 : MooseEnum predictor_type("OLD NEW STAR");
38 164 : params.addRequiredParam<MooseEnum>(
39 : "predictor_type",
40 : predictor_type,
41 : "One of: OLD, NEW, STAR. Indicates which velocity to use in the predictor.");
42 :
43 : // Optional parameters
44 164 : params.addParam<MaterialPropertyName>("mu_name", "mu", "The name of the dynamic viscosity");
45 164 : params.addParam<MaterialPropertyName>("rho_name", "rho", "The name of the density");
46 :
47 82 : return params;
48 82 : }
49 :
50 44 : INSChorinPredictor::INSChorinPredictor(const InputParameters & parameters)
51 : : Kernel(parameters),
52 :
53 : // Current velocities
54 44 : _u_vel(coupledValue("u")),
55 44 : _v_vel(_mesh.dimension() >= 2 ? coupledValue("v") : _zero),
56 44 : _w_vel(_mesh.dimension() == 3 ? coupledValue("w") : _zero),
57 :
58 : // Old velocities
59 44 : _u_vel_old(coupledValueOld("u")),
60 44 : _v_vel_old(_mesh.dimension() >= 2 ? coupledValueOld("v") : _zero),
61 44 : _w_vel_old(_mesh.dimension() == 3 ? coupledValueOld("w") : _zero),
62 :
63 : // Star velocities
64 44 : _u_vel_star(coupledValue("u_star")),
65 44 : _v_vel_star(_mesh.dimension() >= 2 ? coupledValue("v_star") : _zero),
66 44 : _w_vel_star(_mesh.dimension() == 3 ? coupledValue("w_star") : _zero),
67 :
68 : // Velocity Gradients
69 44 : _grad_u_vel(coupledGradient("u")),
70 44 : _grad_v_vel(_mesh.dimension() >= 2 ? coupledGradient("v") : _grad_zero),
71 44 : _grad_w_vel(_mesh.dimension() == 3 ? coupledGradient("w") : _grad_zero),
72 :
73 : // Old Velocity Gradients
74 44 : _grad_u_vel_old(coupledGradientOld("u")),
75 44 : _grad_v_vel_old(_mesh.dimension() >= 2 ? coupledGradientOld("v") : _grad_zero),
76 44 : _grad_w_vel_old(_mesh.dimension() == 3 ? coupledGradientOld("w") : _grad_zero),
77 :
78 : // Star Velocity Gradients
79 44 : _grad_u_vel_star(coupledGradient("u_star")),
80 44 : _grad_v_vel_star(_mesh.dimension() >= 2 ? coupledGradient("v_star") : _grad_zero),
81 44 : _grad_w_vel_star(_mesh.dimension() == 3 ? coupledGradient("w_star") : _grad_zero),
82 :
83 : // Variable numberings
84 44 : _u_vel_var_number(coupled("u")),
85 44 : _v_vel_var_number(_mesh.dimension() >= 2 ? coupled("v") : libMesh::invalid_uint),
86 44 : _w_vel_var_number(_mesh.dimension() == 3 ? coupled("w") : libMesh::invalid_uint),
87 :
88 : // Star velocity numberings
89 44 : _u_vel_star_var_number(coupled("u_star")),
90 44 : _v_vel_star_var_number(_mesh.dimension() >= 2 ? coupled("v_star") : libMesh::invalid_uint),
91 44 : _w_vel_star_var_number(_mesh.dimension() == 3 ? coupled("w_star") : libMesh::invalid_uint),
92 :
93 : // Required parameters
94 88 : _component(getParam<unsigned>("component")),
95 88 : _predictor_enum(getParam<MooseEnum>("predictor_type")),
96 :
97 : // Material properties
98 88 : _mu(getMaterialProperty<Real>("mu_name")),
99 132 : _rho(getMaterialProperty<Real>("rho_name"))
100 : {
101 44 : }
102 :
103 : Real
104 11776000 : INSChorinPredictor::computeQpResidual()
105 : {
106 : // Vector object for test function
107 : RealVectorValue test;
108 11776000 : test(_component) = _test[_i][_qp];
109 :
110 : // Tensor object for test function gradient
111 : RealTensorValue grad_test;
112 47104000 : for (unsigned k = 0; k < 3; ++k)
113 35328000 : grad_test(_component, k) = _grad_test[_i][_qp](k);
114 :
115 : // Decide what velocity vector, gradient to use:
116 : RealVectorValue U;
117 : RealTensorValue grad_U;
118 :
119 11776000 : switch (_predictor_enum)
120 : {
121 0 : case OLD:
122 : {
123 0 : U = RealVectorValue(_u_vel_old[_qp], _v_vel_old[_qp], _w_vel_old[_qp]);
124 0 : grad_U = RealTensorValue(_grad_u_vel_old[_qp], _grad_v_vel_old[_qp], _grad_w_vel_old[_qp]);
125 0 : break;
126 : }
127 11776000 : case NEW:
128 : {
129 11776000 : U = RealVectorValue(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
130 11776000 : grad_U = RealTensorValue(_grad_u_vel[_qp], _grad_v_vel[_qp], _grad_w_vel[_qp]);
131 11776000 : break;
132 : }
133 0 : case STAR:
134 : {
135 : // Note: Donea and Huerta's book says you are supposed to use "star" velocity to make Chorin
136 : // implicit, not U^{n+1}.
137 0 : U = RealVectorValue(_u_vel_star[_qp], _v_vel_star[_qp], _w_vel_star[_qp]);
138 0 : grad_U = RealTensorValue(_grad_u_vel_star[_qp], _grad_v_vel_star[_qp], _grad_w_vel_star[_qp]);
139 0 : break;
140 : }
141 0 : default:
142 0 : mooseError("Unrecognized Chorin predictor type requested.");
143 : }
144 :
145 : //
146 : // Compute the different parts
147 : //
148 :
149 : // Note: _u is the component'th entry of "u_star" in Chorin's method.
150 11776000 : RealVectorValue U_old(_u_vel_old[_qp], _v_vel_old[_qp], _w_vel_old[_qp]);
151 11776000 : Real symmetric_part = (_u[_qp] - U_old(_component)) * _test[_i][_qp];
152 :
153 : // Convective part. Remember to multiply by _dt!
154 11776000 : Real convective_part = _dt * (grad_U * U) * test;
155 :
156 : // Viscous part - we are using the Laplacian form here for simplicity.
157 : // Remember to multiply by _dt!
158 11776000 : Real viscous_part = _dt * (_mu[_qp] / _rho[_qp]) * grad_U.contract(grad_test);
159 :
160 11776000 : return symmetric_part + convective_part + viscous_part;
161 : }
162 :
163 : Real
164 47104000 : INSChorinPredictor::computeQpJacobian()
165 : {
166 : // The mass matrix part is always there.
167 47104000 : Real mass_part = _phi[_j][_qp] * _test[_i][_qp];
168 :
169 : // The on-diagonal Jacobian contribution depends on whether the predictor uses the
170 : // 'new' or 'star' velocity.
171 : Real other_part = 0.;
172 47104000 : switch (_predictor_enum)
173 : {
174 : case OLD:
175 : case NEW:
176 : break;
177 0 : case STAR:
178 : {
179 0 : RealVectorValue U_star(_u_vel_star[_qp], _v_vel_star[_qp], _w_vel_star[_qp]);
180 : Real convective_part =
181 0 : _dt * ((U_star * _grad_phi[_j][_qp]) + _phi[_j][_qp] * _grad_u[_qp](_component)) *
182 0 : _test[_i][_qp];
183 : Real viscous_part =
184 0 : _dt * ((_mu[_qp] / _rho[_qp]) * (_grad_phi[_j][_qp] * _grad_test[_i][_qp]));
185 0 : other_part = convective_part + viscous_part;
186 : break;
187 : }
188 0 : default:
189 0 : mooseError("Unrecognized Chorin predictor type requested.");
190 : }
191 :
192 47104000 : return mass_part + other_part;
193 : }
194 :
195 : Real
196 150732800 : INSChorinPredictor::computeQpOffDiagJacobian(unsigned jvar)
197 : {
198 150732800 : switch (_predictor_enum)
199 : {
200 : case OLD:
201 : {
202 : return 0.;
203 : }
204 :
205 150732800 : case NEW:
206 : {
207 150732800 : if ((jvar == _u_vel_var_number) || (jvar == _v_vel_var_number) || (jvar == _w_vel_var_number))
208 : {
209 : // Derivative of grad_U wrt the velocity component
210 : RealTensorValue dgrad_U;
211 :
212 : // Initialize to invalid value, then determine correct value.
213 : unsigned vel_index = 99;
214 :
215 : // Map jvar into the indices (0,1,2)
216 75366400 : if (jvar == _u_vel_var_number)
217 : vel_index = 0;
218 :
219 37683200 : else if (jvar == _v_vel_var_number)
220 : vel_index = 1;
221 :
222 0 : else if (jvar == _w_vel_var_number)
223 : vel_index = 2;
224 :
225 : // Fill in the vel_index'th row of dgrad_U with _grad_phi[_j][_qp]
226 301465600 : for (unsigned k = 0; k < 3; ++k)
227 226099200 : dgrad_U(vel_index, k) = _grad_phi[_j][_qp](k);
228 :
229 : // Vector object for test function
230 : RealVectorValue test;
231 75366400 : test(_component) = _test[_i][_qp];
232 :
233 : // Vector object for U
234 75366400 : RealVectorValue U(_u_vel[_qp], _v_vel[_qp], _w_vel[_qp]);
235 :
236 : // Tensor object for test function gradient
237 : RealTensorValue grad_test;
238 301465600 : for (unsigned k = 0; k < 3; ++k)
239 226099200 : grad_test(_component, k) = _grad_test[_i][_qp](k);
240 :
241 : // Compute the convective part
242 : RealVectorValue convective_jac =
243 75366400 : _phi[_j][_qp] * RealVectorValue(_grad_u_vel[_qp](vel_index),
244 75366400 : _grad_v_vel[_qp](vel_index),
245 75366400 : _grad_w_vel[_qp](vel_index));
246 :
247 : // Extra contribution in vel_index component
248 75366400 : convective_jac(vel_index) += U * _grad_phi[_j][_qp];
249 :
250 : // Be sure to scale by _dt!
251 75366400 : Real convective_part = _dt * (convective_jac * test);
252 :
253 : // Compute the viscous part, be sure to scale by _dt. Note: the contracted
254 : // value should be zero unless vel_index and _component match.
255 75366400 : Real viscous_part = _dt * (_mu[_qp] / _rho[_qp]) * dgrad_U.contract(grad_test);
256 :
257 : // Return the result
258 75366400 : return convective_part + viscous_part;
259 : }
260 : else
261 : return 0;
262 : }
263 :
264 0 : case STAR:
265 : {
266 0 : if (jvar == _u_vel_star_var_number)
267 : {
268 0 : return _dt * _phi[_j][_qp] * _grad_u[_qp](0) * _test[_i][_qp];
269 : }
270 :
271 0 : else if (jvar == _v_vel_star_var_number)
272 : {
273 0 : return _dt * _phi[_j][_qp] * _grad_u[_qp](1) * _test[_i][_qp];
274 : }
275 :
276 0 : else if (jvar == _w_vel_star_var_number)
277 : {
278 0 : return _dt * _phi[_j][_qp] * _grad_u[_qp](2) * _test[_i][_qp];
279 : }
280 :
281 : else
282 : return 0;
283 : }
284 :
285 0 : default:
286 0 : mooseError("Unrecognized Chorin predictor type requested.");
287 : }
288 : }
|
