wetting with free surface or large density ratio

[grasingerm]

Hello all,

I’ve been working with the free energy model as I’m very interested in studying the influence of surface energies and wetting on the filling of molds, cavities, etc. I feel comfortable with using setFreeEnergyWallBoundary for this and have been able to reproduce some standard contact angle results, as well as simulate some other interesting problems. However, the model assumes equal densities of all fluids–which, for my case, means equal densities between the fluid filling the mold and the fluid being displaced. Unfortunately, this isn’t realistic for my interests as the fluid being displaced is ‘air’.

I’ve also had some success modelling mold/cavity filling with the free surface model; however, I wasn’t sure how to include wetting effects in those models.

How would you recommend I model wetting at solid boundaries with a free surface flow model (or, alternatively, a multi-component model at large density ratio)? Can setFreeEnergyWallBoundary be made to work with the free surface model? Can setFreeEnergyWallBoundary be used with a large density ratio Shan and Chen model? Or would it be best to try something else entirely?

Warm regards,
– Matt

[savis]

Hello Matt,

Unfortunately, the current implementation of the free energy model only considers equal densities. A non-equal model is something that we are interested in implementing in the future, but is not currently included.

It is worth considering if the problem actually requires a high density ratio. If the Reynolds number of the system is small (Re<<1) then an equal density model may be able to accurately simulate the system so long as the correct viscosities are assigned. Although, if inertia is important (Re>=1) then the density ratio must be accurate.

Regarding the setFreeEnergyWallBoundary function, it can only be used in conjunction with the free energy model. I don’t believe there are any similar wetting conditions implemented yet for the free surface or Shan-Chen model.

For Shan-Chen, a wetting condition could be implemented using an additional solid-liquid interaction (for example, see this paper: https://doi.org/10.1103/PhysRevE.76.066701). As for the density ratio, I believe there is a working prototype with density ratios of up to 1000, but this is not fully implemented yet.

And of course, all OpenLB models can be adapted and extended. In the framework of a common project we are happy to help with this!

Best wishes,
Sam

[grasingerm]

Hi Sam,

Thanks for the reply.

Assuming I understand correctly, the free energy wall boundary condition works by 1) formulating a contribution to the free energy functional which depends on the order parameter at the boundary and 2) taking the negative gradient of the free energy as a kind of thermodynamic force. It seems to me that something very similar could be done for the free surface model where instead (of the order parameter) the “thermodynamic force” depends on the local mass.

I would be happy to try to implement this, test it using an analogous contact angle simulation, and share my code if successful. That being said, I would also welcome any advice on how to best get started. I understand that what I most likely need to implement is a “postProcessor”. I’ve had some difficulty tracing back through the call stack of setFreeEnergyWallBoundary and understanding exactly what is going on there. Are there any additional examples, documentation, references, etc., that you think may help me in getting up to speed on postProcessors and how to use one to simulate wetting in a free surface model? I understand C++ and LBM theory well enough but I am new to OpenLB. I’m glad the community here is consistently responsive and helpful.

Warm regards,
– Matt

[savis]

Hi Matt,

Yes, that just about sums up the free energy boundary condition. Crucially, the thermodynamic force is already calculated on the bulk fluid nodes from gradients of the order parameter. So, all that is needed to implement the force at the solid surface is to modify the order parameter slightly on the boundary nodes. I am not that familiar with the free surface model, so I am not sure if the same approach would work or if you would need to assign a force itself.

To give you an idea how to implement a boundary condition I will summarise how setFreeEnergyWallBoundary works, as that’s what I’m most familiar with.

Firstly, there are a variety of different ‘setFreeEnergyWallBoundary’ functions that can be called depending on different parameters, but ultimately they call this function (for 2D): https://www.openlb.net/DoxyGen/html/db/d1b/namespaceolb.html#a385301cd3914ac9a54cf782e9299046f

This calls a different ‘setFreeEnergyWallBoundary’ function on each processor / block:
https://www.openlb.net/DoxyGen/html/db/d1b/namespaceolb.html#a1441af06e059821ec0dee3a9639a249b
In here, for each solid boundary node, the normal vector to the boundary is found (discreteNormal), the dynamics are set (bounce-back, etc), and then the post-processor generators are assigned. For the free-energy boundary there are two post-processors, but you would probably just need one. Also, the ‘addend’ and ‘latticeNumber’ parameters are specific to this boundary condition, so you would not need them.

The post-processor generator is largely a wrapper for the post-processor. It has been a while since I have looked at it, so I can’t quite remember its purpose: https://www.openlb.net/DoxyGen/html/d3/d21/classolb_1_1FreeEnergyWallProcessorGenerator2D.html

Finally, the post processor class contains the main functionality inside the ‘subprocess’ function: https://www.openlb.net/DoxyGen/html/d3/dee/classolb_1_1FreeEnergyWallProcessor2D.html
To clarify what is going on in this example, it takes the density / order parameter (somewhat confusingly called rho) on the adjacent node and adds the modification to get the desired value. The 0th distribution function is then modified on the boundary node so that the sum of the distribution functions gives this value.

Much of this code is boilerplate so to make a new boundary condition not much would need to be changed. The main things to modify are the second function which you would need the correct dynamics and postprocessor. As for the postprocessor, you would need to modify the subprocess function.

Let me know if any of this is unclear or if you need any further help.

Best,
Sam

[grasingerm]

Thanks, Sam. This is indeed very helpful.

[GASciences]

Hello,

First of all, thank you for both of you for your discussion, it made me progress a lot on understanding how to customize boundary conditions.

I have a couple of question though:

– @Matt, have you progressed since your last message?
– @Sam, you provide a reference in your first reply that suggests to add and solid/fluid interaction term to the Shan-Chen fluid/fluid interaction force. But wouldn’t doing so via a similar approach as in the ‘setFreeEnergyWallBoundary’ function necessitate some kind of “coupling” ? (just like the Shan-Chen forced postProcessor in the phase separation example in 2D)
– @both of you guies, don’t you think that, according to the way this additional force term is formulated in the cited Huang et al paper, adding this solid/fluid interaction force term directly into the bulk resolution? (by applying this force to every bulk node and just use a “solid neighbor indicator function” as a factor as suggested in the paper). In that case I would start from the ShanChenForcedSingleComponentPostProcessor (or any similar one depending on the issue to solve) and modify it directly. (this approach seems safer to me, a pure roockie of OpenLB, since I am not sure yet if one can modify the FORCE field in two speparate postProcessor (I am not sure I understand yet the difference between EXTERNAL_FORCE and FORCE fields))

Cheers,

Gilles

[mathias]

FORCE is basically for the surface tension and EXTERNAL_FORCE is for any other forceslike gravity, magnatic,..

[GASciences]

Thank you for your answer @mathias. I am tempted to ask another question then: Let’s take again the example of “ShanChenForcedSingleComponentPostProcessor” if we look at the end of this function it seems that the content of the FORCE field is purely overridden by the EXTERNAL_FORCE field + the additional Shan-Chen model’s term :

auto force = cell.template getField<descriptors::EXTERNAL_FORCE>();
cell.template setField<descriptors::FORCE>(force – g*rhoContribution/rho);

But what if somewhere else on the application the FORCE field is used? I assume that the application developper will have to re-write this portion of code then.

Also, if I look at the forcing.h file, all the forcing schemes seems to take the FORCE field into account but none takes the EXTERNAL_FORCE. Which leads me to think that if an application developer wants to use the EXTERNAL_FORCE field he/she either has to override the FORCE field the same way “ShanChenForcedSingleComponentPostProcessor” does or to re-write and adapt a forcing scheme, is that correct?

Sorry for this basic questions but the concept of EXTERNAL_FORCE Vs FORCE is puzzling. For example, the user manual only mentions the FORCE field and in quite general examples (no mention of the surface tension I believe).

Thanks again, I believe this discussion already provided me enough understanding to implement the Huang et al liquid/fluid interaction force.

Gilles

[Adrian]

The FORCE field describes the local force field to be applied per-cell and is read by the local (forced) dynamics. I.e. for simple cases one can set the FORCE field directly while more complex cases as Shan Chen require a separation between a user-defined EXTERNAL_FORCE and the computed FORCE.

As all fields, including (EXTERNAL_)FORCE are local to each cell a code modification should only required if multiple models apply to the same cell (in which case there likely are further model problems that have to be considered). => You are correct in your assumptions but feel free to ask further questions.

As a general note I would recommend to implement any new post processors and coupling operators in the new style introduced by OpenLB 1.5. This will require much less boilerplate than the legacy generator / processSubDomain approach while also transparently working on GPUs and improving performance also on CPU targets. The recent user guide for 1.6 contains FAQ entries explaining on how to write new operators in this style.

[Amirmansour Jafari]

Hi @savis. Thank you for your thorough explanation of high density and viscosity ratios. I am interested in checking that the density of two fluids is considered the same and the problem of simulating immiscible two-phase fluids is solved with the difference in high viscosity ratio. I would be grateful if you could introduce a reference in this field.

Best,
Amirmansour

[TimBingert]

Hi Amirmansour,

for a model with viscosity ratios of the two phases one has to connect the relaxation times responsible for macroscopic viscosities to the order parameter (in case of the OLB FreeEnergy, in case of the ShanChen to the density), as can be seen in https://doi.org/10.1103/PhysRevE.87.043301. Here, the authors suggest multiple schemes for viscosity interpolation between the phases through the order parameter. This macroscopic interpolation can then be mapped on the relaxation times, so they have to be updated each timestep.

Hope, this helps.
Kind regards,
Tim

[Author]

Posts