[vrohde]

Concerning the VOF method:

1) As I understand it, the surface tension force is modeled as a volumetric source term in the Navier-Stokes equations. This implies that each interface cell is associated with a single force at a given time. However, considering Newton’s third law, where is the corresponding reaction force? Shouldn’t there be a force exerted by one fluid and an equal and opposite response from the other across the interface — i.e., two forces instead of one? Am I misunderstanding something here, or is this a known limitation of the model?

2) While working with the VOF method, I noticed that the velocity field is often not tangent to the interface. As a result, streamlines can cross it, which seems physically questionable. Is this behavior to be expected?

Thank you,
VR

[jd113]

Originally Posted by vrohde

Concerning the VOF method:

1) As I understand it, the surface tension force is modeled as a volumetric source term in the Navier-Stokes equations. This implies that each interface cell is associated with a single force at a given time. However, considering Newton’s third law, where is the corresponding reaction force? Shouldn’t there be a force exerted by one fluid and an equal and opposite response from the other across the interface — i.e., two forces instead of one? Am I misunderstanding something here, or is this a known limitation of the model?

2) While working with the VOF method, I noticed that the velocity field is often not tangent to the interface. As a result, streamlines can cross it, which seems physically questionable. Is this behavior to be expected?

Thank you,
VR

You might have integrated the traction on the cylinder in flow past a cylinder to calculate lift and drag. So interface just represents a traction jump and that difference is the sigma*curvature*unit normal term that's used in the continuum surface force model in VOF or level set. I don't do FVM so I can't tell you about how VOF works in there but in FEM this part is easy to understand. You integrate by parts the divergence of stress term so a traction term gets thrown to the boundary (evaluation of definite integral which is one of the 2 terms you get after IBP). This term is part of the force vector or RHS which also contains the source terms. At the interface theres a jump of traction on both sides. So (traction1-traction2).n = surface tension x curvature.n. So you are saying that
the difference in tractions acting in normal to interface direction (ie normal stresses) are balancing the surface tension which is what's causing the interface to curve.

About the second question. Again no idea about FVM but in FEM traction= zero exists throughout the geometry on each element to element boundary . Rather the opposite normals cancel each other. And traction is always zero at inter element boundaries or else there is no flow. Similar to how zero Neumann or zero flux exists across the entire geomtry at interelement boundaries for a heat flow problem. Without this the transport cannot occur. So when you have an oil water interface they refuse to mix so the transport between them ceases that's why a difference in tractions will be finite here.

[vrohde]

Originally Posted by jd113

You might have integrated the traction on the cylinder in flow past a cylinder to calculate lift and drag. So interface just represents a traction jump and that difference is the sigma*curvature*unit normal term that's used in the continuum surface force model in VOF or level set. I don't do FVM so I can't tell you about how VOF works in there but in FEM this part is easy to understand. You integrate by parts the divergence of stress term so a traction term gets thrown to the boundary (evaluation of definite integral which is one of the 2 terms you get after IBP). This term is part of the force vector or RHS which also contains the source terms. At the interface theres a jump of traction on both sides. So (traction1-traction2).n = surface tension x curvature.n. So you are saying that
the difference in tractions acting in normal to interface direction (ie normal stresses) are balancing the surface tension which is what's causing the interface to curve.

About the second question. Again no idea about FVM but in FEM traction= zero exists throughout the geometry on each element to element boundary . Rather the opposite normals cancel each other. And traction is always zero at inter element boundaries or else there is no flow. Similar to how zero Neumann or zero flux exists across the entire geomtry at interelement boundaries for a heat flow problem. Without this the transport cannot occur. So when you have an oil water interface they refuse to mix so the transport between them ceases that's why a difference in tractions will be finite here.

Hey, JD

I do work with FVM rather than FEM, but your insights are very helpful.

Q1) OK, just to check if I got it right. Physically, we can think of 2 forces: one, F1, of the interface acting on fluid 1, and the other, F2, of the interface acting on fluid 2, so the net force at the interface is F = F2 - F1, which is nonzero for a curved interface and corresponds to the surface tension force. Now, in terms of the VOF-CSF model, these individual forces aren’t applied explicitly to the fluids. Instead, the interaction is handled implicitly by adding a term to the momentum equation that accounts for the net force F. That’s it, right?

Q2) I get your point about stress continuity in FEM and how zero-flux is implicitly enforced at inter-element boundaries. But in FVM/VOF, I'm still wondering — does the velocity field used to advect the interface allow for components normal to it? As far as I know, there’s no explicit enforcement of zero flux across the interface in VOF, so I'm still a bit unsure how this translates into interface behavior in dynamic cases. Also, there’s the known issue of spurious/parasitic currents. But when streamlines cross the interface, is that necessarily the reason? Take, for instance, a VOF simulation of a spray jet: it’s hard to imagine streamlines staying perfectly confined within each droplet, right? …

Again, thank you – your answer is much appreciated!

[jd113]

Originally Posted by vrohde

Hey, JD

I do work with FVM rather than FEM, but your insights are very helpful.

Q1) OK, just to check if I got it right. Physically, we can think of 2 forces: one, F1, of the interface acting on fluid 1, and the other, F2, of the interface acting on fluid 2, so the net force at the interface is F = F2 - F1, which is nonzero for a curved interface and corresponds to the surface tension force. Now, in terms of the VOF-CSF model, these individual forces aren’t applied explicitly to the fluids. Instead, the interaction is handled implicitly by adding a term to the momentum equation that accounts for the net force F. That’s it, right?

Q2) I get your point about stress continuity in FEM and how zero-flux is implicitly enforced at inter-element boundaries. But in FVM/VOF, I'm still wondering — does the velocity field used to advect the interface allow for components normal to it? As far as I know, there’s no explicit enforcement of zero flux across the interface in VOF, so I'm still a bit unsure how this translates into interface behavior in dynamic cases. Also, there’s the known issue of spurious/parasitic currents. But when streamlines cross the interface, is that necessarily the reason? Take, for instance, a VOF simulation of a spray jet: it’s hard to imagine streamlines staying perfectly confined within each droplet, right? …

Again, thank you – your answer is much appreciated!

1) Yes what you said is what is my understanding too. I haven't done FVM but from what I understand: Stokes flow would be volume integral of divergence of force/stress/traction= circular integral of force on cell faces = 0. So at interface that would be divergence of force = circular force integral on cell faces = non zero source (surface tension). But unlike FEM where this comes from weak form's integration by parts which allows you to set this force balance on any boundary of however curved or complex(so you consider the interface as internal boundary and model it with moving mesh)...here you smear out the interface and solve the force integral on face = source.
2) But in laminar flow the equation would become quite different circular integral of inertial flux+ force = source. Because you have the divergence of( rho u u) sort of term there. So guess that's an added thing to solve here. Maybe the flux reconstruction for this inertial flux is the issue. Or maybe it's not I don't know because I never coded this. In FEM that inertial term stays as volume integral and doesn't go to boundary. But it's the same equation so I can't imagine why streamlines cross the boundary in your case. My experience is not that high in CFD. I am mainly an experimentalist so I really don't understand beyond the very basics

[janesonn]

You're absolutely right in recognizing that surface tension is modeled as a volumetric force in VOF — typically via the Continuum Surface Force (CSF) model introduced by Brackbill et al. This approach smears the interface over a few cells and distributes the surface tension force volumetrically, proportional to the curvature and gradient of the volume fraction field.

[vrohde]

Originally Posted by janesonn

You're absolutely right in recognizing that surface tension is modeled as a volumetric force in VOF — typically via the Continuum Surface Force (CSF) model introduced by Brackbill et al. This approach smears the interface over a few cells and distributes the surface tension force volumetrically, proportional to the curvature and gradient of the volume fraction field.

Hi Janesonn,

Thank you for your reply!

Would you happen to have any insight on the second question — whether streamlines can cross the interface, given that VOF is a "one-fluid formulation"?

[jd113]

Originally Posted by vrohde

Hi Janesonn,

Thank you for your reply!

Would you happen to have any insight on the second question — whether streamlines can cross the interface, given that VOF is a "one-fluid formulation"?

Sorry this is a probably stupid to ask but have you minused the interface velocity before calculating streamlines. Because without that in droplet I used to get streamlines crossing at times. I don't know your test case

[vrohde]

Originally Posted by jd113

Sorry this is a probably stupid to ask but have you minused the interface velocity before calculating streamlines. Because without that in droplet I used to get streamlines crossing at times. I don't know your test case

I use a visualization software to generate streamlines from the raw velocity field.
(I work with bubbles/droplets, but I would like to hear from your experience.)

So, you're saying that you manually subtracted the interface velocity from the velocity field in order to calculate streamlines in the reference frame of the interface? (Please correct me if I’m wrong.)

How did you handle this, considering that each interface cell has a different velocity? Or, instead, did you make the velocity components normal to the interface equal to zero?

