OilfieldToolsNet/LBPM
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

modfiy relperm.csv output structure

Browse Source
This commit is contained in:
James E McClure
2019-04-24 14:16:09 -04:00
parent 2efa0c8fc3
commit a95d6b63a7
2 changed files with 44 additions and 20 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
analysis/SubPhase.cpp
View File

@@ -251,8 +251,8 @@ void SubPhase::Basic(){
double fractional_flow= water_flow_rate / total_flow_rate;
double h = Dm->voxel_length;
double krn = h*h*rho_n*nu_n*not_water_flow_rate / force_mag ;
double krw = h*h*rho_w*nu_w*water_flow_rate / force_mag;
double krn = h*h*nu_n*not_water_flow_rate / force_mag ;
double krw = h*h*nu_w*water_flow_rate / force_mag;
//printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow);
fprintf(TIMELOG,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",saturation,krw,krn,h*h*h*water_flow_rate,h*h*h*not_water_flow_rate, gwb.p, gnb.p);
fflush(TIMELOG);

60
models/ColorModel.cpp
View File

@@ -611,6 +611,8 @@ void ScaLBL_ColorModel::Run(){
if ( morph_timesteps > morph_interval ){
double volB = Averages->gwb.V;
double volA = Averages->gnb.V;
volA /= double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
volB /= double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
double vA_x = Averages->gnb.Px/Averages->gnb.M;
double vA_y = Averages->gnb.Py/Averages->gnb.M;
double vA_z = Averages->gnb.Pz/Averages->gnb.M;
@@ -624,12 +626,18 @@ void ScaLBL_ColorModel::Run(){
double dir_x = Fx/force_mag;
double dir_y = Fy/force_mag;
double dir_z = Fz/force_mag;
double flow_rate_A = (vA_x*dir_x + vA_y*dir_y + vA_z*dir_z);
double flow_rate_B = (vB_x*dir_x + vB_y*dir_y + vB_z*dir_z);
if (force_mag == 0.0){
// default to z direction
dir_x = 0.0;
dir_y = 0.0;
dir_z = 1.0;
force_mag = 1.0;
}
double current_saturation = volB/(volA+volB);
double Ca = fabs(volA*muA*flow_rate_A + volB*muB*flow_rate_B)/(5.796*alpha*double((Nx-2)*(Ny-2)*(Nz-2)*nprocs));
double force_magnitude = sqrt(Fx*Fx + Fy*Fy + Fz*Fz);
double flow_rate_A = volA*(vA_x*dir_x + vA_y*dir_y + vA_z*dir_z);
double flow_rate_B = volB*(vB_x*dir_x + vB_y*dir_y + vB_z*dir_z);
double Ca = fabs(muA*flow_rate_A + muB*flow_rate_B)/(5.796*alpha);
bool isSteady = false;
if ( (fabs((Ca - Ca_previous)/Ca) < tolerance && CURRENT_STEADY_TIMESTEPS > MIN_STEADY_TIMESTEPS))
@@ -649,11 +657,27 @@ void ScaLBL_ColorModel::Run(){
if (rank==0){
printf("** WRITE STEADY POINT *** ");
printf("Ca = %f, (previous = %f) \n",Ca,Ca_previous);
volA /= double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
volB /= double((Nx-2)*(Ny-2)*(Nz-2)*nprocs);
FILE * kr_log_file = fopen("relperm.csv","a");
fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g ",CURRENT_STEADY_TIMESTEPS,muA,muB,5.796*alpha,Fx,Fy,Fz);
fprintf(kr_log_file,"%.5g %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",volA,volB,vA_x,vA_y,vA_z,vB_x,vB_y,vB_z);
double pA = Averages->gnb.p;
double pB = Averages->gwb.p;
double h = Dm->voxel_length;
double kAeff = h*h*muA*flow_rate_A/(rhoA*force_mag);
double kBeff = h*h*muB*flow_rate_B/(rhoB*force_mag);
double pAB = (pA-pB)/(h*5.796*alpha);
double viscous_pressure_drop = (rhoA*volA + rhoB*volB)*force_mag;
double Mobility = muA/muB;
bool WriteHeader=false;
FILE * kr_log_file = fopen("relperm.csv","r");
if (kr_log_file != NULL)
fclose(kr_log_file);
else
WriteHeader=true;
kr_log_file = fopen("relperm.csv","a");
if (WriteHeader)
fprintf(kr_log_file,"timesteps sat.water eff.perm.oil eff.perm.water cap.pressure pressure.drop Ca M\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,pAB,viscous_pressure_drop,Ca,Mobility);
fprintf(kr_log_file,"%i %.5g %.5g %.5g %.5g %.5g %.5g %.5g\n",CURRENT_STEADY_TIMESTEPS,current_saturation,kAeff,kBeff,pAB,viscous_pressure_drop,Ca,Mobility);
fclose(kr_log_file);
printf(" Measured capillary number %f \n ",Ca);
@@ -664,15 +688,15 @@ void ScaLBL_ColorModel::Run(){
Fy *= capillary_number / Ca;
Fz *= capillary_number / Ca;
if (force_magnitude > 1e-3){
Fx *= 1e-3/force_magnitude; // impose ceiling for stability
Fy *= 1e-3/force_magnitude;
Fz *= 1e-3/force_magnitude;
if (force_mag > 1e-3){
Fx *= 1e-3/force_mag; // impose ceiling for stability
Fy *= 1e-3/force_mag;
Fz *= 1e-3/force_mag;
}
if (force_magnitude < 1e-7){
Fx *= 1e-7/force_magnitude; // impose floor
Fy *= 1e-7/force_magnitude;
Fz *= 1e-7/force_magnitude;
if (force_mag < 1e-7){
Fx *= 1e-7/force_mag; // impose floor
Fy *= 1e-7/force_mag;
Fz *= 1e-7/force_mag;
}
if (rank == 0) printf(" -- adjust force by factor %f \n ",capillary_number / Ca);
Averages->SetParams(rhoA,rhoB,tauA,tauB,Fx,Fy,Fz,alpha,beta);