SubPhase.cpp : "pressure_drop" - "force_mag"

[Ricardoleite]

Hello LBPM community,

I would like to bring to your attention some possible errors that I have identified in the "pressure_drop" and "force_mag" calculation within the Core Flooding routine.

While reviewing the code in the file "SubPhase.cpp," I have noticed that at line 419, the "pressure_drop" parameter is calculated by equation:

419         double pressure_drop = (Pressure(Nx * Ny + Nx + 1) - 1.0) / 3.0;

In the equation above "/3.0" is outside of the parenthesis dividing the "Pressure(Nx * Ny + Nx + 1)". The correct form should be represented by:

419         double pressure_drop = (Pressure(Nx * Ny + Nx + 1) - 1.0/ 3.0);

where "1.0/ 3.0" represents the pressure of the outlet constant pressure boundary condition in axis-Z.

Additionally, the inlet pressure in the "pressure_drop" calculation is represented by a pressure in a fixed position of "Nx * Ny + Nx + 1" ("Nx * Ny + Nx + 1" correspond to a 3d cartesian coordenate of "(x,y,z)=(1,1,1)"). The issue arise when this point corresponde to a solid position, resulting in "Pressure(Nx * Ny + Nx + 1)=0". Consequently, the "force_mag" given by:

421         force_mag -= pressure_drop / length;

used to calculate the relative permeability will be calculated wrong.

To illustrate the case of the inlet pressure to be in a solid point, I conducted a simulation of a viscous fingering problem within a 3D channel, having walls only on the upper and side boundaries (this problem can also be represented for a 2D fluid displacement in a channel). The simulation setup is given by:

Color {
   protocol = "core flooding"
   capillary_number = 0.238            // capillary number for the displacement
   timestepMax = 4000              // maximum timtestep
   alpha = 0.015561779                     // controls interfacial tension
   rhoA = 1.0                         // controls the density of fluid A
   rhoB = 1.0                         // controls the density of fluid B
   tauA = 1.5                         // controls the viscosity of fluid A
   tauB = 1.5                         // controls the viscosity of fluid B
   F = 0, 0, 0                        // body force
   WettingConvention = "SCAL"         // convention for sign of wetting affinity
   ComponentLabels = 0, -1, -2        // image labels for solid voxels
   ComponentAffinity = 0.0, 1.0, 0.6  // controls the wetting affinity for each label
   Restart = false
}
Domain {
   Filename = "/home/ricardo/LBPM/Results/ChannelDisplace/Channel-Displace-400-68-5.raw"
   ReadType = "8bit"              // data type
   N = 5, 68, 400             // size of original image
   nproc = 1, 1, 1                // process grid
   n = 5, 68, 400              // sub-domain size
   offset = 0, 0, 0         // offset to read sub-domain
   voxel_length = 1.0            // voxel length (in microns)
   ReadValues = -2, -1, 0, 1, 2   // labels within the original image
   WriteValues = -2, -1, 0, 2, 1  // associated labels to be used by LBPM
   BC = 4                         // boundary condition type (0 for periodic)
}
Analysis {
   analysis_interval = 500           // logging interval for timelog.csv
   subphase_analysis_interval = 500  // loggging interval for subphase.csv
   visualization_interval = 500    // interval to write visualization files
   N_threads = 1                      // number of analysis threads (GPU version only)
   restart_interval = 1000000         // interval to write restart file
   restart_file = "Restart"           // base name of restart file
}
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
}
FlowAdaptor {
}        

Due to the extra communication layers of the MPI interface, the point "(x,y,z)=(1,1,1)" in the simulation corresponds to the point "(x,y,z)=(0,0,0)" in the image. Notably, this point represents a solid point in the 3D channel. Consequently, in the output file "timelog.csv", the "force_mag" value is reported as "force" by the "analysis_interval", being equal to 0.00083333333, i.e.,

$$ force_mag=-\frac{\underbrace{P_{in}}_{=0}-1/3}{N_{z}}=-\frac{-1/3}{400}=0.00083333333 $$

There are numerous methods to address the issue of locating a fluid node at the inlet interface. I propose running a "for" around the inlet interface section of "rank=0" (the rank responsible for counting), to check for pressure values higher than 0. Here is an example of the mentioned code:

n=Nx * Ny + Nx + 1;
double pressure_drop = (Pressure(n) - 1.0/ 3.0);
for (i = 0; i < (Nx - 2) * (Ny - 2); i++) {
    if(Pressure(n + i)>0.0){
        pressure_drop = (Pressure( n + i) - 1.0/ 3.0);
        break;
    }
}

[JamesEMcClure]

I believe that the LBPM Domain object will remove layers at the inlet if an external boundary condition is applied. See line 1407 from common/Domain.cpp

https://github.com/OPM/LBPM/blob/master/common/Domain.cpp

This should also be visible in the output data. The reason this is done is that the presence of solid can influence the value of the pressure and/or flux that would need to be set (since the interfacial forces would also need to be accounted for).

It is possible that a pressure / flux BC will "work" without removing these layers, but I expect that it would reduce the range of stability for the core flooding protocol.

There are internal data structures that store the pore sites along a face, since these are used in the communication structures. These correspond to the values that need to be sent / recieved by the MPI communications, and could be used within this context if necessary.

[Ricardoleite]

Thank you for your attention, Professor James McClure,

Regarding our discussion about the inlet position of "Pressure (Nx * Ny + Nx + 1)", I would like to emphasize the importance of correcting line 419. It is crucial to accurately account for the relative permeability in this correction.

419         double pressure_drop = (Pressure(Nx * Ny + Nx + 1) - 1.0) / 3.0;

In the equation above "/3.0" is outside of the parenthesis dividing the "Pressure(Nx * Ny + Nx + 1)". The correct form should be represented by:

419         double pressure_drop = (Pressure(Nx * Ny + Nx + 1) - 1.0/ 3.0);

Also, if possible, please take a look at the other issue I submitted (#85). You might have missed it