From aaeddba4669147c1b40e7ca0f3cdc62477499b7b Mon Sep 17 00:00:00 2001 From: Kutalmis Bercin Date: Wed, 28 Jul 2021 10:38:56 +0100 Subject: [PATCH] ENH: electricPotential: new solver function object --- src/functionObjects/solvers/Make/files | 1 + .../electricPotential/electricPotential.C | 423 ++++++++++++++++++ .../electricPotential/electricPotential.H | 275 ++++++++++++ 3 files changed, 699 insertions(+) create mode 100644 src/functionObjects/solvers/electricPotential/electricPotential.C create mode 100644 src/functionObjects/solvers/electricPotential/electricPotential.H diff --git a/src/functionObjects/solvers/Make/files b/src/functionObjects/solvers/Make/files index ce5306b1c8..b37d8bafb0 100644 --- a/src/functionObjects/solvers/Make/files +++ b/src/functionObjects/solvers/Make/files @@ -1,4 +1,5 @@ scalarTransport/scalarTransport.C energyTransport/energyTransport.C +electricPotential/electricPotential.C LIB = $(FOAM_LIBBIN)/libsolverFunctionObjects diff --git a/src/functionObjects/solvers/electricPotential/electricPotential.C b/src/functionObjects/solvers/electricPotential/electricPotential.C new file mode 100644 index 0000000000..92d6b8913d --- /dev/null +++ b/src/functionObjects/solvers/electricPotential/electricPotential.C @@ -0,0 +1,423 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "electricPotential.H" +#include "fvc.H" +#include "fvm.H" +#include "calculatedFvPatchField.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ +namespace functionObjects +{ + defineTypeNameAndDebug(electricPotential, 0); + addToRunTimeSelectionTable(functionObject, electricPotential, dictionary); +} +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +Foam::volScalarField& +Foam::functionObjects::electricPotential::operandField() +{ + if (!foundObject(fieldName_)) + { + auto tfldPtr = tmp::New + ( + IOobject + ( + fieldName_, + mesh_.time().timeName(), + mesh_, + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh_ + ); + store(fieldName_, tfldPtr); + } + + return lookupObjectRef(fieldName_); +} + + +Foam::tmp +Foam::functionObjects::electricPotential::sigma() const +{ + const IOobject sigmaIO + ( + IOobject::scopedName(typeName, "sigma"), + mesh_.time().timeName(), + mesh_.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ); + + if (phases_.size()) + { + tmp tsigma = phases_[0]*sigmas_[0]; + + for (label i = 1; i < phases_.size(); ++i) + { + tsigma.ref() += phases_[i]*sigmas_[i]; + } + + return tmp::New + ( + sigmaIO, + tsigma, + calculatedFvPatchField::typeName + ); + } + + return tmp::New + ( + sigmaIO, + mesh_, + sigma_, + calculatedFvPatchField::typeName + ); +} + + +Foam::tmp +Foam::functionObjects::electricPotential::epsilonm() const +{ + // Vacuum permittivity (aka the electric constant) [A^2 s^4/(kg m^3)] + const dimensionedScalar epsilon0 + ( + sqr(dimCurrent)*pow4(dimTime)/(dimMass*pow3(dimLength)), + 8.8541878128e-12 // CODATA value + ); + + const IOobject epsilonrIO + ( + IOobject::scopedName(typeName, "epsilonr"), + mesh_.time().timeName(), + mesh_.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ); + + if (phases_.size()) + { + tmp tepsilonr = phases_[0]*epsilonrs_[0]; + + for (label i = 1; i < phases_.size(); ++i) + { + tepsilonr.ref() += phases_[i]*epsilonrs_[i]; + } + + return tmp::New + ( + epsilonrIO, + epsilon0*tepsilonr, + calculatedFvPatchField::typeName + ); + } + + return tmp::New + ( + epsilonrIO, + mesh_, + epsilon0*epsilonr_, + calculatedFvPatchField::typeName + ); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::functionObjects::electricPotential::electricPotential +( + const word& name, + const Time& runTime, + const dictionary& dict +) +: + fvMeshFunctionObject(name, runTime, dict), + phasesDict_(dict.subOrEmptyDict("phases")), + phaseNames_(), + phases_(), + sigmas_(), + sigma_ + ( + dimensionedScalar + ( + sqr(dimCurrent)*pow3(dimTime)/(dimMass*pow3(dimLength)), + dict.getCheckOrDefault + ( + "sigma", + scalar(1), + scalarMinMax::ge(SMALL) + ) + ) + ), + epsilonrs_(), + epsilonr_ + ( + dimensionedScalar + ( + dimless, + dict.getCheckOrDefault + ( + "epsilonr", + scalar(1), + scalarMinMax::ge(SMALL) + ) + ) + ), + fieldName_ + ( + dict.getOrDefault + ( + "field", + IOobject::scopedName(typeName, "V") + ) + ), + nCorr_(1), + writeDerivedFields_(false) +{ + read(dict); + + // Force creation of transported field so any BCs using it can + // look it up + volScalarField& eV = operandField(); + eV.correctBoundaryConditions(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::functionObjects::electricPotential::read(const dictionary& dict) +{ + if (fvMeshFunctionObject::read(dict)) + { + Log << type() << " read: " << name() << endl; + + dict.readIfPresent("sigma", sigma_); + dict.readIfPresent("epsilonr", epsilonr_); + dict.readIfPresent("nCorr", nCorr_); + dict.readIfPresent("writeDerivedFields", writeDerivedFields_); + + // If flow is multiphase + if (!phasesDict_.empty()) + { + phaseNames_.setSize(phasesDict_.size()); + phases_.setSize(phasesDict_.size()); + sigmas_.setSize(phasesDict_.size()); + + if (writeDerivedFields_) + { + epsilonrs_.setSize(phasesDict_.size()); + } + + label phasei = 0; + forAllConstIters(phasesDict_, iter) + { + const word& key = iter().keyword(); + + if (!phasesDict_.isDict(key)) + { + FatalErrorInFunction + << "Found non-dictionary entry " << iter() + << " in top-level dictionary " << phasesDict_ + << exit(FatalError); + } + + const dictionary& subDict = phasesDict_.subDict(key); + + phaseNames_[phasei] = key; + + sigmas_.set + ( + phasei, + new dimensionedScalar + ( + sqr(dimCurrent)*pow3(dimTime)/(dimMass*pow3(dimLength)), + subDict.getCheck + ( + "sigma", + scalarMinMax::ge(SMALL) + ) + ) + ); + + if (writeDerivedFields_) + { + epsilonrs_.set + ( + phasei, + new dimensionedScalar + ( + dimless, + subDict.getCheck + ( + "epsilonr", + scalarMinMax::ge(SMALL) + ) + ) + ); + } + + ++phasei; + } + + forAll(phaseNames_, i) + { + phases_.set + ( + i, + mesh_.getObjectPtr(phaseNames_[i]) + ); + } + } + + return true; + } + + return false; +} + + +bool Foam::functionObjects::electricPotential::execute() +{ + Log << type() << " execute: " << name() << endl; + + tmp tsigma = this->sigma(); + const volScalarField& sigma = tsigma(); + + volScalarField& eV = operandField(); + + for (label i = 1; i <= nCorr_; ++i) + { + fvScalarMatrix eVEqn + ( + - fvm::laplacian(sigma, eV) + ); + + eVEqn.relax(); + + eVEqn.solve(); + } + + Log << endl; + + return true; +} + + +bool Foam::functionObjects::electricPotential::write() +{ + Log << type() << " write: " << name() << nl + << tab << fieldName_ + << endl; + + volScalarField& eV = operandField(); + + if (writeDerivedFields_) + { + // Write the electric field + const volVectorField E + ( + IOobject + ( + IOobject::scopedName(typeName, "E"), + mesh_.time().timeName(), + mesh_.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + -fvc::grad(eV), + calculatedFvPatchField::typeName + ); + + Log << tab << E.name() << endl; + + E.write(); + + + // Write the current density field + tmp tsigma = this->sigma(); + + auto eJ = tmp::New + ( + IOobject + ( + IOobject::scopedName(typeName, "J"), + mesh_.time().timeName(), + mesh_.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + -tsigma*fvc::grad(eV), + calculatedFvPatchField::typeName + ); + + Log << tab << eJ().name() << endl; + + eJ->write(); + + + // Write the free-charge density field + tmp tepsilonm = this->epsilonm(); + + auto erho = tmp::New + ( + IOobject + ( + IOobject::scopedName(typeName, "rho"), + mesh_.time().timeName(), + mesh_.time(), + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + fvc::div(tepsilonm*E), + calculatedFvPatchField::typeName + ); + + Log << tab << erho().name() << endl; + + erho->write(); + } + + eV.write(); + + return true; +} + + +// ************************************************************************* // diff --git a/src/functionObjects/solvers/electricPotential/electricPotential.H b/src/functionObjects/solvers/electricPotential/electricPotential.H new file mode 100644 index 0000000000..83fdfbd7b2 --- /dev/null +++ b/src/functionObjects/solvers/electricPotential/electricPotential.H @@ -0,0 +1,275 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 OpenCFD Ltd. +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Class + Foam::functionObjects::electricPotential + +Group + grpSolversFunctionObjects + +Description + Computes the steady-state equation of charge conservation to obtain + the electric potential by strictly assuming a quasi-static electrostatic + field for single-phase and multiphase applications. + + The steady-state equation of the charge conservation: + + \f[ + \nabla \cdot \left( \sigma \nabla V \right) = 0 + \f] + + where + \vartable + V | Electric potential [volt = kg m^2/(A s^3)] + \sigma | Isotropic conductivity of mixture [S/m = A^2 s^3/(kg m^3)] + \endvartable + + Optionally, electric field, current density and free-charge + density fields can be written out by using the following equations: + + \f[ + \vec{E} = - \nabla V + \f] + + \f[ + \vec{J} = \sigma \vec{E} = - \sigma \nabla V + \f] + + \f[ + \rho_E = \nabla \cdot \left(\epsilon_m \vec{E} \right) + = \nabla \cdot \left(\epsilon_0 \epsilon_r \vec{E} \right) + \f] + + where + \vartable + \vec{E} | Electric field [m kg/(s^3 A)] + \vec{J} | Current density [A/m^2] + \rho_E | Volume charge density [C/m^3 = A s/m^3] + \epsilon_m | Isotropic permittivity of mixture [F/m = A^2 s^4/(kg m^3)] + \epsilon_0 | Isotropic vacuum permittivity [F/m = A^2 s^4/(kg m^3)] + \epsilon_r | Isotropic relative permittivity of mixture [-] + \endvartable + + For multiphase applications, \c sigma and \c epsilonr are blended + (to consider their interface values) by using the simple weighted + arithmetic mean interpolation, for example: + + \f[ + \sigma = \alpha_1 \sigma_1 + \alpha_2 \sigma_2 + = \alpha_1 \sigma_1 + (1 - \alpha_1) \sigma_2 + \f] + +Usage + Minimal example by using \c system/controlDict.functions: + \verbatim + electricPotential1 + { + // Mandatory entries + type electricPotential; + libs (solverFunctionObjects); + + // Conditional entries + + // Option-1: single-phase + sigma ; + epsilonr ; + + // Option-2: multiphase + phases + { + alpha.air + { + sigma ; + epsilonr ; + } + alpha.water + { + sigma ; + epsilonr ; + } + alpha.mercury + { + sigma ; + epsilonr ; + } + ... + } + + // Optional entries + nCorr