Improved entrainment model Medina in terms of stability

src/avalanche/entrainment/Medina/Medina.C

@@ -60,6 +60,7 @@ Foam::entrainmentModels::Medina::Medina
     entrainmentModel(type(), entrainmentProperties, Us, h, hentrain, pb, tau),
     tauc_("tauc", coeffDict_),
     mu_("mu", coeffDict_),
+    relaxationFactor_(coeffDict_.lookupOrDefault<scalar>("relax", 1.)),
     gs_(Us_.db().lookupObject<areaVectorField>("gs")),
     gn_(Us_.db().lookupObject<areaScalarField>("gn"))
 {
@@ -73,12 +74,13 @@ Foam::entrainmentModels::Medina::Medina
 const Foam::areaScalarField&
 Foam::entrainmentModels::Medina::Sm() const
 {
-    areaScalarField taures(tauc_ + pb_*mu_);
+    areaScalarField taures(tauc_ + pb_.oldTime()*mu_);
     //no entrainment in "dry" areas:
-    areaScalarField eflag(pos(h_*rho_*gn_-tauc_/10.));
-    areaScalarField hent(eflag*(mag(tau_) - taures)/(rho_*(gn_*mu_-mag(gs_))));
+    areaScalarField eflag(pos(h_-dimensionedScalar("small", dimLength, 1e-2)));
+    areaScalarField eflag2(pos(mag(tau_.oldTime()) - taures));
+    areaScalarField hent(eflag*eflag2*(mag(tau_.oldTime()) - taures)/(rho_*max(gn_*mu_-mag(gs_), dimensionedScalar("small", dimAcceleration, 1e-2))));
 
-    Sm_ = hent/Us_.db().time().deltaT();
+    Sm_ = hent/Us_.db().time().deltaT() * relaxationFactor_;
 
     Sm_ = max(Sm_, dimensionedScalar(dimVelocity));
     Sm_ = min(Sm_, hentrain_/Us_.db().time().deltaT());

src/avalanche/entrainment/Medina/Medina.H

@@ -79,6 +79,8 @@ class Medina
         //- Friction coefficient tan(phi) of the basal material
         dimensionedScalar mu_;
 
+        scalar relaxationFactor_;
+
         //- Reference to the gravitation fields
         const areaVectorField& gs_;