Bidirectional CFD-DEM variant of the particle transport

channel-transport-particles/fluid-nutils/README.md

@@ -0,0 +1,26 @@
+# Hints and TODOs
+
+## General hints
+
+- we already have a unidirectional mesh-particle tracing [tutorial](https://precice.org/tutorials-channel-transport-particles.html). This scenario is also available in the [tutorials repository](https://github.com/precice/tutorials/tree/develop/channel-transport-particles) and couples Nutils with MercuryDPM
+- the unidirectional tutorial only couples via the fluid velocity and should run out-of-the-box, use `setup-mercurydpm.sh` to run and compile the MercuryDPM participant
+- the tutorial uses OpenFOAM or Nutils, and MercuryDPM
+- all Nutils simulations use Nutils version 7
+- the bidirectional scenario builds upon this tutorial case
+
+## Before starting
+
+- make sure you have the latest develop of preCICE
+
+## The bidirectional coupling
+
+- the bidirectional coupling is pull request [706](https://github.com/precice/tutorials/pull/706)
+- from the [tutorials](https://github.com/precice/tutorials/tree/develop/channel-bidirectional), checkout the branch `channel-bidirectional` and navigate to `~/tutorials/channel-transport-particles/`
+- run the `setup-mercurydpm.sh` script again (it checks out another branch for the solver)
+- the overall setup is no too much physics-motivated, using the same boundary conditions leads to the following
+
+![snapshot-final](../images/tutorials-channel-transport-particles-t-final.png)
+
+- after you completed the case, I would validate it with a more complex scenario
+- while the simulation here is 2D, the final validation case will be 3D
+- I put all further comments directly into the `fluid.py` code

channel-transport-particles/fluid-nutils/fluid.py

@@ -13,9 +13,9 @@
 
 def main():
 
+    # COMMENT: nothing to change here
     print("Running Nutils")
 
-    # define the Nutils mesh
     nx = 48
     ny = 16
     step_start = nx // 3
@@ -28,13 +28,19 @@ def main():
         step_start:step_end, :step_hight
     ].withboundary(wall="left,top,right")
 
-    # cloud of Gauss points
     gauss = domain.sample("gauss", degree=4)
 
-    # Nutils namespace
     ns = function.Namespace()
     ns.x = geom
 
+    # COMMENT: looks like the weak is being defined here: the weak form now needs to incorporate the volume fraction
+    # Note that preCICE/MercuryDPM provides a field called "Alpha", which is the solid volume fraction. The formula
+    # 7.1 and 7.2 from the PDF use the fluid volume frection (epsilon_f), which is epsilon_f = (1 - Alpha).
+    # Note that reading from preCICE makes only sense after the initialize call.
+    # Moreover, we need to incorporate gravity (ε_f*ρ_f*g) in 7.2 and the drag force. The drag force is here zero either
+    # Note that the drag force is called "ExplicitMomentum" in preCICE. I
+    # would maybe later still try to implement the explicit-implicit split,
+    # that's the only reason
     ns.ubasis = domain.basis("std", degree=2).vector(2)
     ns.pbasis = domain.basis("std", degree=1)
     ns.u_i = "ubasis_ni ?u_n"  # solution
@@ -48,7 +54,7 @@ def main():
     ures = gauss.integral("ubasis_ni,j σ_ij d:x" @ ns)
     pres = gauss.integral("pbasis_n u_k,k d:x" @ ns)
 
-    # define Dirichlet boundary condition
+    # Dirichlet boundary condition
     sqr = domain.boundary["inflow"].integral("(u_0 - uin)^2 d:x" @ ns, degree=2)
     sqr += domain.boundary["inflow,outflow"].integral("u_1^2 d:x" @ ns, degree=2)
     sqr += domain.boundary["wall"].integral("u_k u_k d:x" @ ns, degree=2)
@@ -57,13 +63,16 @@ def main():
     # preCICE setup
     participant = precice.Participant("Fluid", "../precice-config.xml", 0, 1)
 
-    # define coupling mesh
+    # define coupling mesh (nothing to change)
     mesh_name = "Fluid-Mesh"
     vertices = gauss.eval(ns.x)
     vertex_ids = participant.set_mesh_vertices(mesh_name, vertices)
 
-    # coupling data
+    # COMMENT: names of the new coupling data fields in the preCICE configuration. The velocity was already there before
+    solid_fraction_name = "Alpha"
+    drag_force_name = "ExplicitMomentum"
     data_name = "Velocity"
+    pressure_gradient_name = "PressureGradientFull"
 
     participant.initialize()
 
@@ -90,13 +99,24 @@ def main():
         # potentially adjust non-matching timestep sizes
         dt = min(solver_dt, precice_dt)
 
+        # COMMENT: here we would need to read the solid volume fraction and the drag force from preCICE and incorporate them into the weak form. The pressure gradient would also need to be written to preCICE, but let's start with the velocity first
+        # Note that the read_data call takes a dt argument, which is the relative read time. For Euler backwards, it would simply be dt
+        # solid_fraction_values = participant.read_data(mesh_name, solid_fraction_name, vertex_ids, ...)
+        # drag_force_name = participant.read_data(mesh_name, drag_force_name, vertex_ids, ...)
+        # Checkpointing for implicit coupling is generally not required
+
         # solve Nutils timestep
         state["u0"] = state["u"]
         state["dt"] = dt
         state = solver.newton(("u", "p"), (ures, pres), constrain=cons, arguments=state).solve(1e-10)
 
         velocity_values = gauss.eval(ns.u, **state)
+        # COMMENT: here we would also need to write the pressure gradient to preCICE
+        # Important note: OpenFOAM scales the pressue (and thereby also its gradient) by the fluid density, so we would
+        # need to do the same scaling here. I set the fluid density to 1.2
         participant.write_data(mesh_name, data_name, vertex_ids, velocity_values)
+        # rho_f = 1.2
+        # participant.write_data(mesh_name, pressure_gradient_name, vertex_ids, pressure_gradient_values)
 
         # do the coupling
         participant.advance(dt)

channel-transport-particles/fluid-openfoam/constant/transportProperties

@@ -8,4 +8,7 @@ FoamFile
 
 transportModel  Newtonian;
 
-nu              nu [ 0 2 -1 0 0 0 0 ] 1;
+nu              [0 2 -1 0 0 0 0] 1.5e-05; // = air at 20°C m^2/s
+rho             1.2; // kg/m^3
+g               (0 -9.81 0);
+

channel-transport-particles/fluid-openfoam/system/controlDict

@@ -6,7 +6,7 @@ FoamFile
     object      controlDict;
 }
 
-application     pimpleFoam;
+application     AndersonJacksonFoam;
 
 startFrom       startTime;
 

channel-transport-particles/fluid-openfoam/system/fvSchemes

@@ -18,7 +18,7 @@ gradSchemes
 
 divSchemes
 {
-    default         none;
+    default         Gauss linear;
     div(phi,U)      bounded Gauss upwind;
     div((nuEff*dev2(T(grad(U))))) Gauss linear;
 }

channel-transport-particles/fluid-openfoam/system/preciceDict

@@ -12,6 +12,14 @@ participant Fluid;
 
 modules (FF);
 
+FF
+{
+  nameAlpha alphaSolid;
+  nameExplicitMomentum ExplicitMomentum;
+  nameImplicitMomentum ImplicitMomentum;
+};
+
+
 interfaces
 {
   Interface1
@@ -22,11 +30,15 @@ interfaces
 
     readData
     (
+        ImplicitMomentum
+        ExplicitMomentum
+        Alpha
     );
 
     writeData
     (
         Velocity
+        PressureGradientFull
     );
   };
-};
+};

channel-transport-particles/particles-mercurydpm/setup-mercurydpm.sh

@@ -4,7 +4,7 @@ set -euo pipefail
 # ---- Config (you can override via env vars) ----
 REPO_HTTPS_DEFAULT="https://bitbucket.org/davidscn/mercurydpm.git"
 REPO_URL="${REPO_URL:-$REPO_HTTPS_DEFAULT}"          # or git@bitbucket.org:davidscn/mercurydpm.git
-BRANCH="${BRANCH:-channel-transport-tutorial}"
+BRANCH="${BRANCH:-channel-bidirectional}"
 SRC_DIR="${SRC_DIR:-$PWD/mercurydpm}"
 BUILD_DIR="${BUILD_DIR:-$SRC_DIR/build}"
 

channel-transport-particles/precice-config.xml

@@ -1,43 +1,84 @@
 <?xml version="1.0" encoding="UTF-8" ?>
-<precice-configuration experimental="true">
+<precice-configuration experimental="True">
   <log>
     <sink
-      filter="%Severity% > debug and %Rank% = 0"
+      type="stream"
+      output="stdout"
       format="---[precice] %ColorizedSeverity% %Message%"
-      enabled="true" />
+      filter="(%Severity% > debug) and not ((%Rank% != 0))" />
   </log>
 
-  <profiling mode="fundamental" synchronize="false" />
+  <profiling mode="off" />
 
-  <data:vector name="Velocity" />
+  <data:scalar name="Alpha" waveform-degree="2" />
+  <data:vector name="ExplicitMomentum" waveform-degree="2" />
+  <data:scalar name="ImplicitMomentum" waveform-degree="2" />
+  <data:vector name="Velocity" waveform-degree="0" />
+  <data:vector name="PressureGradientFull" waveform-degree="0" />
 
   <mesh name="Fluid-Mesh" dimensions="2">
     <use-data name="Velocity" />
+    <use-data name="PressureGradientFull" />
+    <use-data name="Alpha" />
+    <use-data name="ExplicitMomentum" />
+    <use-data name="ImplicitMomentum" />
   </mesh>
 
   <participant name="Fluid">
     <provide-mesh name="Fluid-Mesh" />
     <write-data name="Velocity" mesh="Fluid-Mesh" />
+    <write-data name="PressureGradientFull" mesh="Fluid-Mesh" />
+    <read-data name="Alpha" mesh="Fluid-Mesh" />
+    <read-data name="ExplicitMomentum" mesh="Fluid-Mesh" />
+    <read-data name="ImplicitMomentum" mesh="Fluid-Mesh" />
   </participant>
 
   <participant name="Particles">
-    <receive-mesh name="Fluid-Mesh" from="Fluid" api-access="true" />
+    <receive-mesh name="Fluid-Mesh" from="Fluid" safety-factor="0.5" api-access="true" />
     <read-data name="Velocity" mesh="Fluid-Mesh" />
+    <read-data name="PressureGradientFull" mesh="Fluid-Mesh" />
+    <write-data name="Alpha" mesh="Fluid-Mesh" />
+    <write-data name="ExplicitMomentum" mesh="Fluid-Mesh" />
+    <write-data name="ImplicitMomentum" mesh="Fluid-Mesh" />
     <mapping:rbf-pum-direct
       direction="read"
       from="Fluid-Mesh"
       constraint="consistent"
       vertices-per-cluster="30">
-      <basis-function:compact-polynomial-c4 support-radius="1.5" />
+      <basis-function:compact-polynomial-c4 support-radius="0.5" />
     </mapping:rbf-pum-direct>
+    <mapping:coarse-graining
+      direction="write"
+      to="Fluid-Mesh"
+      constraint="conservative"
+      radius="0.05" />
   </participant>
 
   <m2n:sockets acceptor="Fluid" connector="Particles" exchange-directory=".." />
 
-  <coupling-scheme:parallel-explicit>
+  <coupling-scheme:serial-explicit>
+    <participants second="Fluid" first="Particles" />
     <time-window-size value="0.005" />
     <max-time value="0.25" />
-    <participants first="Fluid" second="Particles" />
-    <exchange data="Velocity" mesh="Fluid-Mesh" from="Fluid" to="Particles" />
-  </coupling-scheme:parallel-explicit>
+    <exchange data="Velocity" mesh="Fluid-Mesh" from="Fluid" to="Particles" substeps="false" />
+    <exchange
+      data="PressureGradientFull"
+      mesh="Fluid-Mesh"
+      from="Fluid"
+      to="Particles"
+      substeps="false" />
+    <exchange data="Alpha" mesh="Fluid-Mesh" from="Particles" to="Fluid" substeps="true" />
+    <exchange
+      data="ExplicitMomentum"
+      mesh="Fluid-Mesh"
+      from="Particles"
+      to="Fluid"
+      substeps="true" />
+    <exchange
+      data="ImplicitMomentum"
+      mesh="Fluid-Mesh"
+      from="Particles"
+      to="Fluid"
+      substeps="true" />
+  </coupling-scheme:serial-explicit>
 </precice-configuration>