ENH: interfaceProperties: add smoother for interfacial curvatures (closes #2531)

The new algorithm introduces an optional entry 'nAlphaSmoothCurvature',
and aims to smooth interface curvatures to reduce spurious currents.

src/transportModels/interfaceProperties/interfaceProperties.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,6 +33,7 @@ License
 #include "fvcDiv.H"
 #include "fvcGrad.H"
 #include "fvcSnGrad.H"
+#include "fvcAverage.H"
 #include "unitConversion.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
@@ -109,10 +111,27 @@ void Foam::interfaceProperties::calculateK()
     const surfaceVectorField& Sf = mesh.Sf();
 
     // Cell gradient of alpha
-    const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat"));
+    tmp<volVectorField> tgradAlpha;
+    if (nAlphaSmoothCurvature_ < 1)
+    {
+        tgradAlpha = fvc::grad(alpha1_, "nHat");
+    }
+    else
+    {
+        // Smooth interface curvature to reduce spurious currents
+        auto talpha1L = tmp<volScalarField>::New(alpha1_);
+        auto& alpha1L = talpha1L.ref();
+
+        for (int i = 0; i < nAlphaSmoothCurvature_; ++i)
+        {
+            alpha1L = fvc::average(fvc::interpolate(alpha1L));
+        }
+
+        tgradAlpha = fvc::grad(talpha1L, "nHat");
+    }
 
     // Interpolated face-gradient of alpha
-    surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
+    surfaceVectorField gradAlphaf(fvc::interpolate(tgradAlpha));
 
     //gradAlphaf -=
     //    (mesh.Sf()/mesh.magSf())
@@ -157,6 +176,11 @@ Foam::interfaceProperties::interfaceProperties
 )
 :
     transportPropertiesDict_(dict),
+    nAlphaSmoothCurvature_
+    (
+        alpha1.mesh().solverDict(alpha1.name()).
+            getOrDefault<int>("nAlphaSmoothCurvature", 0)
+    ),
     cAlpha_
     (
         alpha1.mesh().solverDict(alpha1.name()).get<scalar>("cAlpha")

src/transportModels/interfaceProperties/interfaceProperties.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,6 +34,15 @@ Description
     -# Correct the alpha boundary condition for dynamic contact angle.
     -# Calculate interface curvature.
 
+    References:
+    \verbatim
+        Smoother algorithm for interface curvatures (tag:YK):
+            Yamamoto, T. & Komarov, S. V. (2022).
+            Evaluation on different volume of fluid methods in
+            unstructured solver under the optimized condition. Preprint.
+            DOI:10.13140/RG.2.2.17223.37281
+    \endverbatim
+
 SourceFiles
     interfaceProperties.C
 
@@ -62,6 +72,9 @@ class interfaceProperties
         //- Keep a reference to the transportProperties dictionary
         const dictionary& transportPropertiesDict_;
 
+        //- Number of iterations to smooth interfacial curvatures
+        int nAlphaSmoothCurvature_;
+
         //- Compression coefficient
         scalar cAlpha_;