Agenda

📚 Thermal porous convection (implicit and steady-state solution)

💻 Julia's`Project`

environment

🚧 Exercises:

Porous convection in 2D

Rayleigh's number

The diffusion equation

Spatial discretisation: 1D and 2D

Finite-differences and staggered grids

Accelerated pseudo-transient method

Convection is a fluid flow driven by any instability arising from the interaction between the fluid properties such as density, and external forces such as gravity. If a layer of denser fluid lays on top of a layer of fluid with lower density, they will eventually mix and swap.

An example of such fluids would be oil and water. In thermal convection, the density difference is caused by the thermal expansion of the fluid, i.e., the dependence of density on temperature. Usually, higher temperatures correspond to lower densities.

Fluid flows in porous materials such as rocks and soil could also be a result of convection.

In this course, we only consider porous convection since it build on the already acquired knowledge of steady-state and transient diffusion processes.

Porous convection often arises in nature and geoengineering. For example, water circulation in hydrothermal systems is caused by thermal convection, and mixing of CO$_2$ with saline water during geological storage results from chemical convection.

In the following, we will introduce the equation governing the thermal porous convection, demonstrate that in the simple cases these equations reduce to the already familiar steady-state and transient diffusion equations.

Consider a layer of porous material of size $l_x \times l_y$. We assume that this layer is saturated with fluid, i.e., the pore space is completely filled by fluid. We introduce the *porosity* $\phi$, the volume fraction of material taken by pore space. The conservation of mass for the fluid requires:

Here $\rho$ is the density of the fluid, and $\boldsymbol{v}$ is the fluid velocity. If the porous material is undeformable, i.e., $\phi = \mathrm{const}$ and the fluid is incompressible, i.e., $\mathrm{d}\rho/\mathrm{d}t = \partial\rho/\partial t + \boldsymbol{v}\cdot\nabla\rho = 0$, the conservation of mass reduces to the following:

$\nabla\cdot(\phi\boldsymbol{v}) = 0$We define the quantity $\boldsymbol{q_D} = \phi\boldsymbol{v}$, which is called the *Darcy flux* or *Darcy velocity* and is the volumetric flow rate per unit area of the porous media. Standard approximation is the linear dependence between $\boldsymbol{q_D}$ and the pressure gradient, known as the *Darcy's law*:

where $k$ is the proportionality coefficient called *permeability*, $\eta$ is the fluid viscosity, $p$ is pressure, $\boldsymbol{g}$ is the gravity vector.

Substituting the Darcy's law into the mass conservation law for incompressible fluid, we obtain the steady-state diffusion equation for the pressure $p$:

$\nabla\cdot\left[\frac{k}{\eta}(\nabla p - \rho\boldsymbol{g})\right] = 0$In the following, we assume for simplicity that the porosity is constant: $\phi=\mathrm{const}$.

Conservation of energy in the fluid results in the following equation for the temperature $T$:

$\rho c_p \frac{\partial T}{\partial t} + \rho c_p\boldsymbol{v}\cdot\nabla T + \nabla\cdot\boldsymbol{q_F} = 0$Here, $c_p$ is the specific heat capacity of the fluid, and $\boldsymbol{q_F}$ is the conductive heat flux. Note that the time $t$ here is the physical time.

In a very similar manner to the Darcy's law, we relate the heat flux to the gradient of temperature (Fourier's law):

$\boldsymbol{q_F} = -\lambda\nabla T$where $\lambda$ is the *thermal conductivity*. Assuming $\lambda=\mathrm{const}$, substituting the Fourier's law into the energy equation, and using the definition of the Darcy flux we obtain:

This is the transient advection-diffusion equation. The temperature is advected with the fluid velocity $\boldsymbol{q_D}/\phi$ and diffused with the diffusion coefficient $\lambda/(\rho c_p)$

Now the transient pressure diffusion and steady pressure diffusion are coupled in a one-way fashion through the Darcy flux in the temperature equation. To model convection, we need to incorporate the dependency of the density on temperature in the equations. The simplest way is to introduce the linear dependency:

$\rho = \rho_0\left[1-\alpha (T-T_0)\right]$where $\alpha$ is the thermal expansion coefficient.

However, the equations were formulated for the incompressible case!

In convection problems, the gravity term $\rho\boldsymbol{g}$ makes the dominant contribution to the force balance. Therefore, variations in density due to thermal expansion are often accounted for only in the $\rho\boldsymbol{g}$ term and are neglected in other places. This is called the *Boussinesq approximation*.

Using the Boussinesq approximation, the incompressible equations remain valid, and the modified densities only appear in the definition of the Darcy flux:

$\boldsymbol{q_D} = -\frac{k}{\eta}(\nabla P - \rho_0\left[1-\alpha (T-T_0)\right]\boldsymbol{g})$We already discussed how the steady-state and transient equations could be solved efficiently by adding the pseudo-transient terms to the governing equations.

Let's apply this strategy to solve the thermal porous convection!

The thermal porous convection is a coupled system of equations describing the evolution of pressure and temperature. To simplify the equations, we solve not for the absolute values for pressure and temperature, but for the deviation of pressure from hydrostatic gradient $\int_{z}\rho g\,dz$, and deviation of temperature from the reference value $T_0$.

In addition, to reduce the number of independent variables in the code, instead of using the Fourier heat flux $\boldsymbol{q_F}$ we use the temperature diffusion flux $\boldsymbol{q_T}=\boldsymbol{q_F}/(\rho_0 c_p)$.

With these reformulations in mind, the full system of equations to solve is:

$\boldsymbol{q_D} = -\frac{k}{\eta}(\nabla p - \rho_0\alpha\boldsymbol{g}T)$ $\nabla\cdot\boldsymbol{q_D} = 0$ $\boldsymbol{q_T} = -\frac{\lambda}{\rho_0 c_p}\nabla T$ $\frac{\partial T}{\partial t} + \frac{1}{\phi}\boldsymbol{q_D}\cdot\nabla T + \nabla\cdot\boldsymbol{q_T} = 0$We formulate the pseudo-transient system of equations by augmenting the system with pseudo-physical terms. We add inertial terms to the Darcy and temperature diffusion fluxes:

$\theta_D\frac{\partial \boldsymbol{q_D}}{\partial\tau} + \boldsymbol{q_D} = -\frac{k}{\eta}(\nabla p - \rho_0\alpha\boldsymbol{g}T)$ $\theta_T\frac{\partial \boldsymbol{q_T}}{\partial\tau} + \boldsymbol{q_T} = -\frac{\lambda}{\rho_0 c_p}\nabla T$Here, $\theta_D$ and $\theta_T$ are the characteristic relaxation times for pressure and heat diffusion, respectively, and $\tau$ is the pseudo-time.

Then, we add the pseudo-compressibility to the mass balance equation. For each physical time step we discretise the physical time derivative and add the pseudo-time derivative (dual-time method):

$\beta\frac{\partial p}{\partial\tau} + \nabla\cdot\boldsymbol{q_D} = 0$ $\frac{\partial T}{\partial \tau} + \frac{T-T_\mathrm{old}}{\mathrm{d}t} + \frac{1}{\phi}\boldsymbol{q_D}\cdot\nabla T + \nabla\cdot\boldsymbol{q_T} = 0$Here, $\beta$ is the pseudo-compressibility, $\mathrm{d}t$ is the physical time step, and $T_\mathrm{old}$ is the distribution of temperature at the previous physical time step.

This new system of equations is amendable to the efficient solution by the pseudo-transient method. We'll implement the thermal porous convection solver in 2 stages:

in the first stage, we'll program the efficient elliptic solver for the pressure, leaving the temperature update explicit, and;

in the second stage, we'll make the temperature (advection-diffusion) solver also implicit.

Refer to the exercises for the steps to take to implement the porous convection solver.

Recall that we are using conservative finite-differences and thus a *staggered grid*.

For 2D grids, we will have to handle scalar quantity and two fluxes as depicted below, taking care about correct staggering:

You can use `heatmap()`

function from `Plots.jl`

, to plot e.g. `C`

as function of the spatial coordinates `xc`

and `yc`

:

`heatmap(xc, yc, C')`

*note the transpose '*

Use `display()`

to force the display of the plot, e.g., in the time loop every `nout`

.

More advanced implementation, one can define the plotting options and apply them in the `heatmap()`

call:

```
opts = (aspect_ratio=1, xlims=(xc[1], xc[end]), ylims=(yc[1], yc[end]), clims=(0.0, 1.0), c=:turbo, xlabel="Lx", ylabel="Ly", title="time = $(round(it*dt, sigdigits=3))")
display(heatmap(xc, yc, C'; opts...))
```

On GitHub, make sure to create a new folder for each week's exercises.

Each week's folder should be a Julia project, i.e. contain a `Project.toml`

file.

This can be achieved by typing entering the Pkg mode from the Julia REPL in the target folder

```
julia> ]
(@v1.10) pkg> activate .
(lectureXX) pkg> add Plots
```

and adding at least one package.

In addition, it is recommended to have the following structure and content:

lectureXX

`README.md`

`Project.toml`

`Manifest.toml`

docs/

scripts/

Codes could be placed in the `scripts/`

folder. Output material to be displayed in the `README.md`

could be placed in the `docs/`

folder.

💡 Note

The

`Manifest.toml`

file should be kept local. An automated way of doing so is to add it as entry to a `.gitignore`

file in the root of your repo. Mac users may also add `.DS_Store`

to their `.gitignore`

⚠️ Warning!

Exercises have to be handed in as monolithic Julia scripts (one code per script) and uploaded to your private (shared) GitHub repository, in a **specific folder for each lecture**. The git commit hash (or SHA) of the final push needs to be uploaded on Moodle (more).

👉 See Logistics for submission details.

The goal of this exercise is to:

Implement coupled diffusion equations in 2D

Consolidate the implicit updates and dual timestepping

In this first exercise, you will finalise the thermal porous convection discussed and implemented in class. The following tasks combine the major steps needed to program 2D thermal porous convection starting from 1D steady state diffusion.

Create a new folder named `lecture4`

in your GitHub repository for this week's (lecture 4) exercises. In there, create a new Julia script named `porous_convection_2D.jl`

for this homework. Take the 1D steady diffusion script `l3_steady_diffusion_1D.jl`

as a basis.

Rename variables so that we solve it for the pressure:

`C`

should now be replaced by`Pf`

`qx`

becomes`qDx`

and should be initialised with size`nx+1`

. Modify the rest of the code accordingly such that entire`Pf`

array can be updated having the boundary condition set on the flux.`dc`

becomes`k_ηf`

which is now the permeability over fluid viscosity`dτ./(ρ*dc .+ dτ)`

from flux update becomes`1.0./(1.0 + θ_dτ)`

`dτ`

from pressure update becomes`β_dτ`

. Note that instead of multiplying by`dτ`

one should divide by`β_dτ`

Use following definition for `θ_dτ`

and `β_dτ`

:

```
θ_dτ = lx/re/cfl/dx
β_dτ = (re*k_ηf)/(cfl*dx*lx)
```

and introduce `cfl = 1.0/sqrt(2.1)`

in the numerics section. Also, `ρ`

is no longer needed and we can move `re`

to numerics section as well.

Finally, update residual and error checking calculation as following (initialising `r_Pf`

accordingly)

```
r_Pf .= diff(qDx)./dx
err_Pf = maximum(abs.(r_Pf))
```

renaming `err`

to `err_Pf`

.

Convert this script to 2D. The script should produce a `heatmap()`

plot after the iterations converge. Start by renaming the main function `porous_convection_2D`

.

Use `lx=40.0`

and `ly=20.0`

and make sure the `dx`

and `dy`

cell size to be equal by selecting the appropriate number of grid points `ny`

in $y$ direction.

In the numerical parameters, take `min(dx,dy)`

, `max(lx,ly)`

and `max(nx,ny)`

whenever possible.

Make the initial condition a Gaussian distribution in 2D; you can use following

`@. exp(-(xc-lx/4)^2 -(yc'-ly/4)^2)`

taking care not to omit the transpose `'`

.

Also, initialise the flux in the $y$ direction `qDy`

. Update the residual calculation to account for the 2D setup.

Add a line to print out the convergence results within the `ncheck`

block. You can take inspiration from:

`@printf(" iter/nx=%.1f, err_Pf=%1.3e\n",iter/nx,err_Pf)`

Finally, remove the error saving procedure used to plot convergence.

Wrap the iteration loop into a time loop. Make `nt=10`

time steps. Move visualisation part out of the iteration loop and add also an error monitoring step after the iteration loop as

`@printf("it = %d, iter/nx=%.1f, err_Pf=%1.3e\n",it,iter/nx,err_Pf)`

to verify convergence.

Add new fields for the temperature evolution (advection and diffusion) using a fully explicit update approach outside of the iteration loop, right before visualisation. Don't forget to use an upwind scheme!

Use following time step definition (and compute it right before the temperature update)

```
dt_adv = ϕ*min(dx/maximum(abs.(qDx)), dy/maximum(abs.(qDy)))/2.1
dt = min(dt_diff,dt_adv)
```

