37 #include "libmesh/utility.h"
38 #include "libmesh/string_to_enum.h"
39 #include "libmesh/getpot.h"
40 #include "libmesh/fem_system.h"
41 #include "libmesh/quadrature.h"
46 template<
class Conductivity>
61 template<
class Conductivity>
67 template<
class Conductivity>
71 libMesh::Utility::string_to_enum<GRINSEnums::FEFamily>( input(
"Physics/"+
axisymmetric_heat_transfer+
"/FE_family",
"LAGRANGE") );
100 template<
class Conductivity>
104 this->_dim = system->get_mesh().mesh_dimension();
106 _T_var = system->add_variable( _T_var_name, _T_order, _T_FE_family);
109 _u_r_var = system->add_variable(_u_r_var_name, _V_order, _V_FE_family);
110 _u_z_var = system->add_variable(_u_z_var_name, _V_order, _V_FE_family);
115 template<
class Conductivity>
119 system->time_evolving(_T_var);
123 template<
class Conductivity>
129 context.get_element_fe(_T_var)->get_JxW();
130 context.get_element_fe(_T_var)->get_phi();
131 context.get_element_fe(_T_var)->get_dphi();
132 context.get_element_fe(_T_var)->get_xyz();
134 context.get_side_fe(_T_var)->get_JxW();
135 context.get_side_fe(_T_var)->get_phi();
136 context.get_side_fe(_T_var)->get_dphi();
137 context.get_side_fe(_T_var)->get_xyz();
145 template<
class Conductivity>
150 #ifdef GRINS_USE_GRVY_TIMERS
151 this->_timer->BeginTimer(
"AxisymmetricHeatTransfer::element_time_derivative");
155 const unsigned int n_T_dofs = context.get_dof_indices(_T_var).size();
156 const unsigned int n_u_dofs = context.get_dof_indices(_u_r_var).size();
164 const std::vector<libMesh::Real> &JxW =
165 context.get_element_fe(_T_var)->get_JxW();
168 const std::vector<std::vector<libMesh::Real> >& T_phi =
169 context.get_element_fe(_T_var)->get_phi();
172 const std::vector<std::vector<libMesh::Real> >& vel_phi =
173 context.get_element_fe(_u_r_var)->get_phi();
177 const std::vector<std::vector<libMesh::RealGradient> >& T_gradphi =
178 context.get_element_fe(_T_var)->get_dphi();
181 const std::vector<libMesh::Point>& u_qpoint =
182 context.get_element_fe(_u_r_var)->get_xyz();
185 libMesh::DenseSubVector<libMesh::Number> &FT = context.get_elem_residual(_T_var);
187 libMesh::DenseSubMatrix<libMesh::Number> &KTT = context.get_elem_jacobian(_T_var, _T_var);
189 libMesh::DenseSubMatrix<libMesh::Number> &KTr = context.get_elem_jacobian(_T_var, _u_r_var);
190 libMesh::DenseSubMatrix<libMesh::Number> &KTz = context.get_elem_jacobian(_T_var, _u_z_var);
199 unsigned int n_qpoints = context.get_element_qrule().n_points();
201 for (
unsigned int qp=0; qp != n_qpoints; qp++)
203 const libMesh::Number r = u_qpoint[qp](0);
206 libMesh::Number u_r, u_z;
207 u_r = context.interior_value(_u_r_var, qp);
208 u_z = context.interior_value(_u_z_var, qp);
210 libMesh::Gradient grad_T;
211 grad_T = context.interior_gradient(_T_var, qp);
213 libMesh::NumberVectorValue U (u_r,u_z);
215 libMesh::Number k = this->_k( context, qp );
222 for (
unsigned int i=0; i != n_T_dofs; i++)
225 (-_rho*_Cp*T_phi[i][qp]*(U*grad_T)
226 -k*(T_gradphi[i][qp]*grad_T) );
228 if (compute_jacobian)
230 libmesh_assert (context.get_elem_solution_derivative() == 1.0);
232 for (
unsigned int j=0; j != n_T_dofs; j++)
237 KTT(i,j) += JxW[qp] * context.get_elem_solution_derivative() *r*
238 (-_rho*_Cp*T_phi[i][qp]*(U*T_gradphi[j][qp])
239 -k*(T_gradphi[i][qp]*T_gradphi[j][qp]));
245 for (
unsigned int j=0; j != n_T_dofs; j++)
248 KTT(i,j) -= JxW[qp] * context.get_elem_solution_derivative() *r*( dk_dT*T_phi[j][qp]*T_gradphi[i][qp]*grad_T );
254 for (
unsigned int j=0; j != n_u_dofs; j++)
256 KTr(i,j) += JxW[qp] * context.get_elem_solution_derivative() *r*(-_rho*_Cp*T_phi[i][qp]*(vel_phi[j][qp]*grad_T(0)));
257 KTz(i,j) += JxW[qp] * context.get_elem_solution_derivative() *r*(-_rho*_Cp*T_phi[i][qp]*(vel_phi[j][qp]*grad_T(1)));
265 #ifdef GRINS_USE_GRVY_TIMERS
266 this->_timer->EndTimer(
"AxisymmetricHeatTransfer::element_time_derivative");
272 template<
class Conductivity>
277 #ifdef GRINS_USE_GRVY_TIMERS
278 this->_timer->BeginTimer(
"AxisymmetricHeatTransfer::side_time_derivative");
281 std::vector<BoundaryID> ids = context.side_boundary_ids();
283 for( std::vector<BoundaryID>::const_iterator it = ids.begin();
284 it != ids.end(); it++ )
286 libmesh_assert (*it != libMesh::BoundaryInfo::invalid_id);
288 _bc_handler->apply_neumann_bcs( context, cache, compute_jacobian, *it );
291 #ifdef GRINS_USE_GRVY_TIMERS
292 this->_timer->EndTimer(
"AxisymmetricHeatTransfer::side_time_derivative");
298 template<
class Conductivity>
303 #ifdef GRINS_USE_GRVY_TIMERS
304 this->_timer->BeginTimer(
"AxisymmetricHeatTransfer::mass_residual");
311 const std::vector<libMesh::Real> &JxW =
312 context.get_element_fe(_T_var)->get_JxW();
315 const std::vector<std::vector<libMesh::Real> >& phi =
316 context.get_element_fe(_T_var)->get_phi();
319 const unsigned int n_T_dofs = context.get_dof_indices(_T_var).size();
322 const std::vector<libMesh::Point>& u_qpoint =
323 context.get_element_fe(_u_r_var)->get_xyz();
326 libMesh::DenseSubVector<libMesh::Real> &F =
327 context.get_elem_residual(_T_var);
329 libMesh::DenseSubMatrix<libMesh::Real> &M =
330 context.get_elem_jacobian(_T_var, _T_var);
332 unsigned int n_qpoints = context.get_element_qrule().n_points();
334 for (
unsigned int qp = 0; qp != n_qpoints; ++qp)
336 const libMesh::Number r = u_qpoint[qp](0);
344 context.interior_rate(_T_var, qp, T_dot);
346 for (
unsigned int i = 0; i != n_T_dofs; ++i)
348 F(i) -= JxW[qp]*r*(_rho*_Cp*T_dot*phi[i][qp] );
350 if( compute_jacobian )
352 for (
unsigned int j=0; j != n_T_dofs; j++)
355 M(i,j) -= JxW[qp] * context.get_elem_solution_rate_derivative() *r*_rho*_Cp*phi[j][qp]*phi[i][qp] ;
363 #ifdef GRINS_USE_GRVY_TIMERS
364 this->_timer->EndTimer(
"AxisymmetricHeatTransfer::mass_residual");
370 template<
class Conductivity>
372 (
const std::string & param_name,
377 _k.register_parameter(param_name, param_pointer);
virtual void read_input_options(const GetPot &input)
Read options from GetPot input file.
const PhysicsName incompressible_navier_stokes
GRINS::ICHandlingBase * _ic_handler
virtual void side_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for boundaries of elements on the domain boundary.
GRINS::BCHandlingBase * _bc_handler
INSTANTIATE_HEAT_TRANSFER_SUBCLASS(AxisymmetricHeatTransfer)
Physics abstract base class. Defines API for physics to be added to MultiphysicsSystem.
const PhysicsName axisymmetric_heat_transfer
virtual void register_parameter(const std::string ¶m_name, libMesh::ParameterMultiPointer< libMesh::Number > ¶m_pointer) const
Each subclass will register its copy of an independent.
virtual void init_variables(libMesh::FEMSystem *system)
Initialization AxisymmetricHeatTransfer variables.
Base class for reading and handling initial conditions for physics classes.
virtual void set_time_evolving_vars(libMesh::FEMSystem *system)
Sets velocity variables to be time-evolving.
~AxisymmetricHeatTransfer()
virtual void init_context(AssemblyContext &context)
Initialize context for added physics variables.
const std::string T_var_name_default
temperature
virtual void register_parameter(const std::string ¶m_name, libMesh::ParameterMultiPointer< libMesh::Number > ¶m_pointer) const
Each subclass will register its copy of an independent.
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...
Base class for reading and handling boundary conditions for physics classes.
const std::string u_z_var_name_default
z-velocity (axisymmetric case)
virtual void element_time_derivative(bool compute_jacobian, AssemblyContext &context, CachedValues &cache)
Time dependent part(s) of physics for element interiors.
const std::string u_r_var_name_default
r-velocity (axisymmetric case)
AxisymmetricHeatTransfer()