adaptativeSimpleFoam
Dans cet article, nous proposons un solveur qui s’appuie sur l’estimateur d’erreur décrit dans cet article .
L’idée est de proposer un solveur de type SIMPLE qui, à chaque pas de temps, calcule un estimateur d’erreur et raffine ou dé-raffine le maillage selon l’importance de l’erreur liée à la discrétisation.
- files
adaptativeSimpleFoam.C
EXE = $(FOAM_APPBIN)/adaptativeSimpleFoam
Options
Voir code
EXE_INC = \ -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ -I$(FOAM_DEV)/ressources/errorEstimation/lnInclude \ -I$(LIB_SRC)/transportModels \ -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \ -I$(LIB_SRC)/dynamicFvMesh/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/sampling/lnInclude \ -fpermissive EXE_LIBS = \ -lturbulenceModels \ -lerrorEstimation \ -lincompressibleTurbulenceModels \ -lincompressibleTransportModels \ -ldynamicFvMesh \ -ltopoChangerFvMesh \ -ldynamicMesh \ -lfiniteVolume \ -lmeshTools \ -lfvOptions \ -lsamplingcreateFields.HVoir code
/* * createFileds.h * * Copyright 2018 arep <arep@debian01> * * This program 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. * * This program 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 this program; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, * MA 02110-1301, USA. */ 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 ); Info<< "Reading transportProperties\n" << endl; IOdictionary transportProperties ( IOobject ( "transportProperties", runTime.constant(), mesh, IOobject::MUST_READ_IF_MODIFIED, IOobject::NO_WRITE ) ); dimensionedScalar nu ( "nu", dimViscosity, transportProperties.lookup("nu") ); Info<< "Creating field error\n" << endl; volScalarField error ( IOobject ( "error", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), mesh, dimensionedScalar("error", dimensionSet(0,1,-1,0,0,0,0), Foam::scalar(0)) ); #include "createPhi.H" label pRefCell = 0; scalar pRefValue = 0.0; setRefCell(p, simple.dict(), pRefCell, pRefValue); mesh.setFluxRequired(p.name()); singlePhaseTransportModel laminarTransport(U, phi); autoPtr<incompressible::turbulenceModel> turbulence ( incompressible::turbulenceModel::New(U, phi, laminarTransport) ); #include "createMRF.H"<br>UEqn.HVoir code
// Momentum predictor MRF.correctBoundaryVelocity(U); tmp<fvVectorMatrix> tUEqn ( fvm::div(phi, U) + MRF.DDt(U) + turbulence->divDevReff(U) == fvOptions(U) ); fvVectorMatrix& UEqn = tUEqn.ref(); UEqn.relax(); fvOptions.constrain(UEqn); if (simple.momentumPredictor()) { solve(UEqn == -fvc::grad(p)); fvOptions.correct(U); }pEqn.HVoir code
{ volScalarField rAU(1.0/UEqn.A()); volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA)); MRF.makeRelative(phiHbyA); adjustPhi(phiHbyA, U, p); tmp<volScalarField> rAtU(rAU); if (simple.consistent()) { rAtU = 1.0/(1.0/rAU - UEqn.H1()); phiHbyA += fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf(); HbyA -= (rAU - rAtU())*fvc::grad(p); } tUEqn.clear(); // Update the pressure BCs to ensure flux consistency constrainPressure(p, U, phiHbyA, rAtU(), MRF); // Non-orthogonal pressure corrector loop while (simple.correctNonOrthogonal()) { fvScalarMatrix pEqn ( fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA) ); pEqn.setReference(pRefCell, pRefValue); pEqn.solve(); if (simple.finalNonOrthogonalIter()) { phi = phiHbyA - pEqn.flux(); } } #include "continuityErrs.H" // Explicitly relax pressure for momentum corrector p.relax(); // Momentum corrector U = HbyA - rAtU()*fvc::grad(p); U.correctBoundaryConditions(); fvOptions.correct(U); }eEqn.HVoir code
/* * eEqn * * Copyright 2018 arep <arep@debian01> * * This program 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. * * This program 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 this program; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, * MA 02110-1301, USA. * * */ errorEstimate<vector> eEqn ( resError::div(phi, U) - resError::laplacian(nu, U) == -fvc::grad(p) ); error=mag(eEqn.error());adaptativeSimpleFoam.CVoir code
/*---------------------------------------------------------------------------*\ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License * Copyright 2018 arep <arep@debian01> * * This program 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. * * This program 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 this program; if not, write to the Free Software * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, * MA 02110-1301, USA. Application adaptativeSimpleFoam Description Steady-state solver for incompressible, turbulent flow with remeshing capabilities using the SIMPLE algorithm. \*---------------------------------------------------------------------------*/ #include "fvCFD.H" #include "dynamicFvMesh.H" #include "singlePhaseTransportModel.H" #include "turbulentTransportModel.H" #include "simpleControl.H" #include "fvOptions.H" #include "errorEstimate.H" #include "resError.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "postProcess.H" #include "setRootCase.H" #include "createTime.H" #include "createDynamicFvMesh.H" #include "createControl.H" #include "createFields.H" #include "createFvOptions.H" #include "initContinuityErrs.H" turbulence->validate(); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Calculate absolute flux from the mapped surface velocity //phi = mesh.Sf() & Uf; while (simple.loop()) { Info<< "Time = " << runTime.timeName() << nl << endl; mesh.update(); // --- Pressure-velocity SIMPLE corrector { #include "UEqn.H" #include "pEqn.H" #include "eEqn.H" } laminarTransport.correct(); turbulence->correct(); runTime.write(); Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; } // ************************************************************************* //