Correction of velocity calculation in the lbpm_permeability_simulator

[Ricardoleite]

There are two issues in the fucntion "ScaLBL_D3Q19_Momentum(double *dist, double *vel, int Np) ", which is dedicated to calculating velocity array used to estimate the absolute permeability. I list and explain them below:

• The velocity is currently computed as $V_{\alpha}=\sum_{i}f_{i}c_{\alpha,i}$ (not correct) instead of $V_{\alpha}=\sum_{i}f_{i}c_{\alpha,i}/\rho$ (Correct), of coarse, if we consider a first-order discretization of space-time in the presence of an external force that drving the flow. However, it is well known in the LBM literature that problems with an external force term require a second-order discretization of space–time. In this case, the correct formulation becomes: $V_{\alpha}=\sum_{i}f_{i}c_{\alpha,i}/\rho + \delta_{t}F_{\alpha}/2$.

• As the code is currently structured, the analysis routine is performed after the AA-even time step, when the distribution functions are in a post-collisional state and not yet underwent streaming . To clarify, consider the D2Q9 lattice to illustrate odd and even steps (I am following the odd and even convention presented by LBPM). In the Lattice Boltzmann Method, the macroscopic moments should be computed is after the streaming step. Thus, in LBPM velocity is calculated after the even time step, when the distributions are still post-collisional. This is problematic because, in the presence of an external force, the first-order moment is not a conserved quantity. In other words, pre-collisional and post-collisional moments differ. With the current code structure these can be fixed by using the post-collisional distribution functions to compute the pre-collisional velocity after the even time step.

$$V_{\alpha} = \sum_{i} f_{i}^{*} c_{\alpha,i} / \rho - \delta_{t}F_{\alpha}/2 $$

where $f_{i}^{*}$ is post-collision distribution function.

BCC (Benchmark)

The goal of this correction is to improve accuracy and reduce the artificial dependence of permeability on viscosity, as discussed by Pan et al. (2006) (https://doi.org/10.1016/j.compfluid.2005.03.008).

The figure below compares results before and after applying the correction for a body-centered cubic (BCC) array. The BCC geometry follows the parameters in Pan et al. (2006). The analytical solution used for comparison is given below.

BCC analytical solution

Click to expand code

import warnings
warnings.filterwarnings("ignore", category=UserWarning, module='matplotlib')
import numpy as np
import matplotlib.pyplot as plt

#Coefficients up to the order 25
alpha_s = np.array([1.0,1.575834,2.483254,3.233022,4.022864,4.650320,5.281412,
                    5.826374,6.258376,6.544504,6.878396,7.190839,7.268068,
                    7.304025,7.301217,7.2364410,7.298014,7.369847,7.109497,
                    6.228418,5.235796,4.476874,3.541982,2.939353,3.935484])
a=11.77/32.0 #Dimensionless Ratio
a0=11.77 #Dimensional Ratio
chi=( (64.0*a**(3.))/(np.sqrt(3.0)*3.0*(1.0)**(3.0)) )**(1./3.) #Ratio between solid volume and max volume of solid
#-------------Dimensionless Drag Force Calculation-----------------------------
d=0.0
for i in range (0,len(alpha_s)):
    d=d+alpha_s[i]*chi**(i)
#-------------------------------------------------------------------------------
ks=1.0/(d*6.0*np.pi*a) #Dimensionless Permeability
print("--------------------------------------------------------------------------------")
print("d=",d, "        Dimensionless Drag Force")
print("k=",ks, "     Dimensionless Permeability")
print("chi=",chi, "     Ratio between solid volume and max volume of solid")
print("--------------------------------------------------------------------------------")

Code to generate the BCC Geometry

Click to expand code

import numpy as np
import matplotlib.pylab as plt
cx = 32
cy = 32
cz = 32
Rr=11.77
D=np.ones((cx,cy,cz),dtype="uint8")*1
for i in range(0,cx):
    for j in range (0, cy):
        for k in range (0,cz):
            dist = np.sqrt((i+0.5-cx/2)*(i+0.5-cx/2) + (j+0.5-cy/2)*(j+0.5-cy/2) + (k+0.5-cz/2)*(k+0.5-cz/2))
            if (dist < Rr) :
                D[i,j,k] = 0
            dist = np.sqrt((i+0.5-cx)*(i+0.5-cx) + (j+0.5)*(j+0.5) + (k+0.5)*(k+0.5))
            if (dist < Rr) :
                D[i,j,k] = 0
            dist = np.sqrt((i+0.5)*(i+0.5) + (j+0.5-cy)*(j+0.5-cy) + (k+0.5)*(k+0.5))
            if (dist < Rr) :
                D[i,j,k] = 0
            dist = np.sqrt((i+0.5)*(i+0.5) + (j+0.5)*(j+0.5) + (k+0.5-cz)*(k+0.5-cz))
            if (dist < Rr) :
                D[i,j,k] = 0
            dist = np.sqrt((i+0.5)*(i+0.5) + (j+0.5)*(j+0.5) + (k+0.5)*(k+0.5))
            if (dist < Rr) :
                D[i,j,k] = 0
            dist = np.sqrt((i+0.5)*(i+0.5) + (j+0.5-cy)*(j+0.5-cy) + (k+0.5-cz)*(k+0.5-cz))
            if (dist < Rr) :
                D[i,j,k] = 0
            dist = np.sqrt((i+0.5-cx)*(i+0.5-cx) + (j+0.5)*(j+0.5) + (k+0.5-cz)*(k+0.5-cz))
            if (dist < Rr) :
                D[i,j,k] = 0
            dist = np.sqrt((i+0.5-cx)*(i+0.5-cx) + (j+0.5-cy)*(j+0.5-cy) + (k+0.5)*(k+0.5))
            if (dist < Rr) :
                D[i,j,k] = 0
            dist = np.sqrt((i+0.5-cx)*(i+0.5-cx) + (j+0.5-cy)*(j+0.5-cy) + (k+0.5-cz)*(k+0.5-cz))
            if (dist < Rr) :
                D[i,j,k] = 0

D.tofile("bcc-%d.raw"% (cx)) # Array save in a file .raw 

BCC .db file

Click to expand code

MRT {
   tau = 0.54
   F = 0.0, 0.0, 0.000001
   timestepMax = 6000000
   tolerance = 0.000001
}
Domain {
   Filename = "../../bcc-32.raw"
   ReadType = "8bit"
   nproc = 1, 1, 1
   N = 32, 32, 32
   n = 32, 32, 32
   voxel_length = 1.0
   ReadValues = 0, 1, 2
   WriteValues = 0, 1, 2
   BC = 0
}
Visualization {
   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 = true    // save pressure field within SILO database
   save_velocity = true    // save velocity field within SILO database
}

LV60A (3D Digital Rock)

Applying the same correction to the LV60A sandpack (https://figshare.com/articles/dataset/LV60A_sandpack/1153795/2
) yields the same outcome: a much weaker (closer to negligible) dependence of permeability on viscosity.

LV60A .db file

Click to expand code

MRT {
   tau = 0.56
   F = 0.0, 0.0, 0.000001
   timestepMax = 6000000
   tolerance = 0.000001
}
Domain {
   Filename = "../../lv60a.raw"
   ReadType = "8bit"
   nproc = 1, 1, 1
   N = 450, 450, 450
   n = 450, 450, 450
   voxel_length = 10.002
   ReadValues = 0, 1, 2
   WriteValues = 0, 1, 2
   BC = 0
}
Visualization {
   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 = true    // save pressure field within SILO database
   save_velocity = true    // save velocity field within SILO database
}
}

[JamesEMcClure]

By construction the density is equal to 1. To reduce the effects of compressibility, the first moment is proportional to the pressure and the unit density is implied in the calculations. I believe the original reference is https://doi.org/10.1023/B:JOSS.0000015179.12689.e4 , but the MRT extensions were done elsewhere (based on the same idea)

I agree that the syntax for the function is confusing, but this is the reason for the construction of ScaLBL_D3Q19_Momentum

This is more important in the color model, because the momentum equation is only tracking the pressure and the momentum; the density must be obtained from the mass transport LBEs. You should never use the first moment from the D3Q19 distributions in LBPM as density. My suggestion is to replace vel in the argument list with momentum so as to prevent this confusion.

For the permeability issue, note that the implied wall location is viscosity dependent. When you work with a "stair-step" type of boundary wall location, then you end up with this undesirable behavior as shown above. You can do what Pan et al. suggest and interpolate to get a better bounce-back rule, and this will be more accurate. However, there is a tradeoff due to the CFL condition (i.e. the \delta x that goes into the CFL condition can be arbitrarily close to the fluid lattice size, which then imposes a bad constraint on the timestep and you can run into stability problems that impact robustness). For permeability this doesn't matter so much, but in the color model it will. Our thinking was that the single phase permeability should rely on a similar bounce-back rule to the multiphase model so that the two approaches are consistent.

Note that the effects above will be considerably smaller if you increase the resolution. LBPM will converge to the right permeability if you have a well-resolved system.

[JamesEMcClure]

With respect to the debate on forcing and velocity, please refer to this paper
https://doi.org/10.1016/j.ijheatmasstransfer.2009.11.014
LBPM uses the Luo convention defined from Eq. (1) in this reference. In this case the velocity should not be shifted (see comment after Eq. 14)