Simulate two phase flow (immiscible)

[Jupihawx]

Hi! First thanks for the library. I was wondering if it was possible to simulate a two phase flow, the same way that can be done through OpenFOAM with "interFoam" solver for example? If so, are there any examples available?

Thanks!

[holl-]

Hi, it's most likely possible but will require a small amount of coding. I'm not familiar with interFoam. Do you know how it works or have a reference?

[Jupihawx]

It seems that interfoam solves for two immiscible fluid, with local mesh property are related to the fluid fraction "alpha". It also takes into account the surface tension. The equations solved can be found here : https://openfoamwiki.net/index.php/InterFoam But as OpenFOAM is C++ based, I have a bit of trouble making the equivalence to this library. I will need to dive deeper into it I guess!

[holl-]

Okay, looks like they use a quantity alpha to identify the fraction of each liquid. They have gravity and surface tension acting on the fluid.

The following code should help you get started. It computes all the terms from the equations on grids but possibly in a different way (I havn't looked at the c++ implementation). This requires version 3.0 which you can get with pip install --upgrade git+https://github.com/tum-pbs/PhiFlow@3.0. Note that it does not treat the interface in any special way, so alpha gradually smears out over time.

from phi.flow import *
from phi.jax.flow import *

density1 = 2
density2 = 1
alpha = CenteredGrid(lambda x, y: y > 50, ZERO_GRADIENT, x=64, y=64, bounds=Box(x=100, y=100))
velocity = StaggeredGrid(Noise(), boundary=0, bounds=alpha.bounds, x=64, y=64) * .01

def curvature(a: Field):
    n = a.gradient()
    n = n.with_values(math.vec_normalize(n.values, epsilon=1e-5))  # normalize
    return - n.divergence()

@jit_compile
def step(a: Field, v: Field, dt=1., gravity=vec(x=0, y=-9.81), surface_tension_constant=10):
    a = advect.mac_cormack(a, v, dt)
    surface_tension = surface_tension_constant * curvature(a) * a.gradient()
    density = a * density1 + (1 - a) * density2
    force = resample(density * gravity + surface_tension, to=v)
    v = advect.semi_lagrangian(v, v, dt) + dt * force
    v, p = fluid.make_incompressible(v)
    return a, v

for i in range(100):
    alpha, velocity = step(alpha, velocity)
    if i % 5 == 0:
        show(alpha)

Hope this helps!

[Jupihawx]

Great, thank you so much for the template! I will look deeper into it and try to make it fit my usage case. Your example runs fine, but there is quite a lot smearing like you said, I will look into fixing that aswel!

[holl-]

By the way, when running without JIT-compilation, you can put a show(...) anywhere to quickly plot any quantity (like curvature, force, etc).