adaptativeSimpleFoam

In this section, we propose a solver based on the error estimator described in [this page].(/outils/scripts/openfoam/errorestimate/).

The idea is to propose a SIMPLE solver which, at each time step, calculates an error estimator and refines or de-refines the mesh according to the size of the error linked to the discretisation.

  • 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 \
        -lsampling
    

  • createFields.H

    Show 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.H

    Show 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.H

    Show 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.H

    Show 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.C

    Show 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;
    }
    
    
    // ************************************************************************* //