Skip to content

Commit

Permalink
Merge pull request #724 from Yashmandaokar/SharedTurbulenceModel
Browse files Browse the repository at this point in the history
Implementation of Standard k-ε Turbulence Model within ESPH
  • Loading branch information
Xiangyu-Hu authored Jan 30, 2025
2 parents 8a0759c + a2f0dd4 commit 4934037
Show file tree
Hide file tree
Showing 37 changed files with 142,937 additions and 13 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/ci.yml
Original file line number Diff line number Diff line change
Expand Up @@ -632,4 +632,4 @@ jobs:
if: ${{ steps.fourth-try.outcome == 'failure' }}
run: |
cd build
ctest --rerun-failed --output-on-failure --timeout 1000
ctest --rerun-failed --output-on-failure --timeout 1000
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@ void EulerianIntegration1stHalf<Inner<>, RiemannSolverType>::interaction(size_t
FluidStateIn state_j(rho_[index_j], vel_[index_j], p_[index_j]);
FluidStateOut interface_state = riemann_solver_.InterfaceState(state_i, state_j, e_ij);
Matd convect_flux = interface_state.rho_ * interface_state.vel_ * interface_state.vel_.transpose();

momentum_change_rate -= 2.0 * Vol_[index_i] * (convect_flux + interface_state.p_ * Matd::Identity()) * e_ij * dW_ijV_j;
}
dmom_dt_[index_i] = momentum_change_rate;
Expand Down Expand Up @@ -105,6 +106,7 @@ void EulerianIntegration2ndHalf<Inner<>, RiemannSolverType>::interaction(size_t

FluidStateIn state_j(rho_[index_j], vel_[index_j], p_[index_j]);
FluidStateOut interface_state = riemann_solver_.InterfaceState(state_i, state_j, e_ij);

mass_change_rate -= 2.0 * Vol_[index_i] * (interface_state.rho_ * interface_state.vel_).dot(e_ij) * dW_ijV_j;
}
dmass_dt_[index_i] = mass_change_rate;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -139,7 +139,6 @@ class PressureForceFromFluidInFVM : public BaseForceFromFluidInFVM
}
};
};

} // namespace fluid_dynamics
} // namespace SPH
#endif // COMMON_WEAKLY_COMPRESSIBLE_FVM_CLASSES_H
Original file line number Diff line number Diff line change
@@ -1,9 +1,9 @@
/**
* @file test_3d_incompressible_channel_flow.cpp
* @file test_3d_FVM_incompressible_channel_flow.cpp
* @brief This is the inviscid incompressible channel flow test for the realization of FVM in the SPHinXsys.
* @author Yash Mandaokar, Zhentong Wang and Xiangyu Hu
*/
#include "test_3d_incompressible_channel_flow.h"
#include "test_3d_FVM_incompressible_channel_flow.h"

using namespace SPH;
//----------------------------------------------------------------------
Expand Down Expand Up @@ -35,11 +35,11 @@ int main(int ac, char *av[])
// Define the main numerical methods used in the simulation.
// Note that there may be data dependence on the constructors of these methods.
//----------------------------------------------------------------------
SimpleDynamics<InvCFInitialCondition> initial_condition(air_block);
/** Here we introduce the limiter in the Riemann solver and 0 means the no extra numerical dissipation.
* the value is larger, the numerical dissipation larger. */
InteractionWithUpdate<fluid_dynamics::EulerianIntegration1stHalfInnerRiemann> pressure_relaxation(air_block_inner, 500.0);
InteractionWithUpdate<fluid_dynamics::EulerianIntegration2ndHalfInnerRiemann> density_relaxation(air_block_inner, 8000.0);
SimpleDynamics<InvCFInitialCondition> initial_condition(air_block);
/** Time step size with considering sound wave speed. */
ReduceDynamics<fluid_dynamics::WCAcousticTimeStepSizeInFVM> get_fluid_time_step_size(air_block, read_mesh_data.MinMeshEdge(), 0.6);
/** Boundary conditions set up */
Expand All @@ -57,9 +57,9 @@ int main(int ac, char *av[])
//----------------------------------------------------------------------
Real &physical_time = *sph_system.getSystemVariableDataByName<Real>("PhysicalTime");
size_t number_of_iterations = 0;
int screen_output_interval = 1000;
Real end_time = 30;
Real output_interval = 2; /**< time stamps for output. */
int screen_output_interval = 5000;
Real end_time = 15;
Real output_interval = 3; /**< time stamps for output. */
//----------------------------------------------------------------------
// Statistics for CPU time
//----------------------------------------------------------------------
Expand Down
Original file line number Diff line number Diff line change
@@ -1,11 +1,11 @@
/**
* @file test_3d_incompressible_channel_flow.h
* @file test_3d_FVM_incompressible_channel_flow.h
* @brief This is a test to show inviscid incompressible channel flow using .msh files from ICEM and FLUENT.
* @author Yash Mandaokar, Zhentong Wang and Xiangyu Hu
*/

