[FlorisvdBeek]

Hi folks,

For a project, Im trying to simulate the effects of a body force created by a plasma actuator on the flow. The solver is based on rhoPimpleFoam (OpenFOAM 3.0.1), with an additional term in the momentum equation, added to the UEqn.H file.

I have the body force fields as external files (from measurements), so I would like to read in the force fields from these files. As exactly this is happening at the initialization (createFields.H) in folder 0/, I created a readForceField.H:

// Read in the body force

Info<< "Reading field bodyForce\n" << endl; // bodyForce term
volVectorField bodyForce
(
    IOobject
    (
	"bodyForce",
	runTime.timeName(),
	mesh,
	IOobject::MUST_READ,
	IOobject::NO_WRITE
    ),
    mesh
);

Which is executed right before entering the PIMPLE loop:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 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/>.

Application
    rhoPimpleFoam

Description
    Transient solver for laminar or turbulent flow of compressible fluids
    for HVAC and similar applications.

    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
    pseudo-transient simulations.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "psiThermo.H"
#include "turbulentFluidThermoModel.H"
#include "bound.H"
#include "pimpleControl.H"
#include "fvIOoptionList.H"
#include "localEulerDdtScheme.H"
#include "fvcSmooth.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

    pimpleControl pimple(mesh);

    // initialization	
    scalar currentTime=0;
//     scalar currentForcePhase=0;
//     scalar currentEnergyPhase=0;
    scalar period=0;
    scalar globalBodyForceVal=0;
    scalar globalEnergyTermVal=0;
    
    
    #include "createTimeControls.H"
    #include "createRDeltaT.H"
    #include "initContinuityErrs.H"
    #include "createFields.H"
    #include "createMRF.H"
    #include "createFvOptions.H"
    
    // two extra header files that are needed  mk


    #include "createBodyForce4.H"
    
    scalar auxTime=auxTime0;

    #include "mathematicalConstants.H"

    if (!LTS)
    {
        #include "compressibleCourantNo.H"
        #include "setInitialDeltaT.H"
    }

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    

    
    Info<< "\nStarting time loop\n" << endl;

    while (runTime.run())
    {
        #include "readTimeControls.H"

        if (LTS)
        {
            #include "setRDeltaT.H"
        }
        else
        {
            #include "compressibleCourantNo.H"
            #include "setDeltaT.H"
        }

        runTime++;

        Info<< "Time = " << runTime.timeName() << nl << endl;
	
	currentTime=runTime.timeOutputValue(); 
	
	if (  currentTime == auxTime ) 
	{
	#include "readForceField.H"    
	  auxTime +=auxTime0; 
	}


        if (pimple.nCorrPIMPLE() <= 1)
        {

	  #include "rhoEqn.H"
        }

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"
            #include "EEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                if (pimple.consistent())
                {
                    #include "pcEqn.H"
                }
                else
                {
                    #include "pEqn.H"
                }
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

The files are located in folders which would be otherwise created (i.e. 0/, 0.001/, 0.002/ etc..). The read-in is controlled by the if-statement, which checks whether the currentTime is equal to a multitude of auxTime (which initially has the value of writeInterval, the time step after which the solver should write its output, and is increased once currentTime has exceded it). This is necessary as the timestep between the body force measurements is much larger than the timestep of the solver, so the solver should read in the files only once every "writeInterval" timestep. Now my problems:

• The first is the if-statement: "if ( currentTime == auxTime ) ", which does not seem to work.. Both values are scalars, but the solver is skipping the loop even if I monitor that currentTime is equal to auxTime.
• My second problem (if I remove the if-statement), the solver does not change the values of the body force after 'reading in' the file at the new timestep: nothing is happening.

Any ideas/clues/alternative ways to do what Id like to do?

Thanks in advance !

[fedvasu]

Hi Florisvdbeek,

I am pretty sure, you have solved this issue by now...

but just for documentation.

if ( currentTime == auxTime )

This won't and shouldn't work, because comparing floating points is both pointless and perilious.

either convert the times into strings and make sure that each time that directory has the correct time(string value) or if you really want to compare times, make times as label (int).


I am not sure about the second part.