OPM
SubPhase.cpp : "pressure_drop" - "force_mag"
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;
}
}
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.
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 (https://github.com/OPM/LBPM/issues/85). You might have missed it