📚 Thermal porous convection (implicit and steady-state solution)
Porous convection in 2D
👉 get started with exercises
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 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 . We assume that this layer is saturated with fluid, i.e., the pore space is completely filled by fluid. We introduce the porosity , the volume fraction of material taken by pore space. The conservation of mass for the fluid requires:
Here is the density of the fluid, and is the fluid velocity. If the porous material is undeformable, i.e., and the fluid is incompressible, i.e., , the conservation of mass reduces to the following:
We define the quantity , 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 and the pressure gradient, known as the Darcy's law:
where is the proportionality coefficient called permeability, is the fluid viscosity, is pressure, 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 :
In the following, we assume for simplicity that the porosity is constant: .
Conservation of energy in the fluid results in the following equation for the temperature :
Here, is the specific heat capacity of the fluid, and is the conductive heat flux. Note that the time 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 is the thermal conductivity. Assuming , 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 and diffused with the diffusion coefficient
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 is the thermal expansion coefficient.
However, the equations were formulated for the incompressible case!
In convection problems, the gravity term makes the dominant contribution to the force balance. Therefore, variations in density due to thermal expansion are often accounted for only in the 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:
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 , and deviation of temperature from the reference value .
In addition, to reduce the number of independent variables in the code, instead of using the Fourier heat flux we use the temperature diffusion flux .
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, and are the characteristic relaxation times for pressure and heat diffusion, respectively, and 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, is the pseudo-compressibility, is the physical time step, and 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.
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
heatmap(xc, yc, C')
note the transpose
display() to force the display of the plot, e.g., in the time loop every
More advanced implementation, one can define the plotting options and apply them in the
opts = (aspect_ratio=1, xlims=(xc, xc[end]), ylims=(yc, 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
This can be achieved by typing entering the Pkg mode from the Julia REPL in the target folder
] activate . add Plots
and adding at least one package.
In addition, it is recommended to have the following structure and content:
Codes could be placed in the
scripts/ folder. Output material to be displayed in the
README.md could be placed in the
Manifest.tomlfile should be kept local. An automated way of doing so is to add it as entry to a
.gitignorefile in the root of your repo. Mac users may also add
⤴ back to Content
👉 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
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.
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
Use following definition for
θ_dτ = lx/re/cfl/dx β_dτ = (re*k_ηf)/(cfl*dx*lx)
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 .= diff(qDx)./dx err_Pf = maximum(abs.(r_Pf))
Convert this script to 2D. The script should produce a
heatmap() plot after the iterations converge. Start by renaming the main function
ly=20.0 and make sure the
dy cell size to be equal by selecting the appropriate number of grid points
ny in direction.
In the numerical parameters, take
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 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:
" 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 ass also an error monitoring step after the iteration loop as
"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
yc coordinates such that they span for direction, from
ly, define coordinates such the they span from
-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:
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 qDx 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,xc[end]),ylims=(yc,yc[end]),aspect_ratio=1,c=:turbo) display(quiver!(Xp[:], Yp[:], quiver=(qDx_p[:], qDy_p[:]), lw=0.5, c=:black)) end
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
The final convection animation you produced should be similar to the one displayed hereafter (using the parameters listed above):
⤴ back to Content
👉 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
θ_dτ should be replaced with
β_dτ should be replaced with
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
qTy are different from the arrays for the Darcy fluxes
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_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_Pf = 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_Pf,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)) (" 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
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
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