Flame structure evolution using OpenFOAM fvm

This commit is contained in:
ignis 2017-08-29 18:29:50 +09:00
parent d1fcf79dee
commit 7d8700ab2e
5 changed files with 519 additions and 0 deletions

View file

@ -6,5 +6,6 @@ targetType=libso
wmake $targetType chemistryModel_POSTECH
wmake $targetType combustionModels_POSTECH
wmake $targetType thermos
#------------------------------------------------------------------------------

View file

@ -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)

View file

@ -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 <http://www.gnu.org/licenses/>.
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<labelField> 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<LinearMesh> 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
// ************************************************************************* //

View file

@ -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();
}

View file

@ -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<Field<vector>>(CMCPoints),
CMCFaces.xfer(),
CMCOwners.xfer(),
CMCNeighbors.xfer()
);
if (mfMesh.boundaryMesh().empty())
{
List<polyPatch*> 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<PtrList<volScalarField> > FStruct (nf);
for(label k = 0 ; k < nf ; k++)
{
FStruct.set (k, new PtrList<volScalarField> (Ysize+3));
PtrList<volScalarField> &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));
}
}