Updated Leaky Dielectric Model for Electrohydrodynamics

Allwmake

@@ -1,5 +1,5 @@
 #!/bin/bash
-
+. "$WM_PROJECT_DIR"/wmake/scripts/AllwmakeParseArguments
 wmake libso src/VoF
 wmake all src/boundaryConditions
 wmake all src/thermoDynamics

solver/leakyDielectricEHD/EHD/casesEHD.H

@@ -0,0 +1,36 @@
+Info << "Coupling Maxwell tensor with UEqn: " << coupleMaxwellTensor << endl;
+
+Info << "EHD correction loop: " << ehdCorr <<endl;
+
+if (ionConductivityAvg == "arithmetic")
+{
+    Info << "Using arithmetic averaging for ionConductivityAvg" << endl;
+}
+else if (ionConductivityAvg == "harmonic")
+{
+    Info << "Using harmonic averaging for ionConductivityAvg" << endl;
+}
+else
+{
+    FatalErrorInFunction
+        << "Unknown interfaceAvg for ionConductivityAvg: " << ionConductivityAvg << nl
+        << "Valid options are: arithmetic, harmonic"
+        << exit(FatalError);
+}
+
+if (permittivityAvg == "arithmetic")
+{
+    Info << "Using arithmetic averaging for permittivityAvg" << endl;
+}
+else if (permittivityAvg == "harmonic")
+{
+    Info << "Using harmonic averaging for permittivityAvg" << endl;
+}
+else
+{
+    FatalErrorInFunction
+        << "Unknown interfaceAvg for permittivityAvg: " << permittivityAvg << nl
+        << "Valid options are: arithmetic, harmonic"
+        << exit(FatalError);
+}
+

solver/leakyDielectricEHD/EHD/createEHDFields.H

@@ -0,0 +1,202 @@
+Info<< " Reading electrical properties\n" << nl << endl;
+
+IOdictionary electricalProperties
+(
+    IOobject
+    (
+        "electricalProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+dimensionedScalar permittivity1
+    (
+        "permittivity1",
+        dimensionSet(-1, -3, 4, 0, 0, 2, 0),
+        electricalProperties
+    );
+dimensionedScalar permittivity2
+    (
+        "permittivity2",
+        dimensionSet(-1, -3, 4, 0, 0, 2, 0),
+        electricalProperties
+    );
+    dimensionedScalar ionConductivity1
+    (
+        "ionConductivity1",
+        dimensionSet(-1, -3, 3, 0, 0, 2, 0),
+        electricalProperties
+    );
+ dimensionedScalar ionConductivity2
+    (
+        "ionConductivity2",
+        dimensionSet(-1, -3, 3, 0, 0, 2, 0),
+        electricalProperties
+    );
+
+Info<< "Reading field chargeDensity\n" << endl;
+
+volScalarField chargeDensity
+(
+    IOobject
+    (
+        "chargeDensity",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field potentialDifference\n" << endl;
+volScalarField potentialDifference
+(
+    IOobject
+    (
+        "potentialDifference",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+
+Info<< "Reading field ionConductivity\n" << endl;
+
+
+volScalarField ionConductivity
+(
+    IOobject
+    (
+        "ionConductivity",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    (alpha1*ionConductivity1)+(alpha2*ionConductivity2)
+);
+// ionConductivity.oldTime();
+
+surfaceScalarField ionConductivityF
+(
+    IOobject
+    (
+        "ionConductivityF",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    (fvc::interpolate(alpha1)*ionConductivity1 + fvc::interpolate(alpha2)*ionConductivity2)
+);
+
+Info<< "Reading field permittivity\n" << endl;
+volScalarField permittivity
+(
+    IOobject
+    (
+        "permittivity",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    (alpha1*permittivity1)+(alpha2*permittivity2)
+);
+// permittivity.oldTime();
+
+surfaceScalarField permittivityF
+(
+    IOobject
+    (
+        "permittivityF",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    (fvc::interpolate(alpha1)*permittivity1 + fvc::interpolate(alpha2)*permittivity2)
+);
+
+volVectorField electricField
+(
+    IOobject
+    (
+        "electricField",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedVector("electricField", dimensionSet(1,1,-3,0,0,-1,0), vector::zero)
+);
+
+Info<< "Calculating field maxwellTensor\n" << endl;
+volTensorField maxwellTensor
+(
+    IOobject
+    (
+        "maxwellTensor",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    permittivity*((electricField*electricField)-(double(0.5)*magSqr(electricField)*tensor::I))
+);
+
+volVectorField ehdForce
+(
+    IOobject
+    (
+        "ehdForce",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    dimensionedVector("ehdForce", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+);
+
+
+const dictionary& ehdDict =
+    mesh.solutionDict().subDict("ehdSolution");
+
+Switch coupleConvChargeDensity
+(
+    ehdDict.lookupOrDefault<Switch>("coupleConvChargeDensity", true)
+);
+
+Switch coupleMaxwellTensor
+(
+    ehdDict.lookupOrDefault<Switch>("coupleMaxwellTensor", true)
+);
+
+label ehdCorr
+(
+    ehdDict.lookupOrDefault<label>("ehdCorr", 1)
+);
+
+word ionConductivityAvg
+(
+    ehdDict.lookupOrDefault<word>("ionConductivityAvg", "arithmetic")
+);
+
+word permittivityAvg
+(
+    ehdDict.lookupOrDefault<word>("permittivityAvg", "arithmetic")
+);
+
+#include "casesEHD.H"
+
+
+

solver/leakyDielectricEHD/EHD/ehdEqn.H

@@ -0,0 +1,25 @@
+/*********** Solving Herrera Model ************/
+
+// Inspired from Basilisk method of equation solving
+fvScalarMatrix potentialDifferenceEqn
+(
+    fvm::laplacian(ionConductivityF, potentialDifference)
+);
+// potentialDifferenceEqn.relax();
+potentialDifferenceEqn.solve();
+
+// Calculating electric Field
+electricField = -fvc::grad(potentialDifference);
+
+//(LopezHerrera2011, Roghair2015)
+chargeDensity = -fvc::laplacian(permittivityF, potentialDifference);
+
+permittivity = fvc::average(permittivityF);
+ionConductivity = fvc::average(ionConductivityF);
+
+maxwellTensor = permittivity*((electricField*electricField)-(double(0.5)*magSqr(electricField)*tensor::I));
+
+if (coupleMaxwellTensor)
+{
+    ehdForce =  fvc::div(maxwellTensor);
+}

solver/leakyDielectricEHD/Make/files

@@ -0,0 +1,3 @@
+leakyDielectricEHD.C
+
+EXE = $(FOAM_APPBIN)/leakyDielectricEHD

solver/leakyDielectricEHD/Make/options

@@ -0,0 +1,30 @@
+EXE_INC = \
+    -I./EHD \
+    -I$(LIB_SRC)/transportModels/twoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/immiscibleIncompressibleTwoPhaseMixture/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I../../src/VoF/lnInclude \
+    -I../../src/surfaceForces/lnInclude
+
+EXE_LIBS = \
+    -limmiscibleIncompressibleTwoPhaseMixture \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -ldynamicFvMesh \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -lwaveModels \
+    -L$(FOAM_LIBBIN) \
+    -lVoF \
+    -lsurfaceForces
+

solver/leakyDielectricEHD/UEqn.H

@@ -0,0 +1,39 @@
+    MRF.correctBoundaryVelocity(U);
+
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
+      + MRF.DDt(rho, U)
+      + turbulence->divDevRhoReff(rho, U)
+     ==
+        fvOptions(rho, U)
+    );
+    if (coupleMaxwellTensor)
+    {
+        UEqn -= fvc::div(maxwellTensor);
+    }
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    //mixture.surfaceTensionForce()
+                  surfForces.surfaceTensionForce()
+                  //- ghf*fvc::snGrad(rho)
+                  + surfForces.accelerationForce()
+                  - fvc::snGrad(p_rgh)
+                ) * mesh.magSf()
+            )
+        );
+
+        fvOptions.correct(U);
+    }

solver/leakyDielectricEHD/alphaControls.H

@@ -0,0 +1,3 @@
+const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+
+label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));