#ifndef TEST_3D_INCOMPRESSIBLE_CHANNEL_FLOW_H
#define TEST_3D_INCOMPRESSIBLE_CHANNEL_FLOW_H
#ifndef TEST_3D_FVM_INCOMPRESSIBLE_CHANNEL_FLOW_H
#define TEST_3D_FVM_INCOMPRESSIBLE_CHANNEL_FLOW_H

#include "common_weakly_compressible_FVM_classes.h"

Expand All @@ -24,7 +24,6 @@ BoundingBox system_domain_bounds(Vec3d(-0.3, 0.0, 0.0), Vec3d(0.469846, 0.5, 0.0
//----------------------------------------------------------------------
Real rho0_f = 1.0; /**< Density. */
Real U_f = 1.0; /**< freestream velocity. */
Real P_ref = 101325; /*operating pressure fluent [Pa]*/
Real c_f = 10.0 * U_f;
Real mu_f = 0.0; /**< Dynamics viscosity. */

Expand Down Expand Up @@ -57,7 +56,10 @@ class InvCFInitialCondition
: FluidInitialCondition(sph_body),
rho_(particles_->registerStateVariable<Real>("Density")),
p_(particles_->registerStateVariable<Real>("Pressure")),
vel_(particles_->registerStateVariable<Vecd>("Velocity")){};
vel_(particles_->registerStateVariable<Vecd>("Velocity")),
Vol_(this->particles_->template getVariableDataByName<Real>("VolumetricMeasure")),
mass_(this->particles_->template getVariableDataByName<Real>("Mass")),
mom_(this->particles_->template getVariableDataByName<Vecd>("Momentum")){};

protected:
Real *rho_, *p_;
Expand All @@ -68,10 +70,14 @@ class InvCFInitialCondition
vel_[index_i][0] = 1.0;
vel_[index_i][1] = 0.0;
vel_[index_i][2] = 0.0;
mass_[index_i] = rho_[index_i] * Vol_[index_i];
mom_[index_i] = mass_[index_i] * vel_[index_i];
}

protected:
Vecd *vel_;
Real *Vol_, *mass_;
Vecd *mom_;
};
///----------------------------------------------------------------------
// InvCFBoundaryConditionSetup
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,130 @@
#ifndef TURBULENCEMODEL_CPP
#define TURBULENCEMODEL_CPP
#include "turbulence_model.h"
#include "sphinxsys.h"
namespace SPH
{
namespace fluid_dynamics
{
//=================================================================================================//
BaseTurbulence::BaseTurbulence(BaseInnerRelation& inner_relation, GhostCreationFromMesh& ghost_creator)
: BaseIntegration<DataDelegateInner>(inner_relation),
mom_(this->particles_->template getVariableDataByName<Vecd>("Momentum")),
dmom_dt_(this->particles_->template getVariableDataByName<Vecd>("MomentumChangeRate")),
dmass_dt_(this->particles_->template getVariableDataByName<Real>("MassChangeRate")),
K_prod_p_(this->particles_->template registerStateVariable<Real>("TKEProductionInWallAdjCell")),
K_prod_(this->particles_->template registerStateVariable<Real>("TKEProduction")),
Eps_p_(this->particles_->template registerStateVariable<Real>("DissipationRateInWallAdjCell")),
Tau_wall_(this->particles_->template registerStateVariable<Real>("WallShearStress")),
C_mu_(0.09), sigma_k_(1.0), sigma_eps_(1.3), C1_eps_(1.44), C2_eps_(1.92),
K_(this->particles_->template getVariableDataByName<Real>("TKE")),
Eps_(this->particles_->template getVariableDataByName<Real>("Dissipation")),
mu_t_(this->particles_->template getVariableDataByName<Real>("TurblunetViscosity")),
ghost_creator_(ghost_creator),
viscosity_(DynamicCast<Viscosity>(this, particles_->getBaseMaterial()))
{}
//=================================================================================================//
WallAdjacentCells::WallAdjacentCells(BaseInnerRelation &inner_relation, GhostCreationFromMesh &ghost_creator)
: BaseTurbulence(inner_relation, ghost_creator), wall_normal_(this->particles_->template registerStateVariable<Vecd>("WallNormal")),
wall_adjacent_cell_flag_(this->particles_->template registerStateVariable<Real>("FlagForWallAdjacentCells")),
yp_(this->particles_->template registerStateVariable<Real>("WallNormalDistance")),
ymax_(0.0)
{
walladjacentcellyp();
}
//=================================================================================================//
void WallAdjacentCells::walladjacentcellyp()
{
wall_adjacent_index_.resize(particles_->ParticlesBound());
wall_ghost_index_.resize(particles_->ParticlesBound());
wall_eij_.resize(particles_->ParticlesBound());

for (size_t boundary_type = 0; boundary_type < ghost_creator_.each_boundary_type_with_all_ghosts_index_.size(); ++boundary_type)
{
if (!ghost_creator_.each_boundary_type_with_all_ghosts_index_[boundary_type].empty())
{
for (size_t ghost_number = 0; ghost_number != ghost_creator_.each_boundary_type_with_all_ghosts_index_[boundary_type].size(); ++ghost_number)
{
size_t ghost_index = ghost_creator_.each_boundary_type_with_all_ghosts_index_[boundary_type][ghost_number];
size_t index_real = ghost_creator_.each_boundary_type_contact_real_index_[boundary_type][ghost_number];
wall_eij_[index_real] = ghost_creator_.each_boundary_type_with_all_ghosts_eij_[boundary_type][ghost_number];

if (boundary_type == 3)
{
wall_adjacent_index_[index_real] = index_real;
wall_ghost_index_[index_real] = ghost_index;
wall_adjacent_cell_flag_[index_real] = 1.0;
if ((pos_[index_real] - pos_[ghost_index]).dot(wall_eij_[index_real]) > ymax_)
{
ymax_ = (pos_[index_real] - pos_[ghost_index]).dot(wall_eij_[index_real]);
}
}
}
}
}
}
//=================================================================================================//
void WallAdjacentCells::update(size_t index_i, Real dt)
{
Vecd lower_wall = {pos_[index_i][0], 0.0};
Vecd upper_wall = {pos_[index_i][0], 2.0};
Vecd lower_wall_normal = {0.0, 1.0};
Vecd upper_wall_normal = {0.0, -1.0};

bool lower_wall_condition = ((pos_[index_i] - lower_wall).dot(lower_wall_normal) <= 1.0 * ymax_);
bool upper_wall_condition = ((pos_[index_i] - upper_wall).dot(upper_wall_normal) <= 1.0 * ymax_);

if (lower_wall_condition)
{
yp_[index_i] = (pos_[index_i] - lower_wall).dot(lower_wall_normal);
wall_normal_[index_i] = lower_wall_normal;
}
else if (upper_wall_condition)
{
yp_[index_i] = (pos_[index_i] - upper_wall).dot(upper_wall_normal);
wall_normal_[index_i] = upper_wall_normal;
}
}
//=================================================================================================//
StdWallFunctionFVM::StdWallFunctionFVM(BaseInnerRelation &inner_relation, GhostCreationFromMesh &ghost_creator)
: BaseTurbulence(inner_relation, ghost_creator), y_star_(this->particles_->template registerStateVariable<Real>("Ystar")),
yp_(this->particles_->template getVariableDataByName<Real>("WallNormalDistance")),
wall_normal_(this->particles_->template getVariableDataByName<Vecd>("WallNormal")),
vel_gradient_mat_(this->particles_->template registerStateVariable<Matd>("VelocityGradient")),
von_kar_(0.4187), E_(9.793)
{}
//=================================================================================================//
void StdWallFunctionFVM::nearwallquantities(size_t index_i)
{
y_star_[index_i] = (rho_[index_i] * std::pow(C_mu_, 0.25) * std::pow(K_[index_i], 0.5) * yp_[index_i]) / (viscosity_.ReferenceViscosity());
Real u_star;
Vecd veltangential = (vel_[index_i] - wall_normal_[index_i].dot(vel_[index_i]) * (wall_normal_[index_i]));

if (y_star_[index_i] >= 11.225)
{
u_star = (1.0 / von_kar_) * std::log(E_ * y_star_[index_i]);
mu_t_[index_i] = viscosity_.ReferenceViscosity() * ((y_star_[index_i]) / (1 / von_kar_ * std::log(E_ * y_star_[index_i])) - 1.0);

Tau_wall_[index_i] = (veltangential.norm() * std::pow(C_mu_, 0.25) * std::pow(K_[index_i], 0.5) * rho_[index_i]) / (u_star);
vel_gradient_mat_[index_i] = Matd::Zero();
vel_gradient_mat_[index_i](0, 1) = Tau_wall_[index_i] / (rho_[index_i] * pow(C_mu_, 0.25) * pow(K_[index_i], 0.5) * von_kar_ * yp_[index_i]);

K_prod_p_[index_i] = std::pow(Tau_wall_[index_i], 2.0) / (von_kar_ * rho_[index_i] * std::pow(C_mu_, 0.25) * std::pow(K_[index_i], 0.5) * yp_[index_i]);
Eps_p_[index_i] = (std::pow(C_mu_, 3.0 / 4.0) * std::pow(K_[index_i], 1.5)) / (von_kar_ * yp_[index_i]);

}
else if (y_star_[index_i] < 11.225)
{
u_star = y_star_[index_i];
Tau_wall_[index_i] = viscosity_.ReferenceViscosity() * veltangential.norm() / yp_[index_i];
vel_gradient_mat_[index_i] = Matd::Zero();
vel_gradient_mat_[index_i](0, 1) = Tau_wall_[index_i] / viscosity_.ReferenceViscosity();
K_prod_p_[index_i] = 0.0;
Eps_p_[index_i] = (K_[index_i] * 2.0 * viscosity_.ReferenceViscosity()) / (rho_[index_i] * yp_[index_i] * yp_[index_i]);
}
}
//=================================================================================================//
}// namespace fluid_dynamics

}// namespace SPH
#endif // TURBULENCEMODEL_CPP
Loading

0 comments on commit 4934037

Please sign in to comment.