From 6370bb0f5da643f3726abf983f521961c3de82b1 Mon Sep 17 00:00:00 2001 From: ignis Date: Wed, 14 Dec 2016 20:49:28 +0900 Subject: [PATCH] plasmaReactingFoam - snGrad of E for boundary drift calc / Disable mvConvection RTS / Separate surfaceFields for ion fluxes / Write electron flux and ions fluxes --- .../combustion/plasmaReactingFoam/PhiEqn.H | 1 + .../combustion/plasmaReactingFoam/YEqn.H | 26 +++++---- .../plasmaReactingFoam/createFields.H | 56 +++++++++++++++++-- .../plasmaReactingFoam/plasmaReactingFoam.C | 5 +- 4 files changed, 70 insertions(+), 18 deletions(-) diff --git a/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H b/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H index 6a86a1ea..8a912c7d 100644 --- a/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H +++ b/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H @@ -8,6 +8,7 @@ ); E = -fvc::grad(Phi); + snE = -fvc::snGrad(Phi); tmp tMagE (mag(E)); const volScalarField &magE = tMagE(); diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H index 081bd630..9f23be3b 100644 --- a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H +++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H @@ -1,3 +1,4 @@ +/* tmp > mvConvection ( fv::convectionScheme::New @@ -8,6 +9,7 @@ tmp > mvConvection mesh.divScheme("div(phi,Yi_h)") ) ); +*/ { label inertIndex = -1; @@ -35,9 +37,11 @@ tmp > mvConvection // Adding drift flux to boundary patches forAll (bfIonFlux, pidx) { + Info << "Adding drift flux to boundary patches" << pidx << endl; + bfIonFlux[pidx] += - (E.boundaryField()[pidx] - & sf.boundaryField()[pidx]) + snE.boundaryField()[pidx] + * msf.boundaryField()[pidx] * rho.boundaryField()[pidx] * Di.boundaryField()[pidx] / T.boundaryField()[pidx] @@ -113,22 +117,22 @@ tmp > mvConvection if (nCharge != 0) { - phi_drift = phi; - phi_drift += fvc::interpolate((rho*Di/T*(eCharge*z/kB))*E) & mesh.Sf(); + phis[i] = phi; + phis[i] += fvc::interpolate((rho*Di/T*(eCharge*z/kB))*E) & mesh.Sf(); } if (ions.contains(Y[i].name())) { const label ibc = ions[Y[i].name()]; - // phi_drift updated - phi_drift.boundaryField() = ionFluxBFs[ibc]; + // phis[i] updated + phis[i].boundaryField() = ionFluxBFs[ibc]; } else if (neutrals.contains(Y[i].name())) { const label ibc = neutrals[Y[i].name()]; - // update phi_neutral - phi_neutral.internalField() = phi.internalField(); - phi_neutral.boundaryField() = neutralFluxBFs[ibc]; + // update phis[i] + phis[i].internalField() = phi.internalField(); + phis[i].boundaryField() = neutralFluxBFs[ibc]; } @@ -137,9 +141,9 @@ tmp > mvConvection fvm::ddt(rho, Yi) + ( nCharge != 0 - ? fvm::div(phi_drift, Yi, "div(phi,Yi_h)") + ? fvm::div(phis[i], Yi, "div(phi,Yi_h)") : ( neutrals.contains(Y[i].name()) - ? fvm::div(phi_neutral, Yi, "div(phi,Yi_h)") + ? fvm::div(phis[i], Yi, "div(phi,Yi_h)") : fvm::div(phi, Yi, "div(phi,Yi_h)") ) ) diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H index 439e1c34..d48e25cc 100644 --- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H +++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H @@ -92,6 +92,20 @@ volVectorField E -fvc::grad(Phi) ); +surfaceScalarField snE +( + IOobject + ( + "snE", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + -fvc::snGrad(Phi) +); + + Info<< "Creating reaction model\n" << endl; autoPtr reaction @@ -363,8 +377,8 @@ surfaceScalarField ve "ve", runTime.timeName(), mesh, - IOobject::NO_READ, - IOobject::NO_WRITE + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE ), linearInterpolate(Udrift+Uthermal) & mesh.Sf() ); @@ -377,7 +391,7 @@ surfaceScalarField phi_drift "phi_drift", runTime.timeName(), mesh, - IOobject::NO_READ, + IOobject::READ_IF_PRESENT, IOobject::NO_WRITE ), phi @@ -391,13 +405,42 @@ surfaceScalarField phi_neutral "phi_neutral", runTime.timeName(), mesh, - IOobject::NO_READ, + IOobject::READ_IF_PRESENT, IOobject::NO_WRITE ), phi ); + +PtrList phis (composition.species().size()); + +forAll (composition.species(), isp) +{ + const scalar z(composition.z(isp)); + const label nCharge(z); + + if (nCharge != 0) + { + phis.set(isp, + new surfaceScalarField + ( + IOobject + ( + "phi." + composition.species()[isp], + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + phi + ) + ); + } +} + + // plasmaWallFluxes +Info<< "Reading plasma wall flux bc control\n" << endl; // electron wall flux dictionary wallElectronFlux @@ -418,6 +461,8 @@ forAll (wallPatcheNames, pi) = mesh.boundaryMesh().findPatchID(patchName); } +Info<< "plasma walls are \n" << wallPatcheNames << endl; + // ion wall flux dictionary wallIonFluxes @@ -432,6 +477,8 @@ PtrList targetList (ions.size()); PtrList reflexes (ions.size()); PtrList ionFluxBFs (ions.size()); +Info<< ions.size() << " ions are \n" << ions << endl; + label nTargets = 0; forAll (ions, iidx) @@ -485,7 +532,6 @@ for (label nidx = 0; nidx < nTargets; nidx++) neutrals.append(n); nNeutrals += 1; } - } PtrList neutralFluxBFs (neutrals.size()); diff --git a/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C b/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C index ed1062c6..a85b41c4 100644 --- a/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C +++ b/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C @@ -72,13 +72,14 @@ int main(int argc, char *argv[]) runTime++; Info<< "Time = " << runTime.timeName() << nl << endl; - #include "numberDensity.H" - #include "PhiEqn.H" #include "rhoEqn.H" while (pimple.loop()) { + #include "numberDensity.H" + #include "PhiEqn.H" + #include "UEqn.H" reaction->correct();