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
+ f.inner(D_gc * rho_g * f.grad(phi_g), f.grad(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 * mu_l * f.inner(f.grad(u_l), f.grad(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