Computation of the drag force for a rising bubble

[Giancarlo Esposito]

Hello everyone,

I'm currently simulating the rise of a bubble in a quiescent Newtonian fluid using Basilisk, and I would like to accurately compute the drag force acting on the bubble.

From basic principles, the drag should be obtained by integrating the traction vector n⋅T, where T is the total stress tensor defined as:

T=−pI+μ(∇u+∇uT)

Although the geometric VOF method gives us a sharp interface, I'm running into challenges when trying to compute stresses accurately at the interface, since pressure and velocity gradients are smeared across a few cells and are sensitive to how interfacial points are selected.

So far, I’ve tried two different approaches (in an axisymmetric setup for simplicity):

1. Direct surface integration: I loop over interfacial cells (0<f<1), compute the interface normal and area, then evaluate the stress tensor using Newton's law and integrate to get the total force.

2. Divergence theorem: I convert the surface integral into a volume integral of the divergence of the stress tensor over the bubble volume.

In steady-state, the drag force should balance the buoyancy force. In dimensionless terms (for a spherical bubble of unit radius), I expect the drag to equal 4/3π (meaning the drag force balances the buoyancy force), but my computed values are significantly off and strongly depend on the selected mesh, although the velocity of the centroid is pretty accurate with respect to previous results established in literature.

This leads to my main concerns:

• How can one best define the interface location and extract stress values robustly with VOF in Basilisk?

• Is one method (interface integration vs volume integral) more reliable or recommended in this context?

• Are there established practices, example cases, or built-in tools in Basilisk that might help with this kind of force computation?

Any guidance, suggestions, or code snippets would be greatly appreciated. I’m attaching my current version of the code in case it's helpful.

Thank you very much for your time!

Best regards,
Giancarlo

[Stephane Popinet]

Hi Giancarlo,

As you found out, this is not a simple question. All the options you
suggest are correct, whether they work however depend on the details of
their implementation. In principle "integral" methods using the
divergence theorem should be more robust that "point" differentiation.

One could also consider recovering the forces from the evolution of the
integrated (linear and angular) momentum of the bubble i.e. from the
linear and angular accelerations.

Note that directly integrating the pressure may cause problems as
described here:

http://basilisk.fr/src/navier-stokes/double-projection.h

for the case of a solid object.

Hope this helps somewhat,

Stephane

See also (again for a solid object):

http://basilisk.fr/src/embed.h#surface-force-and-vorticity
http://basilisk.fr/sandbox/ghigo/src/myembed.h#surface-force-and-torque-and-vorticity