WIP: electricPotential: add space charge density hook

src/functionObjects/solvers/electricPotential/electricPotential.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,16 +45,20 @@ namespace functionObjects
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-Foam::volScalarField&
-Foam::functionObjects::electricPotential::operandField()
+Foam::volScalarField& Foam::functionObjects::electricPotential::getOrReadField
+(
+    const word& fieldName
+) const
 {
-    if (!foundObject<volScalarField>(fieldName_))
+    auto* ptr = mesh_.getObjectPtr<volScalarField>(fieldName);
+
+    if (!ptr)
     {
-        auto tfldPtr = tmp<volScalarField>::New
+        ptr = new volScalarField
         (
             IOobject
             (
-                fieldName_,
+                fieldName,
                 mesh_.time().timeName(),
                 mesh_,
                 IOobject::MUST_READ,
@@ -62,10 +66,10 @@ Foam::functionObjects::electricPotential::operandField()
             ),
             mesh_
         );
-        store(fieldName_, tfldPtr);
+        mesh_.objectRegistry::store(ptr);
     }
 
-    return lookupObjectRef<volScalarField>(fieldName_);
+    return *ptr;
 }
 
 
@@ -79,7 +83,7 @@ Foam::functionObjects::electricPotential::sigma() const
         mesh_.time(),
         IOobject::NO_READ,
         IOobject::NO_WRITE,
-        false
+        IOobject::NO_REGISTER
     );
 
     if (phases_.size())
@@ -126,7 +130,7 @@ Foam::functionObjects::electricPotential::epsilonm() const
         mesh_.time(),
         IOobject::NO_READ,
         IOobject::NO_WRITE,
-        false
+        IOobject::NO_REGISTER
     );
 
     if (phases_.size())
@@ -166,6 +170,7 @@ Foam::functionObjects::electricPotential::electricPotential
 )
 :
     fvMeshFunctionObject(name, runTime, dict),
+    fvOptions_(mesh_),
     phasesDict_(dict.subOrEmptyDict("phases")),
     phaseNames_(),
     phases_(),
@@ -205,14 +210,16 @@ Foam::functionObjects::electricPotential::electricPotential
             IOobject::scopedName(typeName, "V")
         )
     ),
+    erhoName_(),
     nCorr_(1),
-    writeDerivedFields_(false)
+    writeDerivedFields_(false),
+    spaceChargeDensity_(false)
 {
     read(dict);
 
     // Force creation of transported field so any BCs using it can
     // look it up
-    volScalarField& eV = operandField();
+    volScalarField& eV = getOrReadField(fieldName_);
     eV.correctBoundaryConditions();
 }
 
@@ -221,88 +228,101 @@ Foam::functionObjects::electricPotential::electricPotential
 
 bool Foam::functionObjects::electricPotential::read(const dictionary& dict)
 {
-    if (fvMeshFunctionObject::read(dict))
+    if (!fvMeshFunctionObject::read(dict))
     {
-        Log << type() << " read: " << name() << endl;
+        return false;
+    }
 
-        dict.readIfPresent("sigma", sigma_);
-        dict.readIfPresent("epsilonr", epsilonr_);
-        dict.readIfPresent("nCorr", nCorr_);
-        dict.readIfPresent("writeDerivedFields", writeDerivedFields_);
+    Log << type() << " read: " << name() << endl;
 
-        // If flow is multiphase
-        if (!phasesDict_.empty())
+    {
+        const dictionary* dictptr = dict.findDict("fvOptions");
+        if (dictptr)
         {
-            phaseNames_.setSize(phasesDict_.size());
-            phases_.setSize(phasesDict_.size());
-            sigmas_.setSize(phasesDict_.size());
+            fvOptions_.reset(*dictptr);
+        }
+    }
+
+    dict.readIfPresent("sigma", sigma_);
+    dict.readIfPresent("epsilonr", epsilonr_);
+    dict.readIfPresent("nCorr", nCorr_);
+    dict.readIfPresent("writeDerivedFields", writeDerivedFields_);
+
+    if (dict.readIfPresent("spaceChargeDensity", spaceChargeDensity_))
+    {
+        dict.readEntry("erho", erhoName_);
+    }
+
+    // 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;
+        for (const entry& dEntry : phasesDict_)
+        {
+            const word& key = dEntry.keyword();
+
+            if (!dEntry.isDict())
+            {
+                FatalIOErrorInFunction(phasesDict_)
+                    << "Entry " << key << " is not a dictionary" << nl
+                    << exit(FatalIOError);
+            }
+
+            const dictionary& subDict = dEntry.dict();
+
+            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_.setSize(phasesDict_.size());
-            }
-
-            label phasei = 0;
-            for (const entry& dEntry : phasesDict_)
-            {
-                const word& key = dEntry.keyword();
-
-                if (!dEntry.isDict())
-                {
-                    FatalIOErrorInFunction(phasesDict_)
-                        << "Entry " << key << " is not a dictionary" << nl
-                        << exit(FatalIOError);
-                }
-
-                const dictionary& subDict = dEntry.dict();
-
-                phaseNames_[phasei] = key;
-
-                sigmas_.set
+                epsilonrs_.set
                 (
                     phasei,
                     new dimensionedScalar
                     (
-                        sqr(dimCurrent)*pow3(dimTime)/(dimMass*pow3(dimLength)),
+                        dimless,
                         subDict.getCheck<scalar>
                         (
-                            "sigma",
+                            "epsilonr",
                             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])
-                );
-            }
+            ++phasei;
         }
 
-        return true;
+        forAll(phaseNames_, i)
+        {
+            phases_.set
+            (
+                i,
+                mesh_.getObjectPtr<volScalarField>(phaseNames_[i])
+            );
+        }
     }
 
     return false;
@@ -316,17 +336,31 @@ bool Foam::functionObjects::electricPotential::execute()
     tmp<volScalarField> tsigma = this->sigma();
     const volScalarField& sigma = tsigma();
 
-    volScalarField& eV = operandField();
+    tmp<volScalarField> tepsilon = this->epsilonm();
+    const auto& epsilon = tepsilon();
+
+    volScalarField& eV = getOrReadField(fieldName_);
 
     for (label i = 1; i <= nCorr_; ++i)
     {
         fvScalarMatrix eVEqn
         (
           - fvm::laplacian(sigma, eV)
+         ==
+            fvOptions_(eV)
         );
 
+        if (spaceChargeDensity_)
+        {
+            volScalarField& erho = getOrReadField(erhoName_);
+
+            eVEqn += erho/epsilon;
+        }
+
         eVEqn.relax();
 
+        fvOptions_.constrain(eVEqn);
+
         eVEqn.solve();
     }
 
@@ -342,7 +376,7 @@ bool Foam::functionObjects::electricPotential::write()
         << tab << fieldName_
         << endl;
 
-    volScalarField& eV = operandField();
+    volScalarField& eV = getOrReadField(fieldName_);
 
     if (writeDerivedFields_)
     {
@@ -356,7 +390,7 @@ bool Foam::functionObjects::electricPotential::write()
                 mesh_.time(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
-                false
+                IOobject::NO_REGISTER
             ),
             -fvc::grad(eV),
             calculatedFvPatchField<vector>::typeName
@@ -379,7 +413,7 @@ bool Foam::functionObjects::electricPotential::write()
                 mesh_.time(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
-                false
+                IOobject::NO_REGISTER
             ),
             -tsigma*fvc::grad(eV),
             calculatedFvPatchField<vector>::typeName
@@ -402,7 +436,7 @@ bool Foam::functionObjects::electricPotential::write()
                 mesh_.time(),
                 IOobject::NO_READ,
                 IOobject::NO_WRITE,
-                false
+                IOobject::NO_REGISTER
             ),
             fvc::div(tepsilonm*E),
             calculatedFvPatchField<scalar>::typeName

src/functionObjects/solvers/electricPotential/electricPotential.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2021-2023 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -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:
@@ -160,6 +165,7 @@ SourceFiles
 
 #include "fvMeshFunctionObject.H"
 #include "volFields.H"
