36 #include "libmesh/quadrature.h"
41 template<
class Mu,
class SH,
class TC>
49 template<
class Mu,
class SH,
class TC>
55 template<
class Mu,
class SH,
class TC>
60 #ifdef GRINS_USE_GRVY_TIMERS
61 this->_timer->BeginTimer(
"LowMachNavierStokesBraackStabilization::element_time_derivative");
64 this->assemble_continuity_time_deriv( compute_jacobian, context );
65 this->assemble_momentum_time_deriv( compute_jacobian, context );
66 this->assemble_energy_time_deriv( compute_jacobian, context );
68 #ifdef GRINS_USE_GRVY_TIMERS
69 this->_timer->EndTimer(
"LowMachNavierStokesBraackStabilization::element_time_derivative");
74 template<
class Mu,
class SH,
class TC>
79 #ifdef GRINS_USE_GRVY_TIMERS
80 this->_timer->BeginTimer(
"LowMachNavierStokesBraackStabilization::mass_residual");
83 this->assemble_continuity_mass_residual( compute_jacobian, context );
84 this->assemble_momentum_mass_residual( compute_jacobian, context );
85 this->assemble_energy_mass_residual( compute_jacobian, context );
87 #ifdef GRINS_USE_GRVY_TIMERS
88 this->_timer->EndTimer(
"LowMachNavierStokesBraackStabilization::mass_residual");
93 template<
class Mu,
class SH,
class TC>
98 const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
101 const std::vector<libMesh::Real> &JxW =
102 context.get_element_fe(this->_u_var)->get_JxW();
105 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
106 context.get_element_fe(this->_p_var)->get_dphi();
108 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_p_var);
110 unsigned int n_qpoints = context.get_element_qrule().n_points();
112 for (
unsigned int qp=0; qp != n_qpoints; qp++)
114 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
116 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
117 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
119 libMesh::Real T = context.interior_value( this->_T_var, qp );
120 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
122 libMesh::Real mu = this->_mu(T);
123 libMesh::Real k = this->_k(T);
124 libMesh::Real cp = this->_cp(T);
126 libMesh::RealGradient U( context.interior_value( this->_u_var, qp ),
127 context.interior_value( this->_v_var, qp ) );
128 if( this->_dim == 3 )
129 U(2) = context.interior_value( this->_w_var, qp );
131 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
132 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
134 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
135 libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
139 for (
unsigned int i=0; i != n_p_dofs; i++)
141 Fp(i) += ( tau_M*RM_s*p_dphi[i][qp]
142 + tau_E*RE_s*(U*p_dphi[i][qp])/T )*JxW[qp];
150 template<
class Mu,
class SH,
class TC>
155 const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
158 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_v_var).size());
160 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_w_var).size());
163 const std::vector<libMesh::Real> &JxW =
164 context.get_element_fe(this->_u_var)->get_JxW();
167 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
168 context.get_element_fe(this->_u_var)->get_dphi();
170 const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
171 context.get_element_fe(this->_u_var)->get_d2phi();
173 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_u_var);
174 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_v_var);
175 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_w_var);
177 unsigned int n_qpoints = context.get_element_qrule().n_points();
179 for (
unsigned int qp=0; qp != n_qpoints; qp++)
181 libMesh::Real T = context.interior_value( this->_T_var, qp );
182 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
184 libMesh::Real mu = this->_mu(T);
186 libMesh::RealGradient U( context.interior_value(this->_u_var, qp),
187 context.interior_value(this->_v_var, qp) );
189 libMesh::RealGradient grad_u = context.interior_gradient(this->_u_var, qp);
190 libMesh::RealGradient grad_v = context.interior_gradient(this->_v_var, qp);
191 libMesh::RealGradient grad_w;
193 if( this->_dim == 3 )
195 U(2) = context.interior_value(this->_w_var, qp);
196 grad_w = context.interior_gradient(this->_w_var, qp);
199 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
201 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
202 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
204 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
205 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
207 libMesh::Real RC_s = this->compute_res_continuity_steady( context, qp );
208 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
210 for (
unsigned int i=0; i != n_u_dofs; i++)
212 Fu(i) += ( tau_C*RC_s*u_gradphi[i][qp](0)
213 + tau_M*RM_s(0)*rho*U*u_gradphi[i][qp]
214 + mu*tau_M*RM_s(0)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1)
215 + u_hessphi[i][qp](0,0) + u_hessphi[i][qp](0,1)
216 - 2.0/3.0*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,0)) )
219 Fv(i) += ( tau_C*RC_s*u_gradphi[i][qp](1)
220 + tau_M*RM_s(1)*rho*U*u_gradphi[i][qp]
221 + mu*tau_M*RM_s(1)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1)
222 + u_hessphi[i][qp](1,0) + u_hessphi[i][qp](1,1)
223 - 2.0/3.0*(u_hessphi[i][qp](0,1) + u_hessphi[i][qp](1,1)) )
226 if( this->_dim == 3 )
228 Fu(i) += mu*tau_M*RM_s(0)*(u_hessphi[i][qp](2,2) + u_hessphi[i][qp](0,2)
229 - 2.0/3.0*u_hessphi[i][qp](2,0))*JxW[qp];
231 Fv(i) += mu*tau_M*RM_s(1)*(u_hessphi[i][qp](2,2) + u_hessphi[i][qp](1,2)
232 - 2.0/3.0*u_hessphi[i][qp](2,1))*JxW[qp];
234 Fw(i) += ( tau_C*RC_s*u_gradphi[i][qp](2)
235 + tau_M*RM_s(2)*rho*U*u_gradphi[i][qp]
236 + mu*tau_M*RM_s(2)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2)
237 + u_hessphi[i][qp](2,0) + u_hessphi[i][qp](2,1) + u_hessphi[i][qp](2,2)
238 - 2.0/3.0*(u_hessphi[i][qp](0,2) + u_hessphi[i][qp](1,2)
239 + u_hessphi[i][qp](2,2))
249 template<
class Mu,
class SH,
class TC>
254 const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
257 const std::vector<libMesh::Real> &JxW =
258 context.get_element_fe(this->_T_var)->get_JxW();
261 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
262 context.get_element_fe(this->_T_var)->get_dphi();
264 const std::vector<std::vector<libMesh::RealTensor> >& T_hessphi =
265 context.get_element_fe(this->_T_var)->get_d2phi();
267 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_T_var);
269 unsigned int n_qpoints = context.get_element_qrule().n_points();
271 for (
unsigned int qp=0; qp != n_qpoints; qp++)
273 libMesh::Number u, v;
274 u = context.interior_value(this->_u_var, qp);
275 v = context.interior_value(this->_v_var, qp);
277 libMesh::Gradient grad_T = context.interior_gradient(this->_T_var, qp);
279 libMesh::NumberVectorValue U(u,v);
281 U(2) = context.interior_value(this->_w_var, qp);
283 libMesh::Real T = context.interior_value( this->_T_var, qp );
284 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
286 libMesh::Real k = this->_k(T);
287 libMesh::Real cp = this->_cp(T);
289 libMesh::Number rho_cp = rho*this->_cp(T);
291 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
293 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
294 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
296 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
298 libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
300 for (
unsigned int i=0; i != n_T_dofs; i++)
302 FT(i) += ( rho_cp*tau_E*RE_s*U*T_gradphi[i][qp]
303 + tau_E*RE_s*k*(T_hessphi[i][qp](0,0) + T_hessphi[i][qp](1,1) + T_hessphi[i][qp](2,2) )
312 template<
class Mu,
class SH,
class TC>
317 const unsigned int n_p_dofs = context.get_dof_indices(this->_p_var).size();
320 const std::vector<libMesh::Real> &JxW =
321 context.get_element_fe(this->_u_var)->get_JxW();
324 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
325 context.get_element_fe(this->_p_var)->get_dphi();
327 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_p_var);
329 unsigned int n_qpoints = context.get_element_qrule().n_points();
331 for (
unsigned int qp=0; qp != n_qpoints; qp++)
333 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
335 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
336 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
338 libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
339 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
341 libMesh::Real mu = this->_mu(T);
342 libMesh::Real k = this->_k(T);
343 libMesh::Real cp = this->_cp(T);
345 libMesh::RealGradient U( context.fixed_interior_value( this->_u_var, qp ),
346 context.fixed_interior_value( this->_v_var, qp ) );
347 if( this->_dim == 3 )
348 U(2) = context.fixed_interior_value( this->_w_var, qp );
350 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
351 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
353 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp,
false );
354 libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
358 for (
unsigned int i=0; i != n_p_dofs; i++)
360 Fp(i) += ( tau_M*RM_t*p_dphi[i][qp]
361 + tau_E*RE_t*(U*p_dphi[i][qp])/T
369 template<
class Mu,
class SH,
class TC>
374 const unsigned int n_u_dofs = context.get_dof_indices(this->_u_var).size();
377 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_v_var).size());
379 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_w_var).size());
382 const std::vector<libMesh::Real> &JxW =
383 context.get_element_fe(this->_u_var)->get_JxW();
386 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
387 context.get_element_fe(this->_u_var)->get_dphi();
389 const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
390 context.get_element_fe(this->_u_var)->get_d2phi();
392 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_u_var);
393 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_v_var);
394 libMesh::DenseSubVector<libMesh::Number> &Fw = context.get_elem_residual(this->_w_var);
396 unsigned int n_qpoints = context.get_element_qrule().n_points();
397 for (
unsigned int qp=0; qp != n_qpoints; qp++)
399 libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
400 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
402 libMesh::Real mu = this->_mu(T);
404 libMesh::RealGradient U( context.fixed_interior_value(this->_u_var, qp),
405 context.fixed_interior_value(this->_v_var, qp) );
407 libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->_u_var, qp);
408 libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->_v_var, qp);
409 libMesh::RealGradient grad_w;
411 if( this->_dim == 3 )
413 U(2) = context.fixed_interior_value(this->_w_var, qp);
414 grad_w = context.fixed_interior_gradient(this->_w_var, qp);
417 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
419 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
420 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
422 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
423 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
425 libMesh::Real RC_t = this->compute_res_continuity_transient( context, qp );
426 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
428 for (
unsigned int i=0; i != n_u_dofs; i++)
430 Fu(i) += ( tau_C*RC_t*u_gradphi[i][qp](0)
431 + tau_M*RM_t(0)*rho*U*u_gradphi[i][qp]
432 + mu*tau_M*RM_t(0)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1)
433 + u_hessphi[i][qp](0,0) + u_hessphi[i][qp](0,1)
434 - 2.0/3.0*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,0)) )
437 Fv(i) += ( tau_C*RC_t*u_gradphi[i][qp](1)
438 + tau_M*RM_t(1)*rho*U*u_gradphi[i][qp]
439 + mu*tau_M*RM_t(1)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1)
440 + u_hessphi[i][qp](1,0) + u_hessphi[i][qp](1,1)
441 - 2.0/3.0*(u_hessphi[i][qp](0,1) + u_hessphi[i][qp](1,1)) )
444 if( this->_dim == 3 )
446 Fw(i) -= mu*tau_M*RM_t(0)*(u_hessphi[i][qp](2,2) + u_hessphi[i][qp](0,2)
447 - 2.0/3.0*u_hessphi[i][qp](2,0))*JxW[qp];
449 Fv(i) -= mu*tau_M*RM_t(1)*(u_hessphi[i][qp](2,2) + u_hessphi[i][qp](1,2)
450 - 2.0/3.0*u_hessphi[i][qp](2,1))*JxW[qp];
452 Fw(i) -= ( tau_C*RC_t*u_gradphi[i][qp](2)
453 + tau_M*RM_t(2)*rho*U*u_gradphi[i][qp]
454 + mu*tau_M*RM_t(2)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2)
455 + u_hessphi[i][qp](2,0) + u_hessphi[i][qp](2,1) + u_hessphi[i][qp](2,2)
456 - 2.0/3.0*(u_hessphi[i][qp](0,2) + u_hessphi[i][qp](1,2)
457 + u_hessphi[i][qp](2,2))
467 template<
class Mu,
class SH,
class TC>
472 const unsigned int n_T_dofs = context.get_dof_indices(this->_T_var).size();
475 const std::vector<libMesh::Real> &JxW =
476 context.get_element_fe(this->_T_var)->get_JxW();
479 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
480 context.get_element_fe(this->_T_var)->get_dphi();
482 const std::vector<std::vector<libMesh::RealTensor> >& T_hessphi =
483 context.get_element_fe(this->_T_var)->get_d2phi();
485 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_T_var);
487 unsigned int n_qpoints = context.get_element_qrule().n_points();
489 for (
unsigned int qp=0; qp != n_qpoints; qp++)
491 libMesh::Number u, v;
492 u = context.fixed_interior_value(this->_u_var, qp);
493 v = context.fixed_interior_value(this->_v_var, qp);
495 libMesh::Gradient grad_T = context.fixed_interior_gradient(this->_T_var, qp);
497 libMesh::NumberVectorValue U(u,v);
499 U(2) = context.fixed_interior_value(this->_w_var, qp);
501 libMesh::Real T = context.fixed_interior_value( this->_T_var, qp );
502 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
504 libMesh::Real k = this->_k(T);
505 libMesh::Real cp = this->_cp(T);
507 libMesh::Number rho_cp = rho*cp;
509 libMesh::FEBase* fe = context.get_element_fe(this->_u_var);
511 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
512 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
514 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp,
false );
516 libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
518 for (
unsigned int i=0; i != n_T_dofs; i++)
520 FT(i) += ( rho_cp*tau_E*RE_t*U*T_gradphi[i][qp]
521 + tau_E*RE_t*k*(T_hessphi[i][qp](0,0) + T_hessphi[i][qp](1,1) + T_hessphi[i][qp](2,2))
void assemble_energy_mass_residual(bool compute_jacobian, AssemblyContext &context)
Adds VMS-based stabilization to LowMachNavierStokes physics class.
void assemble_energy_time_deriv(bool compute_jacobian, AssemblyContext &context)
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
void assemble_momentum_mass_residual(bool compute_jacobian, AssemblyContext &context)
void assemble_continuity_time_deriv(bool compute_jacobian, AssemblyContext &context)
virtual ~LowMachNavierStokesBraackStabilization()
Adds VMS-based stabilization to LowMachNavierStokes physics class.
void assemble_continuity_mass_residual(bool compute_jacobian, AssemblyContext &context)
LowMachNavierStokesBraackStabilization()
void assemble_momentum_time_deriv(bool compute_jacobian, AssemblyContext &context)
virtual void mass_residual(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Mass matrix part(s) for element interiors. All boundary terms lie within the time_derivative part...