Thanks

[jd113]

Originally Posted by vrohde

I use a visualization software to generate streamlines from the raw velocity field.
(I work with bubbles/droplets, but I would like to hear from your experience.)

So, you're saying that you manually subtracted the interface velocity from the velocity field in order to calculate streamlines in the reference frame of the interface? (Please correct me if I’m wrong.)

How did you handle this, considering that each interface cell has a different velocity? Or, instead, did you make the velocity components normal to the interface equal to zero?

Thanks

Well when I worked with droplet/segmented flows I would subtract the average velocity of the droplet . In the droplets stationary frame of reference the streamlines would be perfectly curved at the interface . Otherwise they would cross the interface because the entire segments of both fluids were being advected one after another like a train. So while the fluid was turning backward at the interface the net motion was still ahead. So that was basically dependant on the ratio of how fast the drop moved and how fast recirculation happened within. If the circulation was weak almost always I saw streamlines pass straight through . This was in COMSOL though...

[vrohde]

Originally Posted by jd113

Well when I worked with droplet/segmented flows I would subtract the average velocity of the droplet . In the droplets stationary frame of reference the streamlines would be perfectly curved at the interface . Otherwise they would cross the interface because the entire segments of both fluids were being advected one after another like a train. So while the fluid was turning backward at the interface the net motion was still ahead. So that was basically dependant on the ratio of how fast the drop moved and how fast recirculation happened within. If the circulation was weak almost always I saw streamlines pass straight through . This was in COMSOL though...

OK, that was really helpful. I tried in a test case here, and I think it worked.

So:

(i) It is important to consider the referecence frame in which streamlines are computed. In a dynamic case, where the interface moves relatively to a fixed coordinate system (xOy), it is expected that streamlines (in the xOy frame) will cross the interface, according to the motion of the entire segment – just like the “train of fluid particles” you described.

(ii) You can shift to a reference frame centered at any point P of your domain. To do that, you compute streamlines based on the (u - uP) field. If P is a point on the interface, then streamlines won’t cross the interface at that specific point.

(iii) In the particular case of a bubble or droplet, if the internal recirculation << overall droplet motion, then you can expect streamlines (relative to xOy) to cross the interface.

(iv) If you take CM as the center of mass of the droplet, then (u - uCM) defines a reference frame stationary to the (CM of the) droplet. In that frame, streamlines typically do not cross the interface, since, for any point P on the interface, (uCM – uP) << uCM, meaning the droplet’s deformation is slow compared to its translational motion. (If the droplet does not deform, then streamlines will be perfectly confined within the droplet).

Would you agree with that?

Thanks again.

[jd113]

Originally Posted by vrohde

OK, that was really helpful. I tried in a test case here, and I think it worked.

So:

(i) It is important to consider the referecence frame in which streamlines are computed. In a dynamic case, where the interface moves relatively to a fixed coordinate system (xOy), it is expected that streamlines (in the xOy frame) will cross the interface, according to the motion of the entire segment – just like the “train of fluid particles” you described.

(ii) You can shift to a reference frame centered at any point P of your domain. To do that, you compute streamlines based on the (u - uP) field. If P is a point on the interface, then streamlines won’t cross the interface at that specific point.

(iii) In the particular case of a bubble or droplet, if the internal recirculation << overall droplet motion, then you can expect streamlines (relative to xOy) to cross the interface.

(iv) If you take CM as the center of mass of the droplet, then (u - uCM) defines a reference frame stationary to the (CM of the) droplet. In that frame, streamlines typically do not cross the interface, since, for any point P on the interface, (uCM – uP) << uCM, meaning the droplet’s deformation is slow compared to its translational motion. (If the droplet does not deform, then streamlines will be perfectly confined within the droplet).

Would you agree with that?

Thanks again.

Largely yes I agree with the train of thoughts. The final point is a bit confusing. I can't imagine that. Why would the drop's deformation due to velocity lead to a change in velocities from centre of mass? And could that get confused because the centre mass might change position based on the geomtry of the deformation itself?

Other factors like slip velocity might come in too. If it is zero and deformations are quite less then yes subtracting the centre of mass velocity from the interface might be a model to think of. But generally that's not the case I think for any real problem. One must know the Bond number, the pressure difference on the front and back of the drop, etc a lot of parameters.

Why not just integrate the total volumetric flow rate in the channel and subtract the velocity estimated from that from the interface velocity. I think that's what I did though it's 10 years ago I don't recall well. But this should take care of the slip ( even if phase 1 is slipping oast the phase2 due to gravity or pressure gradients the total volumetric flow rate will record everything)

[vrohde]

Originally Posted by jd113

Largely yes I agree with the train of thoughts. The final point is a bit confusing. I can't imagine that. Why would the drop's deformation due to velocity lead to a change in velocities from centre of mass? And could that get confused because the centre mass might change position based on the geomtry of the deformation itself?

Other factors like slip velocity might come in too. If it is zero and deformations are quite less then yes subtracting the centre of mass velocity from the interface might be a model to think of. But generally that's not the case I think for any real problem. One must know the Bond number, the pressure difference on the front and back of the drop, etc a lot of parameters.

Why not just integrate the total volumetric flow rate in the channel and subtract the velocity estimated from that from the interface velocity. I think that's what I did though it's 10 years ago I don't recall well. But this should take care of the slip ( even if phase 1 is slipping oast the phase2 due to gravity or pressure gradients the total volumetric flow rate will record everything)

Hey, jd113

From what you said, your main point concerns item (iv), which refers to the particular case of a droplet/bubble.

I was trying to give a reasonable (numerical) explanation for why streamlines would tend to remain confined within the droplet (i.e., why would that be “expected”). But you’re right – I didn’t consider the slip velocity. So let me try to generalize the idea a bit.

For a point P on the interface:

uP = uCM + uDef + uSlip

where uDef is the velocity contribution from the drop internal deformation; and uSlip is the slip velocity (relative motion between the two fluids at the interface).

When shifting to the CM frame, the relative velocity (P-CM) is:

(uP – uCM) = uDef + uSlip

So:

If (uDef + uSlip) is small (<< uCM), streamlines tend do stay confined within the droplet in the CM frame.

If (uDef + uSlip) is significant (~ uCM), streamlines can cross the interface, even in the CM frame.

I hadn’t considered the slip velocity because, in the VOF models I’m used to, there is no slip imposed at the interface (though numerical artifacts can ideed introduce some slip).

Regarding the other points your raised: I don’t think I have much to add. I acknowledge that the vol. flow rate balance you mentioned may be useful if one tries to estimate or correct for slip. As for the other parameters (Bond number, pressure difference, etc.): they are certainly important for the particle deformation, path, flow regime, etc. – but I don’t think they directly affect the (numerical) reasoning I was trying to build in item (iv).

Thanks again.

[jd113]

Originally Posted by vrohde

Hey, jd113

From what you said, your main point concerns item (iv), which refers to the particular case of a droplet/bubble.

I was trying to give a reasonable (numerical) explanation for why streamlines would tend to remain confined within the droplet (i.e., why would that be “expected”). But you’re right – I didn’t consider the slip velocity. So let me try to generalize the idea a bit.

For a point P on the interface:

uP = uCM + uDef + uSlip

where uDef is the velocity contribution from the drop internal deformation; and uSlip is the slip velocity (relative motion between the two fluids at the interface).

When shifting to the CM frame, the relative velocity (P-CM) is:

(uP – uCM) = uDef + uSlip

So:

If (uDef + uSlip) is small (<< uCM), streamlines tend do stay confined within the droplet in the CM frame.

If (uDef + uSlip) is significant (~ uCM), streamlines can cross the interface, even in the CM frame.

I hadn’t considered the slip velocity because, in the VOF models I’m used to, there is no slip imposed at the interface (though numerical artifacts can ideed introduce some slip).

Regarding the other points your raised: I don’t think I have much to add. I acknowledge that the vol. flow rate balance you mentioned may be useful if one tries to estimate or correct for slip. As for the other parameters (Bond number, pressure difference, etc.): they are certainly important for the particle deformation, path, flow regime, etc. – but I don’t think they directly affect the (numerical) reasoning I was trying to build in item (iv).

Thanks again.

Much appreciated 👍☺️
Actually I really loved the multiphase flow experiments from my past and I wrote so much because I felt like meeting a kindred spirit. I hope I didn't disturb or annoy you with my posts.

As for the slip, since my work was mainly experiments the interface is hard to be perfectly pure and always contains contaminants. So the tangential stress balance including slip and Marangoni terms was an everyday scene. It was really fantastic to watch the interfacial convection driven by mass transfer but this happened at really high concentration gradients which were nigh impossible to model 12 years back. DG wasn't so famous then (still don't have blackboxes for it ) but maybe that would have modelled these problems well.

In any case, all the best to your simulations and if you need any person to discuss these topics, you know I am always there 👍☺️