From dffa98cb0b658ac43900d0506fe02a5fab96d873 Mon Sep 17 00:00:00 2001 From: James E McClure Date: Wed, 2 Jun 2021 15:09:51 -0400 Subject: [PATCH 1/5] added films to timelog.csv in color --- analysis/SubPhase.cpp | 46 ++++++++++++++++++++++++++++++++++--------- 1 file changed, 37 insertions(+), 9 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 2fd22993..d74d9857 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -93,7 +93,7 @@ SubPhase::SubPhase(std::shared_ptr dm): { // If timelog is empty, write a short header to list the averages //fprintf(TIMELOG,"--------------------------------------------------------------------------------------\n"); - fprintf(TIMELOG,"sw krw krn vw vn force pw pn wet\n"); + fprintf(TIMELOG,"sw krw krn krwf krnf vw vn force pw pn wet\n"); } } } @@ -167,7 +167,7 @@ void SubPhase::Basic(){ if (Dm->inlet_layers_z > 0 && Dm->kproc() == 0) kmin += Dm->inlet_layers_z; if (Dm->outlet_layers_z > 0 && Dm->kproc() == Dm->nprocz()-1) kmax -= Dm->outlet_layers_z; */ - nb.reset(); wb.reset(); + nb.reset(); wb.reset(); iwn.reset(); double count_w = 0.0; double count_n = 0.0; @@ -245,6 +245,28 @@ void SubPhase::Basic(){ wb.p += Pressure(n); count_w += 1.0; } + /* compute the film contribution */ + else if (SDs(i,j,k) < 2.0){ + if ( phi > 0.0 ){ + nA = 1.0; + iwn.V += 1.0; + iwn.M += nA*rho_n; + // velocity + iwn.Px += rho_n*nA*Vel_x(n); + iwn.Py += rho_n*nA*Vel_y(n); + iwn.Pz += rho_n*nA*Vel_z(n); + } + else{ + nB = 1.0; + iwn.M += nB*rho_w; + iwn.V += 1.0; + + // velocity + iwn.Px += rho_w*nB*Vel_x(n); + iwn.Py += rho_w*nB*Vel_y(n); + iwn.Pz += rho_w*nB*Vel_z(n); + } + } } } } @@ -267,12 +289,6 @@ void SubPhase::Basic(){ //printf("wetting interaction = %f, count = %f\n",total_wetting_interaction,count_wetting_interaction); total_wetting_interaction_global=Dm->Comm.sumReduce( total_wetting_interaction); count_wetting_interaction_global=Dm->Comm.sumReduce( count_wetting_interaction); - /* normalize wetting interactions <-- Don't do this if normalizing laplacian (use solid surface area) - if (count_wetting_interaction > 0.0) - total_wetting_interaction /= count_wetting_interaction; - if (count_wetting_interaction_global > 0.0) - total_wetting_interaction_global /= count_wetting_interaction_global; - */ gwb.V=Dm->Comm.sumReduce( wb.V); gnb.V=Dm->Comm.sumReduce( nb.V); @@ -285,6 +301,11 @@ void SubPhase::Basic(){ gnb.Py=Dm->Comm.sumReduce( nb.Py); gnb.Pz=Dm->Comm.sumReduce( nb.Pz); + giwn.M=Dm->Comm.sumReduce( iwn.M); + giwn.Px=Dm->Comm.sumReduce( iwn.Px); + giwn.Py=Dm->Comm.sumReduce( iwn.Py); + giwn.Pz=Dm->Comm.sumReduce( iwn.Pz); + count_w=Dm->Comm.sumReduce( count_w); count_n=Dm->Comm.sumReduce( count_n); if (count_w > 0.0) @@ -341,14 +362,21 @@ void SubPhase::Basic(){ double saturation=gwb.V/(gwb.V + gnb.V); double water_flow_rate=gwb.V*(gwb.Px*dir_x + gwb.Py*dir_y + gwb.Pz*dir_z)/gwb.M / Dm->Volume; double not_water_flow_rate=gnb.V*(gnb.Px*dir_x + gnb.Py*dir_y + gnb.Pz*dir_z)/gnb.M/ Dm->Volume; + + /* contribution from water films */ + double water_film_flow_rate=gwb.V*(giwn.Pwx*dir_x + giwn.Pwy*dir_y + giwn.Pwz*dir_z)/gwb.M / Dm->Volume; + double not_water_film_flow_rate=gnb.V*(giwn.Pnx*dir_x + giwn.Pny*dir_y + giwn.Pnz*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 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; + /* not counting films */ + double krnf = krn - h*h*nu_n*not_water_film_flow_rate / force_mag ; + double krwf = krw - h*h*nu_w*water_film_flow_rate / force_mag; //printf(" water saturation = %f, fractional flow =%f \n",saturation,fractional_flow); - fprintf(TIMELOG,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",saturation,krw,krn,h*water_flow_rate,h*not_water_flow_rate, force_mag, gwb.p, gnb.p, total_wetting_interaction_global); + fprintf(TIMELOG,"%.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g %.8g\n",saturation,krw,krn,krwf,krnf,h*water_flow_rate,h*not_water_flow_rate, force_mag, gwb.p, gnb.p, total_wetting_interaction_global); fflush(TIMELOG); } if (err==true){ From f5fbde02d42312d7057a3317ac8dcb03d3cc4b8b Mon Sep 17 00:00:00 2001 From: James E McClure Date: Thu, 3 Jun 2021 16:19:31 -0400 Subject: [PATCH 2/5] fix issue with film in timelog.csv --- analysis/SubPhase.cpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index d74d9857..6964197c 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -301,10 +301,15 @@ void SubPhase::Basic(){ gnb.Py=Dm->Comm.sumReduce( nb.Py); gnb.Pz=Dm->Comm.sumReduce( nb.Pz); - giwn.M=Dm->Comm.sumReduce( iwn.M); - giwn.Px=Dm->Comm.sumReduce( iwn.Px); - giwn.Py=Dm->Comm.sumReduce( iwn.Py); - giwn.Pz=Dm->Comm.sumReduce( iwn.Pz); + giwn.Mw=Dm->Comm.sumReduce( iwn.Mw); + giwn.Pwx=Dm->Comm.sumReduce( iwn.Pwx); + giwn.Pwy=Dm->Comm.sumReduce( iwn.Pwy); + giwn.Pwz=Dm->Comm.sumReduce( iwn.Pwz); + + giwn.Mn=Dm->Comm.sumReduce( iwn.Mn); + giwn.Pnx=Dm->Comm.sumReduce( iwn.Pnx); + giwn.Pny=Dm->Comm.sumReduce( iwn.Pny); + giwn.Pnz=Dm->Comm.sumReduce( iwn.Pnz); count_w=Dm->Comm.sumReduce( count_w); count_n=Dm->Comm.sumReduce( count_n); From 96ddd52b079feae8675c517ced0b694af351520f Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 4 Jun 2021 09:53:36 -0400 Subject: [PATCH 3/5] fix bug in film --- analysis/SubPhase.cpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 6964197c..5830c4eb 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -252,19 +252,18 @@ void SubPhase::Basic(){ iwn.V += 1.0; iwn.M += nA*rho_n; // velocity - iwn.Px += rho_n*nA*Vel_x(n); - iwn.Py += rho_n*nA*Vel_y(n); - iwn.Pz += rho_n*nA*Vel_z(n); + iwn.Pnx += rho_n*nA*Vel_x(n); + iwn.Pny += rho_n*nA*Vel_y(n); + iwn.Pnz += rho_n*nA*Vel_z(n); } else{ nB = 1.0; iwn.M += nB*rho_w; iwn.V += 1.0; - // velocity - iwn.Px += rho_w*nB*Vel_x(n); - iwn.Py += rho_w*nB*Vel_y(n); - iwn.Pz += rho_w*nB*Vel_z(n); + iwn.Pwx += rho_w*nB*Vel_x(n); + iwn.Pwy += rho_w*nB*Vel_y(n); + iwn.Pwz += rho_w*nB*Vel_z(n); } } } From 5c474b1285a89f22c72d080c97728c869422175c Mon Sep 17 00:00:00 2001 From: James E McClure Date: Fri, 4 Jun 2021 10:10:56 -0400 Subject: [PATCH 4/5] fix bug in film --- analysis/SubPhase.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 5830c4eb..68866a05 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -250,7 +250,7 @@ void SubPhase::Basic(){ if ( phi > 0.0 ){ nA = 1.0; iwn.V += 1.0; - iwn.M += nA*rho_n; + iwn.Mn += nA*rho_n; // velocity iwn.Pnx += rho_n*nA*Vel_x(n); iwn.Pny += rho_n*nA*Vel_y(n); @@ -258,7 +258,7 @@ void SubPhase::Basic(){ } else{ nB = 1.0; - iwn.M += nB*rho_w; + iwn.Mw += nB*rho_w; iwn.V += 1.0; iwn.Pwx += rho_w*nB*Vel_x(n); From 87c80399190601250d2386ec4b2e85a6455fd27a Mon Sep 17 00:00:00 2001 From: James E McClure Date: Sun, 6 Jun 2021 16:33:46 -0400 Subject: [PATCH 5/5] fix flow orientation for unsteady flow --- analysis/SubPhase.cpp | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/analysis/SubPhase.cpp b/analysis/SubPhase.cpp index 68866a05..42197ad8 100644 --- a/analysis/SubPhase.cpp +++ b/analysis/SubPhase.cpp @@ -335,20 +335,21 @@ void SubPhase::Basic(){ if (gnb.Pz != gnb.Pz) err=true; if (Dm->rank() == 0){ + /* align flow direction based on total mass flux */ + double dir_x = gwb.Px + gnb.Px; + double dir_y = gwb.Py + gnb.Py; + double dir_z = gwb.Pz + gnb.Pz; + double flow_magnitude = dir_x*dir_x + dir_y*dir_y + dir_z*dir_z; double force_mag = sqrt(Fx*Fx+Fy*Fy+Fz*Fz); - double dir_x = 0.0; - double dir_y = 0.0; - double dir_z = 0.0; if (force_mag > 0.0){ dir_x = Fx/force_mag; dir_y = Fy/force_mag; dir_z = Fz/force_mag; } else { - // default to z direction - dir_x = 0.0; - dir_y = 0.0; - dir_z = 1.0; + dir_x /= flow_magnitude; + dir_y /= flow_magnitude; + dir_z /= flow_magnitude; } if (Dm->BoundaryCondition == 1 || Dm->BoundaryCondition == 2 || Dm->BoundaryCondition == 3 || Dm->BoundaryCondition == 4 ){ // compute the pressure drop @@ -356,7 +357,7 @@ void SubPhase::Basic(){ double length = ((Nz-2)*Dm->nprocz()); force_mag -= pressure_drop/length; } - if (force_mag == 0.0){ + if (force_mag == 0.0 && flow_magnitude == 0.0){ // default to z direction dir_x = 0.0; dir_y = 0.0;