diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 265222c4..9c2f4529 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -249,11 +249,12 @@ void SubPhase::Basic(){ double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume; double total_flow_rate = water_flow_rate + not_water_flow_rate; double fractional_flow= water_flow_rate / total_flow_rate; - - double krn = nu_n*not_water_flow_rate / force_mag ; - double krw = nu_w*water_flow_rate / force_mag; + + double h = Dm->voxel_length; + 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,water_flow_rate,not_water_flow_rate, gwb.p, gnb.p); + 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); } diff --git a/models/ColorModel.cpp b/models/ColorModel.cpp index 429fddb1..a254982f 100644 --- a/models/ColorModel.cpp +++ b/models/ColorModel.cpp @@ -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); @@ -687,7 +711,7 @@ void ScaLBL_ColorModel::Run(){ morph_timesteps=0; Ca_previous = Ca; } - if (MORPH_ADAPT ){ + if (MORPH_ADAPT && morph_delta != 0.0){ CURRENT_MORPH_TIMESTEPS += analysis_interval; if (rank==0) printf("***Morphological step with target volume change %f ***\n", delta_volume_target); //double delta_volume_target = volB - (volA + volB)*TARGET_SATURATION; // change in volume to A diff --git a/models/MRTModel.cpp b/models/MRTModel.cpp index 1eb789d9..d04009c1 100644 --- a/models/MRTModel.cpp +++ b/models/MRTModel.cpp @@ -264,7 +264,7 @@ void ScaLBL_MRTModel::Run(){ 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",timestep, Fx, Fy, Fz, mu, - Vs,As,Hs,Xs,vax,vay,vaz, absperm); + h*h*h*Vs,h*h*As,h*Hs,Xs,vax,vay,vaz, absperm); fclose(log_file); } }