diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 8c77af3a..840fada0 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -758,7 +758,7 @@ double ScaLBL_ColorModel::Run(int returntime) { double volA = Averages->gnb.V; volA /= Dm->Volume; volB /= Dm->Volume; - ; + //initial_volume = volA*Dm->Volume; double vA_x = Averages->gnb.Px / Averages->gnb.M; double vA_y = Averages->gnb.Py / Averages->gnb.M; @@ -925,6 +925,10 @@ double ScaLBL_ColorModel::Run(int returntime) { double krnf = kAeff - h * h * muA * not_water_film_flow_rate / force_mag; double krwf = kBeff - h * h * muB * water_film_flow_rate / force_mag; + /* not counting films */ + double krnf_low = (1.0 - current_saturation) * krnf; + double krwf_low = current_saturation * krwf; + // Saturation normalized effective permeability to account for decoupled phases and // effective porosity. double kAeff_low = (1.0 - current_saturation) * h * h * @@ -951,6 +955,8 @@ double ScaLBL_ColorModel::Run(int returntime) { "eff.perm.water.upper.bound "); fprintf(kr_log_file, "eff.perm.oil.film " "eff.perm.water.film "); + fprintf(kr_log_file, "eff.perm.oil.film.lower.bound " + "eff.perm.water.film.lower.bound "); fprintf(kr_log_file, "eff.perm.oil.lower.bound " "eff.perm.water.lower.bound "); fprintf(kr_log_file, @@ -969,6 +975,7 @@ double ScaLBL_ColorModel::Run(int returntime) { current_saturation); fprintf(kr_log_file, "%.5g %.5g ", kAeff, kBeff); fprintf(kr_log_file, "%.5g %.5g ", krnf, krwf); + fprintf(kr_log_file, "%.5g %.5g ", krnf_low, krwf_low); fprintf(kr_log_file, "%.5g %.5g ", kAeff_low, kBeff_low); fprintf(kr_log_file, "%.5g %.5g ", kAeff_connected, kBeff_connected); @@ -990,27 +997,39 @@ double ScaLBL_ColorModel::Run(int returntime) { WriteHeader = true; scal_log_file = fopen("SCAL.csv", "a"); if (WriteHeader) { - fprintf(scal_log_file, "timesteps sat.water "); + fprintf(scal_log_file, "timesteps " + "sat.water "); fprintf(scal_log_file, "eff.perm.oil.upper.bound " - "eff.perm.water.upper.bound "); + "eff.perm.water.upper.bound "); fprintf(scal_log_file, - "eff.perm.oil.lower.bound " - "eff.perm.water.lower.bound "); + "eff.perm.oil.lower.bound " + "eff.perm.water.lower.bound "); fprintf(scal_log_file, "eff.perm.oil.disconnected " - "eff.perm.water.disconnected "); + "eff.perm.water.disconnected "); + fprintf(scal_log_file, + "eff.perm.oil.film.lower.bound " + "eff.perm.water.film.lower.bound "); fprintf(scal_log_file, - "cap.pressure cap.pressure.connected " - "Ca eff.pressure\n"); + "cap.pressure " + "cap.pressure.connected " + "Ca " + "eff.pressure\n"); } + + // Adjust the effperms with porosity term to make the nominal + // values closer to measured data fprintf(scal_log_file, "%i %.5g ", CURRENT_TIMESTEP, current_saturation); - fprintf(scal_log_file, "%.5g %.5g ", kAeff_low, kBeff_low); - fprintf(scal_log_file, "%.5g %.5g ", kAeff_connected_low, - kBeff_connected_low); - fprintf(scal_log_file, "%.5g %.5g ", kAeff_disconnected, - kBeff_disconnected); - fprintf(scal_log_file, "%.5g %.5g %.5g ", pAB, - pAB_connected, Ca); + fprintf(scal_log_file, "%.5g %.5g ", kAeff_low * Mask->Porosity() * mDarcy_converter, + kBeff_low * Mask->Porosity() * mDarcy_converter); + fprintf(scal_log_file, "%.5g %.5g ", kAeff_connected_low * Mask->Porosity() * mDarcy_converter, + kBeff_connected_low * Mask->Porosity() * mDarcy_converter); + fprintf(scal_log_file, "%.5g %.5g ", kAeff_disconnected * Mask->Porosity() * mDarcy_converter, + kBeff_disconnected * Mask->Porosity() * mDarcy_converter); + fprintf(scal_log_file, "%.5g %.5g ", krnf_low * Mask->Porosity() * mDarcy_converter, + krwf_low * Mask->Porosity() * mDarcy_converter); + fprintf(scal_log_file, "%.5g %.5g %.5g ", pAB, + pAB_connected, Ca); fprintf(scal_log_file, "%.5g\n", eff_pres); fclose(scal_log_file); diff --git a/models/ColorModel.h b/models/ColorModel.h index ed7bed63..7089c619 100644 --- a/models/ColorModel.h +++ b/models/ColorModel.h @@ -117,6 +117,7 @@ class ScaLBL_ColorModel { double tauA, tauB, rhoA, rhoB, alpha, beta; double Fx, Fy, Fz, flux; double din, dout, inletA, inletB, outletA, outletB; + const double mDarcy_converter = 1013.0; int Nx, Ny, Nz, N, Np; int rank, nprocx, nprocy, nprocz, nprocs; diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index d65475fd..6d5d41b2 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -20,13 +20,18 @@ #include "models/MRTModel.h" #include "analysis/distance.h" #include "common/ReadMicroCT.h" + + ScaLBL_MRTModel::ScaLBL_MRTModel(int RANK, int NP, const Utilities::MPI &COMM) : rank(RANK), nprocs(NP), Restart(0), timestep(0), timestepMax(0), tau(0), Fx(0), Fy(0), Fz(0), flux(0), din(0), dout(0), mu(0), Nx(0), Ny(0), Nz(0), N(0), Np(0), nprocx(0), nprocy(0), nprocz(0), BoundaryCondition(0), Lx(0), Ly(0), Lz(0), comm(COMM) {} + + ScaLBL_MRTModel::~ScaLBL_MRTModel() {} + void ScaLBL_MRTModel::ReadParams(string filename) { // read the input database db = std::make_shared(filename); @@ -254,7 +259,7 @@ void ScaLBL_MRTModel::Run() { if (WriteHeader) { log_file = fopen("Permeability.csv", "a+"); - fprintf(log_file, "time Fx Fy Fz mu Vs As Js Xs vx vy vz k\n"); + fprintf(log_file, "time Fx Fy Fz mu Vs As Js Xs vx vy vz absperm absperm_adjusted\n"); fclose(log_file); } } @@ -384,22 +389,25 @@ void ScaLBL_MRTModel::Run() { double h = Dm->voxel_length; double absperm = h * h * mu * Mask->Porosity() * flow_rate / force_mag; + double absperm_adjusted = + h * h * mu * Mask->Porosity() * Mask->Porosity() * flow_rate / force_mag; + absperm_adjusted *= 1013.0; // Convert to mDarcy + if (rank == 0) { printf(" %f\n", absperm); FILE *log_file = fopen("Permeability.csv", "a"); fprintf(log_file, "%i %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g " - "%.8g %.8g\n", + "%.8g %.8g %.8g\n", timestep, Fx, Fy, Fz, mu, h * h * h * Vs, h * h * As, - h * Hs, Xs, vax, vay, vaz, absperm); + h * Hs, Xs, vax, vay, vaz, absperm, absperm_adjusted); fclose(log_file); } } } //************************************************************************/ if (rank == 0) - printf("---------------------------------------------------------------" - "----\n"); + printf("--------------------------------------------------------\n"); // Compute the walltime per timestep auto t2 = std::chrono::system_clock::now(); double cputime = std::chrono::duration(t2 - t1).count() / timestep;