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 : // Navier-Stokes includes
11 : #include "NavierStokesMaterial.h"
12 : #include "NS.h"
13 :
14 : // FluidProperties includes
15 : #include "IdealGasFluidProperties.h"
16 :
17 : // MOOSE includes
18 : #include "Assembly.h"
19 : #include "MooseMesh.h"
20 :
21 : #include "libmesh/quadrature.h"
22 :
23 : namespace nms = NS;
24 :
25 : InputParameters
26 85 : NavierStokesMaterial::validParams()
27 : {
28 85 : InputParameters params = Material::validParams();
29 :
30 85 : params.addClassDescription(
31 : "This is the base class of all materials should use if you are trying to "
32 : "use the Navier-Stokes Kernels.");
33 85 : params.addRequiredCoupledVar(nms::velocity_x, "x-velocity");
34 85 : params.addCoupledVar(nms::velocity_y, "y-velocity"); // only required in >= 2D
35 85 : params.addCoupledVar(nms::velocity_z, "z-velocity"); // only required in 3D
36 :
37 85 : params.addRequiredCoupledVar(nms::temperature, "temperature");
38 85 : params.addRequiredCoupledVar(nms::specific_total_enthalpy, "specific total enthalpy");
39 :
40 85 : params.addRequiredCoupledVar(nms::density, "density");
41 85 : params.addRequiredCoupledVar(nms::momentum_x, "x-momentum");
42 85 : params.addCoupledVar(nms::momentum_y, "y-momentum"); // only required in >= 2D
43 85 : params.addCoupledVar(nms::momentum_z, "z-momentum"); // only required in 3D
44 85 : params.addRequiredCoupledVar(nms::total_energy_density, "energy");
45 170 : params.addRequiredParam<UserObjectName>("fluid_properties",
46 : "The name of the user object for fluid properties");
47 85 : return params;
48 0 : }
49 :
50 66 : NavierStokesMaterial::NavierStokesMaterial(const InputParameters & parameters)
51 : : Material(parameters),
52 132 : _mesh_dimension(_mesh.dimension()),
53 66 : _grad_u(coupledGradient(nms::velocity_x)),
54 66 : _grad_v(_mesh_dimension >= 2 ? coupledGradient(nms::velocity_y) : _grad_zero),
55 66 : _grad_w(_mesh_dimension == 3 ? coupledGradient(nms::velocity_z) : _grad_zero),
56 :
57 66 : _viscous_stress_tensor(declareProperty<RealTensorValue>("viscous_stress_tensor")),
58 66 : _thermal_conductivity(declareProperty<Real>("thermal_conductivity")),
59 :
60 : // Declared here but _not_ calculated here
61 : // (See e.g. derived class, bighorn/include/materials/FluidTC1.h)
62 66 : _dynamic_viscosity(declareProperty<Real>("dynamic_viscosity")),
63 :
64 : // The momentum components of the inviscid flux Jacobians.
65 66 : _calA(declareProperty<std::vector<RealTensorValue>>("calA")),
66 :
67 : // "Velocity column" matrices
68 66 : _calC(declareProperty<std::vector<RealTensorValue>>("calC")),
69 :
70 : // Energy equation inviscid flux matrices, "cal E_{kl}" in the notes.
71 66 : _calE(declareProperty<std::vector<std::vector<RealTensorValue>>>("calE")),
72 66 : _vel_grads({&_grad_u, &_grad_v, &_grad_w}),
73 :
74 : // Coupled solution values needed for computing SUPG stabilization terms
75 66 : _u_vel(coupledValue(nms::velocity_x)),
76 66 : _v_vel(_mesh.dimension() >= 2 ? coupledValue(nms::velocity_y) : _zero),
77 66 : _w_vel(_mesh.dimension() == 3 ? coupledValue(nms::velocity_z) : _zero),
78 :
79 66 : _temperature(coupledValue(nms::temperature)),
80 66 : _specific_total_enthalpy(coupledValue(nms::specific_total_enthalpy)),
81 :
82 : // Coupled solution values
83 66 : _rho(coupledValue(nms::density)),
84 66 : _rho_u(coupledValue(nms::momentum_x)),
85 66 : _rho_v(_mesh.dimension() >= 2 ? coupledValue(nms::momentum_y) : _zero),
86 66 : _rho_w(_mesh.dimension() == 3 ? coupledValue(nms::momentum_z) : _zero),
87 66 : _rho_et(coupledValue(nms::total_energy_density)),
88 :
89 : // Time derivative values
90 66 : _drho_dt(coupledDot(nms::density)),
91 66 : _drhou_dt(coupledDot(nms::momentum_x)),
92 66 : _drhov_dt(_mesh.dimension() >= 2 ? coupledDot(nms::momentum_y) : _zero),
93 66 : _drhow_dt(_mesh.dimension() == 3 ? coupledDot(nms::momentum_z) : _zero),
94 66 : _drhoE_dt(coupledDot(nms::total_energy_density)),
95 :
96 : // Gradients
97 66 : _grad_rho(coupledGradient(nms::density)),
98 66 : _grad_rho_u(coupledGradient(nms::momentum_x)),
99 66 : _grad_rho_v(_mesh.dimension() >= 2 ? coupledGradient(nms::momentum_y) : _grad_zero),
100 66 : _grad_rho_w(_mesh.dimension() == 3 ? coupledGradient(nms::momentum_z) : _grad_zero),
101 66 : _grad_rho_et(coupledGradient(nms::total_energy_density)),
102 :
103 : // Material properties for stabilization
104 66 : _hsupg(declareProperty<Real>("hsupg")),
105 66 : _tauc(declareProperty<Real>("tauc")),
106 66 : _taum(declareProperty<Real>("taum")),
107 66 : _taue(declareProperty<Real>("taue")),
108 66 : _strong_residuals(declareProperty<std::vector<Real>>("strong_residuals")),
109 132 : _fp(getUserObject<IdealGasFluidProperties>("fluid_properties"))
110 : {
111 66 : }
112 :
113 : /**
114 : * Must be called _after_ the child class computes dynamic_viscosity.
115 : */
116 : void
117 526848 : NavierStokesMaterial::computeProperties()
118 : {
119 2370816 : for (unsigned int qp = 0; qp < _qrule->n_points(); ++qp)
120 : {
121 : /******* Viscous Stress Tensor *******/
122 : // Technically... this _is_ the transpose (since we are loading these by rows)
123 : // But it doesn't matter....
124 1843968 : RealTensorValue grad_outer_u(_grad_u[qp], _grad_v[qp], _grad_w[qp]);
125 :
126 1843968 : grad_outer_u += grad_outer_u.transpose();
127 :
128 : Real div_vel = 0.0;
129 7375872 : for (unsigned int i = 0; i < 3; ++i)
130 5531904 : div_vel += (*_vel_grads[i])[qp](i);
131 :
132 : // Add diagonal terms
133 7375872 : for (unsigned int i = 0; i < 3; ++i)
134 5531904 : grad_outer_u(i, i) -= 2.0 / 3.0 * div_vel;
135 :
136 1843968 : grad_outer_u *= _dynamic_viscosity[qp];
137 :
138 1843968 : _viscous_stress_tensor[qp] = grad_outer_u;
139 :
140 : // Tabulated values of thermal conductivity vs. Temperature for air (k increases slightly with
141 : // T):
142 : // T (K) k (W/m-K)
143 : // 273 0.0243
144 : // 373 0.0314
145 : // 473 0.0386
146 : // 573 0.0454
147 : // 673 0.0515
148 :
149 : // Pr = (mu * cp) / k ==> k = (mu * cp) / Pr = (mu * gamma * cv) / Pr.
150 : // TODO: We are using a fixed value of the Prandtl number which is
151 : // valid for air, it may or may not depend on temperature? Since
152 : // this is a property of the fluid, it could possibly be moved to
153 : // the FluidProperties module...
154 : const Real Pr = 0.71;
155 1843968 : _thermal_conductivity[qp] = (_dynamic_viscosity[qp] * _fp.cp()) / Pr;
156 :
157 : // Compute stabilization parameters:
158 :
159 : // .) Compute SUPG element length scale.
160 1843968 : computeHSUPG(qp);
161 : // Moose::out << "_hsupg[" << qp << "]=" << _hsupg[qp] << std::endl;
162 :
163 : // .) Compute SUPG parameter values. (Must call this after computeHSUPG())
164 1843968 : computeTau(qp);
165 : // Moose::out << "_tauc[" << qp << "]=" << _tauc[qp] << ", ";
166 : // Moose::out << "_taum[" << qp << "]=" << _taum[qp] << ", ";
167 : // Moose::out << "_taue[" << qp << "]=" << _taue[qp] << std::endl;
168 :
169 : // .) Compute strong residual values.
170 1843968 : computeStrongResiduals(qp);
171 : // Moose::out << "_strong_residuals[" << qp << "]=";
172 : // for (unsigned i=0; i<_strong_residuals[qp].size(); ++i)
173 : // Moose::out << _strong_residuals[qp][i] << " ";
174 : // Moose::out << std::endl;
175 : }
176 526848 : }
177 :
178 : void
179 1843968 : NavierStokesMaterial::computeHSUPG(unsigned int qp)
180 : {
181 : // // Grab reference to linear Lagrange finite element object pointer,
182 : // // currently this is always a linear Lagrange element, so this might need to
183 : // // be generalized if we start working with higher-order elements...
184 : // FEBase*& fe(_assembly.getFE(FEType(), _current_elem->dim()));
185 : //
186 : // // Grab references to FE object's mapping data from the _subproblem's FE object
187 : // const std::vector<Real> & dxidx(fe->get_dxidx());
188 : // const std::vector<Real> & dxidy(fe->get_dxidy());
189 : // const std::vector<Real> & dxidz(fe->get_dxidz());
190 : // const std::vector<Real> & detadx(fe->get_detadx());
191 : // const std::vector<Real> & detady(fe->get_detady());
192 : // const std::vector<Real> & detadz(fe->get_detadz());
193 : // const std::vector<Real> & dzetadx(fe->get_dzetadx()); // Empty in 2D
194 : // const std::vector<Real> & dzetady(fe->get_dzetady()); // Empty in 2D
195 : // const std::vector<Real> & dzetadz(fe->get_dzetadz()); // Empty in 2D
196 : //
197 : // // Bounds checking on element data
198 : // mooseAssert(qp < dxidx.size(), "Insufficient data in dxidx array!");
199 : // mooseAssert(qp < dxidy.size(), "Insufficient data in dxidy array!");
200 : // mooseAssert(qp < dxidz.size(), "Insufficient data in dxidz array!");
201 : //
202 : // mooseAssert(qp < detadx.size(), "Insufficient data in detadx array!");
203 : // mooseAssert(qp < detady.size(), "Insufficient data in detady array!");
204 : // mooseAssert(qp < detadz.size(), "Insufficient data in detadz array!");
205 : //
206 : // if (_mesh_dimension == 3)
207 : // {
208 : // mooseAssert(qp < dzetadx.size(), "Insufficient data in dzetadx array!");
209 : // mooseAssert(qp < dzetady.size(), "Insufficient data in dzetady array!");
210 : // mooseAssert(qp < dzetadz.size(), "Insufficient data in dzetadz array!");
211 : // }
212 : //
213 : // // The velocity vector at this quadrature point.
214 : // RealVectorValue U(_u_vel[qp],_v_vel[qp],_w_vel[qp]);
215 : //
216 : // // Pull out element inverse map values at the current qp into a little dense matrix
217 : // Real dxi_dx[3][3] = {{0.,0.,0.}, {0.,0.,0.}, {0.,0.,0.}};
218 : //
219 : // dxi_dx[0][0] = dxidx[qp]; dxi_dx[0][1] = dxidy[qp];
220 : // dxi_dx[1][0] = detadx[qp]; dxi_dx[1][1] = detady[qp];
221 : //
222 : // // OK to access third entries on 2D elements if LIBMESH_DIM==3, though they
223 : // // may be zero...
224 : // if (LIBMESH_DIM == 3)
225 : // {
226 : // /**/ /**/ dxi_dx[0][2] = dxidz[qp];
227 : // /**/ /**/ dxi_dx[1][2] = detadz[qp];
228 : // }
229 : //
230 : // // The last row of entries available only for 3D elements.
231 : // if (_mesh_dimension == 3)
232 : // {
233 : // dxi_dx[2][0] = dzetadx[qp]; dxi_dx[2][1] = dzetady[qp]; dxi_dx[2][2] = dzetadz[qp];
234 : // }
235 : //
236 : // // Construct the g_ij = d(xi_k)/d(x_j) * d(xi_k)/d(x_i) matrix
237 : // // from Ben and Bova's paper by summing over k...
238 : // Real g[3][3] = {{0.,0.,0.}, {0.,0.,0.}, {0.,0.,0.}};
239 : // for (unsigned int i = 0; i < 3; ++i)
240 : // for (unsigned int j = 0; j < 3; ++j)
241 : // for (unsigned int k = 0; k < 3; ++k)
242 : // g[i][j] += dxi_dx[k][j] * dxi_dx[k][i];
243 : //
244 : // // Compute the denominator of the h_supg term: U * (g) * U
245 : // Real denom = 0.;
246 : // for (unsigned int i = 0; i < 3; ++i)
247 : // for (unsigned int j = 0; j < 3; ++j)
248 : // denom += U(j) * g[i][j] * U(i);
249 : //
250 : // // Compute h_supg. Some notes:
251 : // // .) The 2 coefficient in this term should be a 1 if we are using tets/triangles.
252 : // // .) The denominator will be identically zero only if the velocity
253 : // // is identically zero, in which case we can't divide by it.
254 : // if (denom != 0.0)
255 : // _hsupg[qp] = 2.* sqrt( U.norm_sq() / denom );
256 : // else
257 : // _hsupg[qp] = 0.;
258 :
259 : // Simple (and fast) implementation: Just use hmin for the element!
260 1843968 : _hsupg[qp] = _current_elem->hmin();
261 1843968 : }
262 :
263 : void
264 1843968 : NavierStokesMaterial::computeTau(unsigned int qp)
265 : {
266 : Real velmag =
267 1843968 : std::sqrt(_u_vel[qp] * _u_vel[qp] + _v_vel[qp] * _v_vel[qp] + _w_vel[qp] * _w_vel[qp]);
268 :
269 : // Moose::out << "velmag=" << velmag << std::endl;
270 :
271 : // Make sure temperature >= 0 before trying to take sqrt
272 : // if (_temperature[qp] < 0.)
273 : // {
274 : // Moose::err << "Negative temperature "
275 : // << _temperature[qp]
276 : // << " found at quadrature point "
277 : // << qp
278 : // << ", element "
279 : // << _current_elem->id()
280 : // << std::endl;
281 : // mooseError("Can't continue, would be nice to throw an exception here?");
282 : // }
283 :
284 : // The speed of sound for an ideal gas, sqrt(gamma * R * T). Not needed unless
285 : // we want to use a form of Tau that requires it.
286 : // Real soundspeed = _fp.c_from_v_e(_specific_volume[_qp], _internal_energy[_qp]);
287 :
288 : // If velmag == 0, then _hsupg should be zero as well. Then tau
289 : // will have only the time-derivative contribution (or zero, if we
290 : // are not including dt terms in our taus!) Note that using the
291 : // time derivative contribution in this way assumes we are solving
292 : // unsteady, and guarantees *some* stabilization is added even when
293 : // u -> 0 in certain regions of the flow.
294 1843968 : if (velmag == 0.)
295 : {
296 : // 1.) Tau without dt terms
297 : // _tauc[qp] = 0.;
298 : // _taum[qp] = 0.;
299 : // _taue[qp] = 0.;
300 :
301 : // 2.) Tau *with* dt terms
302 0 : _tauc[qp] = _taum[qp] = _taue[qp] = 0.5 * _dt;
303 : }
304 : else
305 : {
306 : // The element length parameter, squared
307 1843968 : Real h2 = _hsupg[qp] * _hsupg[qp];
308 :
309 : // The viscosity-based term
310 1843968 : Real visc_term = _dynamic_viscosity[qp] / _rho[qp] / h2;
311 :
312 : // The thermal conductivity-based term, cp = gamma * cv
313 1843968 : Real k_term = _thermal_conductivity[qp] / _rho[qp] / _fp.cp() / h2;
314 :
315 : // 1a.) Standard compressible flow tau. Does not account for low Mach number
316 : // limit.
317 : // _tauc[qp] = _hsupg[qp] / (velmag + soundspeed);
318 :
319 : // 1b.) Inspired by Hauke, the sum of the compressible and incompressible tauc.
320 : // _tauc[qp] =
321 : // _hsupg[qp] / (velmag + soundspeed) +
322 : // _hsupg[qp] / (velmag);
323 :
324 : // 1c.) From Wong 2001. This tau is O(M^2) for small M. At small M,
325 : // tauc dominates the inverse square sums and basically makes
326 : // taum=taue=tauc. However, all my flows occur at low Mach numbers,
327 : // so there would basically never be any stabilization...
328 : // _tauc[qp] = (_hsupg[qp] * velmag) / (velmag*velmag + soundspeed*soundspeed);
329 :
330 : // For use with option "1",
331 : // (tau_c)^{-2}
332 : // Real taucm2 = 1./_tauc[qp]/_tauc[qp];
333 : // _taum[qp] = 1. / std::sqrt(taucm2 + visc_term*visc_term);
334 : // _taue[qp] = 1. / std::sqrt(taucm2 + k_term*k_term);
335 :
336 : // 2.) Tau with timestep dependence (guarantees stabilization even
337 : // in zero-velocity limit) incorporated via the "r-switch" method,
338 : // with r=2.
339 1843968 : Real sqrt_term = 4. / _dt / _dt + velmag * velmag / h2;
340 :
341 : // For use with option "2", i.e. the option that uses dt in the definition of tau
342 1843968 : _tauc[qp] = 1. / std::sqrt(sqrt_term);
343 1843968 : _taum[qp] = 1. / std::sqrt(sqrt_term + visc_term * visc_term);
344 1843968 : _taue[qp] = 1. / std::sqrt(sqrt_term + k_term * k_term);
345 : }
346 :
347 : // Debugging
348 : // Moose::out << "_tauc[" << qp << "]=" << _tauc[qp] << std::endl;
349 : // Moose::out << "_hsupg[" << qp << "]=" << _hsupg[qp] << std::endl;
350 : // Moose::out << "velmag[" << qp << "]=" << velmag << std::endl;
351 1843968 : }
352 :
353 : void
354 1843968 : NavierStokesMaterial::computeStrongResiduals(unsigned int qp)
355 : {
356 : // Create storage at this qp for the strong residuals of all the equations.
357 : // In 2D, the value for the z-velocity equation will just be zero.
358 1843968 : _strong_residuals[qp].resize(5);
359 :
360 : // The timestep is stored in the Problem object, which can be accessed through
361 : // the parent pointer of the SubProblem. Don't need this if we are not
362 : // approximating time derivatives ourselves.
363 : // Real dt = _subproblem.parent()->dt();
364 : // Moose::out << "dt=" << dt << std::endl;
365 :
366 : // Vector object for the velocity
367 1843968 : RealVectorValue vel(_u_vel[qp], _v_vel[qp], _w_vel[qp]);
368 :
369 : // A VectorValue object containing all zeros. Makes it easier to
370 : // construct type tensor objects
371 : RealVectorValue zero(0., 0., 0.);
372 :
373 : // Velocity vector magnitude squared
374 : Real velmag2 = vel.norm_sq();
375 :
376 : // Debugging: How large are the time derivative parts of the strong residuals?
377 : // Moose::out << "drho_dt=" << _drho_dt
378 : // << ", drhou_dt=" << _drhou_dt
379 : // << ", drhov_dt=" << _drhov_dt
380 : // << ", drhow_dt=" << _drhow_dt
381 : // << ", drhoE_dt=" << _drhoE_dt
382 : // << std::endl;
383 :
384 : // Momentum divergence
385 1843968 : Real divU = _grad_rho_u[qp](0) + _grad_rho_v[qp](1) + _grad_rho_w[qp](2);
386 :
387 : // Enough space to hold three space dimensions of velocity components at each qp,
388 : // regardless of what dimension we are actually running in.
389 1843968 : _calC[qp].resize(3);
390 :
391 : // Explicitly zero the calC
392 7375872 : for (unsigned int i = 0; i < 3; ++i)
393 5531904 : _calC[qp][i].zero();
394 :
395 : // x-column matrix
396 1843968 : _calC[qp][0](0, 0) = _u_vel[qp];
397 1843968 : _calC[qp][0](1, 0) = _v_vel[qp];
398 1843968 : _calC[qp][0](2, 0) = _w_vel[qp];
399 :
400 : // y-column matrix
401 1843968 : _calC[qp][1](0, 1) = _u_vel[qp];
402 1843968 : _calC[qp][1](1, 1) = _v_vel[qp];
403 1843968 : _calC[qp][1](2, 1) = _w_vel[qp];
404 :
405 : // z-column matrix (this assumes LIBMESH_DIM==3!)
406 1843968 : _calC[qp][2](0, 2) = _u_vel[qp];
407 1843968 : _calC[qp][2](1, 2) = _v_vel[qp];
408 1843968 : _calC[qp][2](2, 2) = _w_vel[qp];
409 :
410 : // The matrix S can be computed from any of the calC via calC_1*calC_1^T
411 1843968 : RealTensorValue calS = _calC[qp][0] * _calC[qp][0].transpose();
412 :
413 : // Enough space to hold five (=n_sd + 2) 3*3 calA matrices at this qp, regarless of dimension
414 1843968 : _calA[qp].resize(5);
415 :
416 : // 0.) _calA_0 = diag( (gam - 1)/2*|u|^2 ) - S
417 1843968 : _calA[qp][0].zero(); // zero this calA entry
418 1843968 : _calA[qp][0](0, 0) = _calA[qp][0](1, 1) = _calA[qp][0](2, 2) =
419 1843968 : 0.5 * (_fp.gamma() - 1.0) * velmag2; // set diag. entries
420 : _calA[qp][0] -= calS;
421 :
422 7375872 : for (unsigned int m = 1; m <= 3; ++m)
423 : {
424 : // Use m_local when indexing into matrices and vectors
425 5531904 : unsigned int m_local = m - 1;
426 :
427 : // For m=1,2,3, calA_m = C_m + C_m^T + diag( (1.-gam)*u_m )
428 5531904 : _calA[qp][m].zero(); // zero this calA entry
429 5531904 : _calA[qp][m](0, 0) = _calA[qp][m](1, 1) = _calA[qp][m](2, 2) =
430 5531904 : (1. - _fp.gamma()) * vel(m_local); // set diag. entries
431 5531904 : _calA[qp][m] += _calC[qp][m_local]; // Note: use m_local for indexing into _calC!
432 5531904 : _calA[qp][m] += _calC[qp][m_local].transpose(); // Note: use m_local for indexing into _calC!
433 : }
434 :
435 : // 4.) calA_4 = diag(gam - 1)
436 1843968 : _calA[qp][4].zero(); // zero this calA entry
437 1843968 : _calA[qp][4](0, 0) = _calA[qp][4](1, 1) = _calA[qp][4](2, 2) = (_fp.gamma() - 1.0);
438 :
439 : // Enough space to hold the 3*5 "cal E" matrices which comprise the inviscid flux term
440 : // of the energy equation. See notes for additional details
441 1843968 : _calE[qp].resize(3); // Three rows, 5 entries in each row
442 :
443 7375872 : for (unsigned int k = 0; k < 3; ++k)
444 : {
445 : // Make enough room to store all 5 E matrices for this k
446 5531904 : _calE[qp][k].resize(5);
447 :
448 : // Store and reuse the velocity column transpose matrix for the
449 : // current value of k.
450 5531904 : RealTensorValue Ck_T = _calC[qp][k].transpose();
451 :
452 : // E_{k0} (density gradient term)
453 5531904 : _calE[qp][k][0].zero();
454 5531904 : _calE[qp][k][0] = (0.5 * (_fp.gamma() - 1.0) * velmag2 - _specific_total_enthalpy[qp]) * Ck_T;
455 :
456 22127616 : for (unsigned int m = 1; m <= 3; ++m)
457 : {
458 : // Use m_local when indexing into matrices and vectors
459 16595712 : unsigned int m_local = m - 1;
460 :
461 : // E_{km} (momentum gradient terms)
462 16595712 : _calE[qp][k][m].zero();
463 16595712 : _calE[qp][k][m](k, m_local) = _specific_total_enthalpy[qp]; // H * D_{km}
464 33191424 : _calE[qp][k][m] += (1. - _fp.gamma()) * vel(m_local) * Ck_T; // (1-gam) * u_m * C_k^T
465 : }
466 :
467 : // E_{k4} (energy gradient term)
468 5531904 : _calE[qp][k][4].zero();
469 5531904 : _calE[qp][k][4] = _fp.gamma() * Ck_T;
470 : }
471 :
472 : // Compute the sum over ell of: A_ell grad(U_ell), store in DenseVector or Gradient object?
473 : // The gradient object might be more useful, since we are multiplying by VariableGradient
474 : // (which is a MooseArray of RealGradients) objects?
475 1843968 : RealVectorValue mom_resid = _calA[qp][0] * _grad_rho[qp] + _calA[qp][1] * _grad_rho_u[qp] +
476 1843968 : _calA[qp][2] * _grad_rho_v[qp] + _calA[qp][3] * _grad_rho_w[qp] +
477 1843968 : _calA[qp][4] * _grad_rho_et[qp];
478 :
479 : // No matrices/vectors for the energy residual strong form... just write it out like
480 : // the mass equation residual. See "Momentum SUPG terms prop. to energy residual"
481 : // section of the notes.
482 : Real energy_resid =
483 1843968 : (0.5 * (_fp.gamma() - 1.0) * velmag2 - _specific_total_enthalpy[qp]) * (vel * _grad_rho[qp]) +
484 3687936 : _specific_total_enthalpy[qp] * divU +
485 1843968 : (1. - _fp.gamma()) * (vel(0) * (vel * _grad_rho_u[qp]) + vel(1) * (vel * _grad_rho_v[qp]) +
486 1843968 : vel(2) * (vel * _grad_rho_w[qp])) +
487 1843968 : _fp.gamma() * (vel * _grad_rho_et[qp]);
488 :
489 : // Now for the actual residual values...
490 :
491 : // The density strong-residual
492 1843968 : _strong_residuals[qp][0] = _drho_dt[qp] + divU;
493 :
494 : // The x-momentum strong-residual, viscous terms neglected.
495 : // TODO: If we want to add viscous contributions back in, should this kernel
496 : // not inherit from NSViscousFluxBase so it can get tau values? This would
497 : // also involve shape function second derivative values.
498 1843968 : _strong_residuals[qp][1] = _drhou_dt[qp] + mom_resid(0);
499 :
500 : // The y-momentum strong residual, viscous terms neglected.
501 1843968 : _strong_residuals[qp][2] = _drhov_dt[qp] + mom_resid(1);
502 :
503 : // The z-momentum strong residual, viscous terms neglected.
504 1843968 : if (_mesh_dimension == 3)
505 0 : _strong_residuals[qp][3] = _drhow_dt[qp] + mom_resid(2);
506 : else
507 1843968 : _strong_residuals[qp][3] = 0.;
508 :
509 : // The energy equation strong residual
510 1843968 : _strong_residuals[qp][4] = _drhoE_dt[qp] + energy_resid;
511 1843968 : }
|
