diff --git a/src/functionObjects/solvers/electricPotential/electricPotential.C b/src/functionObjects/solvers/electricPotential/electricPotential.C
index c88bb2686c201e217a5c136ae72f708a8c815886..c25569518b63015e62924848f9a1cacf3354a234 100644
--- a/src/functionObjects/solvers/electricPotential/electricPotential.C
+++ b/src/functionObjects/solvers/electricPotential/electricPotential.C
@@ -210,8 +210,10 @@ Foam::functionObjects::electricPotential::electricPotential
             IOobject::scopedName(typeName, "V")
         )
     ),
+    erhoName_(),
     nCorr_(1),
-    writeDerivedFields_(false)
+    writeDerivedFields_(false),
+    spaceChargeDensity_(false)
 {
     read(dict);
 
@@ -246,6 +248,11 @@ bool Foam::functionObjects::electricPotential::read(const dictionary& dict)
     dict.readIfPresent("nCorr", nCorr_);
     dict.readIfPresent("writeDerivedFields", writeDerivedFields_);
 
+    if (dict.readIfPresent("spaceChargeDensity", spaceChargeDensity_))
+    {
+        dict.readEntry("erho", erhoName_);
+    }
+
     // If flow is multiphase
     if (!phasesDict_.empty())
     {
@@ -329,6 +336,9 @@ bool Foam::functionObjects::electricPotential::execute()
     tmp<volScalarField> tsigma = this->sigma();
     const volScalarField& sigma = tsigma();
 
+    tmp<volScalarField> tepsilon = this->epsilonm();
+    const auto& epsilon = tepsilon();
+
     volScalarField& eV = getOrReadField(fieldName_);
 
     for (label i = 1; i <= nCorr_; ++i)
@@ -340,6 +350,13 @@ bool Foam::functionObjects::electricPotential::execute()
             fvOptions_(eV)
         );
 
+        if (spaceChargeDensity_)
+        {
+            volScalarField& erho = getOrReadField(erhoName_);
+
+            eVEqn += erho/epsilon;
+        }
+
         eVEqn.relax();
 
         fvOptions_.constrain(eVEqn);
diff --git a/src/functionObjects/solvers/electricPotential/electricPotential.H b/src/functionObjects/solvers/electricPotential/electricPotential.H
index 4f015532deb1ce5e55ef27dbfbb12f4eb0627a59..e89066843226787d22fa54f92d8df1d7358a7b7d 100644
--- a/src/functionObjects/solvers/electricPotential/electricPotential.H
+++ b/src/functionObjects/solvers/electricPotential/electricPotential.H
@@ -37,7 +37,7 @@ Description
     The steady-state equation of the charge conservation:
 
     \f[
-        \nabla \cdot \left( \sigma \nabla V \right) = 0
+        \nabla \cdot \left(\sigma \nabla V \right) = 0
     \f]
 
     where
@@ -121,6 +121,8 @@ Usage
         nCorr                 <label>;
         writeDerivedFields    <bool>;
         fieldName             <word>;
+        spaceChargeDensity    <bool>;
+        erhoName              <word>;
 
         // Inherited entries
         ...
@@ -137,6 +139,9 @@ Usage
       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
+      spaceChargeDensity | Flag to include space charge density | bool | no <!--
+                -->                                                 | false
+      erhoName  | Name of space charge density field  | word | no   | -
     \endtable
 
     The inherited entries are elaborated in:
@@ -206,6 +211,9 @@ class electricPotential
         //- Name of the operand field
         word fieldName_;
 
+        //- Name of the space charge density field
+        word erhoName_;
+
         //- Number of corrector iterations
         label nCorr_;
 
@@ -213,6 +221,9 @@ class electricPotential
         //- electric field, current density and free-charge density
         bool writeDerivedFields_;
 
+        //- Flag to add the space charge density to the rhs of the potential eq
+        bool spaceChargeDensity_;
+
 
     // Private Member Functions