Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix/surafce tension #391

Merged
merged 4 commits into from
Jul 30, 2023
Merged

Fix/surafce tension #391

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
6f3d5a0
fix bugs, may be still other bugs, please test it.
Xiangyu-Hu Jul 26, 2023
12f0593
found a bug in multiphase viscous acceleration
Xiangyu-Hu Jul 29, 2023
0100028
revise wetting case too
Xiangyu-Hu Jul 29, 2023
5c15698
naming of velocity
Xiangyu-Hu Jul 29, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 5 additions & 3 deletions src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,8 +59,8 @@ TransportVelocityCorrectionInnerAdaptive::
//=================================================================================================//
AcousticTimeStepSize::AcousticTimeStepSize(SPHBody &sph_body, Real acousticCFL)
: LocalDynamicsReduce<Real, ReduceMax>(sph_body, Real(0)),
FluidDataSimple(sph_body), fluid_(DynamicCast<Fluid>(this, particles_->getBaseMaterial())), rho_(particles_->rho_),
p_(*particles_->getVariableByName<Real>("Pressure")), vel_(particles_->vel_),
FluidDataSimple(sph_body), fluid_(DynamicCast<Fluid>(this, particles_->getBaseMaterial())),
rho_(particles_->rho_), p_(*particles_->getVariableByName<Real>("Pressure")), vel_(particles_->vel_),
smoothing_length_min_(sph_body.sph_adaptation_->MinimumSmoothingLength()),
acousticCFL_(acousticCFL) {}
//=================================================================================================//
Expand Down Expand Up @@ -118,7 +118,9 @@ VorticityInner::VorticityInner(BaseInnerRelation &inner_relation)
BaseIntegration::BaseIntegration(BaseInnerRelation &inner_relation)
: LocalDynamics(inner_relation.getSPHBody()), FluidDataInner(inner_relation),
fluid_(DynamicCast<Fluid>(this, particles_->getBaseMaterial())), rho_(particles_->rho_),
p_(*particles_->getVariableByName<Real>("Pressure")), drho_dt_(*particles_->registerSharedVariable<Real>("DensityChangeRate")), pos_(particles_->pos_), vel_(particles_->vel_),
p_(*particles_->getVariableByName<Real>("Pressure")),
drho_dt_(*particles_->registerSharedVariable<Real>("DensityChangeRate")),
pos_(particles_->pos_), vel_(particles_->vel_),
acc_(particles_->acc_), acc_prior_(particles_->acc_prior_) {}
//=================================================================================================//
Oldroyd_BIntegration1stHalf ::
Expand Down
5 changes: 3 additions & 2 deletions src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,9 @@ ViscousAccelerationMultiPhase::ViscousAccelerationMultiPhase(BaseInnerRelation &

for (size_t k = 0; k != contact_particles_.size(); ++k)
{
contact_fluids_.push_back(DynamicCast<Fluid>(this, &contact_particles_[k]->getBaseMaterial()));
contact_vel_n_.push_back(&(contact_particles_[k]->vel_));
Real mu_k = DynamicCast<Fluid>(this, &contact_particles_[k]->getBaseMaterial())->ReferenceViscosity();
contact_mu_.push_back(Real(2) * (mu_ * mu_k) / (mu_ + mu_k));
contact_vel_.push_back(&(contact_particles_[k]->vel_));
}
}
//=================================================================================================//
Expand Down
6 changes: 3 additions & 3 deletions src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,8 +54,8 @@ class ViscousAccelerationMultiPhase : public ViscousAccelerationInner, public Mu
inline void interaction(size_t index_i, Real dt = 0.0);

protected:
StdVec<Fluid *> contact_fluids_;
StdVec<StdLargeVec<Vecd> *> contact_vel_n_;
StdVec<Real> contact_mu_;
StdVec<StdLargeVec<Vecd> *> contact_vel_;
};
using ViscousAccelerationMultiPhaseWithWall =
BaseViscousAccelerationWithWall<ViscousAccelerationMultiPhase>;
Expand All @@ -75,7 +75,7 @@ class RelaxationMultiPhase : public RelaxationInnerType, public MultiPhaseContac
protected:
StdVec<Fluid *> contact_fluids_;
StdVec<StdLargeVec<Real> *> contact_p_, contact_rho_n_;
StdVec<StdLargeVec<Vecd> *> contact_vel_n_;
StdVec<StdLargeVec<Vecd> *> contact_vel_;
};

