Cylindrical packed bed (#1369)

Description
The steady-state flow withing a cylindrical packed bed example was working fine, but the documentation was outdated. Furthermore, the figures we atrociously bad and were not representative enough. I have updated the example results and I have added new figures which illustrate the behavior of the VANS equation in a better fashion. I have also added a simple python post-processing script to do these calculations.

Co-authored-by: Laura Prieto Saavedra <[email protected]>

doc/source/examples/unresolved-cfd-dem/cylindrical-packed-bed/cylindrical-packed-bed.rst

@@ -2,7 +2,7 @@
 Cylindrical Packed Bed
 ==================================
 
-It is strongly recommended to visit `DEM parameters <../../../parameters/dem/dem.html>`_  and `CFD-DEM parameters <../../../parameters/unresolved-cfd-dem/unresolved-cfd-dem.html>`_ for more detailed information on the concepts and physical meaning of the parameters ind DEM and CFD-DEM.
+It is strongly recommended to visit `DEM parameters <../../../parameters/dem/dem.html>`_  and `CFD-DEM parameters <../../../parameters/unresolved-cfd-dem/unresolved-cfd-dem.html>`_ for more detailed information on the concepts and physical meaning of the parameters in DEM and CFD-DEM.
 
 
 ----------------------------------
@@ -16,7 +16,7 @@ Features
 
 
 ---------------------------
-Files Used in This Example
+Files Used in this Example
 ---------------------------
 
 Both files mentioned below are located in the example's folder (``examples/unresolved-cfd-dem/cylindrical-packed-bed``).
@@ -51,7 +51,7 @@ The ``mesh`` subsection specifies the computational grid:
       set type               = dealii
       set grid type          = subdivided_cylinder
       set grid arguments     = 16:0.01:0.1
-      set initial refinement = 1
+      set initial refinement = 2
     end
 
 Simulation Control
@@ -95,7 +95,6 @@ The section on model parameters is explained in the DEM examples. We show the ch
 
   subsection model parameters
     set contact detection method               = dynamic
-    set contact detection frequency            = 10
     set neighborhood threshold                 = 1.3
     set particle particle contact force method = hertz_mindlin_limit_overlap
     set particle wall contact force method     = nonlinear
@@ -227,9 +226,8 @@ The simulation is run in steady state. The simulation control section is shown:
 .. code-block:: text
 
     subsection simulation control
-      set method            = bdf1
-      set output name       = result
-      set output path       = ./output/
+      set method      = bdf1
+      set output path = ./output/
     end
    
 Physical Properties
@@ -263,7 +261,7 @@ For the initial conditions, we choose zero initial conditions for the velocity.
 Boundary Conditions
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-For the boundary conditions, we choose a slip boundary condition on the walls of the cylinder (ID = 0) and an inlet velocity of 0.2 m/s at the lower face of the cylinder (ID = 1). 
+For the boundary conditions, we set a slip boundary condition on the walls of the cylinder (ID = 0), an inlet velocity of 0.2 m/s at the lower face of the cylinder (ID = 1) and an outlet at the upper face of the cylinder (ID = 2).
 
 .. code-block:: text
 
@@ -286,6 +284,10 @@ For the boundary conditions, we choose a slip boundary condition on the walls of
           set Function expression = 0
         end
       end
+      subsection bc 2
+        set id   = 2
+        set type = outlet
+      end
     end
 
 
@@ -351,12 +353,13 @@ Linear Solver
   subsection linear solver
     subsection fluid dynamics
       set method                                = gmres
-      set max iters                             = 5000
+      set max iters                             = 2000
+      set max krylov vectors                    = 2000
       set relative residual                     = 1e-3
-      set minimum residual                      = 1e-11
+      set minimum residual                      = 1e-10
       set preconditioner                        = ilu
-      set ilu preconditioner fill               = 1
-      set ilu preconditioner absolute tolerance = 1e-14
+      set ilu preconditioner fill               = 4
+      set ilu preconditioner absolute tolerance = 1e-12
       set ilu preconditioner relative tolerance = 1.00
       set verbosity                             = verbose
     end
@@ -372,18 +375,18 @@ The simulation is run using the ``lethe-fluid-vans`` application. Assuming that
 .. code-block:: text
   :class: copy-button
 
-  lethe-fluid-vans parameter_file.prm
+  mpirun -np 8 lethe-fluid-vans flow-in-bed.prm
 
 
 -------------
 Results VANS
 -------------
-The results are shown in the plots below. We visualise the velocity of the fluid, the void fraction calculated using the particles' locations, and the pressure drop resulting from the particle-fluid interactions (drag). The plots to the right show the local distribution of the quantities at the center-line of the cylinder. 
+The results are shown in the plots below. The first plot visualises the velocity of the fluid and the void fraction along the center axis of the cylinder. We see that the fluid rapidly accelerates in as the void fraction decreases to ensure mass conservation. The second plot displays the pressure drop resulting from the particle-fluid interactions. The quasi totality of the pressure drop occurs within the bed of particles.
 
-.. image:: images/packed-bed-vel.png
+.. image:: images/velocity_void_fraction.png
     :alt: velocity and void fraction distribution
     :align: center
 
-.. image:: images/packed-bed-p.png
+.. image:: images/pressure_drop_void_fraction.png
     :alt: pressure drop in packed bed
     :align: center

examples/unresolved-cfd-dem/cylindrical-packed-bed/flow-in-bed.prm

@@ -12,7 +12,6 @@ set dimension = 3
 
 subsection simulation control
   set method      = steady
-  set output name = result
   set output path = ./output/
 end
 
@@ -44,7 +43,7 @@ subsection mesh
   set type               = dealii
   set grid type          = subdivided_cylinder
   set grid arguments     = 16:0.01:0.1
-  set initial refinement = 1
+  set initial refinement = 2
 end
 
 #---------------------------------------------------
@@ -55,7 +54,7 @@ subsection void fraction
   set mode                = pcm
   set read dem            = true
   set dem file name       = dem
-  set l2 smoothing factor = 0.000005
+  set l2 smoothing factor = 0.00001
   set l2 lower bound      = 0
   set l2 upper bound      = 1
 end
@@ -110,19 +109,11 @@ subsection boundary conditions
     end
   end
   subsection bc 2
-    set id = 2
+    set id   = 2
     set type = outlet
   end
 end
 
-#---------------------------------------------------
-# Mesh Adaptation Control
-#---------------------------------------------------
-
-subsection mesh adaptation
-  set type = none
-end
-
 #---------------------------------------------------
 # Timer
 #---------------------------------------------------
@@ -137,7 +128,7 @@ end
 
 subsection non-linear solver
   subsection fluid dynamics
-    set tolerance      = 1e-9
+    set tolerance      = 1e-8
     set max iterations = 20
     set verbosity      = verbose
   end
@@ -150,13 +141,13 @@ end
 subsection linear solver
   subsection fluid dynamics
     set method                                = gmres
-    set max iters                             = 1000
-    set max krylov vectors                    = 1000
+    set max iters                             = 2000
+    set max krylov vectors                    = 2000
     set relative residual                     = 1e-3
-    set minimum residual                      = 1e-11
+    set minimum residual                      = 1e-10
     set preconditioner                        = ilu
-    set ilu preconditioner fill               = 3
-    set ilu preconditioner absolute tolerance = 1e-14
+    set ilu preconditioner fill               = 4
+    set ilu preconditioner absolute tolerance = 1e-12
     set ilu preconditioner relative tolerance = 1.00
     set verbosity                             = verbose
   end

examples/unresolved-cfd-dem/cylindrical-packed-bed/packing-particles.prm

@@ -23,7 +23,7 @@ end
 #---------------------------------------------------
 
 subsection timer
-  set type = none
+  set type = iteration
 end
 
 #---------------------------------------------------
@@ -119,5 +119,5 @@ subsection mesh
   set type               = dealii
   set grid type          = subdivided_cylinder
   set grid arguments     = 16:0.01:0.1
-  set initial refinement = 1
+  set initial refinement = 2
 end

examples/unresolved-cfd-dem/cylindrical-packed-bed/post_process_packed_bed.py

@@ -0,0 +1,95 @@
+# SPDX-FileCopyrightText: Copyright (c) 2022, 2024 The Lethe Authors
+# SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
+
+"""
+Postprocessing code for cylindrical packed bed
+This code extracts the pressure along the x axis at the center of the cylinder
+
+"""
+
+# Modules
+#-------------------------------------------
+import numpy as np
+import matplotlib.pyplot as plt
+import pyvista as pv
+
+# Plot font and colors
+#---------------------
+font = {'weight' : 'normal',
+        'size'   : 13}
+
+plt.rc('font', **font)
+colors=['#1b9e77','#d95f02','#7570b3','#e7298a','#66a61e','#e6ab02']
+
+#--------------------------------------------
+# Main
+#--------------------------------------------
+
+# Load VTU file
+vtu_file="./output/out.00001.00000.vtu"
+sim = pv.read(vtu_file)
+sim.set_active_scalars("pressure")
+
+# Create begin and end point of line
+a = [-0.1, 0, 0]
+b = [0.1, 0, 0]
+
+# Extract all field over the line using pyvista
+sampled_data=sim.sample_over_line(a, b, resolution=1000)
+
+# Get pressure from sampled data
+x = sampled_data["Distance"]
+p = sampled_data["pressure"][:]
+
+
+
+# Get void fraction from sampled data
+sim.set_active_scalars("void_fraction")
+sampled_data=sim.sample_over_line(a, b, resolution=1000)
+x_2 = sampled_data["Distance"]
+void_fraction = sampled_data["void_fraction"][:]
+
+# Create the figure and the first axis
+fig, ax1 = plt.subplots()
+
+# Plot the data for the first y-axis
+ax1.plot(x,p,label="Pressure drop",color=colors[0])
+ax1.set_xlabel('Position')  # Common x-axis label
+ax1.set_ylabel('Pressure', color=colors[0])
+ax1.set_ylim([-1,10])
+ax1.tick_params(axis='y', labelcolor=colors[0])  # Set y-axis tick color
+ax1.set_xticks([0,0.05,0.1,0.15,0.2]) #xticks(np.arange(min(x), max(x)+1, 1.0))
+
+# Create a second y-axis sharing the same x-axis
+ax2 = ax1.twinx()
+ax2.plot(x_2,void_fraction,label="Void Fraction",color=colors[1])
+ax2.set_ylabel('Void fraction', color=colors[1])
+ax2.tick_params(axis='y', labelcolor=colors[1])  # Set y-axis tick color
+plt.savefig("pressure_drop_void_fraction.png",dpi=200)
+plt.show()
+
+
+# Get u component of the velocity from sampled data
+sim.set_active_vectors("velocity")
+y = sampled_data["Distance"]
+u = sampled_data["velocity"][:,0]
+
+# Create the figure and the first axis
+fig, ax1 = plt.subplots()
+
+# Plot the data for the first y-axis
+ax1.plot(x,u,label="Velocity",color=colors[2])
+ax1.set_xlabel('Position')  # Common x-axis label
+ax1.set_ylabel('Velocity', color=colors[2])
+#ax1.set_ylim([-1,10])
+ax1.tick_params(axis='y', labelcolor=colors[2])  # Set y-axis tick color
+ax1.set_xticks([0,0.05,0.1,0.15,0.2]) #xticks(np.arange(min(x), max(x)+1, 1.0))
+
+# Create a second y-axis sharing the same x-axis
+ax2 = ax1.twinx()
+ax2.plot(x_2,void_fraction,label="Void Fraction",color=colors[1])
+ax2.set_ylabel('Void fraction', color=colors[1])
+ax2.tick_params(axis='y', labelcolor=colors[1])  # Set y-axis tick color
+plt.savefig("velocity_void_fraction.png",dpi=200)
+plt.show()
+

examples/unresolved-cfd-dem/gas-solid-spouted-cylinder-bed/packing-particles.prm

@@ -39,9 +39,10 @@ end
 #---------------------------------------------------
 # Post-processing
 #---------------------------------------------------
+
 subsection post-processing
- # Enable output of grid, granular temperature, and particles' average velocity
- set lagrangian post-processing = true
+  # Enable output of grid, granular temperature, and particles' average velocity
+  set lagrangian post-processing = true
 end
 
 #---------------------------------------------------
@@ -75,7 +76,7 @@ subsection lagrangian physical properties
     set size distribution type            = uniform
     set diameter                          = 0.005
     set number of particles               = 100000
-    set density particles                 = 100 
+    set density particles                 = 100
     set young modulus particles           = 1e7
     set poisson ratio particles           = 0.25
     set restitution coefficient particles = 0.97