Rayleigh-Bénard convection¶

Natural convection in a fluid driven by heat differences. This example is inspired by numerical simulations of convection in the Earth mantle (see for example the ASPECT project).

We assume that viscosity is high (more precisely, the limit of infinite Prandtl number) such that we obtain a convection-diffusion equation for temperature coupled with the Stokes equations for fluid flow. Some references:

The unknowns are temperature $T$ , velocity $\mathbf u$ and pressure $p$ , and the problem being solved is

$\begin{aligned} \frac\partial{\partial t} T + \mathbf u \cdot \nabla T &= \Delta T, \\ -\Delta \mathbf u + \nabla p &= -\text{Ra} T \mathbf g, \\ \nabla\cdot \mathbf u &= 0, \end{aligned}$ where

$\text{Ra}$ is the Rayleigh number and

$\mathbf g$ is the vector describing earth acceleration.

The velocity $\mathbf u = \mathbf u(T)$ is fully determined by the temperature in each time step via solving a Stokes problem. Therefore implicit time stepping is performed for $T$ only, and the problem is linearized by freezing $\mathbf u = \mathbf u(T_{\text{old}})$ in each time step.

In [1]:

%pylab inline
import scipy
from IPython.display import HTML

from pyiga import bspline, assemble, geometry, solvers, vis, approx
from tqdm import tqdm

Populating the interactive namespace from numpy and matplotlib

In [2]:

# define semi-annulus geometry
geo = geometry.outer_product(geometry.line_segment(1, 2),
                             geometry.semicircle())

# define isogeometric Taylor-Hood discretization space
p = 3     # spline degree of velocity space
n_el = (30, 150)   # number of elements in y/x direction

# velocity space:    degree p,   continuity p-2
# pressure space:    degree p-1, continuity p-2
# temperature space: degree p,   continuity p-1
kvs_u = tuple(bspline.make_knots(p,   0.0, 1.0, n, mult=2) for n in n_el)
kvs_p = tuple(bspline.make_knots(p-1, 0.0, 1.0, n, mult=1) for n in n_el)
kvs_T = tuple(bspline.make_knots(p,   0.0, 1.0, n, mult=1) for n in n_el)

m_u = tuple(kv.numdofs for kv in kvs_u)
m_p = tuple(kv.numdofs for kv in kvs_p)

In [3]:

# assemble the Stokes system matrix
A_grad = assemble.assemble('inner(grad(u), grad(v)) * dx', kvs_u, bfuns=[('u',2), ('v',2)], geo=geo)
A_div  = assemble.assemble('div(u) * p * dx', (kvs_u, kvs_p), bfuns=[('u',2,0), ('p',1,1)], geo=geo)

A_stokes = scipy.sparse.bmat(
    [[A_grad, A_div.T],
     [A_div,  None   ]], format='csr')

In [4]:

# homogeneous Dirichlet BCs for the velocity
stokes_bcs = assemble.compute_dirichlet_bcs(kvs_u, geo, ('all', lambda *x: (0,0)))
# fix pressure to 0 at an arbitrary (the first) pressure dof
N = np.prod(m_u)
stokes_bcs = (np.append(stokes_bcs[0], 2*N), np.append(stokes_bcs[1], 0.0))

LS_stokes = assemble.RestrictedLinearSystem(A_stokes, 0, stokes_bcs)
stokes_solver = solvers.make_solver(LS_stokes.A)

In [5]:

def solve_stokes_with_rhs(Ra, T):
    # Ra: Rayleigh number
    # T:  temperature field
    rhs = Ra * assemble.assemble('T * inner(x/norm(x), v) * dx', kvs_u, bfuns=[('v',2)], geo=geo, T=T)
    rhs = concatenate((rhs.ravel(), zeros(prod(m_p))))   # pressure rhs = 0

    # solve Stokes system with the given boundary conditions
    LS = assemble.RestrictedLinearSystem(A_stokes, rhs, stokes_bcs)
    u = LS.complete(stokes_solver.dot(LS.b))
    U = moveaxis(u[:2*N].reshape((2,) + m_u), 0, -1)  # bring components to the back
    return geometry.BSplineFunc(kvs_u, U)

In [6]:

T_bcs = assemble.compute_dirichlet_bcs(kvs_T, geo, [('top', 0), ('bottom', 1)])

# slightly randomize the boundary conditions to break unphysical symmetry
T_bcs = (T_bcs[0], T_bcs[1] * exp(randn(len(T_bcs[1])) / 1000))

T_rhs = 0

T_LS = assemble.RestrictedLinearSystem(assemble.stiffness(kvs_T, geo), T_rhs, T_bcs)
T_M = T_LS.restrict_matrix(assemble.mass(kvs_T, geo))

In [7]:

# start with stationary solution of heat equation
T0_coef = T_LS.complete(solvers.make_solver(T_LS.A).dot(T_LS.b))
T0 = bspline.BSplineFunc(kvs_T, T0_coef)

In [8]:

vis.plot_field(T0, geo, cmap='jet')
colorbar();
axis('scaled');

In [9]:

def implicit_euler_step(T_old, M, A, b, tau):
    R1 = solvers.make_solver(M + tau * A)
    return R1 @ (tau * b + M @ T_old)

def crank_nicolson_step(T_old, M, A, b, tau):
    R0 = M - 0.5 * tau * A
    R1 = solvers.make_solver(M + 0.5 * tau * A)
    return R1 @ (tau * b + R0 @ T_old)

In [10]:

Ra = 1e6    # Rayleigh number
tau = 1e-4  # timestep

In [11]:

T = bspline.BSplineFunc(kvs_T, T0.coeffs.copy())
solutions = [T.coeffs.copy()]

In [12]:

time_step = crank_nicolson_step

A_diff = assemble.stiffness(kvs_T, geo=geo)
conv_asm = assemble.Assembler('inner(vel, grad(u)) * v * dx',
        kvs_T, geo=geo, vel=solve_stokes_with_rhs(Ra, T), updatable=['vel'])

for it in tqdm(range(200)):
    vel = solve_stokes_with_rhs(Ra, T)
    A = T_LS.restrict_matrix(A_diff + conv_asm.assemble(vel=vel))
    
    T_new = time_step(T_LS.restrict(solutions[-1].ravel()), T_M, A, T_LS.b, tau)
    
    T.coeffs.flat[:] = T_LS.complete(T_new)
    solutions.append(T.coeffs.copy())

100%|██████████| 200/200 [00:56<00:00,  3.54it/s]

In [13]:

sol_fields = [bspline.BSplineFunc(kvs_T, Tt) for Tt in solutions]

figsize(16,8)
vis.plot_field(sol_fields[-1], geo, cmap='jet', res=400)
colorbar();
axis('scaled');

In [14]:

%%time
figsize(8,5)
HTML(vis.animate_field(sol_fields, geo, cmap='jet', res=200, interval=70,
                      vrange=None, progress=True).to_html5_video())

202it [01:13,  2.74it/s]

CPU times: user 1min 13s, sys: 812 ms, total: 1min 14s
Wall time: 1min 14s

Out[14]:

In [ ]: