Skip to content

Commit

Permalink
recent changes
Browse files Browse the repository at this point in the history
  • Loading branch information
wangfeng_real committed Jul 23, 2023
1 parent 847944d commit 514e295
Show file tree
Hide file tree
Showing 7 changed files with 92 additions and 46 deletions.
33 changes: 24 additions & 9 deletions tests/unit_tests_src/for_2D_build/test_turbulent_module/test_turbulent_module.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ class BaseTurbulentModule
wall_boundary(sph_system, makeShared<WallBoundary>("Wall"))
{
std::cout << "constructor for BaseTurbulentModule" << std::endl;
water_block.defineParticlesAndMaterial<FluidParticles, WeaklyCompressibleFluid>(rho0_f, c_f, mu_f);
water_block.defineParticlesAndMaterial<BaseParticles, WeaklyCompressibleFluid>(rho0_f, c_f, mu_f);
water_block.generateParticles<ParticleGeneratorLattice>();
wall_boundary.defineParticlesAndMaterial<SolidParticles, Solid>();
wall_boundary.generateParticles<ParticleGeneratorLattice>();
Expand All @@ -95,14 +95,16 @@ class TurbulentModule : public :: testing :: Test, public BaseTurbulentModule
InnerRelation water_block_inner;
ComplexRelation water_block_complex_relation;

InteractionWithUpdate<fluid_dynamics::K_TurtbulentModelComplex, SequencedPolicy> k_equation_relaxation;
InteractionWithUpdate<fluid_dynamics::E_TurtbulentModelComplex> epsilon_equation_relaxation;
InteractionDynamics<fluid_dynamics::TKEnergyAccComplex, SequencedPolicy> turbulent_kinetic_energy_acceleration;
InteractionWithUpdate<fluid_dynamics::K_TurtbulentModelInner, SequencedPolicy> k_equation_relaxation;
InteractionWithUpdate<fluid_dynamics::E_TurtbulentModelInner> epsilon_equation_relaxation;
InteractionDynamics<fluid_dynamics::TKEnergyAccInner, SequencedPolicy> turbulent_kinetic_energy_acceleration;
SimpleDynamics<NormalDirectionFromBodyShape> wall_boundary_normal_direction;
InteractionDynamics<fluid_dynamics::StandardWallFunctionCorrection, SequencedPolicy> standard_wall_function_correction;
InteractionDynamics<fluid_dynamics::TurbulentViscousAccelerationWithWall, SequencedPolicy> turbulent_viscous_acceleration;
SimpleDynamics<fluid_dynamics::TurbulentEddyViscosity, SequencedPolicy> update_eddy_viscosity;

InteractionWithUpdate<CorrectedConfigurationInner> correct_configuration;

InteractionWithUpdate<fluid_dynamics::DensitySummationFreeStreamComplex> update_density_by_summation;
BodyStatesRecordingToVtp write_body_states;
size_t number_of_iterations;
Expand All @@ -115,9 +117,11 @@ class TurbulentModule : public :: testing :: Test, public BaseTurbulentModule
water_block_inner(water_block),
water_block_complex_relation(water_block_inner, { &wall_boundary }),

k_equation_relaxation(water_block_complex_relation),
epsilon_equation_relaxation(water_block_complex_relation),
turbulent_kinetic_energy_acceleration(water_block_complex_relation),
correct_configuration(water_block_inner),

k_equation_relaxation(water_block_inner),
epsilon_equation_relaxation(water_block_inner),
turbulent_kinetic_energy_acceleration(water_block_inner),
wall_boundary_normal_direction(wall_boundary),
standard_wall_function_correction(water_block_complex_relation),
turbulent_viscous_acceleration(water_block_complex_relation),
Expand Down Expand Up @@ -202,7 +206,7 @@ class Test_K_Epsilon_Equation : public fluid_dynamics::FluidInitialCondition

std::vector<double> loadInputData(int num_data, int num_file, std::string file_name)
{
std::ifstream file("./MappingData/FVM16_basedOnSPH4_4/" + file_name, std::ios::binary);
std::ifstream file("./MappingData/FromPy10/" + file_name, std::ios::binary);
if (file)
{
std::string line;
Expand Down Expand Up @@ -279,8 +283,18 @@ TEST_F(TurbulentModule, TestTurbulentKineticEnergyEquation)
int num_iter = 0;
while (GlobalStaticVariables::physical_time_<100.0)
{
update_eddy_viscosity.exec();
//** Try to introduce B correction *
correct_configuration.exec();
//update_eddy_viscosity.exec();


//** Test viscous force *
standard_wall_function_correction.exec(); //the wall viscous relies on wall Func values
turbulent_viscous_acceleration.exec();
write_body_states.writeToFile();
std::cout <<"Test viscous" << std::endl;
system("pause");

//turbulent_kinetic_energy_acceleration.exec();
for (int index_i = 0; index_i < num_fluid_particle; ++index_i)
{
Expand All @@ -301,6 +315,7 @@ TEST_F(TurbulentModule, TestTurbulentKineticEnergyEquation)
{
test_k_ep_equation.clear_acc_prior(index_i);
}
system("pause");
}

ASSERT_NEAR(1, 1, 0.02);
Expand Down
16 changes: 11 additions & 5 deletions tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -66,11 +66,11 @@ int main(int ac, char* av[])
InteractionDynamics<fluid_dynamics::StandardWallFunctionCorrection,SequencedPolicy> standard_wall_function_correction(water_block_complex_relation);

/** TurbulentViscous cal. uses friction velocity and Y+ that are defined in WallFunction . */
InteractionDynamics<fluid_dynamics::TurbulentViscousAccelerationWithWall> turbulent_viscous_acceleration(water_block_complex_relation);
InteractionDynamics<fluid_dynamics::TurbulentViscousAccelerationWithWall, SequencedPolicy> turbulent_viscous_acceleration(water_block_complex_relation);
//InteractionDynamics<fluid_dynamics::ViscousAccelerationWithWall> viscous_acceleration(water_block_complex_relation);


InteractionDynamics<fluid_dynamics::TransportVelocityCorrectionComplex> transport_velocity_correction(water_block_complex_relation);
InteractionDynamics<fluid_dynamics::TransportVelocityCorrectionComplex,SequencedPolicy> transport_velocity_correction(water_block_complex_relation);
InteractionWithUpdate<fluid_dynamics::SpatialTemporalFreeSurfaceIdentificationComplex>
inlet_outlet_surface_particle_indicator(water_block_complex_relation);
InteractionWithUpdate<fluid_dynamics::DensitySummationFreeStreamComplex> update_density_by_summation(water_block_complex_relation);
Expand All @@ -91,7 +91,13 @@ int main(int ac, char* av[])
/** Turbulent eddy viscosity calculation needs values of Wall Y start. */
SimpleDynamics<fluid_dynamics::TurbulentEddyViscosity> update_eddy_viscosity(water_block);


/** Try to introduce B correction */
InteractionWithUpdate<CorrectedConfigurationInner> correct_configuration(water_block_inner);




BodyAlignedBoxByParticle emitter(
water_block, makeShared<AlignedBoxShape>(Transform(Vec2d(emitter_translation)), emitter_halfsize));
SimpleDynamics<fluid_dynamics::EmitterInflowInjection> emitter_inflow_injection(emitter, 50, 0);
Expand Down Expand Up @@ -163,14 +169,15 @@ int main(int ac, char* av[])
turbulent_viscous_acceleration.exec();
//viscous_acceleration.exec();
transport_velocity_correction.exec();

//** Try to introduce B correction *
correct_configuration.exec();
/** Dynamics including pressure relaxation. */
Real relaxation_time = 0.0;
while (relaxation_time < Dt)
{
dt = SMIN(get_fluid_time_step_size.exec(), Dt - relaxation_time);

//turbulent_kinetic_energy_acceleration.exec();
turbulent_kinetic_energy_acceleration.exec();

pressure_relaxation.exec(dt);

Expand All @@ -191,7 +198,6 @@ int main(int ac, char* av[])
//write_body_states.writeToFile();

}

if (number_of_iterations % screen_output_interval == 0)
{
std::cout << std::fixed << std::setprecision(9) << "N=" << number_of_iterations << " Time = "
Expand Down
2 changes: 1 addition & 1 deletion tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ using namespace SPH; // Namespace cite here.
//----------------------------------------------------------------------
// Basic geometry parameters and numerical setup.
//----------------------------------------------------------------------
Real DL = 60; /**< Reference length. */
Real DL = 120; /**< Reference length. */
Real DH = 2; /**< Reference and the height of main channel. */
Real resolution_ref = 0.1; /**< Initial reference particle spacing. */
Real BW = resolution_ref * 4; /**< Reference size of the emitter. */
Expand Down
30 changes: 18 additions & 12 deletions tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_complex.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,9 +86,7 @@ namespace SPH
k_production_[index_i] = k_production;

//** for test */
lap_k_[index_i] += k_lap;
lap_k_term_[index_i] = 0.0;
lap_k_term_[index_i] = (mu_ + turbu_mu_i / sigma_k) * lap_k_[index_i];
k_diffusion_[index_i] += k_lap;
}
//=================================================================================================//
void E_TurtbulentModelComplex::
Expand Down Expand Up @@ -179,21 +177,24 @@ namespace SPH
//vel_derivative = 2.0 * (vel_i - vel_ave_k[index_j]) / (r_ij + 0.01 * this->smoothing_length_);
//acceleration += 2.0 * (this->mu_+ turbu_mu_i) * vel_derivative * contact_neighborhood.dW_ijV_j_[n] / rho_i;
//This is to check whether the wall-sub-nearest fluid particles fric, velo. is zero or not
//if (index_i > 2000 && GlobalStaticVariables::physical_time_ > 5. && vel_fric_i.dot(vel_fric_i) <= 0.0+TinyReal)
//{
// system("pause");
// std::cout << index_j << std::endl;
//}
if (index_i > 2000 && GlobalStaticVariables::physical_time_ > 5. && vel_fric_i.dot(vel_fric_i) <= 0.0+TinyReal&& contact_neighborhood.current_size_>2)
{
system("pause");
std::cout << index_j << std::endl;
std::cout << vel_fric_i << std::endl;
std::cout << distance_to_wall << std::endl;
std::cout << contact_neighborhood.current_size_ << std::endl;
}
//vel_derivative = 2.0 * vel_fric_i.dot(vel_fric_i)* direction_vel_fric;
vel_derivative = distance_to_wall * vel_fric_i.dot(vel_fric_i) * direction_vel_fric / (r_ij + 0.01 * this->smoothing_length_);
acceleration += 4.0 * vel_derivative * contact_neighborhood.dW_ijV_j_[n] ;
}
}
//for test
Real wall_viscous_factor = 1.0;
this->visc_acc_wall_[index_i] = wall_viscous_factor * acceleration;

this->visc_acc_wall_[index_i] = acceleration;

this->acc_prior_[index_i] += acceleration;
this->acc_prior_[index_i] += wall_viscous_factor * acceleration;
}
//=================================================================================================//
void StandardWallFunctionCorrection::interaction(size_t index_i, Real dt)
Expand Down Expand Up @@ -254,6 +255,11 @@ namespace SPH
distance_to_wall_[index_i] = r_wall_normal;
index_nearest[index_i] = index_j;
}
if (distance_to_wall_[index_i] <= 0.0 + TinyReal)
{
std::cout << "strange" << std::endl;
system("pause");
}
Vecd n_k_j_nearest = n_k[index_nearest[index_i]];
if (dimension_ == 2)
{
Expand All @@ -270,7 +276,7 @@ namespace SPH
{
is_near_wall_P1_[index_i] = 1;
}
if (r_wall_normal < (cutoff_radius_ - 0.5 * particle_spacing_) &&
if (r_wall_normal < (cutoff_radius_ - 0.5 * particle_spacing_)+TinyReal &&
r_wall_normal > 0.0 * particle_spacing_ + TinyReal)
{
is_near_wall_P2_[index_i] = 10;
Expand Down
26 changes: 18 additions & 8 deletions tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@ namespace SPH
particle_spacing_min_(inner_relation.real_body_->sph_adaptation_->MinimumSpacing()),
rho_(particles_->rho_), vel_(particles_->vel_),
mu_(DynamicCast<Fluid>(this, particles_->getBaseMaterial()).ReferenceViscosity()), dimension_(Vecd(0).size()),
smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()){}
smoothing_length_(sph_body_.sph_adaptation_->ReferenceSmoothingLength()) {}
//=================================================================================================//
K_TurtbulentModelInner::K_TurtbulentModelInner(BaseInnerRelation& inner_relation)
: BaseTurtbulentModelInner(inner_relation)
Expand All @@ -39,14 +39,12 @@ namespace SPH
particles_->registerSortableVariable<Real>("K_Production");
particles_->addVariableToWrite<Real>("K_Production");

//** for test */
particles_->registerVariable(lap_k_, "Lap_K");
particles_->registerSortableVariable<Real>("Lap_K");
particles_->addVariableToWrite<Real>("Lap_K");
particles_->registerVariable(B_, "CorrectionMatrix");

particles_->registerVariable(lap_k_term_, "Lap_K_term");
particles_->registerSortableVariable<Real>("Lap_K_term");
particles_->addVariableToWrite<Real>("Lap_K_term");
//** for test */
particles_->registerVariable(k_diffusion_, "K_Diffusion");
particles_->registerSortableVariable<Real>("K_Diffusion");
particles_->addVariableToWrite<Real>("K_Diffusion");

particles_->addVariableToWrite<Real>("ChangeRateOfTKE");

Expand Down Expand Up @@ -78,6 +76,18 @@ namespace SPH
{
particles_->registerVariable(dE_dt_, "ChangeRateOfTDR");
particles_->registerSortableVariable<Real>("ChangeRateOfTDR");
particles_->addVariableToWrite<Real>("ChangeRateOfTDR");

particles_->registerVariable(ep_production, "Ep_Production");
particles_->registerSortableVariable<Real>("Ep_Production");
particles_->addVariableToWrite<Real>("Ep_Production");
particles_->registerVariable(ep_dissipation_, "Ep_Dissipation_");
particles_->registerSortableVariable<Real>("Ep_Dissipation_");
particles_->addVariableToWrite<Real>("Ep_Dissipation_");
particles_->registerVariable(ep_diffusion_, "Ep_Diffusion_");
particles_->registerSortableVariable<Real>("Ep_Diffusion_");
particles_->addVariableToWrite<Real>("Ep_Diffusion_");

}
//=================================================================================================//
void E_TurtbulentModelInner::update(size_t index_i, Real dt)
Expand Down
7 changes: 3 additions & 4 deletions tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -75,7 +75,6 @@ namespace SPH
StdLargeVec<Real> turbu_mu_;
StdLargeVec<Real> turbu_k_;
StdLargeVec<Real> turbu_epsilon_;

Real smoothing_length_;
Real particle_spacing_min_;
Real mu_;
Expand All @@ -99,11 +98,11 @@ namespace SPH
protected:
StdLargeVec<Real> dk_dt_;
StdLargeVec<Matd> velocity_gradient;

StdLargeVec<Matd> B_;
StdLargeVec<Real> k_production_;

//** for test */
StdLargeVec<Real> lap_k_, lap_k_term_, vel_x_;
StdLargeVec<Real> k_diffusion_, vel_x_;
StdLargeVec<Matd> velocity_gradient_wall;
};

Expand All @@ -120,7 +119,7 @@ namespace SPH
inline void interaction(size_t index_i, Real dt = 0.0);
void update(size_t index_i, Real dt = 0.0);
protected:
StdLargeVec<Real> dE_dt_;
StdLargeVec<Real> dE_dt_, ep_production, ep_dissipation_, ep_diffusion_;
StdLargeVec<Real>& turbu_mu_;
StdLargeVec<Real>& turbu_k_;
StdLargeVec<Real>& turbu_epsilon_;
Expand Down
24 changes: 17 additions & 7 deletions tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,7 +53,9 @@ namespace SPH
for (size_t n = 0; n != inner_neighborhood.current_size_; ++n)
{
size_t index_j = inner_neighborhood.j_[n];
Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n];
//Vecd nablaW_ijV_j = 0.5* inner_neighborhood.dW_ijV_j_[n] * (B_[index_i] + B_[index_j]) * inner_neighborhood.e_ij_[n];
//Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * B_[index_i] * inner_neighborhood.e_ij_[n];
Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n];
velocity_gradient[index_i] += -(vel_i - vel_[index_j]) * nablaW_ijV_j.transpose();

k_derivative = (turbu_k_i - turbu_k_[index_j]) / (inner_neighborhood.r_ij_[n] + 0.01 * smoothing_length_);
Expand All @@ -68,8 +70,7 @@ namespace SPH
dk_dt_[index_i] = k_production_[index_i] - turbu_epsilon_[index_i] + k_lap;

//** for test */
lap_k_[index_i] = k_lap;
lap_k_term_[index_i] = (mu_ + turbu_mu_i / sigma_k) * lap_k_[index_i];
k_diffusion_[index_i] = k_lap;
}
//=================================================================================================//
void E_TurtbulentModelInner::
Expand Down Expand Up @@ -97,6 +98,12 @@ namespace SPH
epsilon_dissipation = C_2 * turbu_epsilon_i * turbu_epsilon_i / turbu_k_i;

dE_dt_[index_i] = epsilon_production - epsilon_dissipation + epsilon_lap;

//** for test */
ep_production[index_i] = epsilon_production;
ep_dissipation_[index_i] = epsilon_dissipation;
ep_diffusion_[index_i] = epsilon_lap;

}
//=================================================================================================//
void TKEnergyAccInner::interaction(size_t index_i, Real dt)
Expand All @@ -105,15 +112,17 @@ namespace SPH
Vecd acceleration = Vecd::Zero();
Vecd k_gradient = Vecd::Zero();

//strong form is used
const 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];
Vecd nablaW_ijV_j = inner_neighborhood.dW_ijV_j_[n] * inner_neighborhood.e_ij_[n];
k_gradient += -1.0*(turbu_k_i - turbu_k_[index_j]) * nablaW_ijV_j;
acceleration -= (2.0 / 3.0) * k_gradient;
//strong form
//k_gradient += -1.0*(turbu_k_i - turbu_k_[index_j]) * nablaW_ijV_j;
//weak form
k_gradient += 1.0 * (turbu_k_i + turbu_k_[index_j]) * nablaW_ijV_j;
}
acceleration = -1.0 * (2.0 / 3.0) * k_gradient;
acc_prior_[index_i] += acceleration;

//for test
Expand All @@ -128,14 +137,15 @@ namespace SPH
Vecd acceleration = Vecd::Zero();
Vecd vel_derivative = Vecd::Zero();
const 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];

// viscous force
vel_derivative = (vel_[index_i] - vel_[index_j]) / (inner_neighborhood.r_ij_[n] + 0.01 * smoothing_length_);
acceleration += 2.0 * (mu_+ turbu_mu_i) * vel_derivative * inner_neighborhood.dW_ijV_j_[n];
//acceleration += 2.0 * (mu_ ) * vel_derivative * inner_neighborhood.dW_ijV_j_[n];

}
//for test
visc_acc_inner_[index_i] = acceleration / rho_[index_i];
Expand Down

0 comments on commit 514e295

Please sign in to comment.