Download this testcase.

2D Rayleigh-Benard instability test case

This test case simulates a 2D Rayleigh-Benard instability using the Finite Element Method (FEM) in a periodic domain.

Keywords

FEM, fluid, convection, 2D, heat transfer

Description

The test demonstrates: - How to perform a coupled fluid and heat transfer simulation. - How to set up boundary conditions for temperature and velocity fields. - How to set up non constant density.

import gmsh
import numpy as np
from migflow import fluid, advdiff
import shutil, os, sys

Output Directory

Create a clean output directory for simulation results.

outputdir = "output" if len(sys.argv) < 2 else sys.argv[1]
shutil.rmtree(outputdir, ignore_errors=True)
os.makedirs(outputdir)

Mesh generation

The domain is created through the python GMSH api is loaded to generate the mesh and extract the boundaries. The refinement is given by the function set_size_call_back. In this example, a constant mesh size is chosen.

w = 1
h = 0.25
lc = 0.005
gmsh.model.add("box")
gmsh.model.occ.add_rectangle(-w / 2, 0, 0, w, h)
gmsh.model.occ.synchronize()


def get_line(x0, x1, eps=1e-6):
    r = gmsh.model.get_entities_in_bounding_box(
        x0[0] - eps, x0[1] - eps, -eps, x1[0] + eps, x1[1] + eps, eps, 1
    )
    return list(tag for _, tag in r)


gmsh.model.add_physical_group(1, get_line([-w / 2, 0], [w / 2, 0]), name="Bottom")
gmsh.model.add_physical_group(1, get_line([-w / 2, h], [w / 2, h]), name="Top")
gmsh.model.add_physical_group(1, get_line([-w / 2, 0], [-w / 2, h]), name="Left")
gmsh.model.add_physical_group(1, get_line([w / 2, 0], [w / 2, h]), name="Right")
gmsh.model.add_physical_group(2, [1], name="domain")
transform = np.eye(4, 4, dtype=np.float64)
transform[0, -1] = -w
gmsh.model.mesh.set_periodic(
    1,
    get_line([-w / 2, 0], [-w / 2, h]),
    get_line([w / 2, 0], [w / 2, h]),
    transform.reshape(-1),
)
gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, mesh_size: lc)
gmsh.model.mesh.generate(2)

Fluid Problem

The fluid is described by its dimension, the external volume force applied, its dynamic viscosity and its density. Pardiso solver is chosen as a direct solver. As we target a natural convection, the density depends on the temperature. To do so, the density must be seen as a field and not as a constant.

g = np.array([0, -9.81])
rho = 1000
nu = 1e-6
mu = rho * nu
beta = 2e-4
f = fluid.FluidProblem(2, g, mu, density_element=b"triangle_p1")
f.load_msh(None)
f.set_wall_boundary("Top", velocity=[0, 0])
f.set_wall_boundary("Bottom", velocity=[0, 0])
f.set_mean_pressure(0)
f.interpolate(
    velocity_y=lambda x: 1e-5
    * np.sin(np.pi * x[:, 0] / w)
    * np.sin(np.pi * x[:, 1] / h)
)

Temperature Problem

The advection-diffusion problem is described by its dimension, its source term, its conductivity and its capacity.

cp = rho * 1.0
k = 1e-2
Tt = -10  # Top temperature
Tb = 10  # Bottom temperature
c = advdiff.AdvDiffProblem(
    2, 0, k, cp, velocity_element=b"triangle_p1"
)
c.load_msh(None)
c.set_dirichlet_boundary("Top", solution=Tt)
c.set_dirichlet_boundary("Bottom", solution=Tb)

Time integration

The numerical parameters is given and the initial conditions are written in the output directory.

dt = 1
tEnd = 500
outf = 1
t = 0
ii = 0


def get_fields(fluid, temp):
    y = fluid.coordinates_fields()[fluid.pressure_index()][:, 1]
    y = y.reshape(-1, 1)
    p1_element = fluid.get_p1_element()
    fields = fluid.get_default_export()
    fields["temperature"] = (temp.solution().reshape(-1, 1), p1_element)
    fields["density"] = (fluid.density().reshape(-1, 1), p1_element)
    return fields
    # return {
    #     "pressure": (fluid.pressure(), p1_element),
    # # "velocity": (fluid.velocity(), p1_element),
    # # "density": (fluid.density(), p1_element),
    # # "temperature": (temp.solution().reshape(-1, 1), p1_element),
    #     # "porosity": (fluid.porosity(), p1_element),
    # }


while t < tEnd:
    # Problem coupling
    # ----------------
    f.density()[:] = rho * (1 - beta * (c.solution().reshape(-1, 1) - (Tt + Tb) / 2))
    c.velocity()[:] = f.velocity()[:]
    if ii % outf == 0:
        print(f"ii : {ii} ----- t : {t}/{tEnd}")
        f.write_mig(outputdir, t, fields=get_fields(f, c))
    # Solve
    # -----
    c.implicit_euler(dt)
    f.implicit_euler(dt)
    t += dt
    ii += 1

Plot

python3 -m migflow.plot.migplot output --actors fluid --fluid-field density --fluid-vmin 998 --fluid-vmax 1002 --element-type triangle_p1