Skip to content

spurious diffusive effect correction#112

Draft
Ricardoleite wants to merge 1 commit intoOPM:devfrom
poro-labcc:wettability-improvement
Draft

spurious diffusive effect correction#112
Ricardoleite wants to merge 1 commit intoOPM:devfrom
poro-labcc:wettability-improvement

Conversation

@Ricardoleite
Copy link
Contributor

The present proposed correction to the lbpm_color_simulator routine aims to eliminate spurious diffusive effects on the phase field occurring across the fluid–solid interface. This diffusive effect is well know in the literature (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)). It arises when the computed color gradients are not properly incorporated into the recoloring step, causing the non-wetting fluid to segregate from the solid phase and form an artificial film along the interface.

Because the phase field term ($\phi$) is defined within the solid region, the segregation term incorrectly interprets the interface between the non-wetting fluid and the solid as a fluid–fluid interface that must be segregated. As a result, unphysical film formation occurs. This behavior is clearly illustrated in Figure bellow, 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 interface with non-wetting fluid is evident. Moreover, the intensity of this film increases as the non-wetting character of the fluid becomes stronger.

Screenshot from 2026-01-16 11-29-35

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 an implementation that preserves the standard LBPM computation of gradients over the entire domain for the momentum balance equation, since this approach is able to correctly capture the pressure response induced by wettability effects (e.g., via the Young equation). However, for the gradients applied in the mass balance distribution functions, we modify the interface normal vectors at locations where the fluid–fluid 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:

$$ n_{\pm} = \left( \cos(\pm \theta) - \frac{\sin(\pm \theta)\cos(\theta')}{\sin(\theta')} \right) n_{w} + \frac{\sin(\pm \theta)}{\sin(\theta')} n_{c} $$

where $n_{c}$ is the unit normal vector of the fluid–fluid interface and $n_{w}$ is the unit normal vector of the solid wall.

Implementation

Click to expand code

To implement this correction, it is necessary to count the number of solid neighbors around the fluid–solid interface. Therefore, we introduce an IDSolid array in ScaLBL_ColorModel::Create():

Allocate and Initialize IDSolid

// Allocate solid-fluid id
ScaLBL_AllocateDeviceMemory((void **)&IDSolid, sizeof(signed char) * Nx * Ny * Nz);

// Initialize solid-fluid id
signed char *TmpID;
TmpID = new signed char[Nx * Ny * Nz];
for (int k = 0; k < Nz; k++) {
    for (int j = 0; j < Ny; j++) {
        for (int i = 0; i < Nx; i++) {
            int idx = k * Nx * Ny + j * Nx + i;
            if (id[idx] <= 0){
                TmpID[idx] = 0;
            }
            else{
                TmpID[idx] = 1;
            }
        }
    }
}
ScaLBL_CopyToDevice(IDSolid, TmpID, sizeof(signed char) * Nx * Ny * Nz);
ScaLBL_Comm->Barrier();
delete[] TmpID;

Within the functions ScaLBL_D3Q19_AAeven_Color and ScaLBL_D3Q19_AAodd_Color, an initial data-reading step is introduced in which the IDSolid values are loaded together with the phase-field terms:

IDSolid data-reading

// Get the 1D index based on regular data layout
ijk = Map[n];
//					COMPUTE THE COLOR GRADIENT
//........................................................................
//.................Read Phase Indicator Values............................
//........................................................................
nn = ijk - 1; // neighbor index (get convention)
m1 = Phi[nn]; // get neighbor for phi - 1
id1 = IDSolid[nn];
//........................................................................
nn = ijk + 1; // neighbor index (get convention)
m2 = Phi[nn]; // get neighbor for phi - 2
id2 = IDSolid[nn];
//........................................................................
nn = ijk - strideY; // neighbor index (get convention)
m3 = Phi[nn];       // get neighbor for phi - 3
id3 = IDSolid[nn];
//........................................................................
nn = ijk + strideY; // neighbor index (get convention)
m4 = Phi[nn];       // get neighbor for phi - 4
id4 = IDSolid[nn];
//........................................................................
nn = ijk - strideZ; // neighbor index (get convention)
m5 = Phi[nn];       // get neighbor for phi - 5
id5 = IDSolid[nn];
//........................................................................
nn = ijk + strideZ; // neighbor index (get convention)
m6 = Phi[nn];       // get neighbor for phi - 6
id6 = IDSolid[nn];
//........................................................................
nn = ijk - strideY - 1; // neighbor index (get convention)
m7 = Phi[nn];           // get neighbor for phi - 7
id7 = IDSolid[nn];
//........................................................................
nn = ijk + strideY + 1; // neighbor index (get convention)
m8 = Phi[nn];           // get neighbor for phi - 8
id8 = IDSolid[nn];
//........................................................................
nn = ijk + strideY - 1; // neighbor index (get convention)
m9 = Phi[nn];           // get neighbor for phi - 9
id9 = IDSolid[nn];
//........................................................................
nn = ijk - strideY + 1; // neighbor index (get convention)
m10 = Phi[nn];          // get neighbor for phi - 10
id10 = IDSolid[nn];
//........................................................................
nn = ijk - strideZ - 1; // neighbor index (get convention)
m11 = Phi[nn];          // get neighbor for phi - 11
id11 = IDSolid[nn];
//........................................................................
nn = ijk + strideZ + 1; // neighbor index (get convention)
m12 = Phi[nn];          // get neighbor for phi - 12
id12 = IDSolid[nn];
//........................................................................
nn = ijk + strideZ - 1; // neighbor index (get convention)
m13 = Phi[nn];          // get neighbor for phi - 13
id13 = IDSolid[nn];
//........................................................................
nn = ijk - strideZ + 1; // neighbor index (get convention)
m14 = Phi[nn];          // get neighbor for phi - 14
id14 = IDSolid[nn];
//........................................................................
nn = ijk - strideZ - strideY; // neighbor index (get convention)
m15 = Phi[nn];                // get neighbor for phi - 15
id15 = IDSolid[nn];
//........................................................................
nn = ijk + strideZ + strideY; // neighbor index (get convention)
m16 = Phi[nn];                // get neighbor for phi - 16
id16 = IDSolid[nn];
//........................................................................
nn = ijk + strideZ - strideY; // neighbor index (get convention)
m17 = Phi[nn];                // get neighbor for phi - 17
id17 = IDSolid[nn];
//........................................................................
nn = ijk - strideZ + strideY; // neighbor index (get convention)
m18 = Phi[nn];                // get neighbor for phi - 18
id18 = IDSolid[nn];
}

The vector "nx,ny,nz" used in the mass-balance distribution function is replaced by a corrected vector "npx,npy,npz", while the original formulation is preserved for the momentum balance. The computation of "npx,npy,npz" is provided in the code below:

Vector - Angle Correction

//...........Correct wettability vector for Mass Balance.................................
if (id1 == 0 || id2 == 0 || id3 == 0 || id4 == 0 || id5 == 0 || id6 == 0 || id7 == 0 || id8 == 0 || id9 == 0 || id10 == 0 || id11 == 0 || id12 == 0 || id13 == 0 || id14 == 0 || id15 == 0 || id16 == 0 || id17 == 0 || id18 == 0){

    int int_nsx = (id1 - id2) * 2 + (id7 - id8 + id9 - id10 + id11 - id12 + id13 - id14);
    int int_nsy = (id3 - id4) * 2 + (id7 - id8 - id9 + id10 + id15 - id16 + id17 - id18);
    int int_nsz = (id5 - id6) * 2 + (id11 - id12 - id13 + id14 + id15 - id16 - id17 + id18);

    double Mag = sqrt(double(int_nsx * int_nsx + int_nsy * int_nsy + int_nsz * int_nsz));
    if (Mag == 0.0)
        Mag = 1.0f;
    nsx = -int_nsx / Mag;
    nsy = -int_nsy / Mag;
    nsz = -int_nsz / Mag;

    int countid =   id1 + id2 + id3 + id4 + id5 + id6 + id7 + id8 + id9 + id10 + id11 + id12 + id13 + id14 + id15 + id16 + id17 + id18;

    double aff = (m1*(id1-1) + m2*(id2-1) + m3*(id3-1) + m4*(id4-1) + m5*(id5-1) + m6*(id6-1) + m7*(id7-1) + m8*(id8-1) +
    m9*(id9-1) + m10*(id10-1) + m11*(id11-1) + m12*(id12-1) + m13*(id13-1) + m14*(id14-1) + m15*(id15-1) +
    m16*(id16-1) + m17*(id17-1) + m18*(id18-1)) / (18.0f - double(countid));

    Mag = 1.0f-(nx*nsx + ny*nsy + nz*nsz)*(nx*nsx + ny*nsy + nz*nsz);
    Mag = (Mag > 0.0f) ? sqrtf(Mag) : 1.0f;

    npx = (nx - nsx*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsx*aff;
    npy = (ny - nsy*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsy*aff;
    npz = (nz - nsz*(nx*nsx + ny*nsy + nz*nsz))*sqrt(1.0f-aff*aff)/Mag + nsz*aff;
}else{
    npx = nx;
    npy = ny;
    npz = nz;
}

@Ricardoleite
Copy link
Contributor Author

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.

Screenshot from 2026-02-23 16-06-55

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
import plotly.graph_objects as go
import plotly.io as pio
pio.renderers.default = 'notebook'  # or 'browser' if you want

# --- 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 = 1
   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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant