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->_press_var.p()).size();
101 const std::vector<libMesh::Real> &JxW =
102 context.get_element_fe(this->_flow_vars.u())->get_JxW();
105 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
106 context.get_element_fe(this->_press_var.p())->get_dphi();
108 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p());
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->_flow_vars.u());
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->_temp_vars.T(), 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->_flow_vars.u(), qp ),
127 context.interior_value( this->_flow_vars.v(), qp ) );
128 if( this->_dim == 3 )
129 U(2) = context.interior_value( this->_flow_vars.w(), 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->_flow_vars.u()).size();
158 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
160 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
163 const std::vector<libMesh::Real> &JxW =
164 context.get_element_fe(this->_flow_vars.u())->get_JxW();
167 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
168 context.get_element_fe(this->_flow_vars.u())->get_dphi();
170 const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
171 context.get_element_fe(this->_flow_vars.u())->get_d2phi();
173 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
174 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
175 libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
177 if( this->_dim == 3 )
179 Fw = &context.get_elem_residual(this->_flow_vars.w());
182 unsigned int n_qpoints = context.get_element_qrule().n_points();
184 for (
unsigned int qp=0; qp != n_qpoints; qp++)
186 libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
187 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
189 libMesh::Real mu = this->_mu(T);
191 libMesh::RealGradient U( context.interior_value(this->_flow_vars.u(), qp),
192 context.interior_value(this->_flow_vars.v(), qp) );
194 libMesh::RealGradient grad_u = context.interior_gradient(this->_flow_vars.u(), qp);
195 libMesh::RealGradient grad_v = context.interior_gradient(this->_flow_vars.v(), qp);
196 libMesh::RealGradient grad_w;
198 if( this->_dim == 3 )
200 U(2) = context.interior_value(this->_flow_vars.w(), qp);
201 grad_w = context.interior_gradient(this->_flow_vars.w(), qp);
204 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
206 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
207 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
209 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu, this->_is_steady );
210 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
212 libMesh::Real RC_s = this->compute_res_continuity_steady( context, qp );
213 libMesh::RealGradient RM_s = this->compute_res_momentum_steady( context, qp );
215 for (
unsigned int i=0; i != n_u_dofs; i++)
217 Fu(i) += ( tau_C*RC_s*u_gradphi[i][qp](0)
218 + tau_M*RM_s(0)*rho*U*u_gradphi[i][qp]
219 + mu*tau_M*RM_s(0)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1)
220 + u_hessphi[i][qp](0,0) + u_hessphi[i][qp](0,1)
221 - 2.0/3.0*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,0)) )
224 Fv(i) += ( tau_C*RC_s*u_gradphi[i][qp](1)
225 + tau_M*RM_s(1)*rho*U*u_gradphi[i][qp]
226 + mu*tau_M*RM_s(1)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1)
227 + u_hessphi[i][qp](1,0) + u_hessphi[i][qp](1,1)
228 - 2.0/3.0*(u_hessphi[i][qp](0,1) + u_hessphi[i][qp](1,1)) )
231 if( this->_dim == 3 )
233 Fu(i) += mu*tau_M*RM_s(0)*(u_hessphi[i][qp](2,2) + u_hessphi[i][qp](0,2)
234 - 2.0/3.0*u_hessphi[i][qp](2,0))*JxW[qp];
236 Fv(i) += mu*tau_M*RM_s(1)*(u_hessphi[i][qp](2,2) + u_hessphi[i][qp](1,2)
237 - 2.0/3.0*u_hessphi[i][qp](2,1))*JxW[qp];
239 (*Fw)(i) += ( tau_C*RC_s*u_gradphi[i][qp](2)
240 + tau_M*RM_s(2)*rho*U*u_gradphi[i][qp]
241 + mu*tau_M*RM_s(2)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2)
242 + u_hessphi[i][qp](2,0) + u_hessphi[i][qp](2,1) + u_hessphi[i][qp](2,2)
243 - 2.0/3.0*(u_hessphi[i][qp](0,2) + u_hessphi[i][qp](1,2)
244 + u_hessphi[i][qp](2,2))
254 template<
class Mu,
class SH,
class TC>
259 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
262 const std::vector<libMesh::Real> &JxW =
263 context.get_element_fe(this->_temp_vars.T())->get_JxW();
266 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
267 context.get_element_fe(this->_temp_vars.T())->get_dphi();
269 const std::vector<std::vector<libMesh::RealTensor> >& T_hessphi =
270 context.get_element_fe(this->_temp_vars.T())->get_d2phi();
272 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T());
274 unsigned int n_qpoints = context.get_element_qrule().n_points();
276 for (
unsigned int qp=0; qp != n_qpoints; qp++)
278 libMesh::Number u, v;
279 u = context.interior_value(this->_flow_vars.u(), qp);
280 v = context.interior_value(this->_flow_vars.v(), qp);
282 libMesh::Gradient grad_T = context.interior_gradient(this->_temp_vars.T(), qp);
284 libMesh::NumberVectorValue U(u,v);
286 U(2) = context.interior_value(this->_flow_vars.w(), qp);
288 libMesh::Real T = context.interior_value( this->_temp_vars.T(), qp );
289 libMesh::Real rho = this->rho( T, this->get_p0_steady( context, qp ) );
291 libMesh::Real k = this->_k(T);
292 libMesh::Real cp = this->_cp(T);
294 libMesh::Number rho_cp = rho*this->_cp(T);
296 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
298 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
299 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
301 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp, this->_is_steady );
303 libMesh::Real RE_s = this->compute_res_energy_steady( context, qp );
305 for (
unsigned int i=0; i != n_T_dofs; i++)
307 FT(i) += ( rho_cp*tau_E*RE_s*U*T_gradphi[i][qp]
308 + tau_E*RE_s*k*(T_hessphi[i][qp](0,0) + T_hessphi[i][qp](1,1) + T_hessphi[i][qp](2,2) )
317 template<
class Mu,
class SH,
class TC>
322 const unsigned int n_p_dofs = context.get_dof_indices(this->_press_var.p()).size();
325 const std::vector<libMesh::Real> &JxW =
326 context.get_element_fe(this->_flow_vars.u())->get_JxW();
329 const std::vector<std::vector<libMesh::RealGradient> >& p_dphi =
330 context.get_element_fe(this->_press_var.p())->get_dphi();
332 libMesh::DenseSubVector<libMesh::Number> &Fp = context.get_elem_residual(this->_press_var.p());
334 unsigned int n_qpoints = context.get_element_qrule().n_points();
336 for (
unsigned int qp=0; qp != n_qpoints; qp++)
338 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
340 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
341 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
343 libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
344 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
346 libMesh::Real mu = this->_mu(T);
347 libMesh::Real k = this->_k(T);
348 libMesh::Real cp = this->_cp(T);
350 libMesh::RealGradient U( context.fixed_interior_value( this->_flow_vars.u(), qp ),
351 context.fixed_interior_value( this->_flow_vars.v(), qp ) );
352 if( this->_dim == 3 )
353 U(2) = context.fixed_interior_value( this->_flow_vars.w(), qp );
355 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
356 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
358 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp,
false );
359 libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
363 for (
unsigned int i=0; i != n_p_dofs; i++)
365 Fp(i) += ( tau_M*RM_t*p_dphi[i][qp]
366 + tau_E*RE_t*(U*p_dphi[i][qp])/T
374 template<
class Mu,
class SH,
class TC>
379 const unsigned int n_u_dofs = context.get_dof_indices(this->_flow_vars.u()).size();
382 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.v()).size());
384 libmesh_assert (n_u_dofs == context.get_dof_indices(this->_flow_vars.w()).size());
387 const std::vector<libMesh::Real> &JxW =
388 context.get_element_fe(this->_flow_vars.u())->get_JxW();
391 const std::vector<std::vector<libMesh::RealGradient> >& u_gradphi =
392 context.get_element_fe(this->_flow_vars.u())->get_dphi();
394 const std::vector<std::vector<libMesh::RealTensor> >& u_hessphi =
395 context.get_element_fe(this->_flow_vars.u())->get_d2phi();
397 libMesh::DenseSubVector<libMesh::Number> &Fu = context.get_elem_residual(this->_flow_vars.u());
398 libMesh::DenseSubVector<libMesh::Number> &Fv = context.get_elem_residual(this->_flow_vars.v());
399 libMesh::DenseSubVector<libMesh::Real>* Fw = NULL;
401 if( this->_dim == 3 )
403 Fw = &context.get_elem_residual(this->_flow_vars.w());
406 unsigned int n_qpoints = context.get_element_qrule().n_points();
407 for (
unsigned int qp=0; qp != n_qpoints; qp++)
409 libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
410 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
412 libMesh::Real mu = this->_mu(T);
414 libMesh::RealGradient U( context.fixed_interior_value(this->_flow_vars.u(), qp),
415 context.fixed_interior_value(this->_flow_vars.v(), qp) );
417 libMesh::RealGradient grad_u = context.fixed_interior_gradient(this->_flow_vars.u(), qp);
418 libMesh::RealGradient grad_v = context.fixed_interior_gradient(this->_flow_vars.v(), qp);
419 libMesh::RealGradient grad_w;
421 if( this->_dim == 3 )
423 U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp);
424 grad_w = context.fixed_interior_gradient(this->_flow_vars.w(), qp);
427 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
429 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
430 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
432 libMesh::Real tau_M = this->_stab_helper.compute_tau_momentum( context, qp, g, G, rho, U, mu,
false );
433 libMesh::Real tau_C = this->_stab_helper.compute_tau_continuity( tau_M, g );
435 libMesh::Real RC_t = this->compute_res_continuity_transient( context, qp );
436 libMesh::RealGradient RM_t = this->compute_res_momentum_transient( context, qp );
438 for (
unsigned int i=0; i != n_u_dofs; i++)
440 Fu(i) += ( tau_C*RC_t*u_gradphi[i][qp](0)
441 + tau_M*RM_t(0)*rho*U*u_gradphi[i][qp]
442 + mu*tau_M*RM_t(0)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1)
443 + u_hessphi[i][qp](0,0) + u_hessphi[i][qp](0,1)
444 - 2.0/3.0*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,0)) )
447 Fv(i) += ( tau_C*RC_t*u_gradphi[i][qp](1)
448 + tau_M*RM_t(1)*rho*U*u_gradphi[i][qp]
449 + mu*tau_M*RM_t(1)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1)
450 + u_hessphi[i][qp](1,0) + u_hessphi[i][qp](1,1)
451 - 2.0/3.0*(u_hessphi[i][qp](0,1) + u_hessphi[i][qp](1,1)) )
454 if( this->_dim == 3 )
456 Fu(i) -= mu*tau_M*RM_t(0)*(u_hessphi[i][qp](2,2) + u_hessphi[i][qp](0,2)
457 - 2.0/3.0*u_hessphi[i][qp](2,0))*JxW[qp];
459 Fv(i) -= mu*tau_M*RM_t(1)*(u_hessphi[i][qp](2,2) + u_hessphi[i][qp](1,2)
460 - 2.0/3.0*u_hessphi[i][qp](2,1))*JxW[qp];
462 (*Fw)(i) -= ( tau_C*RC_t*u_gradphi[i][qp](2)
463 + tau_M*RM_t(2)*rho*U*u_gradphi[i][qp]
464 + mu*tau_M*RM_t(2)*(u_hessphi[i][qp](0,0) + u_hessphi[i][qp](1,1) + u_hessphi[i][qp](2,2)
465 + u_hessphi[i][qp](2,0) + u_hessphi[i][qp](2,1) + u_hessphi[i][qp](2,2)
466 - 2.0/3.0*(u_hessphi[i][qp](0,2) + u_hessphi[i][qp](1,2)
467 + u_hessphi[i][qp](2,2))
477 template<
class Mu,
class SH,
class TC>
482 const unsigned int n_T_dofs = context.get_dof_indices(this->_temp_vars.T()).size();
485 const std::vector<libMesh::Real> &JxW =
486 context.get_element_fe(this->_temp_vars.T())->get_JxW();
489 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
490 context.get_element_fe(this->_temp_vars.T())->get_dphi();
492 const std::vector<std::vector<libMesh::RealTensor> >& T_hessphi =
493 context.get_element_fe(this->_temp_vars.T())->get_d2phi();
495 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(this->_temp_vars.T());
497 unsigned int n_qpoints = context.get_element_qrule().n_points();
499 for (
unsigned int qp=0; qp != n_qpoints; qp++)
501 libMesh::Number u, v;
502 u = context.fixed_interior_value(this->_flow_vars.u(), qp);
503 v = context.fixed_interior_value(this->_flow_vars.v(), qp);
505 libMesh::Gradient grad_T = context.fixed_interior_gradient(this->_temp_vars.T(), qp);
507 libMesh::NumberVectorValue U(u,v);
509 U(2) = context.fixed_interior_value(this->_flow_vars.w(), qp);
511 libMesh::Real T = context.fixed_interior_value( this->_temp_vars.T(), qp );
512 libMesh::Real rho = this->rho( T, this->get_p0_transient( context, qp ) );
514 libMesh::Real k = this->_k(T);
515 libMesh::Real cp = this->_cp(T);
517 libMesh::Number rho_cp = rho*cp;
519 libMesh::FEBase* fe = context.get_element_fe(this->_flow_vars.u());
521 libMesh::RealGradient g = this->_stab_helper.compute_g( fe, context, qp );
522 libMesh::RealTensor G = this->_stab_helper.compute_G( fe, context, qp );
524 libMesh::Real tau_E = this->_stab_helper.compute_tau_energy( context, qp, g, G, rho, U, k, cp,
false );
526 libMesh::Real RE_t = this->compute_res_energy_transient( context, qp );
528 for (
unsigned int i=0; i != n_T_dofs; i++)
530 FT(i) += ( rho_cp*tau_E*RE_t*U*T_gradphi[i][qp]
531 + 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...