From 514e2953b0308f7e977198c50d4d87b077e96ddc Mon Sep 17 00:00:00 2001 From: wangfeng_real Date: Sun, 23 Jul 2023 10:57:45 +0200 Subject: [PATCH] recent changes --- .../test_turbulent_module.cpp | 33 ++++++++++++++----- .../2d_turbulent_channel.cpp | 16 ++++++--- .../2d_turbulent_channel.h | 2 +- .../k-epsilon_turbulent_model_complex.hpp | 30 ++++++++++------- .../k-epsilon_turbulent_model_inner.cpp | 26 ++++++++++----- .../k-epsilon_turbulent_model_inner.h | 7 ++-- .../k-epsilon_turbulent_model_inner.hpp | 24 ++++++++++---- 7 files changed, 92 insertions(+), 46 deletions(-) diff --git a/tests/unit_tests_src/for_2D_build/test_turbulent_module/test_turbulent_module.cpp b/tests/unit_tests_src/for_2D_build/test_turbulent_module/test_turbulent_module.cpp index 16fa7b27b9..c37eafc66d 100644 --- a/tests/unit_tests_src/for_2D_build/test_turbulent_module/test_turbulent_module.cpp +++ b/tests/unit_tests_src/for_2D_build/test_turbulent_module/test_turbulent_module.cpp @@ -81,7 +81,7 @@ class BaseTurbulentModule wall_boundary(sph_system, makeShared("Wall")) { std::cout << "constructor for BaseTurbulentModule" << std::endl; - water_block.defineParticlesAndMaterial(rho0_f, c_f, mu_f); + water_block.defineParticlesAndMaterial(rho0_f, c_f, mu_f); water_block.generateParticles(); wall_boundary.defineParticlesAndMaterial(); wall_boundary.generateParticles(); @@ -95,14 +95,16 @@ class TurbulentModule : public :: testing :: Test, public BaseTurbulentModule InnerRelation water_block_inner; ComplexRelation water_block_complex_relation; - InteractionWithUpdate k_equation_relaxation; - InteractionWithUpdate epsilon_equation_relaxation; - InteractionDynamics turbulent_kinetic_energy_acceleration; + InteractionWithUpdate k_equation_relaxation; + InteractionWithUpdate epsilon_equation_relaxation; + InteractionDynamics turbulent_kinetic_energy_acceleration; SimpleDynamics wall_boundary_normal_direction; InteractionDynamics standard_wall_function_correction; InteractionDynamics turbulent_viscous_acceleration; SimpleDynamics update_eddy_viscosity; + InteractionWithUpdate correct_configuration; + InteractionWithUpdate update_density_by_summation; BodyStatesRecordingToVtp write_body_states; size_t number_of_iterations; @@ -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), @@ -202,7 +206,7 @@ class Test_K_Epsilon_Equation : public fluid_dynamics::FluidInitialCondition std::vector 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; @@ -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) { @@ -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); diff --git a/tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.cpp b/tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.cpp index 6b204f0799..c15fc1c7eb 100644 --- a/tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.cpp +++ b/tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.cpp @@ -66,11 +66,11 @@ int main(int ac, char* av[]) InteractionDynamics standard_wall_function_correction(water_block_complex_relation); /** TurbulentViscous cal. uses friction velocity and Y+ that are defined in WallFunction . */ - InteractionDynamics turbulent_viscous_acceleration(water_block_complex_relation); + InteractionDynamics turbulent_viscous_acceleration(water_block_complex_relation); //InteractionDynamics viscous_acceleration(water_block_complex_relation); - InteractionDynamics transport_velocity_correction(water_block_complex_relation); + InteractionDynamics transport_velocity_correction(water_block_complex_relation); InteractionWithUpdate inlet_outlet_surface_particle_indicator(water_block_complex_relation); InteractionWithUpdate update_density_by_summation(water_block_complex_relation); @@ -91,7 +91,13 @@ int main(int ac, char* av[]) /** Turbulent eddy viscosity calculation needs values of Wall Y start. */ SimpleDynamics update_eddy_viscosity(water_block); + + /** Try to introduce B correction */ + InteractionWithUpdate correct_configuration(water_block_inner); + + + BodyAlignedBoxByParticle emitter( water_block, makeShared(Transform(Vec2d(emitter_translation)), emitter_halfsize)); SimpleDynamics emitter_inflow_injection(emitter, 50, 0); @@ -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); @@ -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 = " diff --git a/tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.h b/tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.h index 31f3b4d149..214d0a2cfc 100644 --- a/tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.h +++ b/tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.h @@ -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. */ diff --git a/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_complex.hpp b/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_complex.hpp index 37910b2c75..fe95a13132 100644 --- a/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_complex.hpp +++ b/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_complex.hpp @@ -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:: @@ -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) @@ -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) { @@ -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; diff --git a/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.cpp b/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.cpp index 72d99b2ff6..b91b7f0464 100644 --- a/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.cpp +++ b/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.cpp @@ -15,7 +15,7 @@ namespace SPH particle_spacing_min_(inner_relation.real_body_->sph_adaptation_->MinimumSpacing()), rho_(particles_->rho_), vel_(particles_->vel_), mu_(DynamicCast(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) @@ -39,14 +39,12 @@ namespace SPH particles_->registerSortableVariable("K_Production"); particles_->addVariableToWrite("K_Production"); - //** for test */ - particles_->registerVariable(lap_k_, "Lap_K"); - particles_->registerSortableVariable("Lap_K"); - particles_->addVariableToWrite("Lap_K"); + particles_->registerVariable(B_, "CorrectionMatrix"); - particles_->registerVariable(lap_k_term_, "Lap_K_term"); - particles_->registerSortableVariable("Lap_K_term"); - particles_->addVariableToWrite("Lap_K_term"); + //** for test */ + particles_->registerVariable(k_diffusion_, "K_Diffusion"); + particles_->registerSortableVariable("K_Diffusion"); + particles_->addVariableToWrite("K_Diffusion"); particles_->addVariableToWrite("ChangeRateOfTKE"); @@ -78,6 +76,18 @@ namespace SPH { particles_->registerVariable(dE_dt_, "ChangeRateOfTDR"); particles_->registerSortableVariable("ChangeRateOfTDR"); + particles_->addVariableToWrite("ChangeRateOfTDR"); + + particles_->registerVariable(ep_production, "Ep_Production"); + particles_->registerSortableVariable("Ep_Production"); + particles_->addVariableToWrite("Ep_Production"); + particles_->registerVariable(ep_dissipation_, "Ep_Dissipation_"); + particles_->registerSortableVariable("Ep_Dissipation_"); + particles_->addVariableToWrite("Ep_Dissipation_"); + particles_->registerVariable(ep_diffusion_, "Ep_Diffusion_"); + particles_->registerSortableVariable("Ep_Diffusion_"); + particles_->addVariableToWrite("Ep_Diffusion_"); + } //=================================================================================================// void E_TurtbulentModelInner::update(size_t index_i, Real dt) diff --git a/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.h b/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.h index 2d1fa66bf9..12eb01e1dc 100644 --- a/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.h +++ b/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.h @@ -75,7 +75,6 @@ namespace SPH StdLargeVec turbu_mu_; StdLargeVec turbu_k_; StdLargeVec turbu_epsilon_; - Real smoothing_length_; Real particle_spacing_min_; Real mu_; @@ -99,11 +98,11 @@ namespace SPH protected: StdLargeVec dk_dt_; StdLargeVec velocity_gradient; - + StdLargeVec B_; StdLargeVec k_production_; //** for test */ - StdLargeVec lap_k_, lap_k_term_, vel_x_; + StdLargeVec k_diffusion_, vel_x_; StdLargeVec velocity_gradient_wall; }; @@ -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 dE_dt_; + StdLargeVec dE_dt_, ep_production, ep_dissipation_, ep_diffusion_; StdLargeVec& turbu_mu_; StdLargeVec& turbu_k_; StdLargeVec& turbu_epsilon_; diff --git a/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.hpp b/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.hpp index fa13361ea5..e88e90824a 100644 --- a/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.hpp +++ b/tests/user_examples/extra_src/for_2D_build/k-epsilon_turbulent_model_inner.hpp @@ -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_); @@ -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:: @@ -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) @@ -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 @@ -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];