From 7d8700ab2ea72008cc93cc13909e954579ab721f Mon Sep 17 00:00:00 2001 From: ignis Date: Tue, 29 Aug 2017 18:29:50 +0900 Subject: [PATCH] Flame structure evolution using OpenFOAM fvm --- libs/Allwmake | 1 + .../LagrangianCMCFoam/LagrangianCMCFoam.C | 7 + solvers_post/LagrangianCMCFoam/LinearMesh.H | 337 ++++++++++++++++++ solvers_post/LagrangianCMCFoam/QiEqn.H | 29 ++ .../LagrangianCMCFoam/createLinearMesh.H | 145 ++++++++ 5 files changed, 519 insertions(+) create mode 100644 solvers_post/LagrangianCMCFoam/LinearMesh.H create mode 100644 solvers_post/LagrangianCMCFoam/QiEqn.H create mode 100644 solvers_post/LagrangianCMCFoam/createLinearMesh.H diff --git a/libs/Allwmake b/libs/Allwmake index a494487..16f1269 100755 --- a/libs/Allwmake +++ b/libs/Allwmake @@ -6,5 +6,6 @@ targetType=libso wmake $targetType chemistryModel_POSTECH wmake $targetType combustionModels_POSTECH +wmake $targetType thermos #------------------------------------------------------------------------------ diff --git a/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C b/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C index 8ed9abd..8a49cbc 100644 --- a/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C +++ b/solvers_post/LagrangianCMCFoam/LagrangianCMCFoam.C @@ -48,6 +48,9 @@ Contact #include "IFstream.H" #include "OFstream.H" #include "Math.H" //Mathematical functions for CMC (AMC, gammaln, TDMA) + +#include "LinearMesh.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) @@ -73,6 +76,8 @@ int main(int argc, char *argv[]) #include "startSummary.H" //Make logSummary file for CMC #include "readCMCProperties.H" //Read and set fields for CMC calculation + #include "createLinearMesh.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nStarting time loop\n" << endl; @@ -87,6 +92,8 @@ int main(int argc, char *argv[]) Info<< "Time = " << runTime.timeName() << nl << endl; + #include "QiEqn.H" + #include "rhoEqn.H" #include "UEqn.H" if(init_start_CMC == true) diff --git a/solvers_post/LagrangianCMCFoam/LinearMesh.H b/solvers_post/LagrangianCMCFoam/LinearMesh.H new file mode 100644 index 0000000..862eae3 --- /dev/null +++ b/solvers_post/LagrangianCMCFoam/LinearMesh.H @@ -0,0 +1,337 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +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::LinearMesh + +Description + +SourceFiles + LinearMeshI.H + LinearMesh.C + LinearMeshIO.C + +\*---------------------------------------------------------------------------*/ + +#ifndef LinearMesh_H +#define LinearMesh_H + +#include "fvMesh.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward declaration of classes +class Istream; +class Ostream; + +// Forward declaration of friend functions and operators +class LinearMesh; +Istream& operator>>(Istream&, LinearMesh&); +Ostream& operator<<(Ostream&, const LinearMesh&); + + +/*---------------------------------------------------------------------------*\ + Class LinearMesh Declaration +\*---------------------------------------------------------------------------*/ + +class LinearMesh +{ + // Private data + + //- Lower end value of 1-D domain + scalar x0_; + + //- Upper end value of 1-D domain + scalar x1_; + + //- Total number of cells + label nCells_; + + //- Total number of cells + pointField points_; + + //- Total number of cells + faceList faces_; + + //- Total number of cells + labelList owner_; + + //- Total number of cells + labelList neighbour_; + + //- Total number of cells + label nPatches_; + + //- Total number of cells + wordList patchNames_; + + //- Total number of cells + labelList patchStartIndices_; + + //- Total number of cells + labelList patchSizes_; + + + // Private Member Functions + + //- Disallow Construct null + LinearMesh(); + + //- Disallow default bitwise copy construct + LinearMesh(const LinearMesh&); + + //- Disallow default bitwise assignment + void operator=(const LinearMesh&); + + //- Generate points and faces + void init(const scalarField &etaValue) + { + for(label i = 0; i < (nCells_+1); i++) + { + points_[4*i+0] = vector(etaValue[i], 0.0, 0.0); + points_[4*i+1] = vector(etaValue[i], 0.1, 0.0); + points_[4*i+2] = vector(etaValue[i], 0.1, 0.1); + points_[4*i+3] = vector(etaValue[i], 0.0, 0.1); + } + + const label nBoundaryFaces = 2; + const label nInternalFaces = nCells_-1; + const label nEmptyFaces = 4*nCells_; + + for(label i = 0; i < nInternalFaces; i++) + { + IStringStream faceFormatStream ("4(4 5 6 7)"); + labelField faceFormat (faceFormatStream); + faces_[i] = face(faceFormat + 4*i); + } + + { + IStringStream faceFormatStream ("4(0 3 2 1)"); + faces_[nInternalFaces] = face(faceFormatStream); + } + + { + IStringStream faceFormatStream ("4(4 5 6 7)"); + labelField faceFormat (faceFormatStream); + faces_[nInternalFaces+1] = face(faceFormat + 4*nInternalFaces); + } + + for(label i = 0; i < nCells_; i++) + { + IStringStream faceFormatStream ("4(4(0 1 5 4)4(1 2 6 5)4(2 3 7 6)4(3 0 4 7))"); + List faceFormats (faceFormatStream); + + label fi = nInternalFaces + nBoundaryFaces + 4*i; + faces_[fi+0] = face(faceFormats[0] + 4*i); + faces_[fi+1] = face(faceFormats[1] + 4*i); + faces_[fi+2] = face(faceFormats[2] + 4*i); + faces_[fi+3] = face(faceFormats[3] + 4*i); + } + + + for(label i = 0; i < nInternalFaces; i++) + { + owner_[i] = i; + neighbour_[i] = i+1; + } + + { + owner_[nInternalFaces] = 0; + } + + { + owner_[nInternalFaces+1] = nCells_ - 1; + } + + for(label i = 0; i < nCells_; i++) + { + label fi = nInternalFaces + nBoundaryFaces + 4*i; + owner_[fi+0] = i; + owner_[fi+1] = i; + owner_[fi+2] = i; + owner_[fi+3] = i; + } + + + nPatches_ = 3; + + IStringStream sInd + ( + word("3(") + + name(nInternalFaces) + " " + + name(nInternalFaces+1) + " " + + name(nInternalFaces+2) + " " + + word(")") + ); + patchStartIndices_ = labelList(sInd); + + IStringStream sSize(word("3(1 1 ")+name(4*nCells_)+word(")")); + patchSizes_ = labelList(sSize); + + IStringStream sName("3(lowerEnd upperEnd Sides)"); + patchNames_ = wordList(sName); + } + +public: + + // Static data members + + //- Static data staticData + // static const dataType staticData; + + + // Constructors + + //- Construct from components + LinearMesh(const scalar x0, const scalar x1, const label nCells) + : + x0_(x0), + x1_(x1), + nCells_(nCells), + points_(4*(nCells+1)), + faces_(5*nCells+1), + owner_(faces_.size()), + neighbour_(nCells-1) + { + + scalarField etaValue(nCells+1, 0.0); + scalar delta = (x1 - x0) / scalar(nCells); + forAll(etaValue, i) + { + etaValue[i] = x0 + delta*i; + } + + + init(etaValue); + + } + + //- Construct from components + LinearMesh(const scalarField &x) + : + x0_(x[0]), + x1_(x.last()), + nCells_(x.size()-1), + points_(4*(nCells_+1)), + faces_(5*nCells_+1), + owner_(faces_.size()), + neighbour_(nCells_-1) + { + + init(x); + + } + + //- Construct from Istream + LinearMesh(Istream&); + + //- Construct as copy + // LinearMesh(const LinearMesh&); + + + // Selectors + + //- Select null constructed + // static autoPtr New(); + + + //- Destructor + ~LinearMesh() {} + + + // Member Functions + + // Access + //- Lower end value of 1-D domain + scalar x0() const {return x0_;}; + + //- Upper end value of 1-D domain + scalar x1() const {return x1_;}; + + //- Total number of cells + label nCells() const {return nCells_;}; + + //- Total number of cells + const pointField &points() const {return points_;}; + + //- Total number of cells + const faceList &faces() const {return faces_;}; + + //- Total number of cells + const labelList &owner() const {return owner_;}; + + //- Total number of cells + const labelList &neighbour() const {return neighbour_;}; + + //- Total number of cells + label nPatches () const {return nPatches_;}; + + //- Total number of cells + word patchName (const label i) const {return patchNames_[i];}; + + //- Total number of cells + label patchStartIndex (const label i) const {return patchStartIndices_[i];}; + + //- Total number of cells + label patchSize (const label i) const {return patchSizes_[i];}; + + + // Check + + // Edit + + // Write + + + // Member Operators + + // void operator=(const LinearMesh&); + + + // Friend Functions + + // Friend Operators + + // IOstream Operators + + friend Istream& operator>>(Istream&, LinearMesh&); + friend Ostream& operator<<(Ostream&, const LinearMesh&); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// #include "LinearMeshI.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/solvers_post/LagrangianCMCFoam/QiEqn.H b/solvers_post/LagrangianCMCFoam/QiEqn.H new file mode 100644 index 0000000..c8e84c8 --- /dev/null +++ b/solvers_post/LagrangianCMCFoam/QiEqn.H @@ -0,0 +1,29 @@ +{ + fvScalarMatrix testSource(testEtaMesh, testEtaMesh.dimensions()*dimVol/dimTime); + + testSource.source() = 1.0; + + fvScalarMatrix QiEqn + ( + fvm::ddt(testEtaMesh) +// + mvConvection->fvmDiv(phi, he) +// + fvc::ddt(rho, K) + fvc::div(phi, K) +// + ( +// he.name() == "ea" +// ? fvc::div +// ( +// fvc::absolute(phi/fvc::interpolate(rho), U), +// p, +// "div(phiv,p)" +// ) +// : -dpdt +// ) + - fvm::laplacian(DQi, testEtaMesh) + == + testSource + ); + + QiEqn.relax(); + + QiEqn.solve(); +} diff --git a/solvers_post/LagrangianCMCFoam/createLinearMesh.H b/solvers_post/LagrangianCMCFoam/createLinearMesh.H new file mode 100644 index 0000000..7749bcf --- /dev/null +++ b/solvers_post/LagrangianCMCFoam/createLinearMesh.H @@ -0,0 +1,145 @@ +Info<< "Constructing single cell mesh" << nl << endl; + +LinearMesh lMesh(0, 1, 100); + +pointField CMCPoints(lMesh.points()); +faceList CMCFaces(lMesh.faces()); +labelList CMCOwners(lMesh.owner()); +labelList CMCNeighbors(lMesh.neighbour()); + +fvMesh mfMesh +( + IOobject + ( + "mfSpace", + runTime.timeName(), + runTime, + IOobject::READ_IF_PRESENT + ), + xferMove>(CMCPoints), + CMCFaces.xfer(), + CMCOwners.xfer(), + CMCNeighbors.xfer() +); + +if (mfMesh.boundaryMesh().empty()) +{ + List mfPatches(3); + + mfPatches[0] = new polyPatch + ( + lMesh.patchName(0), + lMesh.patchSize(0), + lMesh.patchStartIndex(0), + 0, + mfMesh.boundaryMesh(), + polyPatch::typeName + ); + + mfPatches[1] = new polyPatch + ( + lMesh.patchName(1), + lMesh.patchSize(1), + lMesh.patchStartIndex(1), + 1, + mfMesh.boundaryMesh(), + polyPatch::typeName + ); + + mfPatches[2] = new emptyPolyPatch + ( + lMesh.patchName(2), + lMesh.patchSize(2), + lMesh.patchStartIndex(2), + 2, + mfMesh.boundaryMesh(), + emptyPolyPatch::typeName + ); + + mfMesh.addFvPatches(mfPatches); +} + +IStringStream tempStream("3(fixedValue fixedValue empty)"); +wordList testEtaPatchTypeNames(tempStream); + +volScalarField testEtaMesh +( + IOobject + ( + "testEtaMesh", + runTime.timeName(), + mfMesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mfMesh, + dimensionedScalar("testEtaMesh", dimless, 1.0), + testEtaPatchTypeNames +); + +const volVectorField &cEtaMesh(mfMesh.C()); +forAll(cEtaMesh, celli) +{ + const scalar x = cEtaMesh[celli].x(); + if (x > 0.5) + { + testEtaMesh[celli] = 1.0; + } + else + { + testEtaMesh[celli] = 0.0; + } + +} + +volScalarField DQi +( + IOobject + ( + "DQi", + runTime.timeName(), + mfMesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mfMesh, + dimensionedScalar("DQi", dimArea/dimTime, 0.01), + testEtaPatchTypeNames +); + +PtrList > FStruct (nf); +for(label k = 0 ; k < nf ; k++) +{ + FStruct.set (k, new PtrList (Ysize+3)); + + PtrList &Fk(FStruct[k]); + Info << Fk.size() << endl; + for (label i = 0; i < Ysize + 3; i++) + { + Fk.set (i, + new volScalarField + ( + IOobject + ( + "test", + runTime.timeName(), + mfMesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mfMesh, + dimensionedScalar("test", dimless, 0.0), + testEtaPatchTypeNames + ) + ); + + } +} + +for(label k = 0 ; k < nf ; k++) +{ + for (label i = 0; i < Ysize + 3; i++) + { + FStruct[k][i].rename("test."+name(k)+"."+name(i)); + } +}