Merge branch 'dev_doury' of https://gitlab.com/phorgue/porousMultiphaseFoam into dev_doury

Author: phorgue

solvers/impesEvapFoam/CoatsNo.H

@@ -20,9 +20,14 @@
 
     CFLCoats /= mesh.V();
     CFLUse = gMax(CFLCoats);
-    if (runTime.value() < 500.0) 
-    	maxDeltaTFact = 0.005/(CFLUse + SMALL); 
-    else
+    if (runTime.value() < 100000000.0) 
+    	maxDeltaTFact = (maxCo/3.0)/(CFLUse + SMALL); 
+    else if  (runTime.value() < 1000000000.0) 
+    {
+    	scalar max = (maxCo/3.0) + (runTime.value()-100000000.0) * (maxCo - (maxCo/3.0))/(900000000.0); 
+    	maxDeltaTFact = max/(CFLUse + SMALL);
+    }
+    else  
     	maxDeltaTFact = maxCo/(CFLUse + SMALL);
 
     Info<< "Coats CFL Number mean: " << gAverage(CFLCoats) << " max: " << CFLUse << endl;

solvers/impesEvapFoam/SEqn.H

@@ -40,7 +40,7 @@
     }
 
     Info << " Lost Volume of Water" << "=" << lost_water_vol << "m3" << endl;
-    Info << " Lost Volume of Water" << "=" << gain_water_vol << "m3" << endl;
+    Info << " Gain Volume of Water" << "=" << gain_water_vol << "m3" << endl;
     Info << "Saturation b: " << " Min(Sb) = " << gMin(Sb) << " Max(Sb) = " << gMax(Sb) << endl;
 
 }

solvers/impesEvapFoam/TEqn.H

@@ -13,6 +13,8 @@
 	//	);
 	//TEqn.solve();
 	//Info << T << endl;
+	
+    	// Info << (hwa-hwb) * fvc::interpolate(m) << endl;
 	fvScalarMatrix TEqn
 		(
 			rhoa * eps * (fvc::ddt(u_a) - fvc::ddt(u_a,Sb))
@@ -22,5 +24,6 @@
 			==  fvm::laplacian(lambda_pm,T)  + (hwa-hwb) * m - wallSourceTerm 
 		);
 	TEqn.solve();
+	//Info << "gradT" << rhoa * fvc::div(phia,(Xw*hwa+(1-Xw)*hair)) +  rhob * fvc::div(phib*(hwb)) - fvm::laplacian(lambda_pm,T) << endl; 
 
 }

solvers/impesEvapFoam/createFields.H

@@ -133,14 +133,38 @@ dimensionedScalar hair(transportProperties.lookupOrDefault("hair",dimensionedSca
 // non-wetting phase properties
 autoPtr<incompressiblePhase> phasea = incompressiblePhase::New(mesh,transportProperties,"a");
 volVectorField& Ua = phasea->U();
-surfaceScalarField& phia = phasea->phi();
+//surfaceScalarField& phia = phasea->phi();
+surfaceScalarField phia
+(
+    IOobject
+    (
+        "phia",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    phasea->phi()
+);
 const dimensionedScalar& rhoa = phasea->rho();
 const dimensionedScalar& mua = phasea->mu();
 
 // wetting phase properties
 autoPtr<incompressiblePhase> phaseb = incompressiblePhase::New(mesh,transportProperties,"b");
 volVectorField& Ub = phaseb->U();
-surfaceScalarField& phib = phaseb->phi();
+//surfaceScalarField& phib = phaseb->phi();
+surfaceScalarField phib
+(
+    IOobject
+    (
+        "phib",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    phaseb->phi()
+);
 const dimensionedScalar& rhob = phaseb->rho();
 const dimensionedScalar& mub = phaseb->mu();
 
@@ -491,6 +515,7 @@ forAll(mesh.boundaryMesh(),patchi)
        waterMassBalanceCSV << " Xw";
        waterMassBalanceCSV << " p";
     }
+    waterMassBalanceCSV << "Heat Flux(" << T.boundaryField()[patchi].patch().name() << ")";
 }
 waterMassBalanceCSV << " Sb_var" ;
 waterMassBalanceCSV << endl;

solvers/impesEvapFoam/pEqn.H

@@ -1,21 +1,22 @@
 {
+
     fvScalarMatrix pEqn
         (
-            fvm::laplacian(-Mf, p) +  fvc::div(phiG)
+            fvm::laplacian(-Mf, p) + fvc::div(phiG)
             // capillary term
             + fvc::div(phiPc)*activateCapillarity
             ==
             // event source terms
             - sourceTerm
         );
-    //pEqn.setReference(860,0);
+    // pEqn.setReference(860,0);
     pEqn.solve();
 
     phiP = pEqn.flux();
 
-    phi = phiP+phiG+phiPc*activateCapillarity;//
+    phi = phiP+phiG+phiPc*activateCapillarity;
 
-    phib = Fbf*phiP + (Lbf/Lf)*phiG + phiPc*activateCapillarity ; // ;
+    phib = Fbf*phiP + (Lbf/Lf)*phiG + phiPc*activateCapillarity;
     phia = phi - phib;
 
     U = fvc::reconstruct(phi);

solvers/impesEvapFoam/updateSbProperties.H

@@ -1,3 +1,4 @@
+
 //- relative permeability computation
 krModel->correct(Sb);
 kraf = fvc::interpolate(kra,"kra");

solvers/impesEvapFoam/waterMassBalance.H

@@ -8,9 +8,12 @@
             Info << phib.boundaryField()[patchi].patch().name() << " = " <<  gSum(phib.boundaryField()[patchi]) << " ; ";
         }
     }
-
+    
     Info << endl;
 
+    surfaceScalarField heatFlux = fvc::interpolate(-lambda_pm)*(fvc::snGrad(T)) ;//+ (hwa-hwb) * fvc::interpolate(m);
+    
+    
     //- write fields using usual openfoam rules
     runTime.write();
 
@@ -26,6 +29,7 @@
      	   waterMassBalanceCSV << " " << gSum(Xw.boundaryField()[patchi]);
      	   waterMassBalanceCSV << " " << gSum(p.boundaryField()[patchi]);
          }
+        waterMassBalanceCSV << " " << gSum(heatFlux.boundaryField()[patchi]);
       } 
      waterMassBalanceCSV << " " << fvc::domainIntegrate(fvc::ddt(Sb)*eps).value();
      waterMassBalanceCSV << endl;