Minor corrections in the documentation + adding the contribution of e…

…xcess fluid pressure in writeOuput (#199)

* adding excess fluid pressure contribution in writeOutput.H

* corrections in the documentation

doc/tutorials_DyM.doc

@@ -1,6 +1,6 @@
 /**
  * \file        tutorials_DyM.doc
- * \author      Cyrille Bonamy and Julien Chauchat
+ * \author      Eduard Puig Montellà, Cyrille Bonamy and Julien Chauchat
  * \date        June 03, 2021
  * \brief       Tutorials dynamic (doc)
  */
@@ -16,12 +16,12 @@ In this chapter we shall describe in detail the process of setup, simulation and
 
 
 
-Before launching any numerical simulation, make sure overSedDymFoam_rbgh is correctly installed and compiled. In case you are only interested in running dynamic mesh cases, the compilation is done by performing the following steps:
+Before launching any numerical simulation, make sure overSedDymFoam_rbgh is correctly installed and compiled:
 
 <ol>
-<li>First, download the overSedDyMFoam branch:
+<li>First, download the the official SedFoam package:
 
-        git clone -b overSedDyMFoam https://github.com/puigmontella/sedfoam.git
+        git clone --recurse-submodules http://github.com/sedfoam/sedfoam
 
 <li>Compile the program by running Allwmake:
 

doc/tutorials_laminar.doc

@@ -1183,10 +1183,10 @@ and *alphaPlastic=0.520* leads to \f$\phi_0=0.566\f$.
 \subsection meshgen2DAvalanche Mesh generation
 
 
-The numerical domain consists of a rectangle of side length \f$L = 1m\f$ in the x-z plane and height \f$H = 0.07m\f$ in the vertical direction. 
+The numerical domain consists of a rectangle of side length \f$L = 1m\f$ in the x-z plane and height \f$H = 0.07m\f$ in the vertical direction.
 The initial bed depth is  \f$h = 0.0049m\f$. The mesh consists of 1000 elements  in the horizontal direction and 70  elements in  the vertical direction.
 The meshing utility *snappyHexMesh* is  executed to generate a finer mesh in regions of high concentration.
- 
+
 
 The meshing and refinement  process are generated by running:
 
@@ -1413,7 +1413,7 @@ and you should see the following figures:
 
 @image html "tuto2DCollapse_Morphology.png" "Figure 16: Evolution of the morphology during the collapse of initially a) dense and b) loose columns." width=800px
 @image latex "tuto2DCollapse_Morphology.png" "Figure 16: Evolution of the morphology during the collapse of initially a) dense and b) loose columns."
-It is worth mentioning that the numerical cases have been significantly coarsened to obtain the results rapidly and illustrate the potential of the dilatancy model. We refer the reader to Montellà et al. (2023) for more detailed data and a thorough analysis
+It is worth mentioning that the numerical cases have been significantly coarsened to obtain the results rapidly and illustrate the potential of the dilatancy model. We refer the reader to Montellà et al. (2023) for more detailed data and a thorough analysis.
 
 
 */

solver/callGranularStress.H

@@ -36,7 +36,7 @@ if (granularRheology.dilatancy())
      + fvm::div(phia, alphaPlastic, "div(phia,alphaPlastic)")
      - fvm::Sp(fvc::div(phia), alphaPlastic)
      ==
-     -fvm::SuSp(delta*magDtensor, alphaPlastic)    );
+     -fvc::Sp(delta*magDtensor, alphaPlastic)    );
     phi_pl_Eqn.relax();
     phi_pl_Eqn.solve();
 

solver/createFields.H

@@ -637,7 +637,7 @@ volScalarField CohesionDistrb
         runTime.timeName(),
         mesh,
         IOobject::READ_IF_PRESENT,
-        IOobject::NO_WRITE
+        IOobject::AUTO_WRITE
     ),
     mesh,
     dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 0.0)
@@ -921,6 +921,19 @@ volScalarField SUStilde
     dimensionedScalar("one", dimensionSet(0, 0, 0, 0, 0), 1.0)
 );
 
+volVectorField ExcessPorePressureContr
+(
+    IOobject
+    (
+        "ExcessPorePressureContr",
+         runTime.timeName(),
+         mesh,
+         IOobject::NO_READ,
+         IOobject::NO_WRITE
+    ),
+    fvc::grad(p_rbgh)
+);
+
 volVectorField SolidPressureContr
 (
     IOobject

solver/interfacialModels/dragModels/Engelund/Engelund.C

@@ -70,7 +70,8 @@ Foam::tmp<Foam::volScalarField> Foam::Engelund::K
     volScalarField beta(max(scalar(1) - alpha_, scalar(1.0e-6)));
     volScalarField Cds
     (
-      (phasea_.aE()*alpha_*alpha_/beta+ phasea_.bE()*Ur*phasea_.d()
+      (phasea_.aE()*alpha_*max(alpha_, scalar(1.0e-3))/beta
+      + max(alpha_, scalar(1.0e-3))/scalar(0.6)*phasea_.bE()*Ur*phasea_.d()
       /(phaseb_.nu()*beta*beta))
     );
 

solver/overSedDymFoam/createFields.H

@@ -556,7 +556,7 @@ volScalarField CohesionDistrb
         runTime.timeName(),
         mesh,
         IOobject::READ_IF_PRESENT,
-        IOobject::NO_WRITE
+        IOobject::AUTO_WRITE
     ),
     mesh,
     dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 0.0)
@@ -941,6 +941,19 @@ volScalarField SUStilde
     dimensionedScalar("one", dimensionSet(0, 0, 0, 0, 0), 1.0)
 );
 
+volVectorField ExcessPorePressureContr
+(
+    IOobject
+    (
+        "ExcessPorePressureContr",
+         runTime.timeName(),
+         mesh,
+         IOobject::NO_READ,
+         IOobject::NO_WRITE
+    ),
+    fvc::grad(p_rbgh)
+);
+
 volVectorField SolidPressureContr
 (
     IOobject

solver/writeOutput.H

@@ -64,6 +64,7 @@ if (runTime.outputTime())
             /fvc::interpolate(alpha+ alphaSmall)
         );
 
+        ExcessPorePressureContr = fvc::grad(p_rbgh);
         SolidPressureContr = -fvc::grad(pa+pff);
         ViscStressContr = rhoa*alpha*fvc::laplacian(nuEffa, Ua);
         FrictContr = rhoa*fvc::laplacian(nuFra, Ua);
@@ -81,6 +82,8 @@ if (runTime.outputTime())
         SuspContr=-SUS*K*beta*turbulenceb->nut()
                   *((SUS_I*iso-SUS_A*aniso) & fvc::grad(alpha));
 
+        ExcessPorePressureContr.write();
+        SolidPressureContr.write();
         divTauaContr.write();
         gravityContr.write();
         DragContr.write();