Skip to content

Commit

Permalink
Resolves Rayleigh-Bénard Convection issue #1356 (#1403)
Browse files Browse the repository at this point in the history 
Description
The previous version of the Rayleigh-Benard example was a bit outdated. In this new version where we compare our simulation to results of the literature. We also produced another video for the simulation that is more fun to watch. We added a python script for the figures that follow the style of Lethe validation so it could be added to the tests if desired.

Co-authored-by: OGaboriault <[email protected]>
Co-authored-by: mivaia <[email protected]>
Co-authored-by: Bruno Blais <[email protected]>
Co-authored-by: Olivier Guévremont <[email protected]>
Co-authored-by: Amishga Alphonius <[email protected]>
  • Loading branch information
@OresteMarquis @OGaboriault
@mivaia @blaisb @oguevremont @AmishgaAlphonius
6 people authored Jan 15, 2025
1 parent 2f5f4c4 commit c6777e7
Show file tree
Hide file tree
Showing 8 changed files with 1,068 additions and 112 deletions.
Binary file modified doc/source/examples/multiphysics/rayleigh-benard-convection/images/geometry.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
757 changes: 757 additions & 0 deletions ...es/multiphysics/rayleigh-benard-convection/images/schematic_rayleigh-benard.svg
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added ...xamples/multiphysics/rayleigh-benard-convection/images/solution-rayleigh-uy.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Binary file added ...xamples/multiphysics/rayleigh-benard-convection/images/solution-rayleigh-xv.png
View file
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
112 changes: 66 additions & 46 deletions ...examples/multiphysics/rayleigh-benard-convection/rayleigh-benard-convection.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,7 +2,7 @@
Rayleigh-Bénard Convection
==========================

This example simulates two-dimensional Rayleigh–Benard convection [#venturi2010]_ [#mpi2022]_ at Rayleigh numbers of :math:`10^4` and :math:`2.5 \times 10^4` .
This example simulates a two-dimensional Rayleigh–Bénard convection [#ouertatani]_ [#venturi2010]_ [#mpi2022]_ at Rayleigh numbers of :math:`10^4` and :math:`10^6` .


----------------------------------
Expand All @@ -20,78 +20,76 @@ Files Used in This Example

Both files mentioned below are located in the example's folder (``examples/multiphysics/rayleigh-benard-convection``).

- Parameter file for :math:`Ra=10\, 000`: ``rayleigh-benard-convection-Ra10k.prm``
- Parameter file for :math:`Ra=25\, 000`: ``rayleigh-benard-convection-Ra25k.prm``
- Parameter file for :math:`Ra=10^4`: ``rayleigh-benard-convection-Ra10k.prm``
- Parameter file for :math:`Ra=10^6`: ``rayleigh-benard-convection-Ra1M.prm``


-----------------------------
Description of the Case
-----------------------------

In this example, we evaluate the performance of the ``lethe-fluid`` solver in the simulation of the stability of natural convection within a two-dimensional rectangular domain. The following schematic describes the geometry and dimensions of the simulation in the :math:`(x,y)` plane:
In this example, we evaluate the performance of the ``lethe-fluid`` solver in the simulation of the stability of natural convection within a two-dimensional square domain. Our results are tested against the benchmark presented by Ouertatani *et al.* [#ouertatani]_ . The following schematic describes the geometry and dimensions of the simulation:

.. image:: images/geometry.png
:alt: Schematic
:align: center
:width: 400

:alt: Schematic
:align: center
:width: 400

The incompressible Navier-Stokes equations with a Boussinesq approximation for the buoyant force are:
.. math::
\nabla \cdot {\bf{u}} = 0

.. math::
\rho \frac{\partial {\bf{u}}}{\partial t} + \rho ({\bf{u}} \cdot \nabla) {\bf{u}} = -\nabla p + \nabla \cdot {\bf{\tau}} - \rho \beta {\bf{g}} (T - T_0)
.. math::
\nabla \cdot {\bf{u}} &= 0 \\
\rho \frac{\partial {\bf{u}}}{\partial t} + \rho ({\bf{u}} \cdot \nabla) {\bf{u}} &= -\nabla p + \nabla \cdot {\bf{\tau}} - \rho \beta {\bf{g}} (T - T_0)
where :math:`\beta` and :math:`T_0` denote thermal expansion coefficient and a reference temperature, respectively.

A two-dimensional block of fluid is heated from its bottom wall at :math:`t = 0` s. The temperature of the bottom wall is equal to :math:`T_h=50`, the temperature of the top wall is equal to :math:`T_c=0`, and the left and right walls are insulated. By heating the fluid from the bottom wall, the buoyant force (natural convection) creates vortices inside the fluid. The shape and number of these vortices mainly depend on the Rayleigh number [1, 2]:

.. math::
\text{Ra} = \frac{\rho^2 \beta g (T_h - T_c) H^3 c_p}{k \mu}
A two-dimensional block of fluid is heated from its bottom wall at :math:`t = 0` s. The temperature of the bottom wall is equal to :math:`T_\text{h}=10`, the temperature of the top wall is equal to :math:`T_\text{c}=0`, and the left and right walls are insulated. By heating the fluid from the bottom wall, the buoyant force (natural convection) creates vortices inside the fluid. The shape and count of these vortices depend on the Rayleigh number and the Prandtl dimensionless numbers [#ouertatani]_ [#venturi2010]_ :

.. math::
Ra &= \frac{\rho^2 \beta g (T_\text{h} - T_\text{c}) H^3 c_\text{p}}{k \mu} \\
Pr &= \frac{\mu c_\text{p}}{k}
where :math:`\rho` is the fluid density, :math:`g` is the magnitude of gravitational acceleration, :math:`H` denotes the characteristic length, :math:`k` is the thermal conduction coefficient, and :math:`\mu` is the dynamic viscosity, and :math:`c_p` is the specific thermal capacity.
where :math:`\rho` is the fluid density, :math:`g` is the magnitude of gravitational acceleration, :math:`H` denotes the characteristic length, :math:`k` is the thermal conduction coefficient, :math:`\mu` is the dynamic viscosity, and :math:`c_\text{p}` is the specific heat capacity.

In this example, we simulate the Rayleigh-Bénard convection problem at two Rayleigh numbers of 10000 and 25000. According to the literature [1, 2], we should see different numbers (2 and 3 vortices at :math:`Ra=10^4` and :math:`2.5 \times 10^4`, respectively) of vortices in the fluid at these two Rayleigh numbers. The gravity magnitude is set to -10 for both simulations. Additionally, :math:`\rho = 100`, :math:`\beta = 0.0002`, :math:`H = 0.25`, :math:`c_p = 100` and :math:`\mu = 1`. Thus, the Rayleigh number is controlled only by the thermal conduction coefficient for this example. In other words, we change the Rayleigh number by changing the thermal conduction coefficient of the fluid.
In this example, we simulate the Rayleigh-Bénard convection problem at two Rayleigh numbers ( :math:`Ra=10^4` and :math:`Ra=10^6` ) with a Prandtl number of :math:`Pr=0.71` which correspond to air. According to the literature [#ouertatani]_ , we should see one big convective cell at steady-state for both :math:`Ra=10^4` and :math:`Ra=10^6`, but for the latter, there should also be two small vortices in opposite corners rotating in the reverse direction of the big vortex. For simplicity, the gravity magnitude is set to 10 pointing in the opposite direction of the y-axis for both simulations. Additionally, because the two dimensionless numbers above are the only thing that characterize the flow we may choose the remaining parameters as we want. Here we chose to fix :math:`\rho = 1`, :math:`H = 1`, :math:`c_\text{p} = 100` and :math:`\mu = 0.071`. Thus, the Rayleigh number is controlled only by the thermal expansion coefficient (:math:`\beta = 0.71` or :math:`\beta = 71`).

.. note::
All four boundary conditions are ``noslip``, and an external
gravity field of :math:`-10` is applied in the :math:`y` direction.
All four boundary conditions for the velocity are ``noslip``. The sides wall are set to ``convection-radiation-flux`` with a heat transfer coefficient of 0 to have isolated walls. The top and bottom walls are set to a ``temperature`` to force a value of 0 and 10, respectively.


--------------
Parameter File
--------------

For simplicity, we present only the parameter file for :math:`Ra=10^4`.

.. note::
Note that the resolution (256x256) is set to match the one of the article results [#ouertatani]_ . If desired, you can choose to reduced to 7 or 6 the ``initial refinement level`` in the parameter files to reduce the simulation time without loosing too much accuracy.

Simulation Control
~~~~~~~~~~~~~~~~~~

Time integration is handled by a 1st order backward differentiation scheme
`(bdf1)`, for a :math:`10000` s simulation time with an initial
time step of :math:`0.01` second.
The time integration is handled by a 1st order backward differentiation scheme
`(bdf1)`. A simulation time of :math:`24` s has been selected to ensure that the system reaches a quasi-steady state. The initial time-step is set to :math:`0.001` s, but this will be adapted during the simulation.

.. note::
This example uses an adaptive time-stepping method, where the
time-step is modified during the simulation to keep the maximum value of the CFL condition below a given threshold (0.5 here). Using ``output control = iteration``, and ``output frequency = 100`` the simulation results are written every 100 iteration regardless of the time steps.

.. note::
Note that the heating process is slow, and the velocity magnitudes are small inside the fluid. Hence, we expect large time-steps and a long simulation.
time-step is modified during the simulation to keep the maximum value of the CFL condition below a given threshold (0.9 here). Using ``output control = time``, and ``output time frequency = 0.5`` the simulation results are written every 0.5 seconds regardless of the time-steps.

.. code-block:: text
subsection simulation control
set method = bdf1
set time end = 10000
set time end = 24
set time step = 0.01
set adapt = true
set max cfl = 0.5
set max cfl = 0.9
set stop tolerance = 1e-5
set adaptative time step scaling = 1.3
set number mesh adapt = 0
set output name = rayleigh-benard_convection
set output control = iteration
set output frequency = 100
set output control = time
set output time frequency = 0.5
set output path = ./output/
end
Expand Down Expand Up @@ -124,19 +122,17 @@ The ``source term`` subsection defines gravitational acceleration.
Physical Properties
~~~~~~~~~~~~~~~~~~~

The ``physical properties`` subsection defines the physical properties of the fluid. Since we simulate the Rayleigh-Bénard convection at two Rayleigh numbers (:math:`Ra=10^4` and :math:`2.5 \times 10^4`), we use different thermal conductivities to reach mentioned Rayleigh numbers. We change the thermal conductivity of the fluid in the two simulations. Note that any other physical property (that is present in the Rayleigh number equation defined above) can be used instead of thermal conductivity. Both thermal conductivity values (:math:`k=0.15625` for :math:`Ra=10^4`, and :math:`k=0.0625` for :math:`Ra=2.5 \times 10^4`) are added to the parameter handler file. However, only one of them should be uncommented for each simulation.

The ``physical properties`` subsection defines the physical properties of the fluid.

.. code-block:: text
subsection physical properties
set number of fluids = 1
subsection fluid 0
set density = 100
set kinematic viscosity = 0.01
set thermal expansion = 0.0002
set thermal conductivity = 0.15625 # for Ra = 10000
#set thermal conductivity = 0.0625 # for Ra = 25000
set density = 1
set kinematic viscosity = 0.071
set thermal expansion = 0.71
set thermal conductivity = 10
set specific heat = 100
end
end
Expand All @@ -153,40 +149,64 @@ Call the ``lethe-fluid`` by invoking:
mpirun -np 8 lethe-fluid rayleigh-benard-convection-Ra10k.prm
and
or

.. code-block:: text
:class: copy-button
mpirun -np 8 lethe-fluid rayleigh-benard-convection-Ra25k.prm
mpirun -np 8 lethe-fluid rayleigh-benard-convection-Ra1M.prm
to run the simulations using eight CPU cores. Feel free to use more. Note that the first and second commands belong to the simulations at :math:`Ra=10^4` and :math:`Ra=2.5 \times 10^4`, repectively.
to run the simulations using eight CPU cores for the :math:`Ra=10^4` and :math:`Ra=10^6` cases respectively. Feel free to use more CPU if available.


.. warning::
Make sure to compile lethe in `Release` mode and
run in parallel using mpirun. This simulation takes
:math:`\approx` 20 minutes on 8 processes.
run in parallel using mpirun. The first simulation takes
:math:`\approx` 20 minutes on 8 processes and the second at :math:`Ra=10^6` can take a few hours because of the much smaller time-step required to respect the CFL condition.


-------
Results
-------

The following animation shows the results of this simulation:
The following animation shows the evolution of the temperature field with the flow direction for the simulation at :math:`Ra=10^6`:

.. raw:: html

<iframe width="560" height="315" src="https://www.youtube.com/embed/tEg5M-wiCp8" frameborder="0" allowfullscreen></iframe>
<iframe width="640" height="360" src="https://www.youtube.com/embed/NSJJpPauiXo" frameborder="0" allowfullscreen></iframe>

Below, we also present the velocity profiles at steady-state of our simulation compared to the ones presented by Ouertatani *et al.* [#ouertatani]_ as a verfification of the Lethe software.

|fig1| |fig2|

.. |fig1| image:: images/solution-rayleigh-uy.png
:width: 45%

Note that at :math:`Ra=10000`, two vortices exist in the fluid, while an extra (relatively small) vortex appears near the right wall for :math:`Ra=25000`. The velocity magnitude in the vortices is larger at smaller Rayleigh number.
.. |fig2| image:: images/solution-rayleigh-xv.png
:width: 47%

The results can be postprocessed using the provided Python script (``rayleigh-benard-convection.py``). Here is an example of how to call the script:

.. code-block:: text
:class: copy-button
python3 rayleigh-benard-convection.py -f ./output_10k -f ./output_1M -Ra 10k -Ra 1M
This script extracts the velocity in the :math:`x` and :math:`y` directions at the mid-width (:math:`x=0.5`) and mid-height (:math:`y=0.5`) respectively and create the above plots.

.. warning::
The orientation of the vortex rotation obtained with the simulation may differ from the one above due to machine precision that generates the initial instability.

.. important::
You need to ensure that the ``lethe_pyvista_tools`` is working on your machine. Click :doc:`here <../../../tools/postprocessing/postprocessing_pyvista>` for details.


-----------
References
-----------

.. [#ouertatani] \N. Ouertatani, N. Ben Cheikh, B. Ben Beya, T. Lili, "Numerical simulation of two-dimensional Rayleigh-Bénard convection in an enclosure," Comptes Rendus – Mec. 2008;336(5):464–70. `10.1016/j.crme.2008.02.004 <https://comptes-rendus.academie-sciences.fr/mecanique/articles/10.1016/j.crme.2008.02.004/>`_\.
.. [#venturi2010] \D. Venturi, X. Wan, and G. E. Karniadakis, “Stochastic bifurcation analysis of Rayleigh–Bénard convection,” *J. Fluid Mech.*, vol. 650, pp. 391–413, May 2010, doi: `10.1017/S0022112009993685 <https://doi.org/10.1017/S0022112009993685>`_\.
.. [#mpi2022] \“Rayleigh-Bénard Convection” *Max Planck Institute*, Accessed: 17 Jul. 2024, Available: https://archive.ph/XrJXx\.
57 changes: 27 additions & 30 deletions examples/multiphysics/rayleigh-benard-convection/rayleigh-benard-convection-Ra10k.prm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,16 +12,16 @@ set dimension = 2

subsection simulation control
set method = bdf1
set time end = 10000
set time end = 24
set time step = 0.01
set adapt = true
set max cfl = 0.5
set max cfl = 0.9
set stop tolerance = 1e-5
set adaptative time step scaling = 1.3
set number mesh adapt = 0
set output name = rayleigh-benard_convection
set output control = iteration
set output frequency = 100
set output control = time
set output time frequency = 0.5
set output path = ./output/
end

Expand Down Expand Up @@ -61,10 +61,10 @@ end
subsection physical properties
set number of fluids = 1
subsection fluid 0
set density = 100
set kinematic viscosity = 0.01
set thermal expansion = 0.0002
set thermal conductivity = 0.15625 # for Ra = 10000
set density = 1
set kinematic viscosity = 0.071
set thermal expansion = 0.71
set thermal conductivity = 10
set specific heat = 100
end
end
Expand All @@ -75,9 +75,9 @@ end

subsection mesh
set type = dealii
set grid type = subdivided_hyper_rectangle
set grid arguments = 2, 1: 0.5, 0.25 : 0 , 0 : true
set initial refinement = 6
set grid type = hyper_cube
set grid arguments = -0.5 : 0.5 : true
set initial refinement = 8
end

#---------------------------------------------------
Expand Down Expand Up @@ -117,8 +117,8 @@ end
subsection boundary conditions heat transfer
set number = 4
subsection bc 0
set id = 0
set type = convection-radiation-flux
set id = 0
set type = convection-radiation-flux
subsection h
set Function expression = 0
end
Expand All @@ -130,8 +130,8 @@ subsection boundary conditions heat transfer
end
end
subsection bc 1
set id = 1
set type = convection-radiation-flux
set id = 1
set type = convection-radiation-flux
subsection h
set Function expression = 0
end
Expand All @@ -143,15 +143,15 @@ subsection boundary conditions heat transfer
end
end
subsection bc 2
set id = 2
set type = temperature
set id = 2
set type = temperature
subsection value
set Function expression = 50
set Function expression = 10
end
end
subsection bc 3
set id = 3
set type = temperature
set id = 3
set type = temperature
subsection value
set Function expression = 0
end
Expand Down Expand Up @@ -179,9 +179,12 @@ subsection non-linear solver
set max iterations = 20
end
subsection fluid dynamics
set verbosity = verbose
set tolerance = 1e-8
set max iterations = 20
set solver = inexact_newton
set verbosity = verbose
set matrix tolerance = 0.1
set reuse matrix = true
set tolerance = 1e-8
set max iterations = 20
end
end

Expand All @@ -196,10 +199,7 @@ subsection linear solver
set max iters = 1000
set relative residual = 1e-6
set minimum residual = 1e-8
set preconditioner = ilu
set ilu preconditioner fill = 0
set ilu preconditioner absolute tolerance = 1e-14
set ilu preconditioner relative tolerance = 1.00
set preconditioner = amg
end
subsection heat transfer
set verbosity = verbose
Expand All @@ -208,8 +208,5 @@ subsection linear solver
set relative residual = 1e-6
set minimum residual = 1e-8
set preconditioner = ilu
set ilu preconditioner fill = 0
set ilu preconditioner absolute tolerance = 1e-14
set ilu preconditioner relative tolerance = 1.00
end
end
Loading

0 comments on commit c6777e7

Please sign in to comment.