diff --git a/src/functionObjects/solvers/Make/files b/src/functionObjects/solvers/Make/files
index ce5306b1c8..b37d8bafb0 100644
--- a/src/functionObjects/solvers/Make/files
+++ b/src/functionObjects/solvers/Make/files
@@ -1,4 +1,5 @@
 scalarTransport/scalarTransport.C
 energyTransport/energyTransport.C
+electricPotential/electricPotential.C
 
 LIB = $(FOAM_LIBBIN)/libsolverFunctionObjects
diff --git a/src/functionObjects/solvers/electricPotential/electricPotential.C b/src/functionObjects/solvers/electricPotential/electricPotential.C
new file mode 100644
index 0000000000..92d6b8913d
--- /dev/null
+++ b/src/functionObjects/solvers/electricPotential/electricPotential.C
@@ -0,0 +1,423 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "electricPotential.H"
+#include "fvc.H"
+#include "fvm.H"
+#include "calculatedFvPatchField.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+    defineTypeNameAndDebug(electricPotential, 0);
+    addToRunTimeSelectionTable(functionObject, electricPotential, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+Foam::volScalarField&
+Foam::functionObjects::electricPotential::operandField()
+{
+    if (!foundObject<volScalarField>(fieldName_))
+    {
+        auto tfldPtr = tmp<volScalarField>::New
+        (
+            IOobject
+            (
+                fieldName_,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_
+        );
+        store(fieldName_, tfldPtr);
+    }
+
+    return lookupObjectRef<volScalarField>(fieldName_);
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::functionObjects::electricPotential::sigma() const
+{
+    const IOobject sigmaIO
+    (
+        IOobject::scopedName(typeName, "sigma"),
+        mesh_.time().timeName(),
+        mesh_.time(),
+        IOobject::NO_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if (phases_.size())
+    {
+        tmp<volScalarField> tsigma = phases_[0]*sigmas_[0];
+
+        for (label i = 1; i < phases_.size(); ++i)
+        {
+            tsigma.ref() += phases_[i]*sigmas_[i];
+        }
+
+        return tmp<volScalarField>::New
+        (
+            sigmaIO,
+            tsigma,
+            calculatedFvPatchField<scalar>::typeName
+        );
+    }
+
+    return tmp<volScalarField>::New
+    (
+        sigmaIO,
+        mesh_,
+        sigma_,
+        calculatedFvPatchField<scalar>::typeName
+    );
+}
+
+
+Foam::tmp<Foam::volScalarField>
+Foam::functionObjects::electricPotential::epsilonm() const
+{
+    // Vacuum permittivity (aka the electric constant) [A^2 s^4/(kg m^3)]
+    const dimensionedScalar epsilon0
+    (
+        sqr(dimCurrent)*pow4(dimTime)/(dimMass*pow3(dimLength)),
+        8.8541878128e-12    // CODATA value
+    );
+
+    const IOobject epsilonrIO
+    (
+        IOobject::scopedName(typeName, "epsilonr"),
+        mesh_.time().timeName(),
+        mesh_.time(),
+        IOobject::NO_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if (phases_.size())
+    {
+        tmp<volScalarField> tepsilonr = phases_[0]*epsilonrs_[0];
+
+        for (label i = 1; i < phases_.size(); ++i)
+        {
+            tepsilonr.ref() += phases_[i]*epsilonrs_[i];
+        }
+
+        return tmp<volScalarField>::New
+        (
+            epsilonrIO,
+            epsilon0*tepsilonr,
+            calculatedFvPatchField<scalar>::typeName
+        );
+    }
+
+    return tmp<volScalarField>::New
+    (
+        epsilonrIO,
+        mesh_,
+        epsilon0*epsilonr_,
+        calculatedFvPatchField<scalar>::typeName
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::functionObjects::electricPotential::electricPotential
+(
+    const word& name,
+    const Time& runTime,
+    const dictionary& dict
+)
+:
+    fvMeshFunctionObject(name, runTime, dict),
+    phasesDict_(dict.subOrEmptyDict("phases")),
+    phaseNames_(),
+    phases_(),
+    sigmas_(),
+    sigma_
+    (
+        dimensionedScalar
+        (
+            sqr(dimCurrent)*pow3(dimTime)/(dimMass*pow3(dimLength)),
+            dict.getCheckOrDefault<scalar>
+            (
+                "sigma",
+                scalar(1),
+                scalarMinMax::ge(SMALL)
+            )
+        )
+    ),
+    epsilonrs_(),
+    epsilonr_
+    (
+        dimensionedScalar
+        (
+            dimless,
+            dict.getCheckOrDefault<scalar>
+            (
+                "epsilonr",
+                scalar(1),
+                scalarMinMax::ge(SMALL)
+            )
+        )
+    ),
+    fieldName_
+    (
+        dict.getOrDefault<word>
+        (
+            "field",
+            IOobject::scopedName(typeName, "V")
+        )
+    ),
+    nCorr_(1),
+    writeDerivedFields_(false)
+{
+    read(dict);
+
+    // Force creation of transported field so any BCs using it can
+    // look it up
+    volScalarField& eV = operandField();
+    eV.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+bool Foam::functionObjects::electricPotential::read(const dictionary& dict)
+{
+    if (fvMeshFunctionObject::read(dict))
+    {
+        Log << type() << " read: " << name() << endl;
+
+        dict.readIfPresent("sigma", sigma_);
+        dict.readIfPresent("epsilonr", epsilonr_);
+        dict.readIfPresent("nCorr", nCorr_);
+        dict.readIfPresent("writeDerivedFields", writeDerivedFields_);
+
+        // If flow is multiphase
+        if (!phasesDict_.empty())
+        {
+            phaseNames_.setSize(phasesDict_.size());
+            phases_.setSize(phasesDict_.size());
+            sigmas_.setSize(phasesDict_.size());
+
+            if (writeDerivedFields_)
+            {
+                epsilonrs_.setSize(phasesDict_.size());
+            }
+
+            label phasei = 0;
+            forAllConstIters(phasesDict_, iter)
+            {
+                const word& key = iter().keyword();
+
+                if (!phasesDict_.isDict(key))
+                {
+                    FatalErrorInFunction
+                        << "Found non-dictionary entry " << iter()
+                        << " in top-level dictionary " << phasesDict_
+                        << exit(FatalError);
+                }
+
+                const dictionary& subDict = phasesDict_.subDict(key);
+
+                phaseNames_[phasei] = key;
+
+                sigmas_.set
+                (
+                    phasei,
+                    new dimensionedScalar
+                    (
+                        sqr(dimCurrent)*pow3(dimTime)/(dimMass*pow3(dimLength)),
+                        subDict.getCheck<scalar>
+                        (
+                            "sigma",
+                            scalarMinMax::ge(SMALL)
+                        )
+                    )
+                );
+
+                if (writeDerivedFields_)
+                {
+                    epsilonrs_.set
+                    (
+                        phasei,
+                        new dimensionedScalar
+                        (
+                            dimless,
+                            subDict.getCheck<scalar>
+                            (
+                                "epsilonr",
+                                scalarMinMax::ge(SMALL)
+                            )
+                        )
+                    );
+                }
+
+                ++phasei;
+            }
+
+            forAll(phaseNames_, i)
+            {
+                phases_.set
+                (
+                    i,
+                    mesh_.getObjectPtr<volScalarField>(phaseNames_[i])
+                );
+            }
+        }
+
+        return true;
+    }
+
+    return false;
+}
+
+
+bool Foam::functionObjects::electricPotential::execute()
+{
+    Log << type() << " execute: " << name() << endl;
+
+    tmp<volScalarField> tsigma = this->sigma();
+    const volScalarField& sigma = tsigma();
+
+    volScalarField& eV = operandField();
+
+    for (label i = 1; i <= nCorr_; ++i)
+    {
+        fvScalarMatrix eVEqn
+        (
+          - fvm::laplacian(sigma, eV)
+        );
+
+        eVEqn.relax();
+
+        eVEqn.solve();
+    }
+
+    Log << endl;
+
+    return true;
+}
+
+
+bool Foam::functionObjects::electricPotential::write()
+{
+    Log << type() << " write: " << name() << nl
+        << tab << fieldName_
+        << endl;
+
+    volScalarField& eV = operandField();
+
+    if (writeDerivedFields_)
+    {
+        // Write the electric field
+        const volVectorField E
+        (
+            IOobject
+            (
+                IOobject::scopedName(typeName, "E"),
+                mesh_.time().timeName(),
+                mesh_.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            -fvc::grad(eV),
+            calculatedFvPatchField<vector>::typeName
+        );
+
+        Log << tab << E.name() << endl;
+
+        E.write();
+
+
+        // Write the current density field
+        tmp<volScalarField> tsigma = this->sigma();
+
+        auto eJ = tmp<volVectorField>::New
+        (
+            IOobject
+            (
+                IOobject::scopedName(typeName, "J"),
+                mesh_.time().timeName(),
+                mesh_.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            -tsigma*fvc::grad(eV),
+            calculatedFvPatchField<vector>::typeName
+        );
+
+        Log << tab << eJ().name() << endl;
+
+        eJ->write();
+
+
+        // Write the free-charge density field
+        tmp<volScalarField> tepsilonm = this->epsilonm();
+
+        auto erho = tmp<volScalarField>::New
+        (
+            IOobject
+            (
+                IOobject::scopedName(typeName, "rho"),
+                mesh_.time().timeName(),
+                mesh_.time(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            fvc::div(tepsilonm*E),
+            calculatedFvPatchField<scalar>::typeName
+        );
+
+        Log << tab << erho().name() << endl;
+
+        erho->write();
+    }
+
+    eV.write();
+
+    return true;
+}
+
+
+// ************************************************************************* //
diff --git a/src/functionObjects/solvers/electricPotential/electricPotential.H b/src/functionObjects/solvers/electricPotential/electricPotential.H
new file mode 100644
index 0000000000..83fdfbd7b2
--- /dev/null
+++ b/src/functionObjects/solvers/electricPotential/electricPotential.H
@@ -0,0 +1,275 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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/>.
+
+Class
+    Foam::functionObjects::electricPotential
+
+Group
+    grpSolversFunctionObjects
+
+Description
+    Computes the steady-state equation of charge conservation to obtain
+    the electric potential by strictly assuming a quasi-static electrostatic
+    field for single-phase and multiphase applications.
+
+    The steady-state equation of the charge conservation:
+
+    \f[
+        \nabla \cdot \left( \sigma \nabla V \right) = 0
+    \f]
+
+    where
+    \vartable
+      V       | Electric potential                [volt = kg m^2/(A s^3)]
+      \sigma  | Isotropic conductivity of mixture [S/m = A^2 s^3/(kg m^3)]
+    \endvartable
+
+    Optionally, electric field, current density and free-charge
+    density fields can be written out by using the following equations:
+
+    \f[
+        \vec{E} = - \nabla V
+    \f]
+
+    \f[
+        \vec{J} = \sigma \vec{E} = - \sigma \nabla V
+    \f]
+
+    \f[
+        \rho_E = \nabla \cdot \left(\epsilon_m \vec{E} \right)
+               = \nabla \cdot \left(\epsilon_0 \epsilon_r \vec{E} \right)
+    \f]
+
+    where
+    \vartable
+      \vec{E}    | Electric field                       [m kg/(s^3 A)]
+      \vec{J}    | Current density                      [A/m^2]
+      \rho_E     | Volume charge density                [C/m^3 = A s/m^3]
+      \epsilon_m | Isotropic permittivity of mixture    [F/m = A^2 s^4/(kg m^3)]
+      \epsilon_0 | Isotropic vacuum permittivity        [F/m = A^2 s^4/(kg m^3)]
+      \epsilon_r | Isotropic relative permittivity of mixture [-]
+    \endvartable
+
+    For multiphase applications, \c sigma and \c epsilonr are blended
+    (to consider their interface values) by using the simple weighted
+    arithmetic mean interpolation, for example:
+
+    \f[
+        \sigma = \alpha_1 \sigma_1 + \alpha_2 \sigma_2
+               = \alpha_1 \sigma_1 + (1 - \alpha_1) \sigma_2
+    \f]
+
+Usage
+    Minimal example by using \c system/controlDict.functions:
+    \verbatim
+    electricPotential1
+    {
+        // Mandatory entries
+        type            electricPotential;
+        libs            (solverFunctionObjects);
+
+        // Conditional entries
+
+            // Option-1: single-phase
+            sigma       <scalar>;
+            epsilonr    <scalar>;
+
+            // Option-2: multiphase
+            phases
+            {
+                alpha.air
+                {
+                    sigma       <scalar>;
+                    epsilonr    <scalar>;
+                }
+                alpha.water
+                {
+                    sigma       <scalar>;
+                    epsilonr    <scalar>;
+                }
+                alpha.mercury
+                {
+                    sigma       <scalar>;
+                    epsilonr    <scalar>;
+                }
+                ...
+            }
+
+        // Optional entries
+        nCorr                 <label>;
+        writeDerivedFields    <bool>;
+        fieldName             <word>;
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property  | Description                         | Type | Reqd | Deflt
+      type      | Type name: electricPotential        | word | yes  | -
+      libs      | Library name: solverFunctionObjects | word | yes  | -
+      sigma     | Isotropic electrical conductivity of phase | scalar | yes | -
+      epsilonr  | Isotropic relative permittivity of phase   | scalar | no  | -
+      nCorr     | Number of corrector iterations     | label | no   | 1
+      writeDerivedFields | Flag to write extra fields | bool | no   | false
+      fieldName | Name of operand field   | word | no | electricPotential:V
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link functionObject.H \endlink
+
+    Fields written out when the \c writeDerivedFields entry is \c true:
+    \table
+      Operand         | Type             | Location
+      Electric field  | volVectorField   | \<time\>/electricPotential:E
+      Current density | volVectorField   | \<time\>/electricPotential:J
+      Charge density  | volScalarField   | \<time\>/electricPotential:rho
+    \endtable
+
+SourceFiles
+    electricPotential.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef functionObjects_electricPotential_H
+#define functionObjects_electricPotential_H
+
+#include "fvMeshFunctionObject.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace functionObjects
+{
+
+/*---------------------------------------------------------------------------*\
+                       Class electricPotential Declaration
+\*---------------------------------------------------------------------------*/
+
+class electricPotential
+:
+    public fvMeshFunctionObject
+{
+    // Private Data
+
+        //- Dictionary of phase data
+        dictionary phasesDict_;
+
+        //- List of phase names
+        wordList phaseNames_;
+
+        //- Unallocated list of phase fields
+        UPtrList<volScalarField> phases_;
+
+        //- List of isotropic electrical conductivity of phases
+        PtrList<dimensionedScalar> sigmas_;
+
+        //- Isotropic electrical conductivity of a single phase
+        dimensionedScalar sigma_;
+
+        //- List of isotropic relative permittivity of phases
+        PtrList<dimensionedScalar> epsilonrs_;
+
+        //- Isotropic relative permittivity of a single phase
+        dimensionedScalar epsilonr_;
+
+        //- Name of the operand field
+        word fieldName_;
+
+        //- Number of corrector iterations
+        label nCorr_;
+
+        //- Flag to write derived fields of
+        //- electric field, current density and free-charge density
+        bool writeDerivedFields_;
+
+
+    // Private Member Functions
+
+        //- Return reference to the registered operand field
+        volScalarField& operandField();
+
+        //- Return the isotropic electrical conductivity field of the mixture
+        tmp<volScalarField> sigma() const;
+
+        //- Return the isotropic permittivity field of the mixture
+        tmp<volScalarField> epsilonm() const;
+
+
+        //- No copy construct
+        electricPotential(const electricPotential&) = delete;
+
+        //- No copy assignment
+        void operator=(const electricPotential&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("electricPotential");
+
+
+    // Constructors
+
+        //- Construct from Time and dictionary
+        electricPotential
+        (
+            const word& name,
+            const Time& runTime,
+            const dictionary& dict
+        );
+
+
+    //- Destructor
+    virtual ~electricPotential() = default;
+
+
+    // Member Functions
+
+        //- Read the function object data
+        virtual bool read(const dictionary& dict);
+
+        //- Calculate the function object
+        virtual bool execute();
+
+        //- Write the function object output
+        virtual bool write();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace functionObjects
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //