New solver: adjointShapeOptimizationFoam

Uses adjoint method to block regions of the domain causing total pressure loss,
e.g. recirculation zones.

pitzDaily tutorial case supplied.
This commit is contained in:
henry
2010-03-15 17:08:45 +00:00
parent d83d28733d
commit d35773f4fb
24 changed files with 1798 additions and 0 deletions

View File

@ -0,0 +1,5 @@
adjointOutletPressure/adjointOutletPressureFvPatchScalarField.C
adjointOutletVelocity/adjointOutletVelocityFvPatchVectorField.C
adjointShapeOptimizationFoam.C
EXE = $(FOAM_APPBIN)/adjointShapeOptimizationFoam

View File

@ -0,0 +1,11 @@
EXE_INC = \
-I$(LIB_SRC)/turbulenceModels \
-I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-I$(LIB_SRC)/finiteVolume/lnInclude
EXE_LIBS = \
-lincompressibleRASModels \
-lincompressibleTransportModels \
-lfiniteVolume

View File

@ -0,0 +1,47 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Global
continuityErrs
Description
Calculates and prints the continuity errors.
\*---------------------------------------------------------------------------*/
{
scalar sumLocalContErr = runTime.deltaT().value()*
mag(fvc::div(phia))().weightedAverage(mesh.V()).value();
scalar globalContErr = runTime.deltaT().value()*
fvc::div(phia)().weightedAverage(mesh.V()).value();
cumulativeContErr += globalContErr;
Info<< "Adjoint continuity errors : sum local = " << sumLocalContErr
<< ", global = " << globalContErr
<< ", cumulative = " << cumulativeContErr
<< endl;
}
// ************************************************************************* //

View File

