spurious diffusive effect correction

[Ricardoleite]

The proposed correction to the lbpm_color_simulator routine aims to eliminate spurious diffusive effects in the phase field that occur across the fluid–solid interface. This diffusive effect has been reported in the literature by several authors (e.g., Leclaire et al., 2017 (https://doi.org/10.1002/fld.4226); Yu et al., 2017 (https://doi.org/10.1177/0954406217749616); Akai et al., 2018 (https://doi.org/10.1016/j.advwatres.2018.03.014)).

This spurious effect arises when the computed color gradients are not properly incorporated into the recoloring step. Because the phase-field term ($\phi$) is defined within the solid region in the current approach, the segregation term incorrectly interprets the interface between the non-wetting fluid and the solid as a fluid–fluid interface that must be segregated, causing the non-wetting fluid to separate from the solid phase and leading to the formation of an artificial film along the interface.

This behavior is clearly illustrated in the figure below, adapted from McClure et al. (2014) (https://doi.org/10.1016/j.cpc.2014.03.012), where the formation of a fluid film along the solid surface in contact with the non-wetting fluid is evident. Moreover, the intensity of this film increases as the non-wetting character of the fluid becomes stronger.

However, the corrections proposed in the literature to address this problem are generally based on the standard color-gradient formulation, in which two distribution functions are used to discretize the coupled mass and momentum balances of each fluid. Consequently, these corrections cannot be directly applied to the LBPM framework, which discretizes mass and momentum balances separately: a single distribution function is used for the momentum balance of the mixture, while two additional distribution functions are employed for the mass balance of each fluid.

To correct this effect within the LBPM lattice Boltzmann formulation, we propose to retain the current formulation for computing the gradient over the entire domain. However, for the gradients applied in the mass-balance distribution functions, we modify the fluid–fluid interface normal vectors at locations where this interface is in contact with the solid region.

The correction of the unit normal vector follows the same approach proposed by Yu et al. (2017) and Akai et al. (2018). Denoting $\theta$ as the contact angle, $\hat{n}^{\ast}$ as the unit vector pointing in the direction of the color gradient, and $\hat{n}_s$ as the unit normal vector to the solid wall (pointing inwards), these authors proposed enforcing a gradient direction $\hat{n}$ that forms an angle $\theta$ with $\hat{n} _s$ while remaining in the plane defined by $\hat{n}^{\ast}$ and $\hat{n} _s$. This formulation leads to two possible solutions, $\hat{n} _{+}$ and $\hat{n} _{-}$, for the corrected gradient direction:

$$ \hat{n}_\pm = \cos \theta~ \hat{n}_s \pm \frac{\sin \theta}{\sin \theta^\prime} \left[ \hat{n}^{\ast} - \cos \theta^\prime \hat{n}_s \right] $$

where $\cos \theta^\prime = \hat{n}_s \cdot \hat{n}^{\ast}$.

It is worth noting that the correction is only applied if the imposed contact angle is not satisfied by the color gradient, since when $\theta = \theta^\prime$

$$ \hat{\boldsymbol{n}} = \hat{\boldsymbol{n}}^{*}. $$

Implementation

To implement this correction, information about solid neighbors is required for lattice sites located at the fluid–solid interface. For this purpose, an IDSolid array is introduced in ScaLBL_ColorModel::Create(), which is subsequently copied to the device.

https://github.com/poro-labcc/LBPM_fork/blob/474d3c211cea855983709ff0256adacb77f3112b/models/ColorModel.cpp#L509-L527

Within the functions ScaLBL_D3Q19_AAeven_Color and ScaLBL_D3Q19_AAodd_Color, the IDSolid values are loaded together with the phase-field terms:

https://github.com/poro-labcc/LBPM_fork/blob/474d3c211cea855983709ff0256adacb77f3112b/cpu/Color.cpp#L1488-L1558

Using the solid-neighbor information, the unit normal vector to the solid wall is computed and stored in nsx, nsy, and nsz. The fluid–fluid interface normal vector (nx, ny, nz) is then corrected according to the formulation described above, and the corrected vector is stored in npx, npy, and npz. The original formulation is preserved for the momentum-balance calculation.

https://github.com/poro-labcc/LBPM_fork/blob/474d3c211cea855983709ff0256adacb77f3112b/cpu/Color.cpp#L1574-L1604

[Ricardoleite]

Results and Improvements

Droplet-on-a-plane

After applying the proposed correction to the droplet-on-a-plane benchmark, noticeable changes are observed in the description of the contact angle as a function of the component affinity parameter ($\phi_{s}$),which is responsible for setting the phase-field value inside the solid region. In contrast to the results reported in Erfani et al. 2023 (https://doi.org/10.1007/s11242-023-02044-x), where the droplet contact angle ($\theta$) exhibits an approximately linear dependence on $\phi_{s}$, rather than the expected arccosine relationship ($\theta=180-arccos(\phi_{s})$). This behavior is illustrated in Figure A, where the corrected results closely follow the theoretical arccosine curve. In Figure B, we present a comparison between the measured contact angles and the expected arccosine relation for different grid resolutions. The results show a reduction in the average error as the grid spacing is refined.

The input data used to reproduce this analysis are provided below:

Python code used to generate the geometry

Click to expand code

import numpy as np

# --- Domain ---
Nx, Ny, Nz = 100, 100, 100
Lx, Ly, Lz = Nx - 1, Ny - 1, Nz - 1

# --- Flat plane at z = 0 ---
# We'll use value 0 for solid plane; 1 for background; 2 for droplet
npvol = np.ones((Nx, Ny, Nz), dtype='uint8')

# --- Spherical droplet that just touches the plane ---
R = 20
x0, y0 = Nx // 2, Ny // 2
z0 = R-int(R/4)  # center placed so the sphere touches z=0

# Build droplet
x, y, z = np.indices((Nx, Ny, Nz))
mask = (x - x0)**2 + (y - y0)**2 + (z - z0)**2 < R**2
npvol[mask] = 2

npvol[:, :, 0:2] = 0  # single-voxel solid floor at z=0

npvol.tofile(f'droplet3D-{Nx}-{Ny}-{Nz}.raw')

.db file employed to run the simulation

Click to expand code

Domain {
   Filename = "droplet3D-100-100-100.raw"
   ReadType = "8bit"
   N = 100, 100, 100
   nproc = 1, 1, 1
   n = 100, 100, 100
   voxel_length = 1.0
   ReadValues = 0, 1, 2
   WriteValues = 0, 1, 2
   BC = 0
}

Color {
   protocol = "fractional flow"
   capillary_number = 1e-4
   timestepMax = 60000
   alpha = 0.005
   rhoA = 1.0
   rhoB = 1.0
   tauA = 1.7
   tauB = 1.7
   F = 0, 0, 0
   WettingConvention = "SCAL"
   ComponentLabels = 0
   ComponentAffinity = -0.3
   Restart = false
}

Analysis {
   analysis_interval = 200
   subphase_analysis_interval = 100000
   N_threads = 0
   visualization_interval = 200
   restart_interval = 10000000
   restart_file = "Restart"
}

Visualization {
   write_silo = true
   save_8bit_raw = true
   save_phase_field = true
   save_pressure = true
   save_velocity = true
}

FlowAdaptor {
   min_steady_timesteps = 10000000
   max_steady_timesteps = 11000000
   fractional_flow_increment = 0.0
}

Another aspect observed in the droplet-on-a-plane benchmark under high interfacial tension is that, in the current version of LBPM, the droplet tends to spread completely over the surface when high contact angles are prescribed. To illustrate this behavior, we compare in the GIFs below the droplet simulations performed with alpha=0.01 before and after the proposed correction.

In the before-correction case, the droplet in the GIF gradually appears to disappear as spreading progresses. This effect occurs because the phase-field $\phi$ decreases during spreading, and the subsequent segmentation of the .raw data causes interfacial points to vanish from the visualization.

Left-Side: After-Correction and Right-Side: Before-Correction