Lecture 4

Agenda
📚 Thermal porous convection (implicit and steady-state solution)
💻 Julia's Project environment
🚧 Exercises:

• Porous convection in 2D

• Rayleigh's number

Content

👉 get started with exercises

Modelling multi-physics

The goal of lecture 4 is to combine the acquired knowledge of numerical modelling and develop the first practical application: thermal porous convection in 2D.

Building upon:

• The diffusion equation

• Spatial discretisation: 1D and 2D

• Finite-differences and staggered grids

• Accelerated pseudo-transient method

What is convection and why we want to model it

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.

Thermal porous convection: a physical model

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:

Darcy's law

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.

Diffusion equation for pressure

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

Heat convection in porous media

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$:

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):

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)$

Modelling buoyancy: Boussinesq approximation

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:

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:

Solving thermal porous convection using the pseudo-transient method

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:

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:

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):

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.

Useful information

Initialise arrays

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:

2D plotting

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...))

Julia's Project environment

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.

⤴ back to Content

Exercises - lecture 4

Exercise 1 — Thermal porous convection in 2D

👉 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.

Getting started

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.

Task 1

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.

Task 2

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.

Task 3

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.

Task 4

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).

Task 5

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):

⤴ back to Content

Exercise 2 — Thermal porous convection with implicit temperature update

👉 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.

Task 1

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.

Task 2

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.

Task 3

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.

⤴ back to Content

Edit this page on
Last modified: December 01, 2024. Website built with Franklin.jl and the Julia programming language.