[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

[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

[savis]

To my understanding, changing the virtual density would affect the contact angle because the increased / decreased density gradient with the gas increases / decreases the solid-gas surface energy. So the solid-liquid area changes to reduce the total surface energy. This therefore affects the contact angle in a similar way to changing the fluid-solid force, which may be sufficient depending upon what you are doing.

However, I don’t believe it truly models the interaction with the solid because this would require a Neumann boundary condition in the density rather than a Dirichlet condition. There is more detail about this in chapter 9 of ‘The Lattice Boltzmann Method: Principles and Practice’ if you are interested.

Apologies for the delayed response.
Sam

[savis]

Unfortunately it’s not possible in the free energy model that’s currently implemented, and it would be quite a large undertaking to do so. If the effect of different densities is required then it would probably be easiest to implement the fluid-solid interaction in the SC model. I have discussed this before here.

The issue back then was that an indicator for the solid is required but not accessible within the lattice coupling. So a work-around will be required if you want to take this route.

Best,
Sam

[savis]

Hi Bhuttu,

I don’t have a great understanding of the topic myself so I could be mistaken. However, I believe that the contact line singularity is solved because it uses a diffuse interface model that allows matter to diffuse across the boundary. This is detailed in this paper, and it also contains references to other relevant papers.

Best,
Sam

[savis]

I have found the issue, it is a bug in this function. All of the “uPerp = +/- u[0/1/2];” lines have the wrong sign, which causes problems whenever an interface tries to pass the outlet.

Thanks for pointing out this issue. I will fix it for the next release. In the mean time you can also fix it in your own installation.

Best,
Sam

[savis]

Hi Mehrdad,

Would it be possible for you to share your code so that I can take a look to see what may be the cause? There are situations where the outlet boundary condition can become unstable, so that may be the case. However, if it is working in 2D then there may just be an error somewhere. Also, the outlet is not a Dirichlet boundary for phi because it must allow different fluids to pass, so some change is to be expected.

Regarding the uninitialized iC warning, I don’t believe this will be the cause of the issue. The warning seems to be because LatticeCouplingGenerator3D has an iC attribute that is uninitialized in the constructor, whereas the 2D version does not. However, I don’t fully understand the lattice coupling interface, so I simply added a shift function to the example. It would be great if someone else with a better understanding could take a look at this.

Best,
Sam

[savis]

Hi Bhuttu,

The issue is that the densities of the two SuperLattices are not the densities of the fluid components. The density of sLattice1 is the total fluid density (= rho_1 + rho_2), while the “density” of sLattice2 is an order parameter to define whether it is fluid 1 or fluid 2 (= rho_1 – rho_2).

Therefore to plot the total density you only need to plot “density1”. It should be slightly higher inside the droplet due to the Laplace pressure. The reason you currently see a smaller value in the droplet is because “density1+density2” is two times the density of fluid 2, which is zero inside the droplet.

I hope this helps.
Sam

[savis]

Unfortunately external forces can’t currently be used with the multicomponent model. It is something that I have been meaning to implement for a while. The reason that it doesn’t work is because the force is overwritten by the free energy model on each time step. (Line 361 of this method)

If you don’t mind editing the OpenLB source code directly then it would be quite simple to include a gravity force by editing the method linked above.
For example, to include gravity in the negative z-direction with a force density of 0.01 (in lattice units), the following line could be added before line 361:
`forceZ -= 0.01;
Note that this applies the same force to all of the different components. If you need different forces for the different components then it must depend upon the values of phi and psi.

Let me know if anything is unclear.

Sam

[savis]

Sorry, the surface tension (gamma) would be 0.67 so the value of kappa would be 2.

[savis]

Dear nv4dll,

I probably shouldn’t have used gamma for the surface tension, the value of gama in the code is actually unrelated. Also, the values that I used for the contactAngle2d example were not based upon a physical system so I would not read too much into them.

You are right about how to convert the surface tension into lattice units. So for example for the surface tension of water (0.072 N/m) would give kappa=0.67 for those values.

However you need to be careful, because if kappa is larger than about 0.025 then the simulation will be unstable. So if that physical surface tension value is needed then you will need to change the values in the unit converter initialisation or the resolution.

Sam

[savis]

Also, since you only use 2 lattices, you should ensure that you do not include the arguments kappa3 and h3 in addFreeEnergyWallBoundary. This could potentially also cause problems.

[savis]

Dear nv4dll,

Are you by any chance using different values for kappa1 and kappa2? If so then this might be the cause. The reason is that if the value h1/kappa1 + h2/kappa2 != 0 then the density on the solid boundary will be different. This will then interact with the outlet causing it to become unstable.

To find the correct kappa values you can use the formula gamma = alpha/6 * (kappa1 + kappa2), where gamma is the surface tension converted into lattice units. Since you are only using two components you can therefore set kappa1 = kappa2 = 3*gamma/alpha.

Unfortunately by changing the kappa values you will also change the contact angles, so you may need to recompute the values h1 and h2.

If this isn’t the issue then please do let me know.

Best wishes,
Sam

[savis]

Dear Simon,

To my knowledge, there is no three component implementation of the Shan-Chen model in OpenLB. So this approach would likely involve much more work than simply adding a solid interaction to the existing method.

Best wishes,
Sam

[savis]

Dear Simon,

I believe that this is because of boundary nodes that aren’t beside any inner fluid nodes. You should hopefully be able to resolve this simply using “superGeometry.clean();” after the material numbers are assigned.

Best wishes,
Sam

[Author]

Posts