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

[Ricardoleite]

Irregular solid surfaces

Under irregular solid surfaces (i.e., staircase-like walls resulting from voxelization), diffusive effects become more pronounced. To investigate this behavior, we consider the benchmark case of a droplet resting on a solid sphere. Figure below (Source: Wanget al. (2025) (https://doi.org/10.1103/xz95-6bmt)) presents a three-dimensional slice illustrating the geometric configuration of the problem. The Python code used to generate the computational geometry is provided in the toogle bellow.

Click to expand code

import numpy as np

# Domain size
Nx, Ny, Nz = 100, 100, 100
x, y, z = np.indices((Nx, Ny, Nz))

# Initialize domain with label 1 (fluid/air)
npvol = np.ones((Nx, Ny, Nz), dtype='uint8')

# ---- Fluid Droplet (Label 2) aligned on top ----
R_fluid = 20
z1 = z0 + R_solid + R_fluid // 2  # Slightly above the solid
fluid_mask = ((x - x0)**2 + (y - y0)**2 + (z - z1)**2) < R_fluid**2
npvol[fluid_mask] = 2  # Droplet

# ---- Solid Sphere (Label 0) ----
R_solid = 30
x0, y0, z0 = Nx // 2, Ny // 2, Nz // 2 - R_solid // 2  # Center lower in z
solid_mask = ((x - x0)**2 + (y - y0)**2 + (z - z0)**2) < R_solid**2
npvol[solid_mask] = 0  # Solid

The diffusive effect observed in this problem is illustrated in the GIFs below for the case with component affinity $\phi_{s}=0.9$. In the simulation without the proposed correction, the droplet exhibits excessive spreading over the spherical surface, indicating the presence of spurious diffusion. The set up data in the .db is shared in the toogle bellow.

Click to expand code

Domain {
   Filename = "sphere-sol-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 = 0.7
   tauB = 0.7
   F = 0, 0, 0
   WettingConvention = "SCAL"
   ComponentLabels = 0
   ComponentAffinity = 0.9
   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
}

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

It is important to note that the segmented .raw images used to generate the GIFs visually attenuate this effect, making it appear less pronounced than it actually is. A more accurate representation of the physical behavior can be obtained by examining a $x-y$ slice at $z=50$ of the phase-field variable. This slice clearly reveals the artificial spreading mechanism, as illustrated in the Figure below.

Comparing the measured contact angle with the theoretical (expected) contact angle, Figure below illustrates the results obtained for a grid resolution of $100^{3}$. The results show good agreement with the expected arccosine relationship, with small deviations over the range $-0.9 \leq \phi_{s} \leq 0.9$ and and slightly higher errors compared to the droplet-on-a-plane problem.

However, the employed correction does not eliminate the spurious currents characteristic of the implemented modeling approach. These currents arise because the pressure near the solid surface is modified proportionally to the gradient between the local phase field and the local component affinity. Figure below illustrate the velocity magnitude on an $x-y$ slice at $z=50$ for $\phi_{s}=0.1$. To mitigate the spurious currents originating from the solid surface, further improvements in the wettability modeling are required.

[Ricardoleite]

3D Digitla Rock

Comparing the results in a 3D digital rock scenario, we can qualitatively observe in the videos below, in a core-flooding drainage process, the removal of spurious diffusive effects that artificially create droplets throughout the domain. The digital rock image corresponds to Bentheimer sandstone reported by Neumann et al. (2020) (https://www.doi.org/10.17612/f4h1-w124), obtained from the Digital Porous Media Portal. The .db file used to configure the simulation is provided in toogle bellow.

Click to expand code

Color {
   protocol          = "core flooding"
   capillary_number  = 1e-3            // capillary number for the displacement
   timestepMax       = 150000          // maximum timtestep
   alpha             = 0.01            // controls interfacial tension
   beta              = 0.98
   rhoA              = 1.0             // controls the density of fluid A
   rhoB              = 1.0             // controls the density of fluid B
   tauA              = 0.55            // controls the viscosity of fluid A
   tauB              = 0.55            // controls the viscosity of fluid B
   F                 = 0, 0, 0         // body force
   WettingConvention = "SCAL"          // convention for sign of wetting affinity
   ComponentLabels   = 0        // image labels for solid voxels
   ComponentAffinity = 0.5  // controls the wetting affinity for each label
   Restart           = false
}
Domain {
   Filename          =  "domain.raw"
   ReadType          =  "8bit"           // data type
   N                 =  300,200,400      // size of original image
   nproc             =  1, 1, 1          // process grid
   n                 =  300,200,400      // sub-domain size
   offset            =  0, 0, 0          // offset to read sub-domain
   voxel_length      =  16            // voxel length (in microns)
   ReadValues        =  0, 1          // labels within the original image
   WriteValues       =  0, 2          // associated labels to be used by LBPM
   BC                =  4             // boundary condition type (0 for periodic)
}
Analysis {
   analysis_interval          = 200     // logging interval for timelog.csv
   subphase_analysis_interval = 200     // loggging interval for subphase.csv
   visualization_interval     = 200     // interval to write visualization
   N_threads                  = 0       // number of analysis threads
   restart_interval           = 20000000    // interval to write restart file
   restart_file               = "Restart"   // base name of restart file
}
Visualization {
   format                     = "vtk"
   write_silo                 = true     // write SILO databases with assigned variables
   save_8bit_raw              = true      // write labeled 8-bit binary files with phase assignments
   save_phase_field           = true     // save phase field within SILO database
   save_pressure              = false     // save pressure field within SILO database
   save_velocity              = false     // save velocity field within SILO database
}

After-Correction	Before-Correction
drainage-correc.mp4	drainage-wrong.mp4