Model for dam breach - erosion,deposition,entrainment

restart/Restart_time_0000001000.dat

@@ -0,0 +1 @@
+2.417521457   

restart/Restart_time_0000002000.dat

@@ -0,0 +1 @@
+4.783250666   

restart/Restart_time_0000003000.dat

@@ -0,0 +1 @@
+7.132537260   

restart/Restart_time_0000004000.dat

@@ -0,0 +1 @@
+9.475482172   

restart/Restart_time_0000005000.dat

@@ -0,0 +1 @@
+11.810872703   

restart/Restart_time_0000006000.dat

@@ -0,0 +1 @@
+14.126971049   

restart/Restart_time_0000007000.dat

@@ -0,0 +1 @@
+16.459614298   

restart/Restart_time_0000008000.dat

@@ -0,0 +1 @@
+18.799581662   

restart/Restart_time_0000009000.dat

@@ -0,0 +1 @@
+21.017852872   

src/shared/common/sphinxsys_containers.h

@@ -75,6 +75,8 @@ class Contact; /**< Contact interaction: interaction between a body with one or
 
 class Boundary; /**< Interaction with boundary */
 class Wall;     /**< Interaction with wall boundary */
+class Soil;     /**< Interaction with soil*/
+class Fluid;     /**< Interaction with fluid*/
 class Extended; /**< An extened method of an interaction type */
 
 template <typename...>

src/shared/physical_closure/materials/general_continuum.h

@@ -71,22 +71,28 @@ class GeneralContinuum : public WeaklyCompressibleFluid
       public:
         GeneralContinuumKernel(GeneralContinuum &encloser) : E_(encloser.E_), G_(encloser.G_), K_(encloser.K_),
                                                              nu_(encloser.nu_), contact_stiffness_(encloser.contact_stiffness_),
-                                                             rho0_(encloser.rho0_) {};
+                                                             rho0_(encloser.rho0_),p0_(encloser.p0_) {};
 
         inline Real getYoungsModulus() { return E_; };
         inline Real getPoissonRatio() { return nu_; };
         inline Real getDensity() { return rho0_; };
         inline Real getBulkModulus(Real youngs_modulus, Real poisson_ratio);
         inline Real getShearModulus(Real youngs_modulus, Real poisson_ratio);
         inline Real getLambda(Real youngs_modulus, Real poisson_ratio);
+        Real getPressure(Real rho)
+        {
+            return p0_ * (rho / rho0_ - 1.0);
+        };
 
       protected:
         Real E_;                 /* Youngs or tensile modules  */
         Real G_;                 /* shear modules  */
         Real K_;                 /* bulk modules  */
         Real nu_;                /* Poisson ratio  */
         Real contact_stiffness_; /* contact-force stiffness related to bulk modulus*/
-        Real rho0_;              /* contact-force stiffness related to bulk modulus*/
+        Real rho0_;              /* reference density*/
+        Real p0_;                /* reference Pressure  */
+
     };
 };
 
@@ -98,11 +104,12 @@ class PlasticContinuum : public GeneralContinuum
     Real psi_;                          /* dilatancy angle  */
     Real alpha_phi_;                    /* Drucker-Prager's constants */
     Real k_c_;                          /* Drucker-Prager's constants */
+    Real d_s_;                          /* Mean particle diameter */
     const Real stress_dimension_ = 3.0; /* plain strain condition */
   public:
     explicit PlasticContinuum(Real rho0, Real c0, Real youngs_modulus, Real poisson_ratio, Real friction_angle, Real cohesion = 0, Real dilatancy = 0)
         : GeneralContinuum(rho0, c0, youngs_modulus, poisson_ratio),
-          c_(cohesion), phi_(friction_angle), psi_(dilatancy), alpha_phi_(0.0), k_c_(0.0)
+          c_(cohesion), phi_(friction_angle), psi_(dilatancy), d_s_(0.002), alpha_phi_(0.0), k_c_(0.0)
     {
         material_type_name_ = "PlasticContinuum";
         alpha_phi_ = getDPConstantsA(friction_angle);
@@ -113,6 +120,7 @@ class PlasticContinuum : public GeneralContinuum
     Real getDPConstantsA(Real friction_angle);
     Real getDPConstantsK(Real cohesion, Real friction_angle);
     Real getFrictionAngle() { return phi_; };
+    Real getCohesion() { return c_; };
 
     virtual Mat3d ConstitutiveRelation(Mat3d &velocity_gradient, Mat3d &stress_tensor);
     virtual Mat3d ReturnMapping(Mat3d &stress_tensor);
@@ -122,19 +130,35 @@ class PlasticContinuum : public GeneralContinuum
       public:
         PlasticKernel(PlasticContinuum &encloser) : GeneralContinuum::GeneralContinuumKernel(encloser),
                                                     c_(encloser.c_), phi_(encloser.phi_),
-                                                    psi_(encloser.psi_), alpha_phi_(encloser.alpha_phi_), k_c_(encloser.k_c_) {};
-
+                                                    psi_(encloser.psi_), alpha_phi_(encloser.alpha_phi_), k_c_(encloser.k_c_),
+                                                    d_s_(encloser.d_s_) {};
+        struct ReturnMappingResult
+        {
+            Mat3d stress_tensor;
+            bool has_yielded;
+        };
         inline Real getDPConstantsA(Real friction_angle);
+        inline Real getDPConstantsK(Real cohesion, Real friction_angle);
         inline Mat3d ConstitutiveRelation(Mat3d &velocity_gradient, Mat3d &stress_tensor);
+        inline Mat3d ConstitutiveRelationWithReduction(Mat3d &velocity_gradient, Mat3d &stress_tensor, Real alpha_phi_i, Real k_c_i);
         inline Mat3d ReturnMapping(Mat3d &stress_tensor);
+        inline int ReturnMappingWithReduction(Mat3d &inner_stress_tensor, Real alpha_phi_i, Real k_c_i);
         inline Real getFrictionAngle() { return phi_; };
+        inline Real getCohesion() { return c_; };
+        inline Real gerParicleDiameter() {return d_s_;};
+        inline Real getFrictionVelocity(Real uz, Real z);
+        inline Real calculateThetaCr(Real u_star);
+        inline Real ThetaToFrictionVelcoty(Real theta_cr);
+        inline Mat3d computeBinghamViscousStress(const Mat3d &strain_rate, Real tau_y, Real eta);
+        inline Mat3d computeHBPViscousStress(const Mat3d &strain_rate, Real tau_y, Real eta);
 
       protected:
         Real c_;                                                   /* cohesion  */
         Real phi_;                                                 /* friction angle  */
         Real psi_;                                                 /* dilatancy angle  */
         Real alpha_phi_;                                           /* Drucker-Prager's constants */
         Real k_c_;                                                 /* Drucker-Prager's constants */
+        Real d_s_;                                                  /* Mean particle diameter */
         Real stress_dimension_ = 3.0; /* plain strain condition */ // Temporarily cancel const --need to check
     };
 };

src/shared/physical_closure/materials/general_continuum.hpp

@@ -20,12 +20,16 @@ Real GeneralContinuum::GeneralContinuumKernel::getLambda(Real youngs_modulus, Re
     return nu_ * youngs_modulus / (1.0 + poisson_ratio) / (1.0 - 2.0 * poisson_ratio);
 }
 //=================================================================================================//
-//=================================================================================================//
 Real PlasticContinuum::PlasticKernel::getDPConstantsA(Real friction_angle)
 {
     return tan(friction_angle) / sqrt(9.0 + 12.0 * tan(friction_angle) * tan(friction_angle));
 };
 //=================================================================================================//
+Real PlasticContinuum::PlasticKernel::getDPConstantsK(Real cohesion, Real friction_angle)
+{
+   return 3.0 * cohesion / sqrt(9.0 + 12.0 * tan(friction_angle) * tan(friction_angle));
+};
+//=================================================================================================//
 Mat3d PlasticContinuum::PlasticKernel::ConstitutiveRelation(Mat3d &velocity_gradient, Mat3d &stress_tensor)
 {
     Mat3d strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose());
@@ -63,5 +67,166 @@ Mat3d PlasticContinuum::PlasticKernel::ReturnMapping(Mat3d &stress_tensor)
     }
     return stress_tensor;
 }