+#include "fvOptionList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -178,6 +184,9 @@ class electricPotential
 {
     // Private Data
 
+        //- Run-time selectable finite volume options, e.g. sources, constraints
+        fv::optionList fvOptions_;
+
         //- Dictionary of phase data
         dictionary phasesDict_;
 
@@ -202,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_;
 
@@ -209,11 +221,16 @@ 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
 
-        //- Return reference to the registered operand field
-        volScalarField& operandField();
+        //- Return requested field from the object registry
+        //- or read+register the field to the object registry
+        volScalarField& getOrReadField(const word& fieldName) const;
+
 
         //- Return the isotropic electrical conductivity field of the mixture
         tmp<volScalarField> sigma() const;

src/lagrangian/intermediate/Make/files

@@ -9,6 +9,7 @@ $(TEMPLATECLOUDS)/KinematicCloud/KinematicCloudName.C
 $(TEMPLATECLOUDS)/ThermoCloud/ThermoCloudName.C
 $(TEMPLATECLOUDS)/ReactingCloud/ReactingCloudName.C
 $(TEMPLATECLOUDS)/ReactingMultiphaseCloud/ReactingMultiphaseCloudName.C
+$(TEMPLATECLOUDS)/ChargeCloud/ChargeCloudName.C
 $(TEMPLATECLOUDS)/ReactingHeterogeneousCloud/ReactingHeterogeneousCloudName.C
 $(TEMPLATECLOUDS)/CollidingCloud/CollidingCloudName.C
 $(TEMPLATECLOUDS)/MPPICCloud/MPPICCloudName.C
@@ -20,6 +21,7 @@ $(DERIVEDCLOUDS)/kinematicCloud/kinematicCloud.C
 $(DERIVEDCLOUDS)/thermoCloud/thermoCloud.C
 $(DERIVEDCLOUDS)/reactingCloud/reactingCloud.C
 $(DERIVEDCLOUDS)/reactingMultiphaseCloud/reactingMultiphaseCloud.C
+$(DERIVEDCLOUDS)/chargeCloud/chargeCloud.C
 $(DERIVEDCLOUDS)/heterogeneousReactingCloud/heterogeneousReactingCloud.C
 $(DERIVEDCLOUDS)/kinematicCollidingCloud/kinematicCollidingCloud.C
 $(DERIVEDCLOUDS)/kinematicMPPICCloud/kinematicMPPICCloud.C
@@ -56,12 +58,16 @@ $(REACTINGMPPARCEL)/defineReactingMultiphaseParcel.C
 $(REACTINGMPPARCEL)/makeReactingMultiphaseParcelSubmodels.C
 
 
+/* charge parcel sub-models */
+CHARGEPARCEL=$(DERIVEDPARCELS)/chargeParcel
+$(CHARGEPARCEL)/defineChargeParcel.C
+$(CHARGEPARCEL)/makeChargeParcelSubmodels.C
+
 /* heterogeneous reacting parcel sub-models */
 REACTINGHETERMPPARCEL=$(DERIVEDPARCELS)/heterogeneousReactingParcel
 $(REACTINGHETERMPPARCEL)/defineHeterogeneousReactingParcel.C
 $(REACTINGHETERMPPARCEL)/makeHeterogeneousReactingParcelSubmodels.C
 
-
 /* kinematic MPPIC parcel sub-models */
 KINEMATICMPPICPARCEL=$(DERIVEDPARCELS)/kinematicMPPICParcel
 $(KINEMATICMPPICPARCEL)/defineKinematicMPPICParcel.C
@@ -102,6 +108,11 @@ $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphase
 $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphaseParcelInjectionDataIO.C
 $(REACTINGMPINJECTION)/ReactingMultiphaseLookupTableInjection/reactingMultiphaseParcelInjectionDataIOList.C
 
+CHARGEINJECTION=submodels/Charge/InjectionModel
+$(CHARGEINJECTION)/ChargeLookupTableInjection/chargeParcelInjectionData.C
+$(CHARGEINJECTION)/ChargeLookupTableInjection/chargeParcelInjectionDataIO.C
+$(CHARGEINJECTION)/ChargeLookupTableInjection/chargeParcelInjectionDataIOList.C
+
 MPPICPARTICLESTRESS=submodels/MPPIC/ParticleStressModels
 $(MPPICPARTICLESTRESS)/ParticleStressModel/ParticleStressModel.C
 $(MPPICPARTICLESTRESS)/HarrisCrighton/HarrisCrighton.C

src/lagrangian/intermediate/clouds/Templates/ChargeCloud/ChargeCloud.C

@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "ChargeCloud.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::setModels()
+{
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::cloudReset
+(
+    ChargeCloud<CloudType>& c
+)
+{
+    CloudType::cloudReset(c);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ChargeCloud<CloudType>::ChargeCloud
+(
+    const word& cloudName,
+    const dimensionedVector& g,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const SLGThermo& thermo,
+    bool readFields
+)
+:
+    CloudType(cloudName, g, rho, U, thermo, false),
+    cloudCopyPtr_(nullptr),
+    constProps_(this->particleProperties())
+{
+    if (this->solution().active())
+    {
+        setModels();
+
+        if (readFields)
+        {
+            parcelType::readFields(*this, this->composition());
+            this->deleteLostParticles();
+        }
+    }
+
+    if (this->solution().resetSourcesOnStartup())
+    {
+        resetSourceTerms();
+    }
+}
+
+
+template<class CloudType>
+Foam::ChargeCloud<CloudType>::ChargeCloud
+(
+    ChargeCloud<CloudType>& c,
+    const word& name
+)
+:
+    CloudType(c, name),
+    cloudCopyPtr_(nullptr),
+    constProps_(c.constProps_)
+{}
+
+
+template<class CloudType>
+Foam::ChargeCloud<CloudType>::ChargeCloud
+(
+    const fvMesh& mesh,
+    const word& name,
+    const ChargeCloud<CloudType>& c
+)
+:
+    CloudType(mesh, name, c),
+    cloudCopyPtr_(nullptr),
+    constProps_()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::setParcelThermoProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt
+)
+{
+    CloudType::setParcelThermoProperties(parcel, lagrangianDt);
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::checkParcelProperties
+(
+    parcelType& parcel,
+    const scalar lagrangianDt,
+    const bool fullyDescribed
+)
+{
+    CloudType::checkParcelProperties(parcel, lagrangianDt, fullyDescribed);
+
+    if (fullyDescribed)
+    {
+
+    }
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::storeState()
+{
+    cloudCopyPtr_.reset
+    (
+        static_cast<ChargeCloud<CloudType>*>
+        (
+            clone(this->name() + "Copy").ptr()
+        )
+    );
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::restoreState()
+{
+    cloudReset(cloudCopyPtr_());
+    cloudCopyPtr_.clear();
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::resetSourceTerms()
+{
+    CloudType::resetSourceTerms();
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::evolve()
+{
+    if (this->solution().canEvolve())
+    {
+        typename parcelType::trackingData td(*this);
+
+        this->solve(*this, td);
+    }
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::autoMap
+(
+    const mapPolyMesh& mapper
+)
+{
+    Cloud<parcelType>::autoMap(mapper);
+
+    this->updateMesh();
+}
+
+
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::info()
+{
+    CloudType::info();
+}
+
+/*
+template<class CloudType>
+void Foam::ChargeCloud<CloudType>::writeFields() const
+{
+    if (this->compositionModel_)
+    {
+        CloudType::particleType::writeFields(*this, this->composition());
+    }
+}
+*/
+
+// ************************************************************************* //

src/lagrangian/intermediate/clouds/Templates/ChargeCloud/ChargeCloud.H

@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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::ChargeCloud
+
+Group
+    grpLagrangianIntermediateClouds
+
+Description
+    Templated base class for electric-charge cloud
+
+    - Adds to multiphase reacting cloud
+      - electric charge
+
+SourceFiles
+    ChargeCloudI.H
+    ChargeCloud.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_ChargeCloud_H
+#define Foam_ChargeCloud_H
+
+#include "className.H"
+#include "autoPtr.H"
+#include "Cloud.H"
+#include "volFields.H"
+#include "fvMatricesFwd.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+
+class SLGThermo;
+
+TemplateName(ChargeCloud);
+
+/*---------------------------------------------------------------------------*\
+                        Class ChargeCloud Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class ChargeCloud
+:
+    public CloudType,
+    public ChargeCloudName
+{
+public:
+
+    // Public Typedefs
+
+        //- Type of cloud this cloud was instantiated for
+        typedef CloudType cloudType;
+
+        //- Type of parcel the cloud was instantiated for
+        typedef typename CloudType::particleType parcelType;
+
+        //- Convenience typedef for this cloud type
+        typedef ChargeCloud<CloudType> chargeCloudType;
+
+
+private:
+
+    // Private Data
+
+        //- Cloud copy pointer
+        autoPtr<ChargeCloud<CloudType>> cloudCopyPtr_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        ChargeCloud(const ChargeCloud&) = delete;
+
+        //- No copy assignment
+        void operator=(const ChargeCloud&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        //- Parcel constant properties
+        typename parcelType::constantProperties constProps_;
+
+
+    // Protected Member Functions
+
+        // Initialisation
+
+            //- Set cloud sub-models
+            void setModels();
+
+
+        // Cloud Evolution Functions
+
+            //- Reset state of cloud
+            void cloudReset(ChargeCloud<CloudType>& c);
+
+
+public:
+
+    // Constructors
+
+        //- Construct given carrier gas fields
+        ChargeCloud
+        (
+            const word& cloudName,
+            const dimensionedVector& g,
+            const volScalarField& rho,
+            const volVectorField& U,
+            const SLGThermo& thermo,
+            bool readFields = true
+        );
+
+        //- Copy constructor with new name
+        ChargeCloud
+        (
+            ChargeCloud<CloudType>& c,
+            const word& name
+        );
+
+        //- Copy constructor with new name - creates bare cloud
+        ChargeCloud
+        (
+            const fvMesh& mesh,
+            const word& name,
+            const ChargeCloud<CloudType>& c
+        );
+
+        //- Construct and return clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType>> clone(const word& name)
+        {
+            return autoPtr<Cloud<parcelType>>
+            (
+                new ChargeCloud(*this, name)
+            );
+        }
+
+        //- Construct and return bare clone based on (this) with new name
+        virtual autoPtr<Cloud<parcelType>> cloneBare(const word& name) const
+        {
+            return autoPtr<Cloud<parcelType>>
+            (
+                new ChargeCloud(this->mesh(), name, *this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~ChargeCloud() = default;
+
+
+    // Member Functions
+
+        // Access
+
+        //- Return a reference to the cloud copy
+        inline const ChargeCloud& cloudCopy() const;
+
+        //- Return the constant properties
+        inline const typename parcelType::constantProperties&
+            constProps() const;
+
+        //- Return access to the constant properties
+        inline typename parcelType::constantProperties& constProps();
+
+
+        // Sub-models
+
+
+        // Cloud Evolution
+
+        //- Set parcel thermo properties
+        void setParcelThermoProperties
+        (
+            parcelType& parcel,
+            const scalar lagrangianDt
+        );
+
+        //- Check parcel properties
+        void checkParcelProperties
+        (
+            parcelType& parcel,
+            const scalar lagrangianDt,
+            const bool fullyDescribed
+        );
+
+        //- Store the current cloud state
+        void storeState();
+
+        //- Reset the current cloud to the previously stored state
+        void restoreState();
+
+        //- Reset the cloud source terms
+        void resetSourceTerms();
+
+        //- Evolve the cloud
+        void evolve();
+
+
+        // Mapping
+
+        //- Remap the cells of particles corresponding to the
+        //- mesh topology change with a default tracking data object
+        virtual void autoMap(const mapPolyMesh&);
+
+
+        // I-O
+
+        //- Print cloud information
+        void info();
+
+        //- Write the field data for the cloud
+        // virtual void writeFields() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ChargeCloudI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ChargeCloud.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/clouds/Templates/ChargeCloud/ChargeCloudI.H

@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+inline const Foam::ChargeCloud<CloudType>&
+Foam::ChargeCloud<CloudType>::cloudCopy() const
+{
+    return *cloudCopyPtr_;
+}
+
+
+template<class CloudType>
+inline const typename CloudType::particleType::constantProperties&
+Foam::ChargeCloud<CloudType>::constProps() const
+{
+    return constProps_;
+}
+
+
+template<class CloudType>
+inline typename CloudType::particleType::constantProperties&
+Foam::ChargeCloud<CloudType>::constProps()
+{
+    return constProps_;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/clouds/Templates/ChargeCloud/ChargeCloudName.C

@@ -0,0 +1,37 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "ChargeCloud.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+defineTypeNameAndDebug(ChargeCloudName, 0);
+}
+
+// ************************************************************************* //

src/lagrangian/intermediate/clouds/derived/chargeCloud/chargeCloud.C

@@ -0,0 +1,45 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "chargeCloud.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    addNamedToRunTimeSelectionTable
+    (
+        parcelCloudModel,
+        chargeCloud,
+        thermo,
+        chargeCloud
+    );
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/clouds/derived/chargeCloud/chargeCloud.H

@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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::chargeCloud
+
+Description
+    Cloud class to introduce multi-phase reacting charged parcels
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_chargeCloud_H
+#define Foam_chargeCloud_H
+
+#include "BaseCloud.H"
+#include "KinematicCloud.H"
+#include "ThermoCloud.H"
+#include "ReactingCloud.H"
+#include "ReactingMultiphaseCloud.H"
+#include "ChargeCloud.H"
+#include "ParcelCloudModel.H"
+#include "chargeParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef ParcelCloudModel
+    <
+        ChargeCloud
+        <
+            ReactingMultiphaseCloud
+            <
+                ReactingCloud
+                <
+                    ThermoCloud
+                    <
+                        KinematicCloud
+                        <
+                            BaseCloud
+                            <
+                                chargeParcel
+                            >
+                        >
+                    >
+                >
+            >
+        >
+    > chargeCloud;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/ChargeParcel/ChargeParcel.C

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "ChargeParcel.H"
+#include "mathematicalConstants.H"
+
+using namespace Foam::constant::mathematical;
+
+// * * * * * * * * * * *  Protected Member Functions * * * * * * * * * * * * //
+
+template<class ParcelType>
+template<class TrackCloudType>
+void Foam::ChargeParcel<ParcelType>::setCellValues
+(
+    TrackCloudType& cloud,
+    trackingData& td
+)
+{
+    ParcelType::setCellValues(cloud, td);
+}
+
+
+template<class ParcelType>
+template<class TrackCloudType>
+void Foam::ChargeParcel<ParcelType>::cellValueSourceCorrection
+(
+    TrackCloudType& cloud,
+    trackingData& td,
+    const scalar dt
+)
+{
+    // Re-use correction from reacting parcel
+    ParcelType::cellValueSourceCorrection(cloud, td, dt);
+}
+
+
+template<class ParcelType>
+template<class TrackCloudType>
+void Foam::ChargeParcel<ParcelType>::calc
+(
+    TrackCloudType& cloud,
+    trackingData& td,
+    const scalar dt
+)
+{
+    ParcelType::calc(cloud, td, dt);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::ChargeParcel<ParcelType>::ChargeParcel
+(
+    const ChargeParcel<ParcelType>& p
+)
+:
+    ParcelType(p),
+    eq_(p.eq_)
+{}
+
+
+template<class ParcelType>
+Foam::ChargeParcel<ParcelType>::ChargeParcel
+(
+    const ChargeParcel<ParcelType>& p,
+    const polyMesh& mesh
+)
+:
+    ParcelType(p, mesh),
+    eq_(p.eq_)
+{}
+
+
+// * * * * * * * * * * * * * * IOStream operators  * * * * * * * * * * * * * //
+
+#include "ChargeParcelIO.C"
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/ChargeParcel/ChargeParcel.H

@@ -0,0 +1,367 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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::ChargeParcel
+
+Group
+    grpLagrangianIntermediateParcels
+
+Description
+    Thermodynamic parcel class with one/two-way coupling with the continuous
+    phase. Includes reacting multiphase parcel sub-models, plus:
+    - electric charge
+
+SourceFiles
+    ChargeParcelI.H
+    ChargeParcel.C
+    ChargeParcelIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_ChargeParcel_H
+#define Foam_ChargeParcel_H
+
+#include "particle.H"
+#include "SLGThermo.H"
+#include "demandDrivenEntry.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class ParcelType> class ChargeParcel;
+
+template<class ParcelType>
+Ostream& operator<<(Ostream&, const ChargeParcel<ParcelType>&);
+
+/*---------------------------------------------------------------------------*\
+                        Class ChargeParcel Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class ParcelType>
+class ChargeParcel
+:
+    public ParcelType
+{
+public:
+
+    //- Size in bytes of the fields
+    static const std::size_t sizeofFields;
+
+
+    //- Class to hold charge particle constant properties
+    class constantProperties
+    :
+        public ParcelType::constantProperties
+    {
+    public:
+
+        // Constructors
+
+            //- Null constructor
+            constantProperties();
+
+            //- Copy constructor
+            constantProperties(const constantProperties& cp);
+
+            //- Construct from dictionary
+            constantProperties(const dictionary& parentDict);
+    };
+
+
+    //- Use base tracking data
+    typedef typename ParcelType::trackingData trackingData;
+
+
+protected:
+
+    // Protected data
+
+        // Parcel properties
+
+            //- Electric charge [A s] = [C]
+            scalar eq_;
+
+
+public:
+
+    // Static data members
+
+        //- Runtime type information
+        TypeName("ChargeParcel");
+
+        //- String representation of properties
+        AddToPropertyList
+        (
+            ParcelType,
+            " eq"
+        );
+
+    // Constructors
+
+        //- Construct from mesh, position and topology
+        //  Other properties initialised as null
+        inline ChargeParcel
+        (
+            const polyMesh& mesh,
+            const barycentric& coordinates,
+            const label celli,
+            const label tetFacei,
+            const label tetPti
+        );
+
+        //- Construct from a position and a cell,
+        //- searching for the rest of the required topology.
+        //  Other properties are initialised as null.
+        inline ChargeParcel
+        (
+            const polyMesh& mesh,
+            const vector& position,
+            const label celli
+        );
+
+        //- Construct from components
+        inline ChargeParcel
+        (
+            const polyMesh& mesh,
+            const barycentric& coordinates,
+            const label celli,
+            const label tetFacei,
+            const label tetPti,
+            const label typeId,
+            const scalar nParticle0,
+            const scalar d0,
+            const scalar dTarget0,
+            const vector& U0,
+            const vector& f0,
+            const vector& angularMomentum0,
+            const vector& torque0,
+            const scalarField& Y0,
+            const scalarField& YGas0,
+            const scalarField& YLiquid0,
+            const scalarField& YSolid0,
+            const scalar eq0,
+            const constantProperties& constProps
+        );
+
+        //- Construct from Istream
+        ChargeParcel
+        (
+            const polyMesh& mesh,
+            Istream& is,
+            bool readFields = true,
+            bool newFormat = true
+        );
+
+        //- Construct as a copy
+        ChargeParcel(const ChargeParcel& p);
+
+        //- Construct as a copy
+        ChargeParcel
+        (
+            const ChargeParcel& p,
+            const polyMesh& mesh
+        );
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone() const
+        {
+            return autoPtr<particle>(new ChargeParcel(*this));
+        }
+
+        //- Construct and return a (basic particle) clone
+        virtual autoPtr<particle> clone(const polyMesh& mesh) const
+        {
+            return autoPtr<particle>(new ChargeParcel(*this, mesh));
+        }
+
+        //- Factory class to read-construct particles used for
+        //- parallel transfer
+        class iNew
+        {
+            const polyMesh& mesh_;
+
+        public:
+
+            iNew(const polyMesh& mesh)
+            :
+                mesh_(mesh)
+            {}
+
+            autoPtr<ChargeParcel<ParcelType>> operator()
+            (
+                Istream& is
+            ) const
+            {
+                return autoPtr<ChargeParcel<ParcelType>>
+                (
+                    new ChargeParcel<ParcelType>(mesh_, is, true)
+                );
+            }
+        };
+
+
+    // Member Functions
+
+        // Access
+
+            //- Return const access to electric charge
+            inline scalar eq() const;
+
+
+        // Edit
+
+            //- Return access to electric charge
+            inline scalar& eq();
+
+
+        // Main calculation loop
+
+            //- Set cell values
+            template<class TrackCloudType>
+            void setCellValues(TrackCloudType& cloud, trackingData& td);
+
+            //- Correct cell values using latest transfer information
+            template<class TrackCloudType>
+            void cellValueSourceCorrection
+            (
+                TrackCloudType& cloud,
+                trackingData& td,
+                const scalar dt
+            );
+
+            //- Update parcel properties over the time interval
+            template<class TrackCloudType>
+            void calc
+            (
+                TrackCloudType& cloud,
+                trackingData& td,
+                const scalar dt
+            );
+
+
+        // I-O
+
+            //- Read - composition supplied
+            template<class CloudType, class CompositionType>
+            static void readFields
+            (
+                CloudType& c,
+                const CompositionType& compModel
+            );
+
+            //- Read - no composition
+            template<class CloudType>
+            static void readFields(CloudType& c);
+
+            //- Write - composition supplied
+            template<class CloudType, class CompositionType>
+            static void writeFields
+            (
+                const CloudType& c,
+                const CompositionType& compModel
+            );
+
+            //- Read - no composition
+            template<class CloudType>
+            static void writeFields(const CloudType& c);
+
+            //- Write individual parcel properties to stream
+            void writeProperties
+            (
+                Ostream& os,
+                const wordRes& filters,
+                const word& delim,
+                const bool namesOnly = false
+            ) const;
+
+            //- Read particle fields as objects from the obr registry
+            //  - no composition
+            template<class CloudType>
+            static void readObjects
+            (
+                CloudType& c,
+                const objectRegistry& obr
+            );
+
+            //- Read particle fields as objects from the obr registry
+            template<class CloudType, class CompositionType>
+            static void readObjects
+            (
+                CloudType& c,
+                const CompositionType& compModel,
+                const objectRegistry& obr
+            );
+
+            //- Write particle fields as objects into the obr registry
+            //  - no composition
+            template<class CloudType>
+            static void writeObjects
+            (
+                const CloudType& c,
+                objectRegistry& obr
+            );
+
+            //- Write particle fields as objects into the obr registry
+            template<class CloudType, class CompositionType>
+            static void writeObjects
+            (
+                const CloudType& c,
+                const CompositionType& compModel,
+                objectRegistry& obr
+            );
+
+
+    // Ostream Operator
+
+        friend Ostream& operator<< <ParcelType>
+        (
+            Ostream&,
+            const ChargeParcel<ParcelType>&
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "ChargeParcelI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ChargeParcel.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/ChargeParcel/ChargeParcelI.H

@@ -0,0 +1,153 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline Foam::ChargeParcel<ParcelType>::constantProperties::
+constantProperties()
+:
+    ParcelType::constantProperties()
+{}
+
+
+template<class ParcelType>
+inline Foam::ChargeParcel<ParcelType>::constantProperties::
+constantProperties
+(
+    const constantProperties& cp
+)
+:
+    ParcelType::constantProperties(cp)
+{}
+
+
+template<class ParcelType>
+inline Foam::ChargeParcel<ParcelType>::constantProperties::
+constantProperties
+(
+    const dictionary& parentDict
+)
+:
+    ParcelType::constantProperties(parentDict)
+{}
+
+
+template<class ParcelType>
+inline Foam::ChargeParcel<ParcelType>::ChargeParcel
+(
+    const polyMesh& mesh,
+    const barycentric& coordinates,
+    const label celli,
+    const label tetFacei,
+    const label tetPti
+)
+:
+    ParcelType(mesh, coordinates, celli, tetFacei, tetPti),
+    eq_(0.0)
+{}
+
+
+template<class ParcelType>
+inline Foam::ChargeParcel<ParcelType>::ChargeParcel
+(
+    const polyMesh& mesh,
+    const vector& position,
+    const label celli
+)
+:
+    ParcelType(mesh, position, celli),
+    eq_(0.0)
+{}
+
+
+template<class ParcelType>
+inline Foam::ChargeParcel<ParcelType>::ChargeParcel
+(
+    const polyMesh& mesh,
+    const barycentric& coordinates,
+    const label celli,
+    const label tetFacei,
+    const label tetPti,
+    const label typeId,
+    const scalar nParticle0,
+    const scalar d0,
+    const scalar dTarget0,
+    const vector& U0,
+    const vector& f0,
+    const vector& angularMomentum0,
+    const vector& torque0,
+    const scalarField& Y0,
+    const scalarField& YGas0,
+    const scalarField& YLiquid0,
+    const scalarField& YSolid0,
+    const scalar eq0,
+    const constantProperties& constProps
+)
+:
+    ParcelType
+    (
+        mesh,
+        coordinates,
+        celli,
+        tetFacei,
+        tetPti,
+        typeId,
+        nParticle0,
+        d0,
+        dTarget0,
+        U0,
+        f0,
+        angularMomentum0,
+        torque0,
+        Y0,
+        YGas0,
+        YLiquid0,
+        YSolid0,
+        constProps
+    ),
+    eq_(eq0)
+{}
+
+
+// * * * * * * * * * * ChargeParcel Member Functions * * * * * * * * * * * * //
+
+template<class ParcelType>
+inline Foam::scalar Foam::ChargeParcel<ParcelType>::eq() const
+{
+    return eq_;
+}
+
+
+template<class ParcelType>
+inline Foam::scalar& Foam::ChargeParcel<ParcelType>::eq()
+{
+    return eq_;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/Templates/ChargeParcel/ChargeParcelIO.C

@@ -0,0 +1,346 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "ChargeParcel.H"
+#include "IOstreams.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::string Foam::ChargeParcel<ParcelType>::propertyList_ =
+    Foam::ChargeParcel<ParcelType>::propertyList();
+
+
+template<class ParcelType>
+const std::size_t Foam::ChargeParcel<ParcelType>::sizeofFields
+(
+    sizeof(ChargeParcel<ParcelType>)
+  - offsetof(ChargeParcel<ParcelType>, eq_)
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::ChargeParcel<ParcelType>::ChargeParcel
+(
+    const polyMesh& mesh,
+    Istream& is,
+    bool readFields,
+    bool newFormat
+)
+:
+    ParcelType(mesh, is, readFields, newFormat),
+    eq_(0.0)
+{
+    if (readFields)
+    {
+        if (is.format() == IOstreamOption::ASCII)
+        {
+            is  >> eq_;
+        }
+        else if (!is.checkLabelSize<>() || !is.checkScalarSize<>())
+        {
+            // Non-native label or scalar size
+
+            is.beginRawRead();
+
+            readRawScalar(is, &eq_);
+
+            is.endRawRead();
+        }
+        else
+        {
+            is.read(reinterpret_cast<char*>(&eq_), sizeofFields);
+        }
+    }
+
+    is.check(FUNCTION_NAME);
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::ChargeParcel<ParcelType>::readFields(CloudType& c)
+{
+    const bool valid = c.size();
+
+    ParcelType::readFields(c);
+
+    IOField<scalar> eq
+    (
+        c.fieldIOobject("eq", IOobject::MUST_READ),
+        valid
+    );
+    c.checkFieldIOobject(c, eq);
+
+    label i = 0;
+
+    for (ChargeParcel<ParcelType>& p : c)
+    {
+        p.eq_ = eq[i];
+
+        ++i;
+    }
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::ChargeParcel<ParcelType>::readFields
+(
+    CloudType& c,
+    const CompositionType& compModel
+)
+{
+    const bool valid = c.size();
+
+    ParcelType::readFields(c, compModel);
+
+    IOField<scalar> eq
+    (
+        c.fieldIOobject("eq", IOobject::MUST_READ),
+        valid
+    );
+    c.checkFieldIOobject(c, eq);
+
+    label i = 0;
+
+    for (ChargeParcel<ParcelType>& p : c)
+    {
+        p.eq_ = eq[i];
+
+        ++i;
+    }
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::ChargeParcel<ParcelType>::writeFields(const CloudType& c)
+{
+    ParcelType::writeFields(c);
+
+    const label np = c.size();
+    const bool valid = np;
+
+    IOField<scalar> eq(c.fieldIOobject("eq", IOobject::NO_READ), np);
+
+    label i = 0;
+
+    for (const ChargeParcel<ParcelType>& p : c)
+    {
+        eq[i] = p.eq();
+
+        ++i;
+    }
+
+    eq.write(valid);
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::ChargeParcel<ParcelType>::writeFields
+(
+    const CloudType& c,
+    const CompositionType& compModel
+)
+{
+    ParcelType::writeFields(c, compModel);
+
+    const label np = c.size();
+    const bool valid = np;
+
+    IOField<scalar> eq(c.fieldIOobject("eq", IOobject::NO_READ), np);
+
+    label i = 0;
+
+    for (const ChargeParcel<ParcelType>& p : c)
+    {
+        eq[i] = p.eq();
+
+        ++i;
+    }
+
+    eq.write(valid);
+}
+
+
+template<class ParcelType>
+void Foam::ChargeParcel<ParcelType>::writeProperties
+(
+    Ostream& os,
+    const wordRes& filters,
+    const word& delim,
+    const bool namesOnly
+) const
+{
+    ParcelType::writeProperties(os, filters, delim, namesOnly);
+
+    #undef  writeProp
+    #define writeProp(Name, Value)                                            \
+        ParcelType::writeProperty(os, Name, Value, namesOnly, delim, filters)
+
+    writeProp("eq", eq_);
+
+    #undef writeProp
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::ChargeParcel<ParcelType>::readObjects
+(
+    CloudType& c,
+    const objectRegistry& obr
+)
+{
+    ParcelType::readObjects(c, obr);
+
+    if (!c.size()) return;
+
+    const auto& eq = cloud::lookupIOField<scalar>("eq", obr);
+
+    label i = 0;
+
+    for (ChargeParcel<ParcelType>& p : c)
+    {
+        p.eq_ = eq[i];
+
+        ++i;
+    }
+}
+
+
+template<class ParcelType>
+template<class CloudType>
+void Foam::ChargeParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+
+    const label np = c.size();
+
+    auto& eq = cloud::createIOField<scalar>("eq", np, obr);
+
+    label i = 0;
+
+    for (const ChargeParcel<ParcelType>& p : c)
+    {
+        eq[i] = p.eq();
+
+        ++i;
+    }
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::ChargeParcel<ParcelType>::readObjects
+(
+    CloudType& c,
+    const CompositionType& compModel,
+    const objectRegistry& obr
+)
+{
+    ParcelType::readObjects(c, obr);
+
+    if (!c.size()) return;
+
+    const auto& eq = cloud::lookupIOField<scalar>("eq", obr);
+
+    label i = 0;
+
+    for (ChargeParcel<ParcelType>& p : c)
+    {
+        p.eq_ = eq[i];
+
+        ++i;
+    }
+}
+
+
+template<class ParcelType>
+template<class CloudType, class CompositionType>
+void Foam::ChargeParcel<ParcelType>::writeObjects
+(
+    const CloudType& c,
+    const CompositionType& compModel,
+    objectRegistry& obr
+)
+{
+    ParcelType::writeObjects(c, obr);
+
+    const label np = c.size();
+
+    auto& eq = cloud::createIOField<scalar>("eq", np, obr);
+
+    label i = 0;
+
+    for (const ChargeParcel<ParcelType>& p : c)
+    {
+        eq[i] = p.eq();
+
+        ++i;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+template<class ParcelType>
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const ChargeParcel<ParcelType>& p
+)
+{
+    if (os.format() == IOstreamOption::ASCII)
+    {
+        os  << static_cast<const ParcelType&>(p)
+            << token::SPACE << p.eq();
+    }
+    else
+    {
+        os  << static_cast<const ParcelType&>(p);
+        os.write
+        (
+            reinterpret_cast<const char*>(&p.eq_),
+            KinematicParcel<ParcelType>::sizeofFields
+        );
+    }
+
+    os.check(FUNCTION_NAME);
+    return os;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/chargeParcel/chargeParcel.H

@@ -0,0 +1,81 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2019 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::chargeParcel
+
+Description
+    Definition of reacting multiphase parcel
+
+SourceFiles
+    chargeParcel.C
+    chargeParcelIO.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_chargeParcel_H
+#define Foam_chargeParcel_H
+
+#include "contiguous.H"
+#include "particle.H"
+#include "KinematicParcel.H"
+#include "ThermoParcel.H"
+#include "ReactingParcel.H"
+#include "ReactingMultiphaseParcel.H"
+#include "ChargeParcel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef ChargeParcel
+    <
+        ReactingMultiphaseParcel
+        <
+            ReactingParcel
+            <
+                ThermoParcel
+                <
+                    KinematicParcel
+                    <
+                        particle
+                    >
+                >
+            >
+        >
+    > chargeParcel;
+
+    //- Non-contiguous data for chargeParcel
+    template<>
+    struct is_contiguous<chargeParcel> : std::false_type {};
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/chargeParcel/defineChargeParcel.C

@@ -0,0 +1,40 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "chargeParcel.H"
+#include "Cloud.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug(chargeParcel, 0);
+    defineTemplateTypeNameAndDebug(Cloud<chargeParcel>, 0);
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/derived/chargeParcel/makeChargeParcelSubmodels.C

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "chargeCloud.H"
+
+#include "makeChargeParcelCloudFunctionObjects.H"
+
+// Kinematic
+#include "makeChargeParcelForces.H"
+#include "makeParcelDispersionModels.H"
+#include "makeChargeParcelInjectionModels.H"
+#include "makeParcelPatchInteractionModels.H"
+#include "makeReactingMultiphaseParcelStochasticCollisionModels.H" // MP variant
+#include "makeReactingParcelSurfaceFilmModels.H" // Reacting variant
+
+// Thermodynamic
+#include "makeParcelHeatTransferModels.H"
+
+// Reacting
+#include "makeReactingMultiphaseParcelCompositionModels.H" // MP Variant
+#include "makeReactingParcelPhaseChangeModels.H"
+
+// Reacting multiphase
+#include "makeReactingMultiphaseParcelDevolatilisationModels.H"
+#include "makeReactingMultiphaseParcelSurfaceReactionModels.H"
+
+// MPPIC sub-models
+#include "makeMPPICParcelDampingModels.H"
+#include "makeMPPICParcelIsotropyModels.H"
+#include "makeMPPICParcelPackingModels.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makeChargeParcelCloudFunctionObjects(chargeCloud);
+
+// Kinematic sub-models
+makeChargeParcelForces(chargeCloud);
+makeParcelDispersionModels(chargeCloud);
+makeChargeParcelInjectionModels(chargeCloud);
+makeParcelPatchInteractionModels(chargeCloud);
+makeReactingMultiphaseParcelStochasticCollisionModels
+(
+    chargeCloud
+);
+makeReactingParcelSurfaceFilmModels(chargeCloud);
+
+// Thermo sub-models
+makeParcelHeatTransferModels(chargeCloud);
+
+// Reacting sub-models
+makeReactingMultiphaseParcelCompositionModels(chargeCloud);
+makeReactingParcelPhaseChangeModels(chargeCloud);
+
+// Reacting multiphase sub-models
+makeReactingMultiphaseParcelDevolatilisationModels
+(
+    chargeCloud
+);
+makeReactingMultiphaseParcelSurfaceReactionModels
+(
+    chargeCloud
+);
+
+// MPPIC sub-models
+makeMPPICParcelDampingModels(chargeCloud);
+makeMPPICParcelIsotropyModels(chargeCloud);
+makeMPPICParcelPackingModels(chargeCloud);
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/include/makeChargeParcelCloudFunctionObjects.H

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeChargeParcelCloudFunctionObjects_H
+#define makeChargeParcelCloudFunctionObjects_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "FaceInteraction.H"
+#include "FacePostProcessing.H"
+#include "ParticleCollector.H"
+#include "ParticleErosion.H"
+#include "ParticleTracks.H"
+#include "ParticleTrap.H"
+#include "ParticleZoneInfo.H"
+#include "PatchCollisionDensity.H"
+#include "PatchInteractionFields.H"
+#include "PatchPostProcessing.H"
+#include "PatchParticleHistogram.H"
+#include "RemoveParcels.H"
+#include "VoidFraction.H"
+#include "NusseltNumber.H"
+#include "HeatTransferCoeff.H"
+#include "ThermoReynoldsNumber.H"
+#include "WeberNumberReacting.H"
+#include "ParticleDose.H"
+#include "SpaceChargeDensity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeChargeParcelCloudFunctionObjects(CloudType)                        \
+                                                                               \
+    makeCloudFunctionObject(CloudType);                                        \
+                                                                               \
+    makeCloudFunctionObjectType(FaceInteraction, CloudType);                   \
+    makeCloudFunctionObjectType(FacePostProcessing, CloudType);                \
+    makeCloudFunctionObjectType(ParticleCollector, CloudType);                 \
+    makeCloudFunctionObjectType(ParticleErosion, CloudType);                   \
+    makeCloudFunctionObjectType(ParticleTracks, CloudType);                    \
+    makeCloudFunctionObjectType(ParticleTrap, CloudType);                      \
+    makeCloudFunctionObjectType(ParticleZoneInfo, CloudType);                  \
+    makeCloudFunctionObjectType(PatchCollisionDensity, CloudType);             \
+    makeCloudFunctionObjectType(PatchInteractionFields, CloudType);            \
+    makeCloudFunctionObjectType(PatchPostProcessing, CloudType);               \
+    makeCloudFunctionObjectType(PatchParticleHistogram, CloudType);            \
+    makeCloudFunctionObjectType(RemoveParcels, CloudType);                     \
+    makeCloudFunctionObjectType(VoidFraction, CloudType);                      \
+    makeCloudFunctionObjectType(NusseltNumber, CloudType);                     \
+    makeCloudFunctionObjectType(HeatTransferCoeff, CloudType);                 \
+    makeCloudFunctionObjectType(ThermoReynoldsNumber, CloudType);              \
+    makeCloudFunctionObjectType(WeberNumberReacting, CloudType);               \
+    makeCloudFunctionObjectType(ParticleDose, CloudType);                      \
+    makeCloudFunctionObjectType(SpaceChargeDensity, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/include/makeChargeParcelForces.H

@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_makeChargeParcelForces_H
+#define Foam_makeChargeParcelForces_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "SphereDragForce.H"
+#include "NonSphereDragForce.H"
+
+#include "SaffmanMeiLiftForce.H"
+#include "TomiyamaLiftForce.H"
+
+#include "GravityForce.H"
+#include "NonInertialFrameForce.H"
+#include "ParamagneticForce.H"
+#include "PressureGradientForce.H"
+#include "SRFForce.H"
+#include "VirtualMassForce.H"
+
+#include "CoulombForce.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeChargeParcelForces(CloudType)                                      \
+                                                                               \
+    makeParticleForceModel(CloudType);                                         \
+    makeParticleForceModelType(SphereDragForce, CloudType);                    \
+    makeParticleForceModelType(NonSphereDragForce, CloudType);                 \
+    makeParticleForceModelType(SaffmanMeiLiftForce, CloudType);                \
+    makeParticleForceModelType(TomiyamaLiftForce, CloudType);                  \
+    makeParticleForceModelType(GravityForce, CloudType);                       \
+    makeParticleForceModelType(NonInertialFrameForce, CloudType);              \
+    makeParticleForceModelType(ParamagneticForce, CloudType);                  \
+    makeParticleForceModelType(PressureGradientForce, CloudType);              \
+    makeParticleForceModelType(SRFForce, CloudType);                           \
+    makeParticleForceModelType(VirtualMassForce, CloudType);                   \
+    makeParticleForceModelType(CoulombForce, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/parcels/include/makeChargeParcelInjectionModels.H

@@ -0,0 +1,63 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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/>.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef makeChargeParcelInjectionModels_H
+#define makeChargeParcelInjectionModels_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "CellZoneInjection.H"
+#include "ConeInjection.H"
+#include "ConeNozzleInjection.H"
+#include "FieldActivatedInjection.H"
+#include "ManualInjection.H"
+#include "NoInjection.H"
+#include "PatchInjection.H"
+#include "PatchFlowRateInjection.H"
+#include "ChargeLookupTableInjection.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#define makeChargeParcelInjectionModels(CloudType)                 \
+                                                                               \
+    makeInjectionModel(CloudType);                                             \
+    makeInjectionModelType(CellZoneInjection, CloudType);                      \
+    makeInjectionModelType(ConeInjection, CloudType);                          \
+    makeInjectionModelType(ConeNozzleInjection, CloudType);                    \
+    makeInjectionModelType(FieldActivatedInjection, CloudType);                \
+    makeInjectionModelType(ManualInjection, CloudType);                        \
+    makeInjectionModelType(NoInjection, CloudType);                            \
+    makeInjectionModelType(PatchInjection, CloudType);                         \
+    makeInjectionModelType(PatchFlowRateInjection, CloudType);                 \
+    makeInjectionModelType(ChargeLookupTableInjection, CloudType);
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/CloudFunctionObjects/SpaceChargeDensity/SpaceChargeDensity.C

@@ -0,0 +1,244 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "SpaceChargeDensity.H"
+
+// * * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::SpaceChargeDensity<CloudType>::write()
+{
+    if (erhoPtr_)
+    {
+        erhoPtr_->write();
+    }
+    else
+    {
+        FatalErrorInFunction
+            << "erhoPtr not valid" << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::SpaceChargeDensity<CloudType>::SpaceChargeDensity
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& modelName
+)
+:
+    CloudFunctionObject<CloudType>(dict, owner, modelName, typeName),
+    mPtr_(nullptr),
+    qPtr_(nullptr),
+    cPtr_(nullptr),
+    erhoPtr_(nullptr)
+{}
+
+
+template<class CloudType>
+Foam::SpaceChargeDensity<CloudType>::SpaceChargeDensity
+(
+    const SpaceChargeDensity<CloudType>& vf
+)
+:
+    CloudFunctionObject<CloudType>(vf),
+    mPtr_(nullptr),
+    qPtr_(nullptr),
+    cPtr_(nullptr),
+    erhoPtr_(nullptr)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::SpaceChargeDensity<CloudType>::preEvolve
+(
+    const typename parcelType::trackingData& td
+)
+{
+    const fvMesh& mesh = this->owner().mesh();
+
+    if (mPtr_)
+    {
+        mPtr_->primitiveFieldRef() = 0.0;
+    }
+    else
+    {
+        mPtr_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    this->owner().name() + "m",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                ),
+                mesh,
+                dimensionedScalar(dimMass, Zero)
+            )
+        );
+    }
+
+    if (qPtr_)
+    {
+        qPtr_->primitiveFieldRef() = 0.0;
+    }
+    else
+    {
+        qPtr_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    this->owner().name() + "q",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                ),
+                mesh,
+                dimensionedScalar(dimCurrent*dimTime, Zero)
+            )
+        );
+    }
+
+    if (cPtr_)
+    {
+        cPtr_->primitiveFieldRef() = 0.0;
+    }
+    else
+    {
+        cPtr_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    this->owner().name() + "c",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    IOobject::NO_REGISTER
+                ),
+                mesh,
+                dimensionedScalar(dimMass*dimTime, Zero)
+            )
+        );
+    }
+
+    if (erhoPtr_)
+    {
+        erhoPtr_->primitiveFieldRef() = 0.0;
+    }
+    else
+    {
+        erhoPtr_.reset
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    this->owner().name() + "erho",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
+                mesh,
+                dimensionedScalar(dimCurrent*dimTime/dimVolume, Zero)
+            )
+        );
+    }
+}
+
+
+template<class CloudType>
+void Foam::SpaceChargeDensity<CloudType>::postEvolve
+(
+    const typename parcelType::trackingData& td
+)
+{
+    volScalarField& m = mPtr_();
+    volScalarField& q = qPtr_();
+    volScalarField& c = cPtr_();
+    volScalarField& erho = erhoPtr_();
+
+    auto& own = this->owner();
+    const fvMesh& mesh = own.mesh();
+
+    c.primitiveFieldRef() /= mesh.time().deltaTValue()*mesh.V();
+
+    // (YSSD:Eq. 4)
+
+    // sum masses and charges in cells
+    forAllConstIters(own, parcelIter)
+    {
+        const parcelType& p = parcelIter();
+        m[p.cell()] += p.mass()*p.nParticle();
+        q[p.cell()] += p.eq()*p.nParticle();
+    }
+
+    forAllConstIters(own, parcelIter)
+    {
+        const parcelType& p = parcelIter();
+
+        // (YSSD:Eq. 5)
+        erho[p.cell()] = q[p.cell()]/m[p.cell()]*c[p.cell()];
+    }
+
+    CloudFunctionObject<CloudType>::postEvolve(td);
+}
+
+
+template<class CloudType>
+void Foam::SpaceChargeDensity<CloudType>::postMove
+(
+    parcelType& p,
+    const scalar dt,
+    const point&,
+    bool&
+)
+{
+    volScalarField& c = cPtr_();
+
+    // (YSSD:Eq. 3)
+    c[p.cell()] += p.mass()*p.nParticle()*dt;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/CloudFunctionObjects/SpaceChargeDensity/SpaceChargeDensity.H

@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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::SpaceChargeDensity
+
+Group
+    grpLagrangianIntermediateFunctionObjects
+
+Description
+
+    References:
+    \verbatim
+        Governing equations (tag:YSSD):
+            Ye, Q., Steigleder, T., Scheibe, A., & Domnick, J. (2002).
+            Numerical simulation of the electrostatic powder coating process
+            with a corona spray gun. Journal of Electrostatics, 54(2), 189-205.
+            DOI:10.1016/S0304-3886(01)00181-4
+    \endverbatim
+
+SourceFiles
+    SpaceChargeDensity.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_SpaceChargeDensity_H
+#define Foam_SpaceChargeDensity_H
+
+#include "CloudFunctionObject.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class SpaceChargeDensity Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class SpaceChargeDensity
+:
+    public CloudFunctionObject<CloudType>
+{
+    // Private Data
+
+        // Typedefs
+
+            //- Convenience typedef for parcel type
+            typedef typename CloudType::parcelType parcelType;
+
+        //- Particle mass field
+        autoPtr<volScalarField> mPtr_;
+
+        //- Particle electric-charge field
+        autoPtr<volScalarField> qPtr_;
+
+        //- Particle specific concentration field
+        autoPtr<volScalarField> cPtr_;
+
+        //- Space charge density field
+        autoPtr<volScalarField> erhoPtr_;
+
+
+protected:
+
+    // Protected Member Functions
+
+        //- Write post-processing info
+        virtual void write();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("spaceChargeDensity");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        SpaceChargeDensity
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& modelName
+        );
+
+        //- Copy construct
+        SpaceChargeDensity(const SpaceChargeDensity<CloudType>& vf);
+
+        //- Construct and return a clone
+        virtual autoPtr<CloudFunctionObject<CloudType>> clone() const
+        {
+            return autoPtr<CloudFunctionObject<CloudType>>
+            (
+                new SpaceChargeDensity<CloudType>(*this)
+            );
+        }
+
+        //- No copy assignment
+        void operator=(const SpaceChargeDensity<CloudType>&) = delete;
+
+
+    //- Destructor
+    virtual ~SpaceChargeDensity() = default;
+
+
+    // Member Functions
+
+        // Evaluation
+
+        //- Pre-evolve hook
+        virtual void preEvolve
+        (
+            const typename parcelType::trackingData& td
+        );
+
+        //- Post-evolve hook
+        virtual void postEvolve
+        (
+            const typename parcelType::trackingData& td
+        );
+
+        //- Post-move hook
+        virtual void postMove
+        (
+            parcelType& p,
+            const scalar dt,
+            const point& position0,
+            bool& keepParticle
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "SpaceChargeDensity.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/InjectionModel/ChargeLookupTableInjection/ChargeLookupTableInjection.C

@@ -0,0 +1,275 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "ChargeLookupTableInjection.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::ChargeLookupTableInjection<CloudType>::
+ChargeLookupTableInjection
+(
+    const dictionary& dict,
+    CloudType& owner,
+    const word& modelName
+)
+:
+    InjectionModel<CloudType>(dict, owner, modelName, typeName),
+    inputFileName_(this->coeffDict().lookup("inputFile")),
+    duration_(this->coeffDict().getScalar("duration")),
+    parcelsPerSecond_(this->coeffDict().getScalar("parcelsPerSecond")),
+    randomise_(this->coeffDict().getBool("randomise")),
+    injectors_
+    (
+        IOobject
+        (
+            inputFileName_,
+            owner.db().time().constant(),
+            owner.db(),
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    injectorCells_(injectors_.size()),
+    injectorTetFaces_(injectors_.size()),
+    injectorTetPts_(injectors_.size())
+{
+    updateMesh();
+
+    duration_ = owner.db().time().userTimeToTime(duration_);
+
+    // Determine volume of particles to inject
+    this->volumeTotal_ = 0.0;
+    forAll(injectors_, i)
+    {
+        this->volumeTotal_ += injectors_[i].mDot()/injectors_[i].rho();
+    }
+    this->volumeTotal_ *= duration_;
+}
+
+
+template<class CloudType>
+Foam::ChargeLookupTableInjection<CloudType>::
+ChargeLookupTableInjection
+(
+    const ChargeLookupTableInjection<CloudType>& im
+)
+:
+    InjectionModel<CloudType>(im),
+    inputFileName_(im.inputFileName_),
+    duration_(im.duration_),
+    parcelsPerSecond_(im.parcelsPerSecond_),
+    randomise_(im.randomise_),
+    injectors_(im.injectors_),
+    injectorCells_(im.injectorCells_),
+    injectorTetFaces_(im.injectorTetFaces_),
+    injectorTetPts_(im.injectorTetPts_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::ChargeLookupTableInjection<CloudType>::updateMesh()
+{
+    // Set/cache the injector cells
+    bitSet reject(injectors_.size());
+
+    // Set/cache the injector cells
+    forAll(injectors_, i)
+    {
+        if
+        (
+            !this->findCellAtPosition
+            (
+                injectorCells_[i],
+                injectorTetFaces_[i],
+                injectorTetPts_[i],
+                injectors_[i].x(),
+                !this->ignoreOutOfBounds_
+
+            )
+        )
+        {
+            reject.set(i);
+        }
+    }
+
+    label nRejected = reject.count();
+
+    if (nRejected)
+    {
+        reject.flip();
+        inplaceSubset(reject, injectorCells_);
+        inplaceSubset(reject, injectorTetFaces_);
+        inplaceSubset(reject, injectorTetPts_);
+        inplaceSubset(reject, injectors_);
+
+        Info<< "    " << nRejected
+            << " positions rejected, out of bounds" << endl;
+    }
+}
+
+
+template<class CloudType>
+Foam::scalar
+Foam::ChargeLookupTableInjection<CloudType>::timeEnd() const
+{
+    return this->SOI_ + duration_;
+}
+
+
+template<class CloudType>
+Foam::label
+Foam::ChargeLookupTableInjection<CloudType>::parcelsToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        return floor(injectorCells_.size()*(time1 - time0)*parcelsPerSecond_);
+    }
+
+    return 0;
+}
+
+
+template<class CloudType>
+Foam::scalar
+Foam::ChargeLookupTableInjection<CloudType>::volumeToInject
+(
+    const scalar time0,
+    const scalar time1
+)
+{
+    scalar volume = 0.0;
+    if ((time0 >= 0.0) && (time0 < duration_))
+    {
+        forAll(injectors_, i)
+        {
+            volume += injectors_[i].mDot()/injectors_[i].rho()*(time1 - time0);
+        }
+    }
+
+    return volume;
+}
+
+
+template<class CloudType>
+void Foam::ChargeLookupTableInjection<CloudType>::setPositionAndCell
+(
+    const label parcelI,
+    const label nParcels,
+    const scalar time,
+    vector& position,
+    label& cellOwner,
+    label& tetFacei,
+    label& tetPti
+)
+{
+    label injectorI = 0;
+    if (randomise_)
+    {
+        Random& rnd = this->owner().rndGen();
+        injectorI = rnd.position<label>(0, injectorCells_.size() - 1);
+    }
+    else
+    {
+        injectorI = parcelI*injectorCells_.size()/nParcels;
+    }
+
+    position = injectors_[injectorI].x();
+    cellOwner = injectorCells_[injectorI];
+    tetFacei = injectorTetFaces_[injectorI];
+    tetPti = injectorTetPts_[injectorI];
+}
+
+
+template<class CloudType>
+void Foam::ChargeLookupTableInjection<CloudType>::setProperties
+(
+    const label parcelI,
+    const label nParcels,
+    const scalar,
+    typename CloudType::parcelType& parcel
+)
+{
+    label injectorI = parcelI*injectorCells_.size()/nParcels;
+
+    // set particle velocity
+    parcel.U() = injectors_[injectorI].U();
+
+    // set particle diameter
+    parcel.d() = injectors_[injectorI].d();
+
+    // set particle density
+    parcel.rho() = injectors_[injectorI].rho();
+
+    // set particle temperature
+    parcel.T() = injectors_[injectorI].T();
+
+    // set particle specific heat capacity
+    parcel.Cp() = injectors_[injectorI].Cp();
+
+    // set particle component mass fractions
+    parcel.Y() = injectors_[injectorI].Y();
+
+    // set particle gaseous component mass fractions
+    parcel.YGas() = injectors_[injectorI].YGas();
+
+    // set particle liquid component mass fractions
+    parcel.YLiquid() = injectors_[injectorI].YLiquid();
+
+    // set particle solid component mass fractions
+    parcel.YSolid() = injectors_[injectorI].YSolid();
+
+    // set particle charge
+    parcel.eq() = injectors_[injectorI].eq();
+}
+
+
+template<class CloudType>
+bool
+Foam::ChargeLookupTableInjection<CloudType>::fullyDescribed() const
+{
+    return true;
+}
+
+
+template<class CloudType>
+bool Foam::ChargeLookupTableInjection<CloudType>::validInjection
+(
+    const label
+)
+{
+    return true;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/InjectionModel/ChargeLookupTableInjection/ChargeLookupTableInjection.H

@@ -0,0 +1,206 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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::ChargeLookupTableInjection
+
+Group
+    grpLagrangianIntermediateInjectionSubModels
+
+Description
+    Particle injection sources read from look-up table. Each row corresponds to
+    an injection site.
+
+    (
+       (x y z) (u v w) d rho mDot T cp (Y0..Y2) (Yg0..YgN) (Yl0..YlN) (Ys0..YsN) eq
+       (x y z) (u v w) d rho mDot T cp (Y0..Y2) (Yg0..YgN) (Yl0..YlN) (Ys0..YsN) eq
+       ...
+       (x y z) (u v w) d rho mDot T cp (Y0..Y2) (Yg0..YgN) (Yl0..YlN) (Ys0..YsN) eq
+    );
+
+    where:
+        x, y, z  = global cartesian coordinates [m]
+        u, v, w  = global cartesian velocity components [m/s]
+        d        = diameter [m]
+        rho      = density [kg/m3]
+        mDot     = mass flow rate [kg/s]
+        T        = temperature [K]
+        cp       = specific heat capacity [J/kg/K]
+        Y(3)     = total mass fraction of gas (Y0), liquid (Y1), solid (Y3)
+        Yg(Ngas) = mass fractions of gaseous components
+        Yl(Nliq) = mass fractions of liquid components
+        Ys(Nsld) = mass fractions of solid components
+        eq       = Electric charge [C]
+
+SourceFiles
+    ChargeLookupTableInjection.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_ChargeLookupTableInjection_H
+#define Foam_ChargeLookupTableInjection_H
+
+#include "InjectionModel.H"
+#include "chargeParcelInjectionDataIOList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                Class ChargeLookupTableInjection Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class ChargeLookupTableInjection
+:
+    public InjectionModel<CloudType>
+{
+    // Private Data
+
+        //- Name of file containing injector/parcel data
+        const word inputFileName_;
+
+        //- Injection duration - common to all injection sources
+        scalar duration_;
+
+        //- Number of parcels per injector - common to all injection sources
+        const scalar parcelsPerSecond_;
+
+        //- Flag to indicate to randomise injection positions
+        bool randomise_;
+
+        //- List of injectors
+        chargeParcelInjectionDataIOList injectors_;
+
+        //- List of cell labels corresponding to injector positions
+        labelList injectorCells_;
+
+        //- List of tetFace labels corresponding to injector positions
+        labelList injectorTetFaces_;
+
+        //- List of tetPt labels corresponding to injector positions
+        labelList injectorTetPts_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("chargeLookupTableInjection");
+
+
+    // Constructors
+
+        //- Construct from dictionary
+        ChargeLookupTableInjection
+        (
+            const dictionary& dict,
+            CloudType& owner,
+            const word& modelName
+        );
+
+        //- Copy construct
+        ChargeLookupTableInjection
+        (
+            const ChargeLookupTableInjection<CloudType>& im
+        );
+
+        //- Construct and return a clone
+        virtual autoPtr<InjectionModel<CloudType>> clone() const
+        {
+            return autoPtr<InjectionModel<CloudType>>
+            (
+                new ChargeLookupTableInjection<CloudType>(*this)
+            );
+        }
+
+
+    //- Destructor
+    virtual ~ChargeLookupTableInjection() = default;
+
+
+    // Member Functions
+
+        //- Set injector locations when mesh is updated
+        virtual void updateMesh();
+
+        //- Return the end-of-injection time
+        scalar timeEnd() const;
+
+        //- Number of parcels to introduce relative to SOI
+        virtual label parcelsToInject(const scalar time0, const scalar time1);
+
+        //- Volume of parcels to introduce relative to SOI
+        virtual scalar volumeToInject(const scalar time0, const scalar time1);
+
+
+        // Injection geometry
+
+            //- Set the injection position and owner cell, tetFace and tetPt
+            virtual void setPositionAndCell
+            (
+                const label parcelI,
+                const label nParcels,
+                const scalar time,
+                vector& position,
+                label& cellOwner,
+                label& tetFacei,
+                label& tetPti
+            );
+
+            //- Set the parcel properties
+            virtual void setProperties
+            (
+                const label parcelI,
+                const label nParcels,
+                const scalar time,
+                typename CloudType::parcelType& parcel
+            );
+
+            //- Flag to identify whether model fully describes the parcel
+            virtual bool fullyDescribed() const;
+
+            //- Return flag to identify whether or not injection of parcelI is
+            //- permitted
+            virtual bool validInjection(const label parcelI);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "ChargeLookupTableInjection.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/InjectionModel/ChargeLookupTableInjection/chargeParcelInjectionData.C

@@ -0,0 +1,64 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "chargeParcelInjectionData.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTypeNameAndDebug(chargeParcelInjectionData, 0);
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::chargeParcelInjectionData::
+chargeParcelInjectionData()
+:
+    reactingMultiphaseParcelInjectionData(),
+    eq_()
+{}
+
+
+Foam::chargeParcelInjectionData::
+chargeParcelInjectionData
+(
+    const dictionary& dict
+)
+:
+    reactingMultiphaseParcelInjectionData(dict),
+    eq_(dict.get<scalar>("eq"))
+{}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::chargeParcelInjectionData::
+~chargeParcelInjectionData()
+{}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/InjectionModel/ChargeLookupTableInjection/chargeParcelInjectionData.H

@@ -0,0 +1,143 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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::chargeParcelInjectionData
+
+Description
+    Container class to provide injection data for electric charge parcels
+
+SourceFiles
+    chargeParcelInjectionData.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_chargeParcelInjectionData_H
+#define Foam_chargeParcelInjectionData_H
+
+#include "reactingMultiphaseParcelInjectionData.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declaration of classes
+class chargeParcelInjectionData;
+
+// Forward declaration of friend functions
+
+Ostream& operator<<
+(
+    Ostream&,
+    const chargeParcelInjectionData&
+);
+
+Istream& operator>>
+(
+    Istream&,
+    chargeParcelInjectionData&
+);
+
+/*---------------------------------------------------------------------------*\
+            Class chargeParcelInjectionData Declaration
+\*---------------------------------------------------------------------------*/
+
+class chargeParcelInjectionData
+:
+    public reactingMultiphaseParcelInjectionData
+{
+protected:
+
+    // Parcel properties
+
+        //- Electric charge [C]
+        scalar eq_;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("chargeParcelInjectionData");
+
+    // Constructors
+
+        //- Null constructor
+        chargeParcelInjectionData();
+
+        //- Construct from dictionary
+        chargeParcelInjectionData(const dictionary& dict);
+
+        //- Construct from Istream
+        chargeParcelInjectionData(Istream& is);
+
+
+    //- Destructor
+    virtual ~chargeParcelInjectionData();
+
+
+    // Access
+
+        //- Return const access to the electric charge
+        inline scalar eq() const;
+
+
+    // Edit
+
+        //- Return access to the electric charge
+        inline scalar& eq();
+
+
+    // I-O
+
+        //- Ostream operator
+        friend Ostream& operator<<
+        (
+            Ostream& os,
+            const chargeParcelInjectionData& data
+        );
+
+        //- Istream operator
+        friend Istream& operator>>
+        (
+            Istream& is,
+            chargeParcelInjectionData& data
+        );
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#include "chargeParcelInjectionDataI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/InjectionModel/ChargeLookupTableInjection/chargeParcelInjectionDataI.H

@@ -0,0 +1,44 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "chargeParcelInjectionData.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline Foam::scalar Foam::chargeParcelInjectionData::eq() const
+{
+    return eq_;
+}
+
+
+inline Foam::scalar& Foam::chargeParcelInjectionData::eq()
+{
+    return eq_;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/InjectionModel/ChargeLookupTableInjection/chargeParcelInjectionDataIO.C

@@ -0,0 +1,77 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 "chargeParcelInjectionData.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+Foam::chargeParcelInjectionData::
+chargeParcelInjectionData(Istream& is)
+:
+    reactingMultiphaseParcelInjectionData(is)
+{
+    is.check("reading eq");
+    is >> eq_;
+
+    is.check(FUNCTION_NAME);
+}
+
+
+// * * * * * * * * * * * * * * * IOstream Operators  * * * * * * * * * * * * //
+
+Foam::Ostream& Foam::operator<<
+(
+    Ostream& os,
+    const chargeParcelInjectionData& data
+)
+{
+    os << static_cast<const reactingMultiphaseParcelInjectionData&>(data);
+
+    os << data.eq_;
+
+    return os;
+}
+
+
+Foam::Istream& Foam::operator>>
+(
+    Istream& is,
+    chargeParcelInjectionData& data
+)
+{
+    is >> static_cast<reactingMultiphaseParcelInjectionData&>(data);
+
+    is.check("reading eq");
+    is >> data.eq_;
+
+    is.check(FUNCTION_NAME);
+
+    return is;
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/InjectionModel/ChargeLookupTableInjection/chargeParcelInjectionDataIOList.C

@@ -0,0 +1,42 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "chargeParcelInjectionDataIOList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    defineTemplateTypeNameAndDebug
+    (
+        GlobalIOList<chargeParcelInjectionData>,
+        0
+    );
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/InjectionModel/ChargeLookupTableInjection/chargeParcelInjectionDataIOList.H

@@ -0,0 +1,54 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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::chargeParcelInjectionDataIOList
+
+Description
+
+SourceFiles
+    chargeParcelInjectionDataIOList.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_chargeParcelInjectionDataIOList_H
+#define Foam_chargeParcelInjectionDataIOList_H
+
+#include "GlobalIOList.H"
+#include "chargeParcelInjectionData.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+    typedef GlobalIOList<chargeParcelInjectionData>
+        chargeParcelInjectionDataIOList;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/ParticleForces/Coulomb/CoulombForce.C

@@ -0,0 +1,134 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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 "CoulombForce.H"
+#include "demandDrivenData.H"
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::volVectorField& Foam::CoulombForce<CloudType>::getOrReadField
+(
+    const word& fieldName
+) const
+{
+    auto* ptr = this->mesh().template getObjectPtr<volVectorField>(fieldName);
+
+    if (!ptr)
+    {
+        ptr = new volVectorField
+        (
+            IOobject
+            (
+                fieldName,
+                this->mesh().time().timeName(),
+                this->mesh(),
+                IOobject::MUST_READ,
+                IOobject::AUTO_WRITE
+            ),
+            this->mesh()
+        );
+        this->mesh().objectRegistry::store(ptr);
+    }
+
+    return *ptr;
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class CloudType>
+Foam::CoulombForce<CloudType>::CoulombForce
+(
+    CloudType& owner,
+    const fvMesh& mesh,
+    const dictionary& dict
+)
+:
+    ParticleForce<CloudType>(owner, mesh, dict, typeName, true),
+    Ename_(this->coeffs().template getOrDefault<word>("E", "E")),
+    EInterpPtr_(nullptr)
+{}
+
+
+template<class CloudType>
+Foam::CoulombForce<CloudType>::CoulombForce
+(
+    const CoulombForce& pf
+)
+:
+    ParticleForce<CloudType>(pf),
+    Ename_(pf.Ename_),
+    EInterpPtr_(pf.EInterpPtr_)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class CloudType>
+void Foam::CoulombForce<CloudType>::cacheFields(const bool store)
+{
+    if (store)
+    {
+        const volVectorField& E = getOrReadField(Ename_);
+
+        EInterpPtr_ = interpolation<vector>::New
+        (
+            this->owner().solution().interpolationSchemes(),
+            E
+        ).ptr();
+    }
+    else
+    {
+        deleteDemandDrivenData(EInterpPtr_);
+    }
+}
+
+
+template<class CloudType>
+Foam::forceSuSp Foam::CoulombForce<CloudType>::calcNonCoupled
+(
+    const typename CloudType::parcelType& p,
+    const typename CloudType::parcelType::trackingData& td,
+    const scalar dt,
+    const scalar mass,
+    const scalar Re,
+    const scalar muc
+) const
+{
+    const interpolation<vector>& EInterp = *EInterpPtr_;
+
+    // (YSSD:Eq. 6 - left term)
+    return forceSuSp
+    (
+        p.eq()*EInterp.interpolate(p.coordinates(), p.currentTetIndices()),
+        Zero
+    );
+}
+
+
+// ************************************************************************* //

src/lagrangian/intermediate/submodels/Charge/ParticleForces/Coulomb/CoulombForce.H

@@ -0,0 +1,188 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2023 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::CoulombForce
+
+Group
+    grpLagrangianIntermediateForceSubModels
+
+Description
+    Particle-electric force model wherein the Coulomb force is calculated.
+
+    \f[
+        \vec{F}_\mathrm{E} = q \, \vec{E}
+    \f]
+
+    where
+
+    \vartable
+        \vec{F}_\mathrm{E}  | Coulomb force                 [kg m s^{-2}]
+        q                   | Electric charge of particle   [A s]
+        \vec{E}             | Electric field    [kg m^2 A^{-1} s^{-3} m^{-1}]
+    \endvartable
+
+    References:
+    \verbatim
+        Governing equations (tag:YSSD):
+            Ye, Q., Steigleder, T., Scheibe, A., & Domnick, J. (2002).
+            Numerical simulation of the electrostatic powder coating process
+            with a corona spray gun. Journal of Electrostatics, 54(2), 189-205.
+            DOI:10.1016/S0304-3886(01)00181-4
+    \endverbatim
+
+Usage
+    Minimal example by using \c constant/cloudProperties:
+    \verbatim
+    subModels
+    {
+        particleForces
+        {
+            Coulomb;
+        }
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property     | Description                       | Type | Reqd | Deflt
+      type         | Type name: Coulomb                | word | yes  | -
+    \endtable
+
+SourceFiles
+    CoulombForce.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_CoulombForce_H
+#define Foam_CoulombForce_H
+
+#include "ParticleForce.H"
+#include "interpolation.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class fvMesh;
+
+/*---------------------------------------------------------------------------*\
+                       Class CoulombForce Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class CloudType>
+class CoulombForce
+:
+    public ParticleForce<CloudType>
+{
+    // Private Data
+
+        //- Name of electric field
+        const word Ename_;
+
+        //- Electric-field interpolator
+        const interpolation<vector>* EInterpPtr_;
+
+
+    // Private Member Functions
+
+        //- Return requested volVectorField from the object registry
+        //- or read+register the field to the object registry
+        volVectorField& getOrReadField(const word& fieldName) const;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("Coulomb");
+
+
+    // Constructors
+
+        //- Construct from mesh
+        CoulombForce
+        (
+            CloudType& owner,
+            const fvMesh& mesh,
+            const dictionary& dict
+        );
+
+        //- Copy construct
+        CoulombForce(const CoulombForce& gf);
+
+        //- Construct and return a clone
+        virtual autoPtr<ParticleForce<CloudType>> clone() const
+        {
+            return autoPtr<ParticleForce<CloudType>>
+            (
+                new CoulombForce<CloudType>(*this)
+            );
+        }
+
+        //- No copy assignment
+        void operator=(const CoulombForce<CloudType>&) = delete;
+
+
+    //- Destructor
+    virtual ~CoulombForce() = default;
+
+
+    // Member Functions
+
+        // Evaluation
+
+        //- Cache fields
+        virtual void cacheFields(const bool store);
+
+        //- Calculate the non-coupled force
+        virtual forceSuSp calcNonCoupled
+        (
+            const typename CloudType::parcelType& p,
+            const typename CloudType::parcelType::trackingData& td,
+            const scalar dt,
+            const scalar mass,
+            const scalar Re,
+            const scalar muc
+        ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "CoulombForce.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //