Merge pull request #391 from Xiangyu-Hu/fix/surafce_tension

Fix/surafce tension

src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_inner.cpp

@@ -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) {}
 //=================================================================================================//
@@ -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 ::

src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.cpp

@@ -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_));
     }
 }
 //=================================================================================================//

src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.h

@@ -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>;
@@ -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_;
 };
 
 /**

src/shared/particle_dynamics/fluid_dynamics/fluid_dynamics_multi_phase.hpp

@@ -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::
@@ -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_));
     }
 }
 //=================================================================================================//
@@ -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)

src/shared/particle_dynamics/solid_dynamics/contact_dynamics.cpp

@@ -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_));
     }
 }

src/shared/particle_dynamics/solid_dynamics/contact_dynamics.h

@@ -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)
@@ -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_;
 };

src/shared/particle_dynamics/solid_dynamics/fluid_structure_interaction.cpp

@@ -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());
     }

src/shared/particle_dynamics/solid_dynamics/fluid_structure_interaction.h

@@ -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)
             {
@@ -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_;
 };
@@ -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];
@@ -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)
@@ -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]));

tests/2d_examples/test_2d_square_droplet/src/droplet.cpp

@@ -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.
 //----------------------------------------------------------------------
@@ -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);

tests/2d_examples/test_2d_wetting_effects/src/wetting.cpp

@@ -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);
@@ -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();
 

tests/2d_examples/test_2d_wetting_effects/src/wetting.h

@@ -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;
 //----------------------------------------------------------------------

tests/3d_examples/test_3d_roof_parametric_cvt/test_3d_roof_parametric_cvt.cpp

@@ -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