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 f938bfb2a4..1c5a2f155f 100644 --- a/tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.cpp +++ b/tests/user_examples/2d_turbulent_channel/2d_turbulent_channel.cpp @@ -56,10 +56,10 @@ int main(int ac, char* av[]) //Attention! the original one does use Riemann solver for density Dynamics1Level density_relaxation(water_block_complex_relation); - /** Turbulent. */ - InteractionWithUpdate k_equation_relaxation(water_block_complex_relation); - InteractionWithUpdate epsilon_equation_relaxation(water_block_complex_relation); - InteractionDynamics turbulent_kinetic_energy_acceleration(water_block_complex_relation); + /** Turbulent.Note: When use wall function, K Epsilon and TKE gradient calculation only consider inner */ + InteractionWithUpdate k_equation_relaxation(water_block_inner); + InteractionWithUpdate epsilon_equation_relaxation(water_block_inner); + InteractionDynamics turbulent_kinetic_energy_acceleration(water_block_inner); SimpleDynamics wall_boundary_normal_direction(wall_boundary); /** Turbulent standard wall function needs normal vectors of wall. */ 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 3fde01338d..247fcbc923 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 @@ -51,8 +51,8 @@ namespace SPH Matd strain_rate = Matd::Zero(); Matd Re_stress = Matd::Zero(); //Matd velocity_gradient = Matd::Zero(); - velocity_gradient_wall[index_i] = Matd::Zero(); - + //velocity_gradient_wall[index_i] = Matd::Zero(); + /* for (size_t k = 0; k < FluidContactData::contact_configuration_.size(); ++k) { StdLargeVec& vel_ave_k = *(contact_vel_ave_[k]); @@ -68,18 +68,19 @@ namespace SPH velocity_gradient_wall[index_i] += -0.0 * (vel_i - vel_ave_k[index_j]) * nablaW_ijV_j.transpose(); - /** With standard wall function, diffusion of k to wall is zero */ + // With standard wall function, diffusion of k to wall is zero k_derivative = 0.0; k_lap += 0.0; } } + */ strain_rate = 0.5 * (velocity_gradient[index_i].transpose() + velocity_gradient[index_i]); Re_stress = 2.0 * strain_rate * turbu_mu_i / rho_i - (2.0 / 3.0) * turbu_k_i * Matd::Identity(); Matd k_production_matrix = Re_stress.array() * velocity_gradient[index_i].array(); k_production = k_production_matrix.sum(); /** With standard wall function, epilson on wall is zero */ - dk_dt_[index_i] += k_production - 0.0 + k_lap; + dk_dt_[index_i] += k_production - 0.0 + 0.0; /** Store the production of K for the use of Epsilon euqation and output*/ k_production_[index_i] = k_production; 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 6fc78a1d6c..72d99b2ff6 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 @@ -58,9 +58,6 @@ namespace SPH particles_->registerVariable(velocity_gradient_wall, "Velocity_Gradient_Wall"); particles_->registerSortableVariable("Velocity_Gradient_Wall"); particles_->addVariableToWrite("Velocity_Gradient_Wall"); - particles_->registerVariable(velocity_gradient_inner, "Velocity_Gradient_Inner"); - particles_->registerSortableVariable("Velocity_Gradient_Inner"); - particles_->addVariableToWrite("Velocity_Gradient_Inner"); particles_->registerVariable(vel_x_, "Velocity_X"); particles_->registerSortableVariable("Velocity_X"); 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 0b834d27ca..2d1fa66bf9 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 @@ -58,7 +58,6 @@ namespace SPH //closure coefficients for epsilon model Real C_l, C_2; Real sigma_E; - }; /** @@ -105,7 +104,7 @@ namespace SPH //** for test */ StdLargeVec lap_k_, lap_k_term_, vel_x_; - StdLargeVec velocity_gradient_inner, velocity_gradient_wall; + StdLargeVec velocity_gradient_wall; }; /** 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 8e1c492c66..fa13361ea5 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 @@ -35,23 +35,19 @@ namespace SPH namespace fluid_dynamics { //=================================================================================================// - void K_TurtbulentModelInner:: - interaction(size_t index_i, Real dt) + void K_TurtbulentModelInner::interaction(size_t index_i, Real dt) { - Vecd vel_i = vel_[index_i]; Real rho_i = rho_[index_i]; Real turbu_mu_i = turbu_mu_[index_i]; Real turbu_k_i = turbu_k_[index_i]; - dk_dt_[index_i] = 0.0; - Real k_production(0.0); + dk_dt_[index_i] = 0.0;; Real k_derivative(0.0); Real k_lap(0.0); - //Matd strain_rate = Matd::Zero(); - //Matd Re_stress = Matd::Zero(); + Matd strain_rate = Matd::Zero(); + Matd Re_stress = Matd::Zero(); velocity_gradient[index_i] = Matd::Zero(); - velocity_gradient_inner[index_i] = Matd::Zero(); const Neighborhood& inner_neighborhood = inner_configuration_[index_i]; for (size_t n = 0; n != inner_neighborhood.current_size_; ++n) @@ -63,22 +59,17 @@ namespace SPH k_derivative = (turbu_k_i - turbu_k_[index_j]) / (inner_neighborhood.r_ij_[n] + 0.01 * smoothing_length_); k_lap += 2.0 * (mu_+turbu_mu_i/sigma_k)*k_derivative * inner_neighborhood.dW_ijV_j_[n]/ rho_i; } - //for test - velocity_gradient_inner[index_i] = velocity_gradient[index_i]; - - //strain_rate = 0.5 * (velocity_gradient.transpose() + velocity_gradient); - //Re_stress = 2.0 * strain_rate * turbu_mu_i / rho_i - (2.0 / 3.0) * turbu_k_i * Matd::Identity(); - //Matd k_production_matrix = Re_stress.array() * velocity_gradient.array(); - //k_production = k_production_matrix.sum(); + strain_rate = 0.5 * (velocity_gradient[index_i].transpose() + velocity_gradient[index_i]); + Re_stress = 2.0 * strain_rate * turbu_mu_i / rho_i - (2.0 / 3.0) * turbu_k_i * Matd::Identity(); + Matd k_production_matrix = Re_stress.array() * velocity_gradient[index_i].array(); + k_production_[index_i] = k_production_matrix.sum(); - /** k_production will be calculated in the complex part after getting the velocity contribution of wall*/ - dk_dt_[index_i] = 0.0 - turbu_epsilon_[index_i] + k_lap; + dk_dt_[index_i] = k_production_[index_i] - turbu_epsilon_[index_i] + k_lap; //** for test */ - lap_k_[index_i] = 0.0; lap_k_[index_i] = k_lap; - + lap_k_term_[index_i] = (mu_ + turbu_mu_i / sigma_k) * lap_k_[index_i]; } //=================================================================================================// void E_TurtbulentModelInner:: @@ -101,8 +92,8 @@ namespace SPH epsilon_derivative = (turbu_epsilon_i - turbu_epsilon_[index_j]) / (inner_neighborhood.r_ij_[n] + 0.01 * smoothing_length_); epsilon_lap += 2.0 * (mu_ + turbu_mu_i / sigma_E) * epsilon_derivative * inner_neighborhood.dW_ijV_j_[n] / rho_i; } - /** epsilon_production will be calculated in the complex part*/ - epsilon_production = 0.0; + + epsilon_production = C_l * turbu_epsilon_i * k_production_[index_i] / turbu_k_i; epsilon_dissipation = C_2 * turbu_epsilon_i * turbu_epsilon_i / turbu_k_i; dE_dt_[index_i] = epsilon_production - epsilon_dissipation + epsilon_lap; @@ -121,11 +112,10 @@ namespace SPH 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; - //k_gradient += (turbu_k_i - turbu_k_[index_j]) * nablaW_ijV_j; acceleration -= (2.0 / 3.0) * k_gradient; } - //if (surface_indicator_[index_i] == 0 && pos_[index_i][0] <= 5.95)//To prevent kernel truncation near outlet acc_prior_[index_i] += acceleration; + //for test tke_acc_inner_[index_i] = acceleration; test_k_grad_rslt_[index_i] = k_gradient;