/**
Expand Down
28 changes: 9 additions & 19 deletions src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,35 +14,25 @@ namespace SPH
namespace fluid_dynamics
{
//=================================================================================================//
void ViscousAccelerationMultiPhase::
interaction(size_t index_i, Real dt)
void ViscousAccelerationMultiPhase::interaction(size_t index_i, Real dt)
{
ViscousAccelerationInner::interaction(index_i, dt);

Real rho_i = this->rho_[index_i];
const Vecd &vel_i = this->vel_[index_i];

Vecd acceleration = Vecd::Zero();
Vecd vel_derivative = Vecd::Zero();
for (size_t k = 0; k < this->contact_configuration_.size(); ++k)
{
Real mu_j = this->contact_fluids_[k]->ReferenceViscosity();
StdLargeVec<Vecd> &vel_k = *(this->contact_vel_n_[k]);
Real contact_mu_k = this->contact_mu_[k];
StdLargeVec<Vecd> &vel_k = *(this->contact_vel_[k]);
Neighborhood &contact_neighborhood = (*this->contact_configuration_[k])[index_i];
for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
{
size_t index_j = contact_neighborhood.j_[n];
Real r_ij = contact_neighborhood.r_ij_[n];

vel_derivative = 2.0 * (vel_i - vel_k[index_j]) /
(r_ij + 0.01 * this->smoothing_length_);
Real mu_ij = 2.0 * this->mu_ * mu_j / (this->mu_ + mu_j);
acceleration += 2.0 * mu_ij * vel_derivative *
contact_neighborhood.dW_ijV_j_[n] / rho_i;
Vecd vel_derivative = (this->vel_[index_i] - vel_k[index_j]) /
(contact_neighborhood.r_ij_[n] + 0.01 * this->smoothing_length_);
acceleration += 2.0 * contact_mu_k * vel_derivative * contact_neighborhood.dW_ijV_j_[n];
}
}

acc_prior_[index_i] += acceleration;
acc_prior_[index_i] += acceleration / this->rho_[index_i];
}
//=================================================================================================//
void MultiPhaseColorFunctionGradient::
Expand Down Expand Up @@ -90,7 +80,7 @@ RelaxationMultiPhase<RelaxationInnerType>::
contact_fluids_.push_back(DynamicCast<Fluid>(this, &contact_particles_[k]->getBaseMaterial()));
contact_p_.push_back(contact_particles_[k]->template getVariableByName<Real>("Pressure"));
contact_rho_n_.push_back(&(contact_particles_[k]->rho_));
contact_vel_n_.push_back(&(contact_particles_[k]->vel_));
contact_vel_.push_back(&(contact_particles_[k]->vel_));
}
}
//=================================================================================================//
Expand Down Expand Up @@ -190,7 +180,7 @@ void BaseMultiPhaseIntegration2ndHalf<Integration2ndHalfType>::
Vecd p_dissipation = Vecd::Zero();
for (size_t k = 0; k < this->contact_configuration_.size(); ++k)
{
StdLargeVec<Vecd> &vel_k = *(this->contact_vel_n_[k]);
StdLargeVec<Vecd> &vel_k = *(this->contact_vel_[k]);
CurrentRiemannSolver &riemann_solver_k = riemann_solvers_[k];
Neighborhood &contact_neighborhood = (*this->contact_configuration_[k])[index_i];
for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
Expand Down
2 changes: 1 addition & 1 deletion src/shared/particle_dynamics/solid_dynamics/contact_dynamics.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -147,7 +147,7 @@ DynamicContactForceWithWall::
reference_pressure_ = solid_.ReferenceDensity() * solid_.ContactStiffness();
for (size_t k = 0; k != contact_particles_.size(); ++k)
{
contact_vel_n_.push_back(&(contact_particles_[k]->vel_));
contact_vel_.push_back(&(contact_particles_[k]->vel_));
contact_n_.push_back(&(contact_particles_[k]->n_));
}
}
Expand Down
4 changes: 2 additions & 2 deletions src/shared/particle_dynamics/solid_dynamics/contact_dynamics.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -409,7 +409,7 @@ class DynamicContactForceWithWall : public LocalDynamics, public ContactDynamics
particle_spacing_ratio2 *= 0.1 * particle_spacing_ratio2;

StdLargeVec<Vecd> &n_k = *(contact_n_[k]);
StdLargeVec<Vecd> &vel_n_k = *(contact_vel_n_[k]);
StdLargeVec<Vecd> &vel_n_k = *(contact_vel_[k]);

Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
Expand Down Expand Up @@ -437,7 +437,7 @@ class DynamicContactForceWithWall : public LocalDynamics, public ContactDynamics
Solid &solid_;
StdLargeVec<Real> &Vol_, &mass_;
StdLargeVec<Vecd> &vel_, &acc_prior_;
StdVec<StdLargeVec<Vecd> *> contact_vel_n_, contact_n_;
StdVec<StdLargeVec<Vecd> *> contact_vel_, contact_n_;
Real penalty_strength_;
Real impedance_, reference_pressure_;
};
Expand Down
2 changes: 1 addition & 1 deletion src/shared/particle_dynamics/solid_dynamics/fluid_structure_interaction.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ ViscousForceFromFluid::ViscousForceFromFluid(BaseContactRelation &contact_relati
particles_->registerVariable(force_from_fluid_, "ViscousForceFromFluid");
for (size_t k = 0; k != contact_particles_.size(); ++k)
{
contact_vel_n_.push_back(&(contact_particles_[k]->vel_));
contact_vel_.push_back(&(contact_particles_[k]->vel_));
mu_.push_back(contact_fluids_[k]->ReferenceViscosity());
smoothing_length_.push_back(contact_bodies_[k]->sph_adaptation_->ReferenceSmoothingLength());
}
Expand Down
10 changes: 5 additions & 5 deletions src/shared/particle_dynamics/solid_dynamics/fluid_structure_interaction.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -80,7 +80,7 @@ class ViscousForceFromFluid : public BaseForceFromFluid
{
Real mu_k = mu_[k];
Real smoothing_length_k = smoothing_length_[k];
StdLargeVec<Vecd> &vel_n_k = *(contact_vel_n_[k]);
StdLargeVec<Vecd> &vel_n_k = *(contact_vel_[k]);
Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
for (size_t n = 0; n != contact_neighborhood.current_size_; ++n)
{
Expand All @@ -98,7 +98,7 @@ class ViscousForceFromFluid : public BaseForceFromFluid

protected:
StdLargeVec<Vecd> &vel_ave_;
StdVec<StdLargeVec<Vecd> *> contact_vel_n_;
StdVec<StdLargeVec<Vecd> *> contact_vel_;
StdVec<Real> mu_;
StdVec<Real> smoothing_length_;
};
Expand Down Expand Up @@ -128,7 +128,7 @@ class BasePressureForceAccelerationFromFluid : public BaseForceFromFluid
{
StdLargeVec<Real> &rho_n_k = *(contact_rho_n_[k]);
StdLargeVec<Real> &p_k = *(contact_p_[k]);
StdLargeVec<Vecd> &vel_k = *(contact_vel_n_[k]);
StdLargeVec<Vecd> &vel_k = *(contact_vel_[k]);
StdLargeVec<Vecd> &acc_prior_k = *(contact_acc_prior_[k]);
RiemannSolverType &riemann_solvers_k = riemann_solvers_[k];
Neighborhood &contact_neighborhood = (*contact_configuration_[k])[index_i];
Expand All @@ -150,7 +150,7 @@ class BasePressureForceAccelerationFromFluid : public BaseForceFromFluid
protected:
StdLargeVec<Vecd> &vel_ave_, &acc_prior_, &acc_ave_, &n_;
StdVec<StdLargeVec<Real> *> contact_rho_n_, contact_p_;
StdVec<StdLargeVec<Vecd> *> contact_vel_n_, contact_acc_prior_;
StdVec<StdLargeVec<Vecd> *> contact_vel_, contact_acc_prior_;
StdVec<RiemannSolverType> riemann_solvers_;

BasePressureForceAccelerationFromFluid(bool mostDerived, BaseContactRelation &contact_relation)
Expand All @@ -162,7 +162,7 @@ class BasePressureForceAccelerationFromFluid : public BaseForceFromFluid
for (size_t k = 0; k != contact_particles_.size(); ++k)
{
contact_rho_n_.push_back(&(contact_particles_[k]->rho_));
contact_vel_n_.push_back(&(contact_particles_[k]->vel_));
contact_vel_.push_back(&(contact_particles_[k]->vel_));
contact_p_.push_back(contact_particles_[k]->template getVariableByName<Real>("Pressure"));
contact_acc_prior_.push_back(&(contact_particles_[k]->acc_prior_));
riemann_solvers_.push_back(RiemannSolverType(*contact_fluids_[k], *contact_fluids_[k]));
Expand Down
16 changes: 8 additions & 8 deletions tests/2d_examples/test_2d_square_droplet/src/droplet.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,7 @@ Real rho0_a = 0.001; /**< Reference density of air. */
Real U_max = 1.0; /**< Characteristic velocity. */
Real c_f = 10.0 * U_max; /**< Reference sound speed. */
Real mu_f = 0.2; /**< Water viscosity. */
Real mu_a = 0.0002; /**< Air viscosity. */
Real mu_a = 0.002; /**< Air viscosity. */
//----------------------------------------------------------------------
// Geometric shapes used in this case.
//----------------------------------------------------------------------
Expand Down Expand Up @@ -142,23 +142,23 @@ int main()
SimpleDynamics<TimeStepInitialization> initialize_a_air_step(air_block);
SimpleDynamics<NormalDirectionFromBodyShape> wall_boundary_normal_direction(wall_boundary);
/** Evaluation of density by summation approach. */
InteractionWithUpdate<fluid_dynamics::DensitySummationFreeSurfaceComplex>
update_water_density_by_summation(water_wall_contact, water_air_complex.getInnerRelation());
InteractionWithUpdate<fluid_dynamics::DensitySummationComplex>
update_water_density_by_summation(water_wall_contact, water_air_complex);
InteractionWithUpdate<fluid_dynamics::DensitySummationComplex>
update_air_density_by_summation(air_wall_contact, air_water_complex);
InteractionDynamics<fluid_dynamics::TransportVelocityCorrectionComplex>
air_transport_correction(air_wall_contact, air_water_complex);
air_transport_correction(air_wall_contact, air_water_complex, 0.05);
/** Time step size without considering sound wave speed. */
ReduceDynamics<fluid_dynamics::AdvectionTimeStepSize> get_water_advection_time_step_size(water_block, U_max);
ReduceDynamics<fluid_dynamics::AdvectionTimeStepSize> get_air_advection_time_step_size(air_block, U_max);
/** Time step size with considering sound wave speed. */
ReduceDynamics<fluid_dynamics::AcousticTimeStepSize> get_water_time_step_size(water_block);
ReduceDynamics<fluid_dynamics::AcousticTimeStepSize> get_air_time_step_size(air_block);
/** Pressure relaxation for water by using position verlet time stepping. */
Dynamics1Level<fluid_dynamics::Integration1stHalfRiemannWithWall>
water_pressure_relaxation(water_wall_contact, water_air_complex.getInnerRelation());
Dynamics1Level<fluid_dynamics::Integration2ndHalfRiemannWithWall>
water_density_relaxation(water_wall_contact, water_air_complex.getInnerRelation());
Dynamics1Level<fluid_dynamics::MultiPhaseIntegration1stHalfRiemannWithWall>
water_pressure_relaxation(water_wall_contact, water_air_complex);
Dynamics1Level<fluid_dynamics::MultiPhaseIntegration2ndHalfRiemannWithWall>
water_density_relaxation(water_wall_contact, water_air_complex);
/** Extend Pressure relaxation is used for air. */
Dynamics1Level<fluid_dynamics::ExtendMultiPhaseIntegration1stHalfRiemannWithWall>
air_pressure_relaxation(air_wall_contact, air_water_complex, 2.0);
Expand Down
15 changes: 7 additions & 8 deletions tests/2d_examples/test_2d_wetting_effects/src/wetting.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -49,24 +49,24 @@ int main()
SimpleDynamics<TimeStepInitialization> initialize_a_water_step(water_block);
SimpleDynamics<TimeStepInitialization> initialize_a_air_step(air_block);
/** Evaluation of density by summation approach. */
InteractionWithUpdate<fluid_dynamics::DensitySummationFreeSurfaceComplex>
update_water_density_by_summation(water_wall_contact, water_air_complex.getInnerRelation());
InteractionWithUpdate<fluid_dynamics::DensitySummationComplex>
update_water_density_by_summation(water_wall_contact, water_air_complex);
InteractionWithUpdate<fluid_dynamics::DensitySummationComplex>
update_air_density_by_summation(air_wall_contact, air_water_complex);
/** transport formulation for regularizing particle distribution. */
InteractionDynamics<fluid_dynamics::TransportVelocityCorrectionComplex>
air_transport_correction(air_wall_contact, air_water_complex);
air_transport_correction(air_wall_contact, air_water_complex, 0.05);
/** Time step size without considering sound wave speed. */
ReduceDynamics<fluid_dynamics::AdvectionTimeStepSize> get_water_advection_time_step_size(water_block, U_max);
ReduceDynamics<fluid_dynamics::AdvectionTimeStepSize> get_air_advection_time_step_size(air_block, U_max);
/** Time step size with considering sound wave speed. */
ReduceDynamics<fluid_dynamics::AcousticTimeStepSize> get_water_time_step_size(water_block);
ReduceDynamics<fluid_dynamics::AcousticTimeStepSize> get_air_time_step_size(air_block);
/** Pressure relaxation for water by using position verlet time stepping. */
Dynamics1Level<fluid_dynamics::Integration1stHalfRiemannWithWall>
water_pressure_relaxation(water_wall_contact, water_air_complex.getInnerRelation());
Dynamics1Level<fluid_dynamics::Integration2ndHalfRiemannWithWall>
water_density_relaxation(water_wall_contact, water_air_complex.getInnerRelation());
Dynamics1Level<fluid_dynamics::MultiPhaseIntegration1stHalfRiemannWithWall>
water_pressure_relaxation(water_wall_contact, water_air_complex);
Dynamics1Level<fluid_dynamics::MultiPhaseIntegration2ndHalfRiemannWithWall>
water_density_relaxation(water_wall_contact, water_air_complex);
/** Extend Pressure relaxation is used for air. */
Dynamics1Level<fluid_dynamics::ExtendMultiPhaseIntegration1stHalfRiemannWithWall>
air_pressure_relaxation(air_wall_contact, air_water_complex, 2.0);
Expand Down Expand Up @@ -135,7 +135,6 @@ int main()
update_water_density_by_summation.exec();
update_air_density_by_summation.exec();
air_transport_correction.exec();

air_viscous_acceleration.exec();
water_viscous_acceleration.exec();

Expand Down
2 changes: 1 addition & 1 deletion tests/2d_examples/test_2d_wetting_effects/src/wetting.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ Real rho0_a = 1.0e-3; /**< Reference density of air.
Real U_max = 1.0; /**< Characteristic velocity. */
Real c_f = 10.0 * U_max; /**< Reference sound speed. */
Real mu_f = 5.0e-2; /**< Water viscosity. */
Real mu_a = 5.0e-5; /**< Air viscosity. */
Real mu_a = 5.0e-4; /**< Air viscosity. */
Real contact_angle = (150.0 / 180.0) * 3.1415926; /**< Contact angle with Wall. */
Real tension_force = 0.008;
//----------------------------------------------------------------------
Expand Down
2 changes: 1 addition & 1 deletion tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -270,7 +270,7 @@ return_data roof_under_self_weight(Real dp, bool cvt = true, int particle_number
observer_point_shell point_A;
observer_point_shell point_B;
point_A.pos_0 = Vec3d(radius * std::sin(theta_radian), radius * std::cos(theta_radian) - radius, 0);
point_B.pos_0 = Vec3d(0);
point_B.pos_0 = Vec3d::Zero();
// resolution
const int dp_cm = dp * 100;
Real total_area = length * 2 * arc; // accounting for particles being on the edges
Expand Down
Loading