Skip to content

Commit

Permalink
Merge pull request #633 from Xiangyu-Hu/eulerian_supersonic_flow_with…
Browse files Browse the repository at this point in the history
…_new_BC

add a new Eulerian supersonic flow around a cylinder case
  • Loading branch information
ZhentongWang authored Aug 18, 2024
2 parents c1e062a + f4c678b commit 3ed7659
Show file tree
Hide file tree
Showing 12 changed files with 847 additions and 7 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -26,13 +26,16 @@ CompressibleFluidInitialCondition::CompressibleFluidInitialCondition(SPHBody &sp
E_(*particles_->getVariableDataByName<Real>("TotalEnergy")) {}
//=================================================================================================//
EulerianCompressibleAcousticTimeStepSize::
EulerianCompressibleAcousticTimeStepSize(SPHBody &sph_body)
EulerianCompressibleAcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL)
: AcousticTimeStepSize(sph_body),
rho_(*particles_->getVariableDataByName<Real>("Density")),
p_(*particles_->getVariableDataByName<Real>("Pressure")),
vel_(*particles_->getVariableDataByName<Vecd>("Velocity")),
smoothing_length_(sph_body.sph_adaptation_->ReferenceSmoothingLength()),
compressible_fluid_(CompressibleFluid(1.0, 1.4)){};
compressible_fluid_(CompressibleFluid(1.0, 1.4))
{
acousticCFL_ = acousticCFL;
};
//=================================================================================================//
Real EulerianCompressibleAcousticTimeStepSize::reduce(size_t index_i, Real dt)
{
Expand All @@ -41,7 +44,7 @@ Real EulerianCompressibleAcousticTimeStepSize::reduce(size_t index_i, Real dt)
//=================================================================================================//
Real EulerianCompressibleAcousticTimeStepSize::outputResult(Real reduced_value)
{
return 0.6 / Dimensions * smoothing_length_ / (reduced_value + TinyReal);
return acousticCFL_ / Dimensions * smoothing_length_ / (reduced_value + TinyReal);
}
//=================================================================================================//
} // namespace fluid_dynamics
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -103,7 +103,7 @@ class EulerianCompressibleAcousticTimeStepSize : public AcousticTimeStepSize
Real smoothing_length_;

public:
explicit EulerianCompressibleAcousticTimeStepSize(SPHBody &sph_body);
explicit EulerianCompressibleAcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL = 0.6);
virtual ~EulerianCompressibleAcousticTimeStepSize(){};

Real reduce(size_t index_i, Real dt = 0.0);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -80,8 +80,26 @@ CompressibleFluidStarState HLLCWithLimiterRiemannSolver::
{
Real ul = -e_ij.dot(state_i.vel_);
Real ur = -e_ij.dot(state_j.vel_);
Real s_l = ul - compressible_fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_);
Real s_r = ur + compressible_fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_);

//one approach to obtain the smallest and largest wave velocity
//Real s_l = ul - compressible_fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_);
//Real s_r = ur + compressible_fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_);

//another approach to obtain the smallest and largest wave velocity
Vecd vl = state_i.vel_ - ul * (-e_ij);
Vecd vr = state_j.vel_ - ur * (-e_ij);
Real R_lf = state_j.rho_ / state_i.rho_;
Real u_tlide = (ul + ur * R_lf) / (1.0 + R_lf);
Real v_tlide = (vl.norm() + vr.norm() * R_lf) / (1.0 + R_lf);
Real q_tlide = u_tlide;
Real hl = (state_i.E_ + state_i.p_) / state_i.rho_;
Real hr = (state_j.E_ + state_j.p_) / state_j.rho_;
Real h_tlide = (hl + hr * R_lf) / (1.0 + R_lf);
Real sound_tlide = sqrt((1.4 - 1.0) * (h_tlide - 0.5 * (u_tlide * u_tlide + v_tlide * v_tlide)));
Real s_l = SMIN(ul - compressible_fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_), q_tlide - sound_tlide);
Real s_r = SMAX(ur + compressible_fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_), q_tlide + sound_tlide);

//obtain the middle wave velocity
Real rhol_cl = compressible_fluid_i_.getSoundSpeed(state_i.p_, state_i.rho_) * state_i.rho_;
Real rhor_cr = compressible_fluid_j_.getSoundSpeed(state_j.p_, state_j.rho_) * state_j.rho_;
Real clr = (rhol_cl + rhor_cr) / (state_i.rho_ + state_j.rho_);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ class HLLCWithLimiterRiemannSolver
Real limiter_parameter_;

public:
HLLCWithLimiterRiemannSolver(CompressibleFluid &compressible_fluid_i, CompressibleFluid &compressible_fluid_j, Real limiter_parameter = 5.0);
HLLCWithLimiterRiemannSolver(CompressibleFluid &compressible_fluid_i, CompressibleFluid &compressible_fluid_j, Real limiter_parameter = 1.0);
CompressibleFluidStarState getInterfaceState(const CompressibleFluidState &state_i, const CompressibleFluidState &state_j, const Vecd &e_ij);
};
} // namespace SPH
Expand Down
Original file line number Diff line number Diff line change
@@ -0,0 +1,162 @@
#include "eulerian_ghost_boundary.h"

#include "base_particles.hpp"
namespace SPH
{
//=================================================================================================//
GhostCreationInESPH::GhostCreationInESPH(BaseInnerRelation &inner_relation, Ghost<ReserveSizeFactor> &ghost_boundary)
: DataDelegateInner(inner_relation), sph_body_(inner_relation.getSPHBody()),
ghost_boundary_(ghost_boundary), get_inner_neighbor_(&inner_relation.getSPHBody()),
indicator_(*particles_->getVariableDataByName<int>("Indicator")),
Vol_(*particles_->getVariableDataByName<Real>("VolumetricMeasure")),
pos_(*particles_->getVariableDataByName<Vecd>("Position")),
ghost_bound_(ghost_boundary.GhostBound())
{
ghostGenerationAndAddToConfiguration();
};
//=================================================================================================//
void GhostCreationInESPH::ghostGenerationAndAddToConfiguration()
{
ghost_bound_.second = ghost_bound_.first;
RealAndGhostParticleData real_and_ghost_necessary_data;

for (size_t index_i = 0; index_i != particles_->TotalRealParticles(); ++index_i)
{
if (indicator_[index_i] == 1)
{
Vecd gradient_summation = Vecd::Zero();
Neighborhood &inner_neighborhood = inner_configuration_[index_i];
for (size_t n = 0; n != inner_neighborhood.current_size_; ++n)
{
size_t index_j = inner_neighborhood.j_[n];
Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j];
gradient_summation += dW_ijV_j * inner_neighborhood.e_ij_[n];
}
Real ghost_particle_dW_ijV_j = -gradient_summation.norm();
Vecd ghost_particle_eij = -gradient_summation / ghost_particle_dW_ijV_j;
Real distance_to_ghost = fabs(sph_body_.getInitialShape().findSignedDistance(pos_[index_i]));
Vecd diaplacement_to_ghost = distance_to_ghost * ghost_particle_eij;
Vecd ghost_position = pos_[index_i] - diaplacement_to_ghost;
mutex_create_ghost_particle_.lock();
size_t ghost_particle_index = ghost_bound_.second;
ghost_bound_.second++;
ghost_boundary_.checkWithinGhostSize(ghost_bound_);
particles_->updateGhostParticle(ghost_particle_index, index_i);
//Here, we assue the volume as the same the real particle
Vol_[ghost_particle_index] = Vol_[index_i];
pos_[ghost_particle_index] = ghost_position;
mutex_create_ghost_particle_.unlock();
Real double_distance_to_ghost = 2.0 * distance_to_ghost;
Real ghost_particle_dW_ij = ghost_particle_dW_ijV_j / (Vol_[ghost_particle_index] + TinyReal);
get_inner_neighbor_(inner_neighborhood, double_distance_to_ghost, ghost_particle_dW_ij, ghost_particle_eij, ghost_particle_index);

//add all necessary data of real particle and corresponding ghost particle
real_and_ghost_necessary_data.real_index_ = index_i;
real_and_ghost_necessary_data.ghost_index_ = ghost_particle_index;
real_and_ghost_necessary_data.e_ij_ghost_ = ghost_particle_eij;
real_and_ghost_particle_data_.push_back(real_and_ghost_necessary_data);
}
}
}
//=================================================================================================//
GhostBoundaryConditionSetupInESPH::
GhostBoundaryConditionSetupInESPH(BaseInnerRelation &inner_relation, GhostCreationInESPH &ghost_creation)
: DataDelegateInner(inner_relation),
rho_(*particles_->getVariableDataByName<Real>("Density")),
Vol_(*particles_->getVariableDataByName<Real>("VolumetricMeasure")),
mass_(*particles_->getVariableDataByName<Real>("Mass")),
vel_(*particles_->getVariableDataByName<Vecd>("Velocity")),
pos_(*particles_->getVariableDataByName<Vecd>("Position")),
mom_(*particles_->getVariableDataByName<Vecd>("Momentum")),
ghost_bound_(ghost_creation.ghost_bound_),
real_and_ghost_particle_data_(ghost_creation.real_and_ghost_particle_data_),
boundary_type_(*particles_->registerSharedVariable<int>("BoundaryType")),
W0_(sph_body_.sph_adaptation_->getKernel()->W0(ZeroVecd))
{
setupBoundaryTypes();
}
//=================================================================================================//
void GhostBoundaryConditionSetupInESPH::resetBoundaryConditions()
{
for (size_t ghost_number = 0; ghost_number != real_and_ghost_particle_data_.size(); ++ghost_number)
{
size_t index_i = real_and_ghost_particle_data_[ghost_number].real_index_;
size_t ghost_index = real_and_ghost_particle_data_[ghost_number].ghost_index_;
Vecd e_ig = real_and_ghost_particle_data_[ghost_number].e_ij_ghost_;

// Dispatch the appropriate boundary condition
switch (boundary_type_[index_i])
{
case 3: // this refer to the different types of wall boundary conditions
applyNonSlipWallBoundary(ghost_index, index_i);
applyReflectiveWallBoundary(ghost_index, index_i, e_ig);
break;
case 4:
applyTopBoundary(ghost_index, index_i);
break;
case 5:
applyPressureOutletBC(ghost_index, index_i);
break;
case 7:
applySymmetryBoundary(ghost_index, index_i, e_ig);
break;
case 9:
applyFarFieldBoundary(ghost_index, index_i);
break;
case 10:
applyGivenValueInletFlow(ghost_index);
applyVelocityInletFlow(ghost_index, index_i);
break;
case 36:
applyOutletBoundary(ghost_index, index_i);
break;
}
}
}
//=================================================================================================//
GhostKernelGradientUpdate::GhostKernelGradientUpdate(BaseInnerRelation &inner_relation)
: LocalDynamics(inner_relation.getSPHBody()), DataDelegateInner(inner_relation),
Vol_(*particles_->getVariableDataByName<Real>("VolumetricMeasure")),
kernel_gradient_original_summation_(*particles_->registerSharedVariable<Vecd>("KernelGradientOriginalSummation")),
indicator_(*particles_->getVariableDataByName<int>("Indicator")){};
//=================================================================================================//
void GhostKernelGradientUpdate::interaction(size_t index_i, Real dt)
{
if (indicator_[index_i] == 1)
{
Neighborhood &inner_neighborhood = inner_configuration_[index_i];
Vecd gradient_summation = Vecd::Zero();
for (size_t n = 0; n != inner_neighborhood.current_size_; ++n)
{
size_t index_j = inner_neighborhood.j_[n];
Vecd &e_ij = inner_neighborhood.e_ij_[n];
Real dW_ijV_j = inner_neighborhood.dW_ij_[n] * Vol_[index_j];
if (index_j < particles_->TotalRealParticles())
{
gradient_summation += dW_ijV_j * e_ij;
}
}
kernel_gradient_original_summation_[index_i] = gradient_summation;
}
}
//=================================================================================================//
void GhostKernelGradientUpdate::update(size_t index_i, Real dt)
{
if (indicator_[index_i] == 1)
{
Neighborhood &inner_neighborhood = inner_configuration_[index_i];
for (size_t n = 0; n != inner_neighborhood.current_size_; ++n)
{
size_t index_j = inner_neighborhood.j_[n];
if (index_j >= particles_->TotalRealParticles())
{
Vecd ghost_dW_ijV_j_with_e_ij = -kernel_gradient_original_summation_[index_i];
Real ghost_dW_ijV_j = -fabs(ghost_dW_ijV_j_with_e_ij.norm());
inner_neighborhood.dW_ij_[n] = ghost_dW_ijV_j / Vol_[index_j];
inner_neighborhood.e_ij_[n] = ghost_dW_ijV_j_with_e_ij / (ghost_dW_ijV_j + TinyReal);
}
}
}
}
//=================================================================================================//
} // namespace SPH
Original file line number Diff line number Diff line change
@@ -0,0 +1,115 @@
/* ------------------------------------------------------------------------- *
* SPHinXsys *
* ------------------------------------------------------------------------- *
* SPHinXsys (pronunciation: s'finksis) is an acronym from Smoothed Particle *
* Hydrodynamics for industrial compleX systems. It provides C++ APIs for *
* physical accurate simulation and aims to model coupled industrial dynamic *
* systems including fluid, solid, multi-body dynamics and beyond with SPH *
* (smoothed particle hydrodynamics), a meshless computational method using *
* particle discretization. *
* *
* SPHinXsys is partially funded by German Research Foundation *
* (Deutsche Forschungsgemeinschaft) DFG HU1527/6-1, HU1527/10-1, *
* HU1527/12-1 and HU1527/12-4. *
* *
* Portions copyright (c) 2017-2023 Technical University of Munich and *
* the authors' affiliations. *
* *
* Licensed under the Apache License, Version 2.0 (the "License"); you may *
* not use this file except in compliance with the License. You may obtain a *
* copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
* *
* ------------------------------------------------------------------------- */
/**
* @file eulerian_ghost_boundary.h
* @brief This is the particle dynamics for domain bounding
* @author Zhentong Wang, Xiangyu Hu
*/
#ifndef Eulerian_GHOST_BOUNDARY_H
#define Eulerian_GHOST_BOUNDARY_H

#include "base_general_dynamics.h"
#include "particle_reserve.h"
#include "unstructured_mesh.h"
namespace SPH
{
struct RealAndGhostParticleData
{
size_t real_index_, ghost_index_;
Vecd e_ij_ghost_;
};

/**
* @class GhostCreationInESPH
* @brief The extra storage for the boundary particles, similar with the ghost particle,
* is established to obey the zero-order consistency.
*/
class GhostCreationInESPH : public DataDelegateInner
{
public:
explicit GhostCreationInESPH(BaseInnerRelation &inner_relation, Ghost<ReserveSizeFactor> &ghost_boundary);
virtual ~GhostCreationInESPH(){};
std::vector<RealAndGhostParticleData> real_and_ghost_particle_data_;
void ghostGenerationAndAddToConfiguration();

protected:
SPHBody &sph_body_;
Ghost<ReserveSizeFactor> &ghost_boundary_;
std::mutex mutex_create_ghost_particle_; /**< mutex exclusion for memory conflict */
NeighborBuilderInnerInFVM get_inner_neighbor_;
StdLargeVec<int> &indicator_;
StdLargeVec<Real> &Vol_;
StdLargeVec<Vecd> &pos_;

public:
std::pair<size_t, size_t> &ghost_bound_;
};

//----------------------------------------------------------------------
// GhostBoundaryConditionSetupInESPH
//----------------------------------------------------------------------
class GhostBoundaryConditionSetupInESPH : public DataDelegateInner
{
public:
GhostBoundaryConditionSetupInESPH(BaseInnerRelation &inner_relation, GhostCreationInESPH &ghost_creation);
virtual ~GhostBoundaryConditionSetupInESPH(){};
virtual void applyReflectiveWallBoundary(size_t ghost_index, size_t index_i, Vecd e_ij){};
virtual void applyNonSlipWallBoundary(size_t ghost_index, size_t index_i){};
virtual void applyGivenValueInletFlow(size_t ghost_index){};
virtual void applyOutletBoundary(size_t ghost_index, size_t index_i){};
virtual void applyTopBoundary(size_t ghost_index, size_t index_i){};
virtual void applyFarFieldBoundary(size_t ghost_index, size_t index_i){};
virtual void applyPressureOutletBC(size_t ghost_index, size_t index_i){};
virtual void applySymmetryBoundary(size_t ghost_index, size_t index_i, Vecd e_ij){};
virtual void applyVelocityInletFlow(size_t ghost_index, size_t index_i){};
virtual void setupBoundaryTypes(){};
void resetBoundaryConditions();

protected:
StdLargeVec<Real> &rho_, &Vol_, &mass_;
StdLargeVec<Vecd> &vel_, &pos_, &mom_;
std::pair<size_t, size_t> &ghost_bound_;
std::vector<RealAndGhostParticleData> real_and_ghost_particle_data_;
StdLargeVec<int> &boundary_type_;
Real W0_;
};
/**
* @class GhostKernelGradientUpdate
* @brief Here, this class is used to update the boundary particle
after implementing the kernel correction to strictly follow the zero-order consistency
*/
class GhostKernelGradientUpdate : public LocalDynamics, public DataDelegateInner
{
public:
explicit GhostKernelGradientUpdate(BaseInnerRelation &inner_relation);
virtual ~GhostKernelGradientUpdate(){};
void interaction(size_t index_i, Real dt = 0.0);
void update(size_t index_i, Real dt = 0.0);

protected:
StdLargeVec<Real> &Vol_;
StdLargeVec<Vecd> &kernel_gradient_original_summation_;
StdLargeVec<int> &indicator_;
};
} // namespace SPH
#endif // Eulerian_GHOST_BOUNDARY_H
Loading

0 comments on commit 3ed7659

Please sign in to comment.