From 16ba55100a6283291fddf1a1259884ec4b909372 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://openfoam.org>
Date: Thu, 16 Nov 2017 22:41:35 +0000
Subject: [PATCH] surfaceFilmModels::waxSolventEvaporation: Corrected handling
 of impingement

---
 .../waxSolventEvaporation.C                   | 25 ++++++++++++++-----
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.C b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.C
index 8f27a1530..0228c4a54 100644
--- a/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.C
+++ b/src/regionModels/surfaceFilmModels/submodels/thermo/phaseChangeModel/waxSolventEvaporation/waxSolventEvaporation.C
@@ -316,6 +316,15 @@ void waxSolventEvaporation::correctModel
         ROOTVSMALL/dt
     );
 
+    volScalarField::Internal impingementRate
+    (
+        max
+        (
+           -film.rhoSp()(),
+            dimensionedScalar("zero", film.rhoSp().dimensions(), 0)
+        )
+    );
+
     if (filmPresent)
     {
         // Solve for the solvent mass fraction
@@ -325,17 +334,21 @@ void waxSolventEvaporation::correctModel
           + fvm::div(phi, Ysolvent_)
          ==
             deltaRho0Bydt*Ysolvent_()
+
           + evapRateInf
 
             // Include the effect of the impinging droplets
             // added with Ysolvent = Ysolvent0
-          + max
-            (
-                film.rhoSp()(),
-                dimensionedScalar("zero", film.rhoSp().dimensions(), 0)
-            )*Ysolvent0_
+          + impingementRate*Ysolvent0_
 
-          - fvm::Sp(deltaRho0Bydt + evapRateCoeff, Ysolvent_)
+          - fvm::Sp
+            (
+                deltaRho0Bydt
+              + evapRateCoeff
+              + film.rhoSp()()
+              + impingementRate,
+                Ysolvent_
+            )
         );
 
         YsolventEqn.relax();