# Modelling bubbly flow

Hi all,

I’d like to implement the set of equations of COMSOL’s “bubbly flow”.

I want to simulate a tank with bubbles sparging at the bottom. Basically reproducing this example: https://www.comsol.com/blogs/how-to-model-sparging-in-comsol-multiphysics/

The set of equations is available here: The Bubbly Flow Equations

This is the code i have so far:

``````import fenics as f

width = 0.5
height = 1
mesh = f.RectangleMesh(f.Point(0, 0), f.Point(width, height), 40, 100)

V_ele = f.VectorElement("CG", mesh.ufl_cell(), 2)
Q_ele = f.FiniteElement("CG", mesh.ufl_cell(), 1)
W = f.FunctionSpace(mesh, f.MixedElement([V_ele, Q_ele, Q_ele]))

upphi = f.Function(W)
upphi_n = f.Function(W)

vu, vp, vphi = f.TestFunctions(W)

u_l, p, phi_g = f.split(upphi)
u_l_n, p_n, phi_g_n = f.split(upphi_n)

phi_l = 1 - phi_g
phi_l_n = 1 - phi_g_n

D_gc = 0.0001  # drift diffusion coeff

u_slip = f.Constant((0.0, 0.0))
u_drift = -D_gc * f.grad(phi_g) / (phi_g + 1e-16)
u_g = u_l + u_drift + u_slip

rho_l = 10
rho_g = 0.01  # TODO see COMSOL guide for actual expression
mg_l = f.Constant(0.0)  # TODO mass transfer rate
mu_l = 10  # dynamic visosity

g = f.Constant((0.0, 9.81))
F = f.Constant((0.0, 0.0))

dt = f.Constant(10)

continuity = (
rho_l * (phi_l - phi_l_n) + rho_g * (phi_g - phi_g_n)
) / dt * vp * f.dx + f.inner(
f.div(rho_l * phi_l * u_l + rho_g * phi_g * u_g), vp
) * f.dx

gas_phase_transport = (
(rho_g * (phi_g - phi_g_n)) / dt * vphi * f.dx
+ f.div(rho_g * phi_g * (u_l + u_slip)) * vphi * f.dx
+ mg_l * vphi * f.dx
)

momentum = (
phi_l * rho_l * f.inner((u_l - u_l_n) / dt, vu) * f.dx
+ phi_l * rho_l * f.dot(f.dot(f.grad(u_l), u_l), vu) * f.dx
+ p * f.div(vu) * f.dx
- phi_l * rho_l * f.inner(g, vu) * f.dx
- f.inner(F, vu) * f.dx
)

F = momentum + continuity + gas_phase_transport

# Construct facet markers
bndry = f.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
for facet in f.facets(mesh):
mp = facet.midpoint()
if f.near(mp[1], 0.0):  # bottom wall
bndry[facet] = 1
if f.near(mp[1], 0.0) and width * 0.45 < mp[0] < width * 0.55:  # bottom inlet
bndry[facet] = 2
elif f.near(mp[1], height):  # top
bndry[facet] = 3
elif f.near(mp[0], 0.0) or f.near(mp[0], width):  # walls
bndry[facet] = 1

f.XDMFFile("facets.xdmf").write(bndry)

# boundary conditions

non_slip = f.DirichletBC(W.sub(0), f.Constant((0.0, 0.0)), bndry, 1)
pressure_top = f.DirichletBC(W.sub(1), f.Constant(1), bndry, 3)
inlet_gas = f.DirichletBC(W.sub(2), f.Constant(1e-1), bndry, 2)
bcs = [non_slip, pressure_top, inlet_gas]

file_velocity = f.XDMFFile("velocity.xdmf")
file_pressure = f.XDMFFile("pressure.xdmf")
flie_phi = f.XDMFFile("gas_fraction.xdmf")
t = 0
while t < 100:
f.solve(F == 0, upphi, bcs)
u, p, phi_g = upphi.split()
file_velocity.write(u, t)
file_pressure.write(p, t)
flie_phi.write(phi_g, t)
upphi_n.assign(upphi)
t += float(dt)
``````

For some reason, the liquid velocity seems very poorly affected by the bubbles, which are therefore not advected.

I’m almost certain I messed something up in the variational formulation.

Any idea?

Thanks!

Remi

After a bit of investigation, it may be that I’m missing a slip boundary condition at the top surface.

I thought it was the default “do-nothing” condition but I’m probably wrong.