CHNS Free Energy Postprocessing (#1095)

Description of the problem The Cahn-Hilliard equations are energy based equations for solving multiphase flows : the interface is the result of two energies in competition, bulk and interface energy. The system, which starts from a non-equilibrium state, lowers its total energy (bulk + interface) to reach equilibrium. In order to visualize this process I added the energy computation into the postprocessing section. Very minor : added the computation of the volume of each phases. Description of the solution New functions calculate_phase_energy() and write_phase_energy were implemented as well as new parameter entries in the postprocess subsection : compute phase energy (boolean) and phase energy name How Has This Been Tested? A new test was added : cahn_hilliard_navier_stokes_energy_postprocessing Documentation Postprocessing documentation was updated to include this new feature File modified : post_processing.rst
chaos-polymtl · Apr 18, 2024 · 50d0ffe · 50d0ffe
1 parent 79f3d07
commit 50d0ffe
Show file tree

Hide file tree

Showing 9 changed files with 411 additions and 9 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -3,6 +3,12 @@
 All notable changes to the Lethe project will be documented in this file.
 The format is based on [Keep a Changelog](http://keepachangelog.com/).
 
+## [Master] - 2024-04-17
+
+### Added
+
+- MINOR Added Cahn-Hilliard equations energy computation and output. The documentation was updated accordingly. [#1095](https://github.com/lethe-cfd/lethe/pull/1095)
+
 ## [Master] - 2024-04-16
 
 ### Fixed

diff --git a/applications_tests/lethe-fluid/cahn_hilliard_navier_stokes_energy_postprocessing.output b/applications_tests/lethe-fluid/cahn_hilliard_navier_stokes_energy_postprocessing.output
@@ -0,0 +1,58 @@
+Running on 1 MPI rank(s)...
+   Number of active cells:       64
+   Number of degrees of freedom: 243
+   Volume of triangulation:      4
+   Number of Cahn-Hilliard degrees of freedom: 162
+-----------------
+Phase statistics
+-----------------
+Min: -0.994529
+Max: 1
+Average: 0.607051
+Integral: 2.4282
+Volume phase 0: 3.2141
+Volume phase 1: 0.785898
+-------------
+Phase energy
+-------------
+Bulk energy: 0.667619
+Interface energy: 1.74433
+Total energy: 2.41195
+
+*******************************************************************************
+Transient iteration: 1        Time: 0.5      Time step: 0.5      CFL: 0
+*******************************************************************************
+-----------------
+Phase statistics
+-----------------
+Min: -1.01453
+Max: 1.00534
+Average: 0.60706
+Integral: 2.42824
+Volume phase 0: 3.21412
+Volume phase 1: 0.78588
+-------------
+Phase energy
+-------------
+Bulk energy: 0.687751
+Interface energy: 1.66703
+Total energy: 2.35478
+
+*******************************************************************************
+Transient iteration: 2        Time: 1        Time step: 0.5      CFL: 0.967317
+*******************************************************************************
+-----------------
+Phase statistics
+-----------------
+Min: -0.891115
+Max: 1.0132
+Average: 0.61502
+Integral: 2.46008
+Volume phase 0: 3.23004
+Volume phase 1: 0.76996
+-------------
+Phase energy
+-------------
+Bulk energy: 1.08564
+Interface energy: 1.19563
+Total energy: 2.28127
diff --git a/applications_tests/lethe-fluid/cahn_hilliard_navier_stokes_energy_postprocessing.prm b/applications_tests/lethe-fluid/cahn_hilliard_navier_stokes_energy_postprocessing.prm
@@ -0,0 +1,127 @@
+set dimension = 2
+
+subsection simulation control
+  set method           = bdf1
+  set time end         = 1
+  set time step        = 0.5
+  set output name      = out
+  set output frequency = 0
+end
+
+subsection initial conditions
+  set type = L2projection
+  subsection cahn hilliard
+    set Function constants  = nbr_refs=5
+    set Function expression = -tanh(sqrt(2*0.25)*(0.5 - sqrt(x*x+y*y))/(2^(-nbr_refs+1))); 0
+  end
+end
+
+subsection boundary conditions
+  set number         = 4
+  set time dependent = false
+  subsection bc 0
+    set id   = 0
+    set type = slip
+  end
+  subsection bc 1
+    set id   = 1
+    set type = slip
+  end
+  subsection bc 2
+    set id   = 2
+    set type = slip
+  end
+  subsection bc 3
+    set id   = 3
+    set type = slip
+  end
+end
+
+subsection boundary conditions cahn hilliard
+  set number         = 1
+  subsection bc 0
+    set id   = 0
+    set type = dirichlet
+    subsection phi
+      set Function constants  = nbr_refs=5
+      set Function expression = -tanh(sqrt(2*0.25)*(0.5 - sqrt(x*x+y*y))/(2^(-nbr_refs+1)))
+    end
+  end
+end
+
+subsection multiphysics
+  set fluid dynamics = true
+  set cahn hilliard  = true
+end
+
+subsection cahn hilliard
+  set potential smoothing coefficient = 0.0
+
+  subsection epsilon
+    set method = automatic
+  end
+end
+
+subsection physical properties
+  set number of fluids = 2
+  subsection fluid 0
+    set density             = 1000
+    set kinematic viscosity = 0.01
+  end
+  subsection fluid 1
+    set density             = 100
+    set kinematic viscosity = 0.01
+  end
+  set number of material interactions = 1
+  subsection material interaction 0
+    set type = fluid-fluid
+    subsection fluid-fluid interaction
+      set first fluid id              = 0
+      set second fluid id             = 1
+      set surface tension model       = constant
+      set surface tension coefficient = 24.5
+      # Mobility Cahn-Hilliard
+      set cahn hilliard mobility model    = constant
+      set cahn hilliard mobility constant = 1e-6
+    end
+  end
+end
+
+subsection mesh
+  set type               = dealii
+  set grid type          = hyper_cube
+  set grid arguments     = -1 : 1 : true
+  set initial refinement = 3
+end
+
+subsection post-processing
+  set verbosity                  = verbose
+  set calculate phase statistics = true
+  set phase statistics name      = phase_statistics
+  set calculate phase energy     = true
+end
+
+subsection FEM
+  set phase cahn hilliard order     = 1
+  set potential cahn hilliard order = 1
+  set velocity order                = 1
+  set pressure order                = 1
+end
+
+subsection non-linear solver
+  subsection fluid dynamics
+    set verbosity = quiet
+  end
+  subsection cahn hilliard
+    set verbosity = quiet
+  end
+end
+
+subsection linear solver
+  subsection fluid dynamics
+    set verbosity = quiet
+  end
+  subsection cahn hilliard
+    set verbosity = quiet
+  end
+end
diff --git a/applications_tests/lethe-fluid/time_dependent_boundaries_ch.output b/applications_tests/lethe-fluid/time_dependent_boundaries_ch.output
@@ -10,9 +10,11 @@ Min: -0.994529
 Max: 1
 Average: 0.544551
 Integral: 2.1782
+Volume phase 0: 3.0891
+Volume phase 1: 0.910898
 
 *******************************************************************************
-Transient iteration: 1        Time: 0.5      Time step: 0.5      CFL: 0       
+Transient iteration: 1        Time: 0.5      Time step: 0.5      CFL: 0
 *******************************************************************************
 -----------------
 Phase statistics
@@ -21,6 +23,8 @@ Min: -0.786164
 Max: 0.781734
 Average: 0.423636
 Integral: 1.69454
+Volume phase 0: 2.84727
+Volume phase 1: 1.15273
 
 *******************************************************************************
 Transient iteration: 2        Time: 1        Time step: 0.5      CFL: 0.533525
@@ -32,3 +36,5 @@ Min: -0.780241
 Max: 0.794182
 Average: 0.427977
 Integral: 1.71191
+Volume phase 0: 2.85595
+Volume phase 1: 1.14405
diff --git a/doc/source/parameters/cfd/post_processing.rst b/doc/source/parameters/cfd/post_processing.rst
@@ -67,6 +67,8 @@ This subsection controls the post-processing other than the forces and torque on
     # Other Cahn-Hilliard postprocessing
     set calculate phase statistics       = false
     set phase statistics name            = phase_statistics
+    set calculate phase energy           = false
+    set phase energy name                = phase_energy
     
   end
 
@@ -219,8 +221,36 @@ This subsection controls the post-processing other than the forces and torque on
 
 * ``mass conservation name``: name of the output file containing the mass of both fluids for VOF simulations. The default file name is ``mass_conservation_information``.
 
-* ``calculate phase statistics``: outputs phase statistics from the solution of the Cahn-Hilliard equations, including minimum, maximum, average, and standard deviation of the phase order parameter. This works only with the :doc:`cahn_hilliard` solver.
+* ``calculate phase statistics``: outputs Cahn-Hilliard phase statistics, including minimum, maximum, average, integral of the phase order parameter, and the volume of each phase.
+
+  .. warning ::
+
+      ``calculate phase statistics = true`` only works with the :doc:`cahn_hilliard` solver.
+
+* ``phase statistics name``: name of the output file containing phase order parameter statistics from Cahn-Hilliard simulations. The default file name is ``phase_statistics``. It is stored in the output folder with in a  ``.dat`` file.
+
+* ``calculate phase energy``: outputs Cahn-Hilliard phase energies, including bulk energy, interface energy and total energy. The energies are computed as follow:
+
+  .. math::
+
+     E_{bulk} = \int_{\Omega} (1-\phi^2)^2 \mathrm{d}\Omega 
+      
+  .. math::
+
+     E_{interface} = \int_{\Omega} 0.5\epsilon^2|\nabla \phi |^2 \mathrm{d}\Omega 
+      
+  .. math::
+
+     E_{total} = E_{bulk} + E_{interface}  
+    
+  where :math:`\epsilon` is the numerical interface thickness. Note that these energies are not homogeneous to physical energies. Nonetheless, they are a convenient way to track the system's evolution.
+
+  .. warning ::
+
+      ``calculate phase energy = true`` only works with the :doc:`cahn_hilliard` solver.
+
+
+* ``phase energy name``: name of the output file containing phase energies from Cahn-Hilliard simulations. The default file name is ``phase_energy``.
 
-* ``phase statistics name``: name of the output file containing phase order parameter statistics from Cahn-Hilliard simulations. The default file name is ``phase_statistics``.
 
 
diff --git a/include/core/parameters.h b/include/core/parameters.h
@@ -907,6 +907,12 @@ namespace Parameters
     /// Prefix for the VOF mass conservation output
     std::string mass_conservation_output_name;
 
+    /// Enable energies calculation on the domain in Cahn-Hilliard simulations
+    bool calculate_phase_energy;
+
+    /// Prefix for the energy output in Cahn-Hilliard simulations
+    std::string phase_energy_output_name;
+
     static void
     declare_parameters(ParameterHandler &prm);
     void

diff --git a/include/solvers/cahn_hilliard.h b/include/solvers/cahn_hilliard.h
@@ -396,17 +396,30 @@ class CahnHilliard : public AuxiliaryPhysics<dim, GlobalVectorType>
   copy_local_rhs_to_global_rhs(const StabilizedMethodsCopyData &copy_data);
 
   /**
-   * @brief Calculate phase order parameter integral for monitoring purposes
+   * @brief Calculate phase statistics for monitoring purposes
    */
   void
   calculate_phase_statistics();
 
   /**
-   * @brief Writes the phase integral to an output file
+   * @brief Writes the phase statistics to an output file
    */
   void
   write_phase_statistics();
 
+  /**
+   *
+   * @brief Calculate the phase energies :  bulk energy, interface energy and total energy.
+   */
+  void
+  calculate_phase_energy();
+
+  /*
+   * @brief Write the energy to an output file
+   */
+  void
+  write_phase_energy();
+
   /**
    * @brief Calculates the barycenter of fluid 1 and its velocity
    *
@@ -483,6 +496,9 @@ class CahnHilliard : public AuxiliaryPhysics<dim, GlobalVectorType>
   // Phase statistics table
   TableHandler statistics_table;
 
+  // Phase energy table
+  TableHandler phase_energy_table;
+
   // Phase fraction filter
   std::shared_ptr<CahnHilliardFilterBase> filter;
 };

diff --git a/source/core/parameters.cc b/source/core/parameters.cc
@@ -1828,6 +1828,18 @@ namespace Parameters
         "mass_conservation_information",
         Patterns::FileName(),
         "Name of mass conservation output file in VOF simulations");
+
+      prm.declare_entry(
+        "calculate phase energy",
+        "false",
+        Patterns::Bool(),
+        "Enable calculation of phase energies, including: total energy, bulk energy, and interface energy");
+
+      prm.declare_entry(
+        "phase energy name",
+        "phase_energy",
+        Patterns::FileName(),
+        "Name of energy output file in Cahn-Hilliard simulations. The file is stored in the output folder specified in the simulation control subsection");
     }
     prm.leave_subsection();
   }
@@ -1875,6 +1887,8 @@ namespace Parameters
       barycenter_output_name      = prm.get("barycenter name");
       calculate_mass_conservation = prm.get_bool("calculate mass conservation");
       mass_conservation_output_name = prm.get("mass conservation name");
+      calculate_phase_energy        = prm.get_bool("calculate phase energy");
+      phase_energy_output_name      = prm.get("phase energy name");
 
       // Viscous dissipative fluid
       const std::string op_fluid = prm.get("postprocessed fluid");