Download this testcase.

2D Convective cell

This test case simulates a 2D convective cell using the Finite Element Method (FEM).

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.

height = 1  # domain height
width = 1  # domain width
mesh_size = 0.015  # mesh size
origin = [0, 0]  # leftmost bottom corner


def gen_mesh(width, height, mesh_size, origin=np.array([0, 0])):
    origin = np.asarray(origin)
    gmsh.model.add("box")
    gmsh.model.occ.add_rectangle(origin[0], origin[1], 0, width, height)
    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 dim, tag in r)

    h, w = height, width
    gmsh.model.add_physical_group(
        1, get_line(origin + [0, 0], origin + [w, 0]), name="Bottom"
    )
    gmsh.model.add_physical_group(
        1, get_line(origin + [0, h], origin + [w, h]), name="Top"
    )
    gmsh.model.add_physical_group(
        1, get_line(origin + [0, 0], origin + [0, h]), name="Left"
    )
    gmsh.model.add_physical_group(
        1, get_line(origin + [w, 0], origin + [w, h]), name="Right"
    )
    gmsh.model.add_physical_group(2, [1], name="domain")
    gmsh.model.mesh.set_size_callback(lambda dim, tag, x, y, z, lc: mesh_size)
    gmsh.model.mesh.generate(2)


gen_mesh(width, height, mesh_size, origin)

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])  # gravity
rho = 1.184  # density
nu = 1.562e-5  # kinematic viscosity
mu = rho * nu  # dynamic viscosity
beta = 2e-4  # thermal expansion
cp = 1007  # heat capacity
k = 0.02551  # thermal conductivity
f = fluid.FluidProblem(2, g, mu, density_element=b"triangle_p1")
# The mesh can be loaded through the GMSH api if None is given or through a mesh filename.
# The boundary conditions are described by their physical name.
# To describle fully the pressure, its mean value over all the nodal values is impose to zero.
f.load_msh(None)
f.set_wall_boundary("Top")
f.set_wall_boundary("Bottom")
f.set_wall_boundary("Right")
f.set_wall_boundary("Left")
f.set_mean_pressure(0)
f.interpolate(
    velocity_y=lambda x: 1e-3
    * np.sin(np.pi * x[:, 0] / (2 * width))
    * np.sin(np.pi * x[:, 1] / (2 * height))
)

Temperature Problem

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

Tt = 0  # Top temperature
Tb = 1  # Bottom temperature
c = advdiff.AdvDiffProblem(
    2, 0, k, cp, velocity_element=b"triangle_p1"
)
c.load_msh(None)
c.set_weak_boundary("Top", solution=Tt)
c.set_weak_boundary("Bottom", solution=Tb)
c.set_weak_boundary("Left")
c.set_weak_boundary("Right")
x = c.coordinates_fields()
c.solution()[:] = Tb + (Tt - Tb) * (x[:, 1] - origin[1]) / height
Ra = np.abs(g[1] * beta * cp * rho / (nu * k) * (Tt - Tb) * height**3)
Pr = nu * rho * cp / k
print("Rayleigh number : Ra = %g    ---- Prandlt number : Pr = %g" % (Ra, Pr))

Time integration

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

dt = 1
tEnd = 1000
outf = 5
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


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, check_residual_norm=1)
    f.implicit_euler(dt, check_residual_norm=1)
    t += dt
    ii += 1

Plot

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