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 "INSFEMaterial.h"
11 : #include "MooseMesh.h"
12 : #include "Steady.h"
13 : #include "libmesh/quadrature.h"
14 :
15 : registerMooseObject("NavierStokesApp", INSFEMaterial);
16 : registerMooseObjectRenamed("NavierStokesApp", MDFluidMaterial, "02/01/2024 00:00", INSFEMaterial);
17 :
18 : InputParameters
19 510 : INSFEMaterial::validParams()
20 : {
21 510 : InputParameters params = Material::validParams();
22 :
23 510 : params.addClassDescription("Computes generic material properties related to simulation of fluid "
24 : "flow");
25 1020 : params.addRequiredCoupledVar("u", "velocity in x-direction");
26 1020 : params.addCoupledVar("v", "velocity in y-direction"); // required in 2D/3D
27 1020 : params.addCoupledVar("w", "velocity in z-direction"); // required in 3D
28 :
29 1020 : params.addRequiredCoupledVar("temperature", "temperature");
30 1020 : params.addRequiredCoupledVar("pressure", "pressure");
31 1020 : params.addRequiredParam<UserObjectName>("eos", "The name of equation of state object to use.");
32 1020 : params.addParam<Real>("scaling_velocity", "Global scaling velocity");
33 :
34 1020 : params.addParam<bool>(
35 1020 : "compute_turbulence_viscosity", false, "If turbulent viscosity will be computed.");
36 1020 : params.addParam<FunctionName>("mixing_length",
37 : "Prandtl mixing length to compute turbulent viscosity.");
38 1020 : params.addCoupledVar("turbulence_viscosity_var",
39 : "An aux variable turbulence_viscosity will be computed from.");
40 510 : return params;
41 0 : }
42 :
43 396 : INSFEMaterial::INSFEMaterial(const InputParameters & parameters)
44 : : Material(parameters),
45 792 : _mesh_dimension(_mesh.dimension()),
46 396 : _u_vel(coupledValue("u")),
47 396 : _v_vel(_mesh_dimension >= 2 ? coupledValue("v") : _zero),
48 396 : _w_vel(_mesh_dimension == 3 ? coupledValue("w") : _zero),
49 396 : _temperature(coupledValue("temperature")),
50 396 : _pressure(coupledValue("pressure")),
51 :
52 396 : _grad_u(coupledGradient("u")),
53 396 : _grad_v(_mesh_dimension >= 2 ? coupledGradient("v") : _grad_zero),
54 396 : _grad_w(_mesh_dimension == 3 ? coupledGradient("w") : _grad_zero),
55 :
56 396 : _viscous_stress_tensor(declareProperty<RealTensorValue>("viscous_stress_tensor")),
57 396 : _dynamic_viscosity(declareProperty<Real>("dynamic_viscosity")),
58 396 : _turbulence_viscosity(declareProperty<Real>("turbulence_viscosity")),
59 396 : _inertia_resistance_coeff(declareProperty<RealTensorValue>("inertia_resistance_coeff")),
60 396 : _viscous_resistance_coeff(declareProperty<RealTensorValue>("viscous_resistance_coeff")),
61 396 : _k(declareProperty<Real>("k_fluid")),
62 396 : _k_turbulence(declareProperty<Real>("k_turbulence")),
63 396 : _k_elem(declareProperty<Real>("k_fluid_elem")),
64 396 : _cp(declareProperty<Real>("cp_fluid")),
65 396 : _rho(declareProperty<Real>("rho_fluid")),
66 396 : _hsupg(declareProperty<Real>("hsupg")),
67 396 : _tauc(declareProperty<Real>("tauc")),
68 396 : _taum(declareProperty<Real>("taum")),
69 396 : _taue(declareProperty<Real>("taue")),
70 792 : _compute_visc_turbulenc(getParam<bool>("compute_turbulence_viscosity")),
71 792 : _has_turb_visc_auxvar(isParamValid("turbulence_viscosity_var")),
72 396 : _mixing_length(NULL),
73 396 : _turb_visc_auxvar(_has_turb_visc_auxvar ? coupledValue("turbulence_viscosity_var") : _zero),
74 792 : _has_scale_vel(isParamValid("scaling_velocity")),
75 396 : _scaling_velocity(_has_scale_vel ? getParam<Real>("scaling_velocity") : 1.),
76 792 : _eos(getUserObject<SinglePhaseFluidProperties>("eos"))
77 : {
78 396 : Executioner * myExe = _app.getExecutioner();
79 396 : Steady * ss = dynamic_cast<Steady *>(myExe);
80 396 : _bSteady = (ss != NULL);
81 :
82 396 : if (_compute_visc_turbulenc)
83 : {
84 : // If visc_turbulence is needed, either 'mixing_length' or 'turbulence_viscosity_var' is needed,
85 : // but 'mixing_length' and 'turbulence_viscosity_var' cannot be both specified
86 0 : if (isParamValid("mixing_length") == isParamValid("turbulence_viscosity_var"))
87 0 : mooseError("To compute turbulence viscosity, "
88 : "please provide either 'mixing_length' or 'turbulence_viscosity_var'");
89 :
90 : // If 'turbulence_viscosity_var' is not given, 'mixing_length' is used
91 0 : if (!_has_turb_visc_auxvar)
92 0 : _mixing_length = &getFunction("mixing_length");
93 : }
94 396 : }
95 :
96 : void
97 4222388 : INSFEMaterial::computeProperties()
98 : {
99 : // calculating element average velocity
100 4222388 : _u_elem = 0;
101 4222388 : _v_elem = 0;
102 4222388 : _w_elem = 0;
103 : // if(_qrule->n_points()<4) return; // only calculated element-based material properties
104 20263724 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
105 : {
106 16041336 : _u_elem += _u_vel[_qp] / _qrule->n_points();
107 16041336 : _v_elem += _v_vel[_qp] / _qrule->n_points();
108 16041336 : _w_elem += _w_vel[_qp] / _qrule->n_points();
109 : }
110 4222388 : _vel_mag = std::sqrt(_u_elem * _u_elem + _v_elem * _v_elem + _w_elem * _w_elem);
111 :
112 4222388 : _k_elem_val = 0.0;
113 :
114 20263724 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
115 16041336 : computeQpProperties();
116 :
117 20263724 : for (_qp = 0; _qp < _qrule->n_points(); _qp++)
118 16041336 : _k_elem[_qp] = _k_elem_val;
119 4222388 : }
120 :
121 : void
122 16041336 : INSFEMaterial::computeQpProperties()
123 : {
124 16041336 : computeFluidProperties();
125 :
126 : // Viscous Stress Tensor
127 16041336 : RealTensorValue grad_vel(_grad_u[_qp], _grad_v[_qp], _grad_w[_qp]);
128 :
129 16041336 : grad_vel += grad_vel.transpose();
130 :
131 : // Turbulent viscosity
132 16041336 : _turbulence_viscosity[_qp] = 0.0;
133 16041336 : _k_turbulence[_qp] = 0.0;
134 16041336 : if (_compute_visc_turbulenc)
135 : {
136 0 : if (_has_turb_visc_auxvar)
137 0 : _turbulence_viscosity[_qp] = _turb_visc_auxvar[_qp];
138 : else
139 : {
140 : Real strain_rate_mag = 0.0;
141 0 : for (unsigned int i = 0; i < 3; i++)
142 0 : for (unsigned int j = 0; j < 3; j++)
143 0 : strain_rate_mag += 0.5 * grad_vel(i, j) * grad_vel(i, j);
144 :
145 0 : Real L_mix = _mixing_length->value(_t, _q_point[_qp]);
146 :
147 0 : _turbulence_viscosity[_qp] = _rho[_qp] * L_mix * L_mix * std::sqrt(strain_rate_mag);
148 : }
149 0 : _k_turbulence[_qp] = _cp[_qp] * _turbulence_viscosity[_qp] / 0.9;
150 : }
151 :
152 16041336 : _k_elem_val += (_k[_qp] + _k_turbulence[_qp]) / _qrule->n_points();
153 :
154 16041336 : Real div_vel = _grad_u[_qp](0) + _grad_v[_qp](1) + _grad_w[_qp](2);
155 :
156 : // Add diagonal terms
157 64165344 : for (unsigned int i = 0; i < 3; i++)
158 48124008 : grad_vel(i, i) -= 2.0 / 3.0 * div_vel;
159 :
160 16041336 : _viscous_stress_tensor[_qp] = (_dynamic_viscosity[_qp] + _turbulence_viscosity[_qp]) * grad_vel;
161 :
162 : // place holder which will be handled by specific closure correlations
163 16041336 : _inertia_resistance_coeff[_qp] = RealTensorValue(0, 0, 0, 0, 0, 0, 0, 0, 0);
164 16041336 : _viscous_resistance_coeff[_qp] = RealTensorValue(0, 0, 0, 0, 0, 0, 0, 0, 0);
165 :
166 : // Compute stabilization parameters:
167 16041336 : computeHSUPG();
168 16041336 : computeTau();
169 16041336 : }
170 :
171 : void
172 16041336 : INSFEMaterial::computeFluidProperties()
173 : {
174 16041336 : _dynamic_viscosity[_qp] = _eos.mu_from_p_T(_pressure[_qp], _temperature[_qp]);
175 16041336 : _k[_qp] = _eos.k_from_p_T(_pressure[_qp], _temperature[_qp]);
176 16041336 : _cp[_qp] = _eos.cp_from_p_T(_pressure[_qp], _temperature[_qp]);
177 16041336 : _rho[_qp] = _eos.rho_from_p_T(_pressure[_qp], _temperature[_qp]);
178 16041336 : }
179 :
180 : void
181 16041336 : INSFEMaterial::computeHSUPG()
182 : {
183 : // Just use hmin for the element!
184 16041336 : _hsupg[_qp] = _current_elem->hmin();
185 16041336 : }
186 :
187 : void
188 16041336 : INSFEMaterial::computeTau()
189 : {
190 : // element-based stabilization parameters
191 16041336 : Real h2 = _hsupg[_qp] * _hsupg[_qp];
192 16041336 : Real sqrt_term_supg = 4 * _vel_mag * _vel_mag / h2;
193 :
194 : Real sqrt_term_pspg = 0.;
195 16041336 : if (_has_scale_vel)
196 0 : sqrt_term_pspg = 4 * _scaling_velocity * _scaling_velocity / h2;
197 : else
198 : sqrt_term_pspg = sqrt_term_supg;
199 :
200 16041336 : if (!_bSteady)
201 : {
202 16041336 : sqrt_term_supg += 4. / _dt / _dt;
203 16041336 : sqrt_term_pspg += 4. / _dt / _dt;
204 : }
205 :
206 16041336 : Real vu = (_dynamic_viscosity[_qp] + _turbulence_viscosity[_qp]) / _rho[_qp];
207 16041336 : Real visc_term = 4 * vu / h2;
208 :
209 16041336 : Real diffusivity = _k[_qp] / _rho[_qp] / _cp[_qp];
210 16041336 : Real diff_term = 4 * diffusivity / h2;
211 :
212 16041336 : _tauc[_qp] = 1. / std::sqrt(sqrt_term_pspg);
213 16041336 : _taum[_qp] = 1. / std::sqrt(sqrt_term_supg + visc_term * visc_term);
214 16041336 : _taue[_qp] = 1. / std::sqrt(sqrt_term_supg + diff_term * diff_term);
215 16041336 : }
|