@ -0,0 +1,132 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "adjointOutletPressureFvPatchScalarField.H"
#include "addToRunTimeSelectionTable.H"
#include "fvPatchMapper.H"
#include "volFields.H"
#include "surfaceFields.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::adjointOutletPressureFvPatchScalarField::
adjointOutletPressureFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchScalarField(p, iF)
{}
Foam::adjointOutletPressureFvPatchScalarField::
adjointOutletPressureFvPatchScalarField
(
const adjointOutletPressureFvPatchScalarField& ptf,
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchScalarField(ptf, p, iF, mapper)
{}
Foam::adjointOutletPressureFvPatchScalarField::
adjointOutletPressureFvPatchScalarField
(
const fvPatch& p,
const DimensionedField<scalar, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchScalarField(p, iF)
{
fvPatchField<scalar>::operator=
(
scalarField("value", dict, p.size())
);
}
Foam::adjointOutletPressureFvPatchScalarField::
adjointOutletPressureFvPatchScalarField
(
const adjointOutletPressureFvPatchScalarField& tppsf,
const DimensionedField<scalar, volMesh>& iF
)
:
fixedValueFvPatchScalarField(tppsf, iF)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::adjointOutletPressureFvPatchScalarField::updateCoeffs()
{
if (updated())
{
return;
}
const fvsPatchField<scalar>& phip =
patch().lookupPatchField<surfaceScalarField, scalar>("phi");
const fvsPatchField<scalar>& phiap =
patch().lookupPatchField<surfaceScalarField, scalar>("phia");
const fvPatchField<vector>& Up =
patch().lookupPatchField<volVectorField, vector>("U");
const fvPatchField<vector>& Uap =
patch().lookupPatchField<volVectorField, vector>("Ua");
operator==((phiap/patch().magSf() - 1.0)*phip/patch().magSf() + (Up & Uap));
fixedValueFvPatchScalarField::updateCoeffs();
}
void Foam::adjointOutletPressureFvPatchScalarField::write(Ostream& os) const
{
fvPatchScalarField::write(os);
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
makePatchTypeField
(
fvPatchScalarField,
adjointOutletPressureFvPatchScalarField
);
}
// ************************************************************************* //

View File

@ -0,0 +1,137 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
adjointOutletPressureFvPatchScalarField
Description
SourceFiles
adjointOutletPressureFvPatchScalarField.C
\*---------------------------------------------------------------------------*/
#ifndef adjointOutletPressureFvPatchScalarField_H
#define adjointOutletPressureFvPatchScalarField_H
#include "fixedValueFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class adjointOutletPressureFvPatch Declaration
\*---------------------------------------------------------------------------*/
class adjointOutletPressureFvPatchScalarField
:
public fixedValueFvPatchScalarField
{
public:
//- Runtime type information
TypeName("adjointOutletPressure");
// Constructors
//- Construct from patch and internal field
adjointOutletPressureFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&
);
//- Construct from patch, internal field and dictionary
adjointOutletPressureFvPatchScalarField
(
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const dictionary&
);
//- Construct by mapping given adjointOutletPressureFvPatchScalarField
// onto a new patch
adjointOutletPressureFvPatchScalarField
(
const adjointOutletPressureFvPatchScalarField&,
const fvPatch&,
const DimensionedField<scalar, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct and return a clone
virtual tmp<fvPatchScalarField> clone() const
{
return tmp<fvPatchScalarField>
(
new adjointOutletPressureFvPatchScalarField(*this)
);
}
//- Construct as copy setting internal field reference
adjointOutletPressureFvPatchScalarField
(
const adjointOutletPressureFvPatchScalarField&,
const DimensionedField<scalar, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchScalarField> clone
(
const DimensionedField<scalar, volMesh>& iF
) const
{
return tmp<fvPatchScalarField>
(
new adjointOutletPressureFvPatchScalarField(*this, iF)
);
}
// Member functions
// Evaluation functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,131 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "adjointOutletVelocityFvPatchVectorField.H"
#include "volFields.H"
#include "addToRunTimeSelectionTable.H"
#include "surfaceFields.H"
#include "fvPatchFieldMapper.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::adjointOutletVelocityFvPatchVectorField::
adjointOutletVelocityFvPatchVectorField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchVectorField(p, iF)
{}
Foam::adjointOutletVelocityFvPatchVectorField::
adjointOutletVelocityFvPatchVectorField
(
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const dictionary& dict
)
:
fixedValueFvPatchVectorField(p, iF)
{
fvPatchVectorField::operator=(vectorField("value", dict, p.size()));
}
Foam::adjointOutletVelocityFvPatchVectorField::
adjointOutletVelocityFvPatchVectorField
(
const adjointOutletVelocityFvPatchVectorField& ptf,
const fvPatch& p,
const DimensionedField<vector, volMesh>& iF,
const fvPatchFieldMapper& mapper
)
:
fixedValueFvPatchVectorField(ptf, p, iF, mapper)
{}
Foam::adjointOutletVelocityFvPatchVectorField::
adjointOutletVelocityFvPatchVectorField
(
const adjointOutletVelocityFvPatchVectorField& pivpvf,
const DimensionedField<vector, volMesh>& iF
)
:
fixedValueFvPatchVectorField(pivpvf, iF)
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Update the coefficients associated with the patch field
void Foam::adjointOutletVelocityFvPatchVectorField::updateCoeffs()
{
if (updated())
{
return;
}
const fvsPatchField<scalar>& phiap =
patch().lookupPatchField<surfaceScalarField, scalar>("phia");
const fvPatchField<vector>& Up =
patch().lookupPatchField<volVectorField, vector>("U");
scalarField Un = mag(patch().nf() & Up);
vectorField UtHat = (Up - patch().nf()*Un)/(Un + SMALL);
vectorField Uan = patch().nf()*(patch().nf() & patchInternalField());
vectorField::operator=(phiap*patch().Sf()/sqr(patch().magSf()) + UtHat);
//vectorField::operator=(Uan + UtHat);
fixedValueFvPatchVectorField::updateCoeffs();
}
void Foam::adjointOutletVelocityFvPatchVectorField::write(Ostream& os) const
{
fvPatchVectorField::write(os);
writeEntry("value", os);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
makePatchTypeField
(
fvPatchVectorField,
adjointOutletVelocityFvPatchVectorField
);
}
// ************************************************************************* //

View File

@ -0,0 +1,131 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2007 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Class
adjointOutletVelocityFvPatchVectorField
Description
SourceFiles
adjointOutletVelocityFvPatchVectorField.C
\*---------------------------------------------------------------------------*/
#ifndef adjointOutletVelocityFvPatchVectorField_H
#define adjointOutletVelocityFvPatchVectorField_H
#include "fixedValueFvPatchFields.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class adjointOutletVelocityFvPatch Declaration
\*---------------------------------------------------------------------------*/
class adjointOutletVelocityFvPatchVectorField
:
public fixedValueFvPatchVectorField
{
public:
//- Runtime type information
TypeName("adjointOutletVelocity");
// Constructors
//- Construct from patch and internal field
adjointOutletVelocityFvPatchVectorField
(
const fvPatch&,
const DimensionedField<vector, volMesh>&
);
//- Construct from patch, internal field and dictionary
adjointOutletVelocityFvPatchVectorField
(
const fvPatch&,
const DimensionedField<vector, volMesh>&,
const dictionary&
);
//- Construct by mapping given adjointOutletVelocityFvPatchVectorField
// onto a new patch
adjointOutletVelocityFvPatchVectorField
(
const adjointOutletVelocityFvPatchVectorField&,
const fvPatch&,
const DimensionedField<vector, volMesh>&,
const fvPatchFieldMapper&
);
//- Construct and return a clone
virtual tmp<fvPatchVectorField> clone() const
{
return tmp<fvPatchVectorField>
(
new adjointOutletVelocityFvPatchVectorField(*this)
);
}
//- Construct as copy setting internal field reference
adjointOutletVelocityFvPatchVectorField
(
const adjointOutletVelocityFvPatchVectorField&,
const DimensionedField<vector, volMesh>&
);
//- Construct and return a clone setting internal field reference
virtual tmp<fvPatchVectorField> clone(const DimensionedField<vector, volMesh>& iF) const
{
return tmp<fvPatchVectorField>
(
new adjointOutletVelocityFvPatchVectorField(*this, iF)
);
}
// Member functions
//- Update the coefficients associated with the patch field
virtual void updateCoeffs();
//- Write
virtual void write(Ostream&) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,226 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
ajointShapeOptimizationFoam
Description
Steady-state solver for incompressible, turbulent flow of non-Newtonian
fluids with optimisation of duct shape by applying "blockage" in regions
causing pressure loss as estimated using an adjoint formulation.
References:
@verbatim
"Implementation of a continuous adjoint for topology optimization of
ducted flows"
C. Othmer,
E. de Villiers,
H.G. Weller
AIAA-2007-3947
http://pdf.aiaa.org/preview/CDReadyMCFD07_1379/PV2007_3947.pdf
@endverbatim
Note that this solver optimises for total pressure loss whereas the
above paper describes the method for optimising power-loss.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "RASModel.H"
template<class Type>
void zeroCells
(
GeometricField<Type, fvPatchField, volMesh>& vf,
const labelList& cells
)
{
forAll(cells, i)
{
vf[cells[i]] = pTraits<Type>::zero;
}
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "initContinuityErrs.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "\nStarting time loop\n" << endl;
for (runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "readSIMPLEControls.H"
p.storePrevIter();
laminarTransport.lookup("lambda") >> lambda;
//alpha +=
// mesh.relaxationFactor("alpha")
// *(lambda*max(Ua & U, zeroSensitivity) - alpha);
alpha +=
mesh.relaxationFactor("alpha")
*(min(max(alpha + lambda*(Ua & U), zeroAlpha), alphaMax) - alpha);
zeroCells(alpha, inletCells);
//zeroCells(alpha, outletCells);
// Pressure-velocity SIMPLE corrector
{
// Momentum predictor
tmp<fvVectorMatrix> UEqn
(
fvm::div(phi, U)
+ turbulence->divDevReff(U)
+ fvm::Sp(alpha, U)
);
UEqn().relax();
solve(UEqn() == -fvc::grad(p));
p.boundaryField().updateCoeffs();
volScalarField rAU = 1.0/UEqn().A();
U = rAU*UEqn().H();
UEqn.clear();
phi = fvc::interpolate(U) & mesh.Sf();
adjustPhi(phi, U, p);
// Non-orthogonal pressure corrector loop
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phi)
);
pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}
#include "continuityErrs.H"
// Explicitly relax pressure for momentum corrector
p.relax();
// Momentum corrector
U -= rAU*fvc::grad(p);
U.correctBoundaryConditions();
}
pa.storePrevIter();
// Adjoint Pressure-velocity SIMPLE corrector
{
// Adjoint Momentum predictor
volVectorField adjointTransposeConvection = (fvc::grad(Ua) & U);
//volVectorField adjointTransposeConvection = fvc::reconstruct
//(
// mesh.magSf()*(fvc::snGrad(Ua) & fvc::interpolate(U))
//);
zeroCells(adjointTransposeConvection, inletCells);
tmp<fvVectorMatrix> UaEqn
(
fvm::div(-phi, Ua)
- adjointTransposeConvection
+ turbulence->divDevReff(Ua)
+ fvm::Sp(alpha, Ua)
);
UaEqn().relax();
solve(UaEqn() == -fvc::grad(pa));
pa.boundaryField().updateCoeffs();
volScalarField rAUa = 1.0/UaEqn().A();
Ua = rAUa*UaEqn().H();
UaEqn.clear();
phia = fvc::interpolate(Ua) & mesh.Sf();
adjustPhi(phia, Ua, pa);
// Non-orthogonal pressure corrector loop
for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix paEqn
(
fvm::laplacian(rAUa, pa) == fvc::div(phia)
);
paEqn.setReference(paRefCell, paRefValue);
paEqn.solve();
if (nonOrth == nNonOrthCorr)
{
phia -= paEqn.flux();
}
}
#include "adjointContinuityErrs.H"
// Explicitly relax pressure for adjoint momentum corrector
pa.relax();
// Adjoint momentum corrector
Ua -= rAUa*fvc::grad(pa);
Ua.correctBoundaryConditions();
}
turbulence->correct();
runTime.write();
Info<< "ExecutionTime = "
<< runTime.elapsedCpuTime()
<< " s\n\n" << endl;
}
Info<< "End\n" << endl;
return(0);
}
// ************************************************************************* //

View File

@ -0,0 +1,105 @@
Info << "Reading field p\n" << endl;
volScalarField p
(
IOobject
(
"p",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info << "Reading field U\n" << endl;
volVectorField U
(
IOobject
(
"U",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhi.H"
label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("SIMPLE"), pRefCell, pRefValue);
Info << "Reading field pa\n" << endl;
volScalarField pa
(
IOobject
(
"pa",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
Info << "Reading field Ua\n" << endl;
volVectorField Ua
(
IOobject
(
"Ua",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);
#include "createPhia.H"
label paRefCell = 0;
scalar paRefValue = 0.0;
setRefCell(pa, mesh.solutionDict().subDict("SIMPLE"), paRefCell, paRefValue);
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::RASModel> turbulence
(
incompressible::RASModel::New(U, phi, laminarTransport)
);
dimensionedScalar zeroSensitivity("0", dimVelocity*dimVelocity, 0.0);
dimensionedScalar zeroAlpha("0", dimless/dimTime, 0.0);
dimensionedScalar lambda(laminarTransport.lookup("lambda"));
dimensionedScalar alphaMax(laminarTransport.lookup("alphaMax"));
const labelList& inletCells =
mesh.boundary()[mesh.boundaryMesh().findPatchID("inlet")].faceCells();
//const labelList& outletCells =
// mesh.boundary()[mesh.boundaryMesh().findPatchID("outlet")].faceCells();
volScalarField alpha
(
IOobject
(
"alpha",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
lambda*max(Ua & U, zeroSensitivity)
);
zeroCells(alpha, inletCells);
//zeroCells(alpha, outletCells);

View File

@ -0,0 +1,57 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
\\/ 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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Global
createPhia
Description
Creates and initialises the face-flux field phia.
\*---------------------------------------------------------------------------*/
#ifndef createPhia_H
#define createPhia_H
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< "Reading/calculating face flux field phia\n" << endl;
surfaceScalarField phia
(
IOobject
(
"phia",
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(Ua) & mesh.Sf()
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //