SubPhase.cpp - Flow Rate Calculation

[Ricardoleite]

Hello LBPM community,

I would like to bring to your attention some possible errors that I have identified in the flow rate calculation within the Core Flooding routine.

While reviewing the code in the file "SubPhase.cpp," I have noticed that at line 405, the "flow_magnitude" parameter is calculated using this equation:

405         double flow_magnitude = dir_x * dir_x + dir_y * dir_y + dir_z * dir_z;

To correct compute the "flow_magnitude" it is necessary to insert square root, i.e.,

405         double flow_magnitude = sqrt(dir_x * dir_x + dir_y * dir_y + dir_z * dir_z);

Consequently, it's necessary to correct the "flow_magnitude" for calculating the normalized vectors of flow direction in lines 412, 413, and 414:

412         dir_x /= flow_magnitude;
413         dir_y /= flow_magnitude;
414         dir_z /= flow_magnitude;

To illustrate the impact of this correction, 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 {
}        

This simulation setup should provide values of flow rate around $u_{z}=0.066666666$. The value can be obtained through capillary equation:

$$ Ca=\frac{\mu u_{z}}{\gamma} \quad \rightarrow \quad u_{z}=\frac{\gamma Ca}{\mu}=0.0666666 \quad \rightarrow \quad Q_{z}=\epsilon C_{xy}u_{z}=0.970588 \times 68 \times 5 \times u_{z}=21,99999 $$

where $\mu$ is the dynamic viscosity, $u_{z}$ is the fluid (usually water) velocity, $\gamma$ is the interfacial tension, $Q_{z}$ is the volumetric flux, $C_{xy}$ is the cross-sectional area, and $\epsilon$ is the porosity.

The figure below illustrate the results of $\phi$ field evolution.

The table below shows the values of 'water_flow_rate' (WFR) and 'not_water_flow_rate' (NWFR) for the case without any correction (Case 1) and with correction of line 405 (Case 2).

	Case 1	Case 1	Case 2	Case 2
	WFR	NWFR	WFR	NWFR
t=500	6.0874668e-06	1.2654743e-06	0.034233833	0.0071165953
t=1000	6.023764e-06	1.3291772e-06	0.057217132	0.01262528
t=1500	5.2948992e-06	2.058042e-06	0.06354778	0.024699999
t=2000	4.5222899e-06	2.8306513e-06	0.046822044	0.029307471
t=2500	3.3200053e-06	4.0329359e-06	0.028894388	0.035099105
t=3000	2.7212275e-06	4.6317137e-06	0.020831444	0.03545653
t=3500	2.1033173e-06	5.2496239e-06	0.016925438	0.042243833
t=4000	1.6095557e-06	5.7433854e-06	0.014022473	0.050036457