channel_flow.h

#ifndef __CHANNEL_FLOW_H__
#define __CHANNEL_FLOW_H__

#include "periodic_turbulence.h"
#include "dg/dg_base.hpp"

namespace PHiLiP {
namespace FlowSolver {

#if PHILIP_DIM==1
using Triangulation = dealii::Triangulation<PHILIP_DIM>;
#else
using Triangulation = dealii::parallel::distributed::Triangulation<PHILIP_DIM>;
#endif

template <int dim, int nspecies, int nstate>
class ChannelFlow : public PeriodicTurbulence<dim,nspecies,nstate>
{
public:
    ChannelFlow(const Parameters::AllParameters *const parameters_input);

    ~ChannelFlow() {};

    std::shared_ptr<Triangulation> generate_grid() const override;

    void set_higher_order_grid(std::shared_ptr <DGBase<dim, nspecies, double>> dg) const override;

protected:
    const double channel_height; 
    const double half_channel_height; 
    const double channel_friction_velocity_reynolds_number; 
    const int number_of_cells_x_direction; 
    const int number_of_cells_y_direction; 
    const int number_of_cells_z_direction; 
    const double pi_val; 
    const double domain_length_x; 
    const double domain_length_y; 
    const double domain_length_z; 
    const double domain_volume; 

    const double channel_bulk_velocity_reynolds_number;

    const double channel_centerline_velocity_reynolds_number;

    double minimum_approximate_grid_spacing; 

    void display_additional_flow_case_specific_parameters() const override;

    void display_grid_parameters() const override;

    std::vector<double> get_mesh_step_size_y_direction() const;

    std::vector<double> get_mesh_step_size_y_direction_HOPW() const;

    std::vector<double> get_mesh_step_size_y_direction_Gullbrand() const;

    std::vector<double> get_mesh_step_size_y_direction_carton_de_wiart_et_al() const;

    dealii::Tensor<2,dim,double> zero_tensor; 

    std::shared_ptr< Physics::NavierStokes_ChannelFlowConstantSourceTerm_WallModel<dim,nspecies,dim+2,double> > navier_stokes_channel_flow_constant_source_term_wall_model_physics;

public:
    double get_adaptive_time_step(std::shared_ptr<DGBase<dim,nspecies,double>> dg) const override;

    double get_adaptive_time_step_initial(std::shared_ptr<DGBase<dim,nspecies,double>> dg) override;

    void compute_unsteady_data_and_write_to_table(
            const std::shared_ptr <ODE::ODESolverBase<dim, nspecies, double>> ode_solver,
            const std::shared_ptr <DGBase<dim, nspecies, double>> dg,
            const std::shared_ptr<dealii::TableHandler> unsteady_data_table,
            const bool do_write_unsteady_data_table_file) override;

    unsigned int get_number_of_degrees_of_freedom_per_state_from_poly_degree(const unsigned int poly_degree_input) const override;

    double get_average_wall_shear_stress(DGBase<dim, nspecies, double> &dg) const;

    double get_average_wall_shear_stress_from_wall_model(DGBase<dim, nspecies, double> &dg) const;

    double get_bulk_density() const; 
    double get_bulk_velocity() const; 
    double get_bulk_mass_flow_rate() const; 

    double get_skin_friction_coefficient_from_average_wall_shear_stress(const double avg_wall_shear_stress) const;

    void set_bulk_flow_quantities(DGBase<dim, nspecies, double> &dg);
private:
    double get_stretched_mesh_size(const int i) const;

    double bulk_density; 
    double bulk_mass_flow_rate; 
    double bulk_velocity; 
};

} // FlowSolver namespace
} // PHiLiP namespace
#endif