+//=================================================================================================//
+Mat3d PlasticContinuum::PlasticKernel::ConstitutiveRelationWithReduction(Mat3d &velocity_gradient, Mat3d &stress_tensor, Real alpha_phi_i, Real k_c_i)
+{
+    Mat3d strain_rate = 0.5 * (velocity_gradient + velocity_gradient.transpose());
+    Mat3d spin_rate = 0.5 * (velocity_gradient - velocity_gradient.transpose());
+    Mat3d deviatoric_strain_rate = strain_rate - (1.0 / stress_dimension_) * strain_rate.trace() * Mat3d::Identity();
+    Mat3d stress_rate_elastic = 2.0 * G_ * deviatoric_strain_rate + K_ * strain_rate.trace() * Mat3d::Identity() + stress_tensor * (spin_rate.transpose()) + spin_rate * stress_tensor;
+    Mat3d deviatoric_stress_tensor = stress_tensor - (1.0 / stress_dimension_) * stress_tensor.trace() * Mat3d::Identity();
+    Real stress_tensor_J2 = 0.5 * (deviatoric_stress_tensor.cwiseProduct(deviatoric_stress_tensor.transpose())).sum();
+    Real f = sqrt(stress_tensor_J2) + alpha_phi_i * stress_tensor.trace() - k_c_i;
+    Real lambda_dot_ = 0;
+    Mat3d g = Mat3d::Zero();
+    if (f >= TinyReal)
+    {
+        Real deviatoric_stress_times_strain_rate = (deviatoric_stress_tensor.cwiseProduct(strain_rate)).sum();
+        // non-associate flow rule
+        lambda_dot_ = (3.0 * alpha_phi_i * K_ * strain_rate.trace() + (G_ / (sqrt(stress_tensor_J2)+TinyReal)) * deviatoric_stress_times_strain_rate) / (9.0 * alpha_phi_i * K_ * getDPConstantsA(psi_) + G_);
+        g = lambda_dot_ * (3.0 * K_ * getDPConstantsA(psi_) * Mat3d::Identity() + G_ * deviatoric_stress_tensor / (sqrt(stress_tensor_J2+ TinyReal)));
+    }
+    Mat3d stress_rate_temp = stress_rate_elastic - g;
+    return stress_rate_temp;
+}; 
+//=================================================================================================//
+int PlasticContinuum::PlasticKernel::ReturnMappingWithReduction(Mat3d &inner_stress_tensor, Real alpha_phi_i, Real k_c_i)
+{
+    // First: check if internal stress alone causes yielding (cohesion-dominated)
+    Real stress_I1 = inner_stress_tensor.trace();
+    if (-alpha_phi_i * stress_I1 + k_c_i < 0.0)
+    {
+        // Project to yield surface
+        inner_stress_tensor -= (1.0 / stress_dimension_) * 
+            (stress_I1 - k_c_i / alpha_phi_i) * Mat3d::Identity();
+        return 1;
+    }
+
+    // Compute deviatoric part and J2 for internal stress
+    Mat3d deviatoric_inner = inner_stress_tensor - 
+        (1.0 / stress_dimension_) * stress_I1 * Mat3d::Identity();
+    Real J2_inner = 0.5 * (deviatoric_inner.cwiseProduct(deviatoric_inner.transpose())).sum();
+
+    // Check yield condition based on J2
+    if (-alpha_phi_i * stress_I1 + k_c_i < sqrt(J2_inner))
+    {
+        Real r = (-alpha_phi_i * stress_I1 + k_c_i) / (sqrt(J2_inner) + TinyReal);
+        inner_stress_tensor = r * deviatoric_inner + 
+            (1.0 / stress_dimension_) * stress_I1 * Mat3d::Identity();
+        return 1;
+    }
+
+    return 0; // no yielding
+}
+//=================================================================================================//
+Real PlasticContinuum::PlasticKernel::getFrictionVelocity(Real uz, Real z)
+{
+    const Real nu = 1.0e-6;        // Kinematic viscosity (m^2/s)
+    const Real d50 = d_s_;         // Median grain diameter (m)
+    const Real kappa = 0.41;       // von Kármán constant
+    const Real ks = 2.0 * d50;     // Equivalent roughness height (m)
+    const Real tol = 0.01;         // Relative error tolerance
+
+    // Initial estimate of shear velocity u_star assuming laminar sublayer
+    Real u_star_1 = std::sqrt(uz * nu / z);
+    Real delta_v = 11.6 * nu / u_star_1;
+
+    // Check smooth flow condition
+    if ((ks * u_star_1 / nu < 5.0) && (z < delta_v)) {
+        return u_star_1;
+    }
+
+    // Iterative solution for u_star in rough/turbulent flow
+    Real u_star_old = u_star_1;
+    Real u_star_new = 0.0;
+    int iter = 0;
+
+    while (true) {
+        // Compute roughness length z0 based on u_star
+        Real ksu_nu = ks * u_star_old / nu;
+        Real z0;
+
+        if (ksu_nu < 5.0) {
+            z0 = 0.11 * nu / u_star_old;
+        } else if (ksu_nu > 70.0) {
+            z0 = 0.033 * ks;
+        } else {
+            z0 = 0.11 * nu / u_star_old + 0.033 * ks;
+        }
+
+        // Update u_star using the logarithmic velocity profile
+        u_star_new = (kappa * uz) / std::log(z / z0);
+
+        // Check convergence
+        Real rel_err = std::abs(u_star_new - u_star_old) / u_star_old;
+        if (rel_err < tol || iter > 100) break;
+
+        u_star_old = u_star_new;
+        ++iter;
+    }
+
+    return u_star_new;
+}
+//=================================================================================================//
+Real PlasticContinuum::PlasticKernel::calculateThetaCr(Real u_star)
+{
+    const Real rho_w = 1000;  // 
+    const Real mu_w = 0.001;  // 
+
+    Real Re = rho_w * u_star * d_s_ / mu_w;
+    if (Re <= 500 && Re > 0.0) 
+    {
+        return 0.010595 * log(Re) + 0.110476 / Re + 0.0027197;
+    } else 
+    {
+        return 0.068;
+    } 
+}
+//=================================================================================================//
+Real PlasticContinuum::PlasticKernel::ThetaToFrictionVelcoty(Real theta_cr)
+{
+    const Real rho_w = 1000.0;        // Water density (kg/m³)
+    const Real gravity = 9.8;         // Gravitational acceleration (m/s²)
+    const Real Sn = rho0_ / rho_w;    // Relative density = particle density / water density
+
+    // Compute critical shear velocity using:
+    // u*_cr = sqrt[ θ_cr * ( (ρ_s - ρ_w) * g * d ) / ρ_w ]
+    Real friction_velocity = sqrt(theta_cr * ((Sn - 1.0) * gravity * d_s_) / 1.0);
+
+    return friction_velocity;
+}
+//=================================================================================================//
+Mat3d PlasticContinuum::PlasticKernel::computeBinghamViscousStress(const Mat3d &strain_rate, Real tau_y, Real eta)
+{
+
+    Real shear_rate_norm = sqrt(0.5 * (strain_rate.cwiseProduct(strain_rate)).sum() + TinyReal);
+    Mat3d strain_dir = strain_rate / shear_rate_norm;
+    Mat3d viscous_stress = tau_y * strain_dir + eta * strain_rate;
+
+    return viscous_stress;
+}
+//=================================================================================================//
+Mat3d PlasticContinuum::PlasticKernel::computeHBPViscousStress(const Mat3d &strain_rate, Real tau_y, Real eta)
+{
+    const Real m = 100.0;  // Papanastasiou regularization parameter
+    const Real n = 1.2;    // Flow behavior index
+
+    // Compute shear rate magnitude: γ_dot = sqrt(0.5 * D_ij * D_ij)
+    Real shear_rate = std::sqrt(0.5 * (strain_rate.cwiseProduct(strain_rate)).sum());
+
+    if (shear_rate < TinyReal)
+        return Mat3d::Zero();  // Near-zero strain: no viscous stress
+
+    // Normalize the strain rate tensor to get the direction
+    Mat3d strain_dir = strain_rate / shear_rate;
+
+    // Use max(shear_rate, TinyReal) to ensure numerical stability in pow()
+    Real tau_apparent = tau_y * (1.0 - std::exp(-m * shear_rate)) + eta * std::pow(std::max(shear_rate, TinyReal), n);
+
+    // Compute viscous stress tensor
+    return tau_apparent * strain_dir;
+}
+
+
 }// namespace SPH
 #endif //GENERAL_CONTINUUM_HPP

src/shared/shared_ck/particle_dynamics/continum_dynamics/all_continum_dynamics.h

@@ -36,4 +36,10 @@
 #include "stress_diffusion_ck.h"
 #include "stress_diffusion_ck.hpp"
 #include "initilization_dynamics_ck.h"
-#include "initilization_dynamics_ck.hpp"
+#include "initilization_dynamics_ck.hpp"
+#include "erosion_dynamics_1st_ck.h"
+#include "erosion_dynamics_1st_ck.hpp"
+#include "erosion_dynamics_2nd_ck.h"
+#include "erosion_dynamics_2nd_ck.hpp"
+#include "erosion_interaction_ck.h"
+#include "erosion_interaction_ck.hpp"

src/shared/shared_ck/particle_dynamics/continum_dynamics/continuum_integration_1st_ck.h

@@ -33,6 +33,7 @@
 #include "constraint_dynamics.h"
 #include "general_continuum.h"
 #include "general_continuum.hpp"
+#include "erosion_interaction_ck.hpp"
 namespace SPH
 {
 namespace continuum_dynamics
@@ -49,7 +50,6 @@ class PlasticAcousticStep : public fluid_dynamics::AcousticStep<BaseInteractionT
 
   protected:
     PlasticContinuum &plastic_continuum_;
-    // DiscreteVariable<Vecd> *dv_pos_;
     DiscreteVariable<Mat3d> *dv_stress_tensor_3D_, *dv_strain_tensor_3D_, *dv_stress_rate_3D_, *dv_strain_rate_3D_;
     DiscreteVariable<Matd> *dv_velocity_gradient_;
 };

src/shared/shared_ck/particle_dynamics/continum_dynamics/continuum_integration_1st_ck.hpp

@@ -23,6 +23,7 @@ PlasticAcousticStep<BaseInteractionType>::PlasticAcousticStep(DynamicsIdentifier
     this->particles_->template addEvolvingVariable<Mat3d>("StrainTensor3D");
     this->particles_->template addEvolvingVariable<Mat3d>("StressRate3D");
     this->particles_->template addEvolvingVariable<Mat3d>("StrainRate3D");
+    this->particles_->template addEvolvingVariable<Matd>("VelocityGradient");
 }
 //=================================================================================================//
 template <class RiemannSolverType, class KernelCorrectionType, typename... Parameters>

src/shared/shared_ck/particle_dynamics/continum_dynamics/continuum_integration_2nd_ck.h

@@ -88,11 +88,15 @@ class PlasticAcousticStep2ndHalf<Inner<OneLevel, RiemannSolverType, KernelCorrec
         Real *rho_, *drho_dt_;
         Matd *velocity_gradient_;
         Mat3d *stress_tensor_3D_, *strain_tensor_3D_, *stress_rate_3D_, *strain_rate_3D_;
+        Vecd *shear_vel_;
+        Real *friction_angle_, *cohesion_, *reduction_para_;
+        Real particle_spacing_;
     };
 
   protected:
     KernelCorrectionType correction_;
     RiemannSolverType riemann_solver_;
+    Real dv_particle_spacing_;
 };
 
 template <class RiemannSolverType, class KernelCorrectionType, typename... Parameters>
@@ -127,6 +131,38 @@ class PlasticAcousticStep2ndHalf<Contact<Wall, RiemannSolverType, KernelCorrecti
     RiemannSolverType riemann_solver_;
 };
 
+template <class RiemannSolverType, class KernelCorrectionType, typename... Parameters>
+class PlasticAcousticStep2ndHalf<Contact<Fluid, RiemannSolverType, KernelCorrectionType, Parameters...>>
+    : public PlasticAcousticStep<Interaction<Contact<Parameters...>>>, public Interaction<Fluid>
+{
+    using BaseInteraction = PlasticAcousticStep<Interaction<Contact<Parameters...>>>;
+
+  public:
+    explicit PlasticAcousticStep2ndHalf(Contact<Parameters...> &fluid_contact_relation);
+    virtual ~PlasticAcousticStep2ndHalf() {};
+
+    class InteractKernel : public BaseInteraction::InteractKernel
+    {
+      public:
+        template <class ExecutionPolicy, class EncloserType>
+        InteractKernel(const ExecutionPolicy &ex_policy, EncloserType &encloser, UnsignedInt contact_index);
+        void interact(size_t index_i, Real dt = 0.0);
+
+      protected:
+        KernelCorrectionType correction_;
+        RiemannSolverType riemann_solver_;
+        Real *Vol_, *rho_, *drho_dt_;
+        Vecd *vel_, *force_;
+        Vecd *shear_vel_; /*For erosion*/
+        Real *fluid_Vol_, *fluid_p_;
+        Vecd *fluid_vel_;
+    };
+
+  protected:
+    KernelCorrectionType correction_;
+    RiemannSolverType riemann_solver_;
+};
+
 using PlasticAcousticStep2ndHalfWithWallRiemannCK =
     PlasticAcousticStep2ndHalf<Inner<OneLevel, AcousticRiemannSolverCK, NoKernelCorrection>,
                                Contact<Wall, AcousticRiemannSolverCK, NoKernelCorrection>>;