The temperature update part could contain one update for the diffusion-related process and one for advection:

```
T[2:end-1,2:end-1] .+= dt.*λ_ρCp.*(...)
T[2:end-1,2:end-1] .-= dt./ϕ.*(...)
```

Make sure to add following to the initialisation

```
# physics
αρgx,αρgy = 0.0,1.0
αρg = sqrt(αρgx^2+αρgy^2)
ΔT = 200.0
ϕ = 0.1
Ra = 100
λ_ρCp = 1/Ra*(αρg*k_ηf*ΔT*ly/ϕ) # Ra = αρg*k_ηf*ΔT*ly/λ_ρCp/ϕ
# numerics
dt_diff = min(dx,dy)^2/λ_ρCp/4.1
```

Implement initial and boundary conditions; Remove the fluid pressure perturbation and initialise temperature array `T`

as following taking care of setting upper and lower boundary initial conditions as well; heating from below and cooling from above.

```
T = @. ΔT*exp(-xc^2 - (yc'+ly/2)^2)
T[:,1] .= ΔT/2; T[:,end] .= -ΔT/2
```

Left and right boundaries of the temperature field should have an adiabatic condition, i.e.,

`T[[1,end],:] .= T[[2,end-1],:]`

In addition, also centre the `xc`

and `yc`

coordinates such that they span for $x$ direction, from `-lx/2+dx/2`

to `lx/2-dx/2`

. For `ly`

, define coordinates such the they span from `-ly+dy/2`

to `-dy/2`

(having the 0 at the upper surface).

Add two-way coupling using the Boussinesq approximation, i.e., add the dependence of density on temperature in the Darcy flux. Produce the animation displaying the temperature evolution including arrows (quiver plot) for velocities.

Change parameters as following: `nt=500`

, `nx=127`

and `Ra=1000.0`

.

For visualisation, embed the plotting into a `nvis`

statement setting `nvis=5`

. Below you'll find a sample code for visualisation:

```
# visualisation
if it % nvis == 0
qDxc .= # average qDx in x
qDyc .= # average qDy in y
qDmag .= sqrt.(qDxc.^2 .+ qDyc.^2)
qDxc ./= qDmag
qDyc ./= qDmag
qDx_p = qDxc[1:st:end,1:st:end]
qDy_p = qDyc[1:st:end,1:st:end]
heatmap(xc,yc,T';xlims=(xc[1],xc[end]),ylims=(yc[1],yc[end]),aspect_ratio=1,c=:turbo)
display(quiver!(Xp[:], Yp[:], quiver=(qDx_p[:], qDy_p[:]), lw=0.5, c=:black))
end
```

where `Xp`

and `Yp`

could be initialised as following before the time loop:

```
# visualisation init
st = ceil(Int,nx/25)
Xc, Yc = [x for x=xc, y=yc], [y for x=xc,y=yc]
Xp, Yp = Xc[1:st:end,1:st:end], Yc[1:st:end,1:st:end]
```

Well done 🚀 - you made it. Add the produced gif or animation to the `README.md`

within your `lecture4`

folder.

The final convection animation you produced should be similar to the one displayed hereafter (using the parameters listed above):

👉 See Logistics for submission details.

The goal of this exercise is to:

Implement implicit advection-diffusion and dual-timestepping

Build your intuition about convection and the Rayleigh number

In this exercise you will implement the fully implicit and fully coupled solver for the thermal porous convection problem. Starting from the working solver that uses the explicit update for the temperature, you will introduce the pseudo-transient parameters for the implicit transient diffusion problem, and move the temperature update to the iteration loop. Then you will introduce the Rayleigh number that characterises the intensity of the buoyancy-driven convection, and verify that your numerical code confirms the analytical predictions for the critical Rayleigh number separating the heat diffusion- and advection-dominated flow regimes.

Copy the file `porous_convection_2D.jl`

and name it `porous_convection_implicit_2D.jl`

. Rename the pseudo-transient variables for fluid pressure diffusion to avoid conflicts with the variables for temperature:

`re`

should be replaced with`re_D`

`θ_dτ`

should be replaced with`θ_dτ_D`

`β_dτ`

should be replaced with`β_dτ_D`

Adjust the value of `re_D`

since the physics is now fully coupled:

`re_D = 4π`

Move the time step definition into the beginning of the time loop. For the first time step, use different definition to avoid division by 0. For the time steps > 1, choose among the minimum between scale or flux based definition:

```
for it = 1:nt
T_old .= T
# time step
dt = if it == 1
0.1*min(dx,dy)/(αρg*ΔT*k_ηf)
else
min(5.0*min(dx,dy)/(αρg*ΔT*k_ηf),ϕ*min(dx/maximum(abs.(qDx)), dy/maximum(abs.(qDy)))/2.1)
end
...
end
```

Introduce the pseudo-transient parameters for the temperature update. Recall that the temperature evolution equation is equivalent to the diffusion-reaction equation with advection. Now, the physical timestep `dt`

is determined from the CFL condition and changes every iteration of the time loop. Thus, the pseudo-transient parameters should also be updated every time step:

```
# time step
# dt = ...
re_T = π + sqrt(π^2 + ly^2/λ_ρCp/dt)
θ_dτ_T = max(lx,ly)/re_T/cfl/min(dx,dy)
β_dτ_T = (re_T*λ_ρCp)/(cfl*min(dx,dy)*max(lx,ly))
...
```

Add new arrays to the `# array initialisation`

section to store the physical time derivative of temperature, temperature equation residual, and the temperature diffusion fluxes:

```
dTdt = zeros(nx-2,ny-2)
r_T = zeros(nx-2,ny-2)
qTx = zeros(nx-1,ny-2)
qTy = zeros(nx-2,ny-1)
```

Note that the sizes of the arrays `qTx`

and `qTy`

are different from the arrays for the Darcy fluxes `qDx`

and `qDy`

. The reason for this is that we use the different boundary conditions for the temperature, and don't want to update the temperature at the domain boundaries.

Move the temperature update into the iteration loop. Rename the variable `err_Pf`

to `err_D`

to avoid confusion. Introduce the new variable `err_T`

to store the residual for the temperature evolution equation and modify the exit criteria to break iterations when both errors are less than tolerance:

```
# iteration loop
iter = 1; err_D = 2ϵtol; err_T = 2ϵtol
while max(err_D, err_T) >= ϵtol && iter <= maxiter
...
end
```

Annotate the Darcy fluxes and pressure update with a comment, and introduce the new section for temperature update:

```
while max(err_D, err_T) >= ϵtol && iter <= maxiter
# fluid pressure update
qDx[2:end-1,:] .-= ...
qDy[:,2:end-1] .-= ...
Pf .-= ...
# temperature update
...
end
```

Add the temperature diffusion flux update analogous the the Darcy flux update, but using the different iteration parameters:

```
# temperature update
qTx .-= ...
qTy .-= ...
```

Compute the material physical time derivative as a combination of partial derivative `(T - T_old)./dt`

and upwind advection:

`dTdt .= (T[2:end-1,2:end-1] .- T_old[2:end-1,2:end-1])./dt .+ (...)./ϕ`

The upwind advection part could be simply copied from the previous explicit version, ignoring the `dt`

factor. Finally, compute the temperature update and move the boundary conditions to the iteration loop:

```
T[2:end-1,2:end-1] .-= (dTdt .+ ...)./(1.0/dt + β_dτ_T)
T[[1,end],:] .= T[[2,end-1],:]
```

Add the residual calculation for the temperature evolution equation and the iteration progress reporting:

```
if iter % ncheck == 0
r_Pf .= ...
r_T .= dTdt .+ ...
err_D = maximum(abs.(r_Pf))
err_T = maximum(abs.(r_T))
@printf(" iter/nx=%.1f, err_D=%1.3e, err_T=%1.3e\n",iter/nx,err_D,err_T)
end
```

Run the code, make sure that it works as expected, produce the animation and add it to the `README.md`

within your `lecture4`

folder. Well done! 🔥

Did the number of iterations required for convergence change compared to the version with the explicit temperature update? Try to come up with the explanation for why the number of iterations changed the way it changed and write a sentence about your thoughts on the topic.

Using the newly developed implicit code, realise a numerical experiment varying the Rayleigh number. Theoretical critical value of `Ra`

above which there is convection is approximately `40`

. Confirm that `Ra < 40`

results in no convection. Confirm that the values of `Ra > 40`

result in the development of convection. Try the following range of values for `Ra`

: `10`

, `40`

, `100`

, `1000`

. Produce the animation or the final figure after `nt=100`

timesteps for each value. Add the produced gif or animation to the `README.md`

within your `lecture4`

folder.

**Question:** What is the difference in the results for the different values of `Ra`

, is there an observable trend? Write a comment explaining your observations.

Last modified: October 08, 2024. Website built with Franklin.jl and the Julia programming language.