diff --git a/Allwmake b/Allwmake
index 92cf0e71..51866810 100755
--- a/Allwmake
+++ b/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
diff --git a/solver/leakyDielectricEHD/EHD/casesEHD.H b/solver/leakyDielectricEHD/EHD/casesEHD.H
new file mode 100644
index 00000000..67129b0f
--- /dev/null
+++ b/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);
+}
+
diff --git a/solver/leakyDielectricEHD/EHD/createEHDFields.H b/solver/leakyDielectricEHD/EHD/createEHDFields.H
new file mode 100644
index 00000000..aa67c051
--- /dev/null
+++ b/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"
+
+
+
diff --git a/solver/leakyDielectricEHD/EHD/ehdEqn.H b/solver/leakyDielectricEHD/EHD/ehdEqn.H
new file mode 100644
index 00000000..3ced93f2
--- /dev/null
+++ b/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);
+}
diff --git a/solver/leakyDielectricEHD/Make/files b/solver/leakyDielectricEHD/Make/files
new file mode 100644
index 00000000..5c75d3c0
--- /dev/null
+++ b/solver/leakyDielectricEHD/Make/files
@@ -0,0 +1,3 @@
+leakyDielectricEHD.C
+
+EXE = $(FOAM_APPBIN)/leakyDielectricEHD
diff --git a/solver/leakyDielectricEHD/Make/options b/solver/leakyDielectricEHD/Make/options
new file mode 100644
index 00000000..50eba11b
--- /dev/null
+++ b/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
+    
\ No newline at end of file
diff --git a/solver/leakyDielectricEHD/UEqn.H b/solver/leakyDielectricEHD/UEqn.H
new file mode 100644
index 00000000..17cdc3c0
--- /dev/null
+++ b/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);
+    }
diff --git a/solver/leakyDielectricEHD/alphaControls.H b/solver/leakyDielectricEHD/alphaControls.H
new file mode 100644
index 00000000..2e94f34d
--- /dev/null
+++ b/solver/leakyDielectricEHD/alphaControls.H
@@ -0,0 +1,3 @@
+const dictionary& alphaControls = mesh.solverDict(alpha1.name());
+
+label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
diff --git a/solver/leakyDielectricEHD/alphaCourantNo.H b/solver/leakyDielectricEHD/alphaCourantNo.H
new file mode 100644
index 00000000..4430b7f8
--- /dev/null
+++ b/solver/leakyDielectricEHD/alphaCourantNo.H
@@ -0,0 +1,68 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    alphaCourantNo
+
+Description
+    Calculates and outputs the mean and maximum Courant Numbers.
+
+\*---------------------------------------------------------------------------*/
+
+scalar maxAlphaCo
+(
+    runTime.controlDict().get<scalar>("maxAlphaCo")
+);
+
+scalar maxCapillaryNum
+(
+    runTime.controlDict().get<scalar>("maxCapillaryNum")
+);
+
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+scalar capillaryNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        mixture.nearInterface()().primitiveField()
+       *fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+    
+    capillaryNum = runTime.deltaTValue()/surfForces.capillaryDt(rho1,rho2);
+
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+Info<< "Capillary Number: " << capillaryNum << endl;
+
+// ************************************************************************* //
diff --git a/solver/leakyDielectricEHD/alphaEqn.H b/solver/leakyDielectricEHD/alphaEqn.H
new file mode 100644
index 00000000..f8327610
--- /dev/null
+++ b/solver/leakyDielectricEHD/alphaEqn.H
@@ -0,0 +1,18 @@
+{
+// Update alpha1
+#include "alphaSuSp.H"
+advector->advect(Sp,(Su + divU*min(alpha1(), scalar(1)))());
+// advector->advect();
+#include "rhofs.H"
+rhoPhi = advector->getRhoPhi(rho1f, rho2f);
+
+alpha2 = 1.0 - alpha1;
+// mixture.correct();
+
+}
+
+Info<< "Phase-1 volume fraction = "
+    << alpha1.weightedAverage(mesh.Vsc()).value()
+    << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+    << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
+    << endl;
diff --git a/solver/leakyDielectricEHD/alphaEqnSubCycle.H b/solver/leakyDielectricEHD/alphaEqnSubCycle.H
new file mode 100644
index 00000000..4321b8b7
--- /dev/null
+++ b/solver/leakyDielectricEHD/alphaEqnSubCycle.H
@@ -0,0 +1,113 @@
+if (pimple.nCorrPIMPLE() > 1)
+{
+    // If nOuterCorrectors > 1 then for all but the first loop the advection
+    // of alpha is done using an average, 0.5*phi+0.5*phiNew where phi is
+    // the flux at the beginning of the time step and phiNew is the flux
+    // estimate at the end of the time step from the previous outer
+    // iteration. Similarly we use 0.5*U + 0.5*UNew in later iterations.
+    if (pimple.firstIter())
+    {
+        // To recalculate the alpha1 update in subsequent iterations, we
+        // must store its current value before overwriting with the new
+        // value
+        alpha1.storePrevIter();
+        // Storing initial phi and U for use in later outer iterations.
+        phi.storePrevIter();
+        U.storePrevIter();
+    }
+    else
+    {
+        // Resetting alpha1 to value before advection in first PIMPLE
+        // iteration.
+        alpha1 = alpha1.prevIter();
+
+        // Setting U and phi with which to advect interface.
+        U = 0.5*U.prevIter() + 0.5*U;
+        phi = 0.5*phi.prevIter() + 0.5*phi;
+    }
+}
+
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum
+    (
+        IOobject
+        (
+            "rhoPhiSum",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar(rhoPhi.dimensions(), Zero)
+    );
+
+    tmp<volScalarField> trSubDeltaT;
+
+    if (LTS)
+    {
+        trSubDeltaT =
+            fv::localEulerDdt::localRSubDeltaT(mesh, nAlphaSubCycles);
+    }
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+        #include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "alphaEqn.H"
+}
+
+if (!pimple.firstIter())
+{
+    // Resetting U and phi to value at latest iteration.
+    U = 2.0*U - U.prevIter();
+    phi = 2.0*phi - phi.prevIter();
+}
+
+rho == alpha1*rho1 + alpha2*rho2;
+
+if (permittivityAvg == "arithmetic")
+{
+    permittivityF == fvc::interpolate(alpha1)*permittivity1 + fvc::interpolate(alpha2)*permittivity2;
+}
+else if (permittivityAvg == "harmonic")
+{
+    permittivityF == permittivity1*permittivity2/(fvc::interpolate(alpha1)*permittivity2 + fvc::interpolate(alpha2)*permittivity1);
+}
+
+if (ionConductivityAvg == "arithmetic")
+{
+    ionConductivityF == fvc::interpolate(alpha1)*ionConductivity1 + fvc::interpolate(alpha2)*ionConductivity2;
+}
+else if (ionConductivityAvg == "harmonic")
+{
+    // if(
+    //     max(ionConductivity1.value(),ionConductivity2.value()) / min(ionConductivity1.value(),min(ionConductivity2.value(),SMALL)) > 1E10
+    // )
+    // {
+    //     forAll( alpha1, celli)
+    //     {
+    //     ionConductivity[celli] = (alpha1[celli] < 0.01) ? 
+    //         ionConductivity2.value() : 
+    //         ((alpha1[celli] > 0.99) ? 
+    //             ionConductivity1.value() : 
+    //             (scalar(1)/((alpha1[celli]/ionConductivity1.value())+((scalar(1)-alpha1[celli])/ionConductivity2.value()))));
+    //     }
+    // }
+    // else
+    // {
+    //     ionConductivity == ionConductivity1*ionConductivity2/(alpha1*ionConductivity2 + alpha2*ionConductivity1);
+    // }
+    ionConductivityF == ionConductivity1*ionConductivity2/(fvc::interpolate(alpha1)*ionConductivity2 + fvc::interpolate(alpha2)*ionConductivity1);
+}
+
diff --git a/solver/leakyDielectricEHD/alphaSuSp.H b/solver/leakyDielectricEHD/alphaSuSp.H
new file mode 100644
index 00000000..9145f394
--- /dev/null
+++ b/solver/leakyDielectricEHD/alphaSuSp.H
@@ -0,0 +1,30 @@
+volScalarField::Internal Sp
+(
+    IOobject
+    (
+        "Sp",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar(dimless/dimTime, Zero)
+);
+
+volScalarField::Internal Su
+(
+    IOobject
+    (
+        "Su",
+        runTime.timeName(),
+        mesh
+    ),
+    mesh,
+    dimensionedScalar(dimless/dimTime, Zero)
+);
+
+volScalarField::Internal divU
+(
+    mesh.moving()
+  ? fvc::div(phi + mesh.phi())
+  : fvc::div(phi)
+);
diff --git a/solver/leakyDielectricEHD/correctPhi.H b/solver/leakyDielectricEHD/correctPhi.H
new file mode 100644
index 00000000..fcb50205
--- /dev/null
+++ b/solver/leakyDielectricEHD/correctPhi.H
@@ -0,0 +1,11 @@
+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+    geometricZeroField(),
+    pimple
+);
+
+#include "continuityErrs.H"
diff --git a/solver/leakyDielectricEHD/createFields.H b/solver/leakyDielectricEHD/createFields.H
new file mode 100644
index 00000000..46d9eab5
--- /dev/null
+++ b/solver/leakyDielectricEHD/createFields.H
@@ -0,0 +1,148 @@
+#include "createRDeltaT.H"
+
+Info<< "Reading field p_rgh\n" << endl;
+volScalarField p_rgh
+(
+    IOobject
+    (
+        "p_rgh",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field U\n" << endl;
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "createPhi.H"
+
+
+Info<< "Reading transportProperties\n" << endl;
+immiscibleIncompressibleTwoPhaseMixture mixture(U, phi);
+
+volScalarField& alpha1(mixture.alpha1());
+volScalarField& alpha2(mixture.alpha2());
+
+const dimensionedScalar& rho1 = mixture.rho1();
+const dimensionedScalar& rho2 = mixture.rho2();
+
+// Need to store rho for ddt(rho, U)
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT
+    ),
+    alpha1*rho1 + alpha2*rho2
+);
+rho.oldTime();
+
+// Mass flux
+surfaceScalarField rhoPhi
+(
+    IOobject
+    (
+        "rhoPhi",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    fvc::interpolate(rho)*phi
+);
+
+// Construct incompressible turbulence model
+autoPtr<incompressible::turbulenceModel> turbulence
+(
+    incompressible::turbulenceModel::New(U, phi, mixture)
+);
+
+IOdictionary transportProperties
+(
+    IOobject
+    (
+        "transportProperties",
+        U.time().constant(),
+        U.db(),
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+surfaceForces surfForces(alpha1,phi,U,transportProperties);
+
+
+//#include "readGravitationalAcceleration.H"
+//#include "readhRef.H"
+//#include "gh.H"
+const volScalarField& gh = surfForces.acc();
+
+
+
+volScalarField p
+(
+    IOobject
+    (
+        "p",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    p_rgh + rho*gh
+);
+
+label pRefCell = 0;
+scalar pRefValue = 0.0;
+setRefCell
+(
+    p,
+    p_rgh,
+    pimple.dict(),
+    pRefCell,
+    pRefValue
+);
+
+if (p_rgh.needReference())
+{
+    p += dimensionedScalar
+    (
+        "p",
+        p.dimensions(),
+        pRefValue - getRefCellValue(p, pRefCell)
+    );
+    p_rgh = p - rho*gh;
+}
+
+mesh.setFluxRequired(p_rgh.name());
+mesh.setFluxRequired(alpha1.name());
+
+#include "createMRF.H"
+#include "createFvOptions.H"
+
+autoPtr<advectionSchemes> advector
+(
+    advectionSchemes::New(alpha1,phi,U)
+);
+
+#include "createEHDFields.H"
+
+
+
diff --git a/solver/leakyDielectricEHD/createIsoAdvection.H b/solver/leakyDielectricEHD/createIsoAdvection.H
new file mode 100644
index 00000000..3360cdae
--- /dev/null
+++ b/solver/leakyDielectricEHD/createIsoAdvection.H
@@ -0,0 +1,2 @@
+// Defining isoAdvection
+isoAdvection advector(alpha1, phi, U);
diff --git a/solver/leakyDielectricEHD/initCorrectPhi.H b/solver/leakyDielectricEHD/initCorrectPhi.H
new file mode 100644
index 00000000..03eb7fd7
--- /dev/null
+++ b/solver/leakyDielectricEHD/initCorrectPhi.H
@@ -0,0 +1,34 @@
+tmp<volScalarField> rAU;
+
+if (correctPhi)
+{
+    rAU = new volScalarField
+    (
+        IOobject
+        (
+            "rAU",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rAU", dimTime/dimDensity, 1)
+    );
+
+    #include "correctPhi.H"
+}
+else
+{
+    CorrectPhi
+    (
+        U,
+        phi,
+        p_rgh,
+        dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+        geometricZeroField(),
+        pimple
+    );
+
+    #include "continuityErrs.H"
+}
diff --git a/solver/leakyDielectricEHD/leakyDielectricEHD.C b/solver/leakyDielectricEHD/leakyDielectricEHD.C
new file mode 100644
index 00000000..fe816fcc
--- /dev/null
+++ b/solver/leakyDielectricEHD/leakyDielectricEHD.C
@@ -0,0 +1,166 @@
+#include "advectionSchemes.H"
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "isoAdvection.H"
+#include "EulerDdtScheme.H"
+#include "localEulerDdtScheme.H"
+#include "CrankNicolsonDdtScheme.H"
+#include "subCycle.H"
+#include "immiscibleIncompressibleTwoPhaseMixture.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "CorrectPhi.H"
+#include "fvcSmooth.H"
+#include "surfaceForces.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Solver for two incompressible, isothermal immiscible fluids"
+        " using isoAdvector phase-fraction based interface capturing.\n"
+        "With optional mesh motion and mesh topology changes including"
+        " adaptive re-meshing.\n"
+        "The solver is derived from interFoam"
+    );
+
+    Foam::argList::addBoolOption
+    (
+        "overwrite",
+        "Update and overwrite the existing mesh useful for adaptive mesh refinement"
+    );
+
+    #include "postProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "initContinuityErrs.H"
+    #include "createDyMControls.H"
+    #include "createFields.H"
+    #include "initCorrectPhi.H"
+    #include "createUfIfPresent.H"
+
+    const bool overwrite = args.found("overwrite");
+
+    turbulence->validate();
+
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readDyMControls.H"
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+        #include "setDeltaT.H"
+
+        ++runTime;
+
+        if(overwrite)
+        {
+            runTime.setTime(runTime.value() - runTime.deltaTValue(), 1);
+            runTime.writeAndEnd();
+        }
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            if (pimple.firstIter() || moveMeshOuterCorrectors)
+            {
+                advector->surf().reconstruct();
+
+                mesh.update();
+
+                if (mesh.changing())
+                {
+                    // gets recompute by surfaces forces
+                    // gh = (g & mesh.C()) - ghRef;
+                    // ghf = (g & mesh.Cf()) - ghRef;
+                    advector->surf().mapAlphaField();
+                    alpha2 = 1.0 - alpha1;
+                    alpha2.correctBoundaryConditions();
+                    rho == alpha1*rho1 + alpha2*rho2;
+                    rho.correctBoundaryConditions();
+                    rho.oldTime() = rho;
+                    alpha2.oldTime() = alpha2;
+
+                    MRF.update();
+
+                    if (correctPhi)
+                    {
+                        // Calculate absolute flux
+                        // from the mapped surface velocity
+                        phi = mesh.Sf() & Uf();
+
+                        #include "correctPhi.H"
+
+                        // Make the flux relative to the mesh motion
+                        fvc::makeRelative(phi, U);
+
+                        mixture.correct();
+                    }
+
+                    if (checkMeshCourantNo)
+                    {
+                        #include "meshCourantNo.H"
+                    }
+                }
+            }
+
+            if(overwrite)
+            {
+                continue;
+            }
+
+            #include "alphaControls.H"
+            #include "alphaEqnSubCycle.H"
+
+            mixture.correct();
+
+            surfForces.correct();
+
+            if (pimple.frozenFlow())
+            {
+                continue;
+            }
+            for (int corrEHD=0; corrEHD<ehdCorr; corrEHD++)
+            {
+                #include "ehdEqn.H"
+            }
+
+            #include "UEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //
\ No newline at end of file
diff --git a/solver/leakyDielectricEHD/pEqn.H b/solver/leakyDielectricEHD/pEqn.H
new file mode 100644
index 00000000..59db4825
--- /dev/null
+++ b/solver/leakyDielectricEHD/pEqn.H
@@ -0,0 +1,91 @@
+{
+    if (correctPhi)
+    {
+        rAU.ref() = 1.0/UEqn.A();
+    }
+    else
+    {
+        rAU = 1.0/UEqn.A();
+    }
+
+    surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
+    volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::flux(HbyA)
+      + MRF.zeroFilter(fvc::interpolate(rho*rAU())*fvc::ddtCorr(U, phi, Uf))
+    );
+    MRF.makeRelative(phiHbyA);
+
+    if (p_rgh.needReference())
+    {
+        fvc::makeRelative(phiHbyA, U);
+        adjustPhi(phiHbyA, U, p_rgh);
+        fvc::makeAbsolute(phiHbyA, U);
+    }
+
+    surfaceScalarField phig
+    (
+        (
+          //  mixture.surfaceTensionForce()
+          surfForces.surfaceTensionForce()
+          //- ghf*fvc::snGrad(rho)
+          + surfForces.accelerationForce()
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            p_rgh.relax();
+
+            U = HbyA + rAU()*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
+            fvOptions.correct(U);
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    // Correct Uf if the mesh is moving
+    fvc::correctUf(Uf, U, phi);
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+
+    if (!correctPhi)
+    {
+        rAU.clear();
+    }
+}
diff --git a/solver/leakyDielectricEHD/rhofs.H b/solver/leakyDielectricEHD/rhofs.H
new file mode 100644
index 00000000..5949bf1e
--- /dev/null
+++ b/solver/leakyDielectricEHD/rhofs.H
@@ -0,0 +1,2 @@
+const dimensionedScalar& rho1f(rho1);
+const dimensionedScalar& rho2f(rho2);
diff --git a/solver/leakyDielectricEHD/setDeltaT.H b/solver/leakyDielectricEHD/setDeltaT.H
new file mode 100644
index 00000000..b815a296
--- /dev/null
+++ b/solver/leakyDielectricEHD/setDeltaT.H
@@ -0,0 +1,57 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    setDeltaT
+
+Description
+    Reset the timestep to maintain a constant maximum courant Number.
+    Reduction of time-step is immediate, but increase is damped to avoid
+    unstable oscillations.
+
+\*---------------------------------------------------------------------------*/
+
+if (adjustTimeStep)
+{
+    scalar maxDeltaTFact =  
+        min
+        (
+            maxCo/(CoNum + SMALL),
+            min(maxAlphaCo/(alphaCoNum + SMALL), maxCapillaryNum/(capillaryNum + SMALL))
+        );
+
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+
+    runTime.setDeltaT
+    (
+        min
+        (
+            deltaTFact*runTime.deltaTValue(),
+            maxDeltaT
+        )
+    );
+
+    Info<< "deltaT = " <<  runTime.deltaTValue() << endl;
+}
+
+// ************************************************************************* //