Skip to content

My first EJS simulation: Stopwatch

August 10, 2011

I created my first EJS Stopwatch simulation.

My project manager asked me to do a stopwatch simulation using Easy Java Simulation. And he gave me some constraints. That are

1. Stopwatch should be resizeable

2. It should be placed on position what we given

3. It should have 300 ticks. Each 5th tick is bigger, and each 25th tick is bigger than 5th ticks. We give the numbers to 25th ticks (5, 10, 25 …..)

4. It should be real analog stopwatch. It must match with the second speed of system time.

Design:

my stop will looks like below

The Elements which I have used:

2D drawables –  segment set. text set, arrow, 2d shape

Model:

Variables:

size, position, tick size, tick position, minute hand size, second hand size, text , text angle, degree, degree increment, radian

Initialization:

double StopwatchTempAngle = 2 * Math.PI / 300 ;         //Angle between two ticks. This angle will be used to find position.
int StopwatchIndex ;                                   //Iteration index

//Draw Ticks for Second hand
for(StopwatchIndex = 0 ; StopwatchIndex < 300; StopwatchIndex++) {    StopwatchTickSecondX[StopwatchIndex] = StopwatchPositionX + StopwatchSize * Math.sin(StopwatchIndex * StopwatchTempAngle) ;        //X position of Second hand in circle       StopwatchTickSecondY[StopwatchIndex] = StopwatchPositionY + StopwatchSize * Math.cos(StopwatchIndex * StopwatchTempAngle) ;        //Y position of Second hand in circle      if(StopwatchIndex % 25 == 0)  {      StopwatchTickSecondSizeX[StopwatchIndex]  =   - (StopwatchSize/8) * Math.sin(StopwatchIndex * StopwatchTempAngle) ;    //Second hand tick size X   StopwatchTickSecondSizeY[StopwatchIndex]  =   - (StopwatchSize/8) * Math.cos(StopwatchIndex * StopwatchTempAngle) ;    //Second hand tick size Y           StopwatchTextSecondX[StopwatchIndex/25] = StopwatchTickSecondX[StopwatchIndex] + 1.2 * StopwatchTickSecondSizeX[StopwatchIndex] ;           //Second text X   StopwatchTextSecondY[StopwatchIndex/25] = StopwatchTickSecondY[StopwatchIndex] + 1.2 * StopwatchTickSecondSizeY[StopwatchIndex] ;           //Second text Y      if(StopwatchIndex > 0) StopwatchText[StopwatchIndex/25] = StopwatchIndex/5 + "";                    //Text in each 25th position
  else StopwatchText[0] = "0";                                                                    //Text in 0th position

  StopwatchTextAngle[StopwatchIndex/25] = - StopwatchIndex * StopwatchTempAngle;                  //Text transformation angle
 }

 else if(StopwatchIndex % 5 == 0)
 {

  StopwatchTickSecondSizeX[StopwatchIndex]  =  - (StopwatchSize/10) * Math.sin(StopwatchIndex * StopwatchTempAngle) ;    //Second hand tick size X
  StopwatchTickSecondSizeY[StopwatchIndex]  =  - (StopwatchSize/10) * Math.cos(StopwatchIndex * StopwatchTempAngle) ;    //Second hand tick size Y
  }

 else
 {
  StopwatchTickSecondSizeX[StopwatchIndex]  =   - (StopwatchSize/14) * Math.sin(StopwatchIndex * StopwatchTempAngle) ;    //Second hand tick size X
  StopwatchTickSecondSizeY[StopwatchIndex]  =   - (StopwatchSize/14) * Math.cos(StopwatchIndex * StopwatchTempAngle) ;    //Second hand tick size Y

 }
}

//Draw Ticks for Minute hand
StopwatchTempAngle = 2 * Math.PI / 60 ;         //Angle between two ticks. This angle will be used to find position.

for(StopwatchIndex = 0 ; StopwatchIndex < 60; StopwatchIndex++) {    StopwatchTickMinuteX[StopwatchIndex] = StopwatchPositionX + (StopwatchSize/2) * Math.sin(StopwatchIndex * StopwatchTempAngle) ;        //X position of Minute hand in circle       StopwatchTickMinuteY[StopwatchIndex] = (StopwatchPositionY+(StopwatchSize/4)) + (StopwatchSize/2) * Math.cos(StopwatchIndex * StopwatchTempAngle) ;        //Y position of Minute hand in circle      if(StopwatchIndex % 5 == 0)  {      StopwatchTickMinuteSizeX[StopwatchIndex]  =   - (StopwatchSize/12) * Math.sin(StopwatchIndex * StopwatchTempAngle) ;    //Minute hand tick size X   StopwatchTickMinuteSizeY[StopwatchIndex]  =   - (StopwatchSize/12) * Math.cos(StopwatchIndex * StopwatchTempAngle) ;    //Minute hand tick size Y      StopwatchTextMinuteX[StopwatchIndex/5] = StopwatchTickMinuteX[StopwatchIndex] + 1.2 * StopwatchTickMinuteSizeX[StopwatchIndex] ;           //Minute text X   StopwatchTextMinuteY[StopwatchIndex/5] = StopwatchTickMinuteY[StopwatchIndex] + 1.2 * StopwatchTickMinuteSizeY[StopwatchIndex] ;           //Minute text Y      if(StopwatchIndex > 0) StopwatchText[StopwatchIndex/5] = StopwatchIndex + "";                    //Text in each 5th position
  else StopwatchText[0] = "0";                                                                    //Text in 0th position

  StopwatchTextAngle[StopwatchIndex/5] = - StopwatchIndex * StopwatchTempAngle;                  //Text transformation angle
 }

 else
 {
  StopwatchTickMinuteSizeX[StopwatchIndex]  =   - (StopwatchSize/25) * Math.sin(StopwatchIndex * StopwatchTempAngle) ;    //Minute hand tick size X
  StopwatchTickMinuteSizeY[StopwatchIndex]  =   - (StopwatchSize/25) * Math.cos(StopwatchIndex * StopwatchTempAngle) ;    //Minute hand tick size Y

 }
}

//initialize second hand degree to 0
StopwatchDegreeSecond = 0 ;

Evolution:

FPS (Frame Per Second) – 10

SPD(Steps per Display) – 1

Code:

StopwatchDegreeSecond = StopwatchDegreeSecond + StopwatchDegreeIncrement ;

Fixed Relation:

//Compute Second hand angle
StopwatchRadianSecond = StopwatchDegreeSecond * Math.PI / 180 ; 

//Compute Secon hand position
StopwatchSecondHandX = StopwatchSize * Math.sin(StopwatchRadianSecond) ;
StopwatchSecondHandY = StopwatchSize * Math.cos(StopwatchRadianSecond) ;

//Compute Minute Hand position

StopwatchMinuteHandX = (StopwatchSize/2) * (Math.sin(StopwatchRadianSecond/60)) ;
StopwatchMinuteHandY = (StopwatchSize/2) * (Math.cos(StopwatchRadianSecond/60)) ;

Download Jar file and EJS source file from here

How to run:

Open you terminal and type

java -jar Stopwatch.jar

There are two buttons start/pause, reset
To start stopwatch click on start button.

Download EJS source file from here:

How to reuse this stopwatch in other simulation :

Copy all the page from stopwatch EJS and paste in other simulation. OR manually create variables.

Future Work:

This simulation works perfectly. But currently it has computations. which means that lots loops are used , math function are used. So it taking more computational cost. So I will be work on that, will be post updated stopwatch soon………

How to Install EJS (Easy Java Simulation) in Ubuntu 10.10

July 13, 2011

Introduction:

An important component of software’s success is an inviting and intuitive interface. Graphics and animation often engage a user. However, such features can be the most time-consuming and difficult portions to add to programs. In this article, we introduce and review Easy Java Simulations software, often abbreviated EJS. We will discuss its striking benefits and also its limitations. Most of all, we will see how EJS can do the heavy-lifting of adding graphical user interfaces (GUIs) to your computer programs.

What is Easy Java Simulation?

EJS is a software tool created to help science teachers and students with little or no programming background easily visualize scientific phenomena. Although the freeware application caters to scientific simulations, EJS can be applied to a wide variety of concepts.

Download and Install:

Easy Java Simulations has been written completely in Java language and, as such, can be run on any platform that supports Java.

System Requirements:

Java Runtime Environment (JRE) 1.5 or later.

Linux system.

In my case I used Ubuntu 10.10.

Download the current version from here : http://www.um.es/fem/EjsWiki/uploads/Download/EJS_4.3.3.1_110414.zip

Install:

Unzip the distribution file. This will create a new folder

Start EJS by typing

java -jar EjsConsole.jar

Thats it. EJS console will open to start simulation programming.

Floating Object simulation using OpenFOAM, PythonFLU

June 2, 2011

Introduction:

At any depth in a fluid there is an upward force due to the effect of gravity on the fluid. This results in a pressure applied over an area. If the density of an object in the fluid is greater than the density of the fluid, the object will sink. If the density is less than that of the fluid, the object will float upward due to the buoyancy from the fluid. An object of lower density will float to the top and only be submerged by an amount according to the ratio of the densities.

Problem Setup:

1. Specify mesh, material properties and initial+boundary flow conditions

2. Dynamic mesh type: sixDofMotion. Mesh holds floatingBody object

3. A floating body holds 6 -DOF parameters: mass, moment of inertia, support, forces

4. Variable diffusivity Laplacian motion solver parameters

The geometry is illustrated below

How to run:

I hope you have install required packages (OpenFoam 1.7.1, Paraview 3.8.1, PythonFlu)

Download case file and solver from here

1. Extract the floatingObject.zip in home directory

2. Open your terminal and type following commands

./Allrun

3. This shell script will run mesh generation and solver
4. To view the output. Type the following command in terminal

paraFoam

The script as follows

#!/usr/bin/env python

#----------------------------------------------------------------------------
def _createFields( runTime, mesh ):
    from Foam.OpenFOAM import ext_Info, nl
    from Foam.OpenFOAM import IOdictionary, IOobject, word, fileName
    from Foam.finiteVolume import volScalarField
        
    ext_Info() << "Reading field p_rgh\n" << nl
    p_rgh = volScalarField( IOobject( word( "p_rgh" ),
                                  fileName( runTime.timeName() ),
                                  mesh,
                                  IOobject.MUST_READ,
                                  IOobject.AUTO_WRITE ),
                        mesh )
    
    ext_Info() << "Reading field alpha1\n" << nl
    alpha1 = volScalarField( IOobject( word( "alpha1" ),
                                  fileName( runTime.timeName() ),
                                  mesh,
                                  IOobject.MUST_READ,
                                  IOobject.AUTO_WRITE ),
                                mesh )
    
    ext_Info() << "Reading field U\n" << nl
    from Foam.finiteVolume import volVectorField
    U = volVectorField( IOobject( word( "U" ),
                                  fileName( runTime.timeName() ),
                                  mesh,
                                  IOobject.MUST_READ,
                                  IOobject.AUTO_WRITE ),
                        mesh )
                        
    from Foam.finiteVolume.cfdTools.incompressible import createPhi
    phi = createPhi( runTime, mesh, U )
    
    
    ext_Info() << "Reading transportProperties\n" << nl
    from Foam.transportModels import twoPhaseMixture
    twoPhaseProperties = twoPhaseMixture(U, phi)
    
    rho1 = twoPhaseProperties.rho1()
    rho2 = twoPhaseProperties.rho2()
    
    # Need to store rho for ddt(rho, U)
    rho = volScalarField( IOobject( word( "rho" ),
                                    fileName( runTime.timeName() ),
                                    mesh,
                                    IOobject.READ_IF_PRESENT ),
                          alpha1 * rho1 + ( 1.0 - alpha1 ) * rho2,
                          alpha1.ext_boundaryField().types() )
    rho.oldTime()
    
    # Mass flux
    # Initialisation does not matter because rhoPhi is reset after the
    # alpha1 solution before it is used in the U equation.
    from Foam.finiteVolume import surfaceScalarField
    rhoPhi = surfaceScalarField( IOobject( word( "rho*phi" ),
                                           fileName( runTime.timeName() ),
                                           mesh,
                                           IOobject.NO_READ,
                                           IOobject.NO_WRITE ),
                                 rho1 * phi )
    
    # Construct interface from alpha1 distribution
    from Foam.transportModels import interfaceProperties
    interface = interfaceProperties( alpha1, U, twoPhaseProperties )


    # Construct incompressible turbulence model
    from Foam import incompressible
    turbulence = incompressible.turbulenceModel.New( U, phi, twoPhaseProperties )
    
    from Foam.finiteVolume.cfdTools.general.include import readGravitationalAcceleration
    g = readGravitationalAcceleration( runTime, mesh)

    
    #dimensionedVector g0(g);

    #Read the data file and initialise the interpolation table
    #interpolationTable timeSeriesAcceleration( runTime.path()/runTime.caseConstant()/"acceleration.dat" );
    
    ext_Info() << "Calculating field g.h\n" << nl
    gh = volScalarField( word( "gh" ), g & mesh.C() )
    ghf = surfaceScalarField( word( "ghf" ), g & mesh.Cf() )
    
    p = volScalarField( IOobject( word( "p" ),
                        fileName( runTime.timeName() ),
                        mesh,
                        IOobject.NO_READ,
                        IOobject.AUTO_WRITE ),
                        p_rgh + rho * gh )

    pRefCell = 0
    pRefValue = 0.0
    
    from Foam.finiteVolume import setRefCell, getRefCellValue
    pRefCell, pRefValue = setRefCell( p, mesh.solutionDict().subDict( word( "PISO" ) ), pRefCell, pRefValue )
    
    if p_rgh.needReference():
       p.ext_assign( p + dimensionedScalar( word( "p" ),
                                           p.dimensions(),
                                           pRefValue - getRefCellValue(p, pRefCell) ) )
       p_rgh.ext_assign( p - rho * gh )
       pass

    return p_rgh, p, alpha1, U, phi, rho1, rho2, rho, rhoPhi, twoPhaseProperties, pRefCell, pRefValue, interface, turbulence, g, gh, ghf
    

#--------------------------------------------------------------------------------------
def readControls( runTime, mesh):
    
    from Foam.finiteVolume.cfdTools.general.include import readTimeControls
    adjustTimeStep, maxCo, maxDeltaT = readTimeControls( runTime )

    from Foam.finiteVolume.cfdTools.general.include import readPISOControls
    piso, nCorr, nNonOrthCorr, momentumPredictor, transonic, nOuterCorr = readPISOControls( mesh )
    
    correctPhi = True
    from Foam.OpenFOAM import word, Switch
    if piso.found( word( "correctPhi" ) ):
        correctPhi = Switch(piso.lookup( word( "correctPhi" ) ) )
        pass

    checkMeshCourantNo = False
    if piso.found( word( "checkMeshCourantNo" ) ):
        checkMeshCourantNo = Switch( piso.lookup( word( "checkMeshCourantNo" ) ) )
        pass
    
    return adjustTimeStep, maxCo, maxDeltaT, piso, nCorr, nNonOrthCorr, momentumPredictor, transonic, nOuterCorr, correctPhi, checkMeshCourantNo


#--------------------------------------------------------------------------------------
def setDeltaT( runTime, adjustTimeStep, maxCo, CoNum, maxAlphaCo, alphaCoNum, maxDeltaT ):
    from Foam.OpenFOAM import SMALL
    if adjustTimeStep:
       maxDeltaTFact = min( maxCo / ( CoNum + SMALL ), maxAlphaCo / ( alphaCoNum + SMALL ) )
       deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
       
       runTime.setDeltaT( min( deltaTFact * runTime.deltaT().value(), maxDeltaT ) )
       from Foam.OpenFOAM import ext_Info, nl
       ext_Info() << "deltaT = " << runTime.deltaT().value() << nl
       pass
    return runTime


#--------------------------------------------------------------------------------------
def alphaCourantNo( runTime, mesh, alpha1, phi ):
    from Foam.OpenFOAM import readScalar, word
    maxAlphaCo = readScalar( runTime.controlDict().lookup(word( "maxAlphaCo" ) ) )
    
    alphaCoNum = 0.0
    meanAlphaCoNum = 0.0
    
    if mesh.nInternalFaces():
       from Foam import fvc
       alphaf = fvc.interpolate(alpha1)

       SfUfbyDelta = (alphaf - 0.01).pos() * ( 0.99 - alphaf ).pos() * mesh.deltaCoeffs() * phi.mag()
       alphaCoNum = ( SfUfbyDelta/mesh.magSf() ).ext_max().value() * runTime.deltaT().value()
       meanAlphaCoNum = ( SfUfbyDelta.sum() / mesh.magSf().sum() ).value() * runTime.deltaT().value()
       pass
    
    from Foam.OpenFOAM import ext_Info, nl
    ext_Info() << "Interface Courant Number mean: " << meanAlphaCoNum << " max: " << alphaCoNum << nl
    
    return maxAlphaCo, alphaCoNum, meanAlphaCoNum


#----------------------------------------------------------------------------------------
def fun_correctPhi( runTime, mesh, phi, p, p_rgh, rho, U, cumulativeContErr, nNonOrthCorr, pRefCell, pRefValue ):
    
    from Foam.finiteVolume.cfdTools.incompressible import continuityErrs
    cumulativeContErr = continuityErrs( mesh, phi, runTime, cumulativeContErr )
    from Foam.OpenFOAM import wordList
    from Foam.finiteVolume import zeroGradientFvPatchScalarField
    pcorrTypes = wordList( p_rgh.ext_boundaryField().size(), zeroGradientFvPatchScalarField.typeName )
    
    from Foam.finiteVolume import fixedValueFvPatchScalarField
    for i in range( p.ext_boundaryField().size() ):
       if p_rgh.ext_boundaryField()[i].fixesValue():
          pcorrTypes[i] = fixedValueFvPatchScalarField.typeName
          pass
       pass
    
    from Foam.OpenFOAM import IOdictionary, IOobject, word, fileName, dimensionedScalar
    from Foam.finiteVolume import volScalarField
    pcorr = volScalarField( IOobject( word( "pcorr" ),
                                      fileName( runTime.timeName() ),
                                      mesh,
                                      IOobject.NO_READ,
                                      IOobject.NO_WRITE ),
                            mesh,
                            dimensionedScalar( word( "pcorr" ), p_rgh.dimensions(), 0.0 ),
                            pcorrTypes )
    
    from Foam.OpenFOAM import dimTime
    rAUf = dimensionedScalar( word( "(1|A(U))" ), dimTime / rho.dimensions(), 1.0)
    
    from Foam.finiteVolume import adjustPhi
    adjustPhi(phi, U, pcorr)
    
    from Foam import fvc, fvm
    for nonOrth in range( nNonOrthCorr + 1 ):
        pcorrEqn = fvm.laplacian( rAUf, pcorr ) == fvc.div( phi )

        pcorrEqn.setReference(pRefCell, pRefValue)
        pcorrEqn.solve()

        if nonOrth == nNonOrthCorr:
           phi.ext_assign( phi - pcorrEqn.flux() )
           pass
    
    from Foam.finiteVolume.cfdTools.incompressible import continuityErrs
    cumulativeContErr = continuityErrs( mesh, phi, runTime, cumulativeContErr )
    
    return cumulativeContErr
    

#--------------------------------------------------------------------------------------
def alphaEqn( mesh, phi, alpha1, rhoPhi, rho1, rho2, interface, nAlphaCorr ):
    from Foam.OpenFOAM import word
    alphaScheme = word( "div(phi,alpha)" )
    alpharScheme = word( "div(phirb,alpha)" )
    
    from Foam.finiteVolume import surfaceScalarField
    phic = surfaceScalarField( ( phi / mesh.magSf() ).mag() )
    phic.ext_assign( ( interface.cAlpha() * phic ).ext_min( phic.ext_max() ) )
    phir = phic * interface.nHatf()
    
    from Foam import fvc
    from Foam import MULES

    for aCorr in range( nAlphaCorr ):
       phiAlpha = fvc.flux( phi, alpha1, alphaScheme ) + fvc.flux( -fvc.flux( -phir, 1.0 - alpha1, alpharScheme ), alpha1, alpharScheme )
       MULES.explicitSolve( alpha1, phi, phiAlpha, 1.0, 0.0 )
       
       rhoPhi.ext_assign( phiAlpha * ( rho1 - rho2 ) + phi * rho2 )

       pass
    from Foam.OpenFOAM import ext_Info, nl
    ext_Info() << "Liquid phase volume fraction = " << alpha1.weightedAverage( mesh.V() ).value() \
               << " Min(alpha1) = " << alpha1.ext_min().value() \
               << " Max(alpha1) = " << alpha1.ext_max().value() < 1):
       totalDeltaT = runTime.deltaT()
       rhoPhiSum = 0.0 * rhoPhi
       from Foam.finiteVolume import subCycle_volScalarField
       alphaSubCycle = subCycle_volScalarField(alpha1, nAlphaSubCycles)
       for item in alphaSubCycle:
           alphaEqn( mesh, phi, alpha1, rhoPhi, rho1, rho2, interface, nAlphaCorr )
           rhoPhiSum.ext_assign( rhoPhiSum + ( runTime.deltaT() / totalDeltaT ) * rhoPhi )
           pass
       # To make sure that variable in the local scope will be destroyed
       # - during destruction of this variable it performs some important actions
       # - there is a difference between C++ and Python memory management, namely
       # if C++ automatically destroys stack variables when they exit the scope,
       # Python relay its memory management of some "garbage collection" algorithm
       # that do not provide predictable behavior on the "exit of scope"
       del alphaSubCycle
       
       rhoPhi.ext_assign( rhoPhiSum )
    else:
       alphaEqn( mesh, phi, alpha1, rhoPhi, rho1, rho2, interface, nAlphaCorr )
       pass
    interface.correct()
    
    rho == alpha1 * rho1 + ( 1.0 - alpha1 ) * rho2

    pass


#--------------------------------------------------------------------------------------
def _UEqn( mesh, alpha1, U, p, p_rgh, ghf, rho, rhoPhi, turbulence, g, twoPhaseProperties, interface, momentumPredictor ):
    from Foam.OpenFOAM import word
    from Foam.finiteVolume import surfaceScalarField
    from Foam import fvc
    muEff = surfaceScalarField( word( "muEff" ),
                                twoPhaseProperties.muf() + fvc.interpolate( rho * turbulence.ext_nut() ) )
    from Foam import fvm

    UEqn = fvm.ddt( rho, U ) + fvm.div( rhoPhi, U ) - fvm.laplacian( muEff, U ) - ( fvc.grad( U ) & fvc.grad( muEff ) )
    
    UEqn.relax()

    if momentumPredictor:
       from Foam.finiteVolume import solve
       solve( UEqn == \
                   fvc.reconstruct( ( fvc.interpolate( interface.sigmaK() ) * fvc.snGrad( alpha1 ) - ghf * fvc.snGrad( rho ) \
                                                                                                 - fvc.snGrad( p_rgh ) ) * mesh.magSf() ) )
       pass
    
    return UEqn


#--------------------------------------------------------------------------------------
def _pEqn( runTime, mesh, UEqn, U, p, p_rgh, gh, ghf, phi, alpha1, rho, g, interface, corr, nCorr, nNonOrthCorr, pRefCell, pRefValue, cumulativeContErr ):
    rAU = 1.0/UEqn.A()
     
    from Foam import fvc
    rAUf = fvc.interpolate( rAU )
    
    U.ext_assign( rAU * UEqn.H() )
    
    from Foam.finiteVolume import surfaceScalarField
    from Foam.OpenFOAM import word
    phiU = surfaceScalarField( word( "phiU" ), fvc.interpolate( U ) & mesh.Sf() )
    
    if p_rgh.needReference():
        fvc.makeRelative( phiU, U )
        from Foam.finiteVolume import adjustPhi
        adjustPhi( phiU, U, p )
        fvc.makeAbsolute( phiU, U )
        pass
    
    phi.ext_assign( phiU + ( fvc.interpolate( interface.sigmaK() ) * fvc.snGrad( alpha1 ) - ghf * fvc.snGrad( rho ) )*rAUf*mesh.magSf() )

    from Foam import fvm
    for nonOrth in range( nNonOrthCorr + 1 ):
        p_rghEqn = fvm.laplacian( rAUf, p_rgh ) == fvc.div( phi )
        p_rghEqn.setReference( pRefCell, pRefValue )

        p_rghEqn.solve( mesh.solver( p_rgh.select(corr == nCorr-1 and nonOrth == nNonOrthCorr) ) )
        
        if nonOrth == nNonOrthCorr:
           phi.ext_assign( phi - p_rghEqn.flux() )
           pass
        pass
    
    U.ext_assign( U + rAU * fvc.reconstruct( ( phi - phiU ) / rAUf ) )
    U.correctBoundaryConditions()

    from Foam.finiteVolume.cfdTools.incompressible import continuityErrs
    cumulativeContErr = continuityErrs( mesh, phi, runTime, cumulativeContErr )
    
    # Make the fluxes relative to the mesh motion
    fvc.makeRelative( phi, U )
    
    p == p_rgh + rho * gh

    if p_rgh.needReference():
       from Foam.OpenFOAM import pRefValue
       p.ext_assign( p + dimensionedScalar( word( "p" ),
                                            p.dimensions(),
                                            pRefValue - getRefCellValue(p, pRefCell) ) )
       p_rgh.ext_assign( p - rho * gh )
       pass
    
    return cumulativeContErr
    
    
#--------------------------------------------------------------------------------------
def main_standalone( argc, argv ):

    from Foam.OpenFOAM.include import setRootCase
    args = setRootCase( argc, argv )

    from Foam.OpenFOAM.include import createTime
    runTime = createTime( args )

    from Foam.dynamicFvMesh import createDynamicFvMesh
    mesh = createDynamicFvMesh( runTime )
    
    from Foam.finiteVolume.cfdTools.general.include import readPISOControls
    piso, nCorr, nNonOrthCorr, momentumPredictor, transonic, nOuterCorr = readPISOControls( mesh() )
    
    from Foam.finiteVolume.cfdTools.general.include import initContinuityErrs
    cumulativeContErr = initContinuityErrs()
    
    p_rgh, p, alpha1, U, phi, rho1, rho2, rho, rhoPhi, twoPhaseProperties, pRefCell, \
                                    pRefValue, interface, turbulence, g, gh, ghf = _createFields( runTime, mesh() )

    from Foam.finiteVolume.cfdTools.general.include import readTimeControls
    adjustTimeStep, maxCo, maxDeltaT = readTimeControls( runTime )
    
    cumulativeContErr = fun_correctPhi( runTime, mesh(), phi, p, p_rgh, rho, U, cumulativeContErr, nNonOrthCorr, pRefCell, pRefValue)
    
    from Foam.finiteVolume.cfdTools.incompressible import CourantNo
    CoNum, meanCoNum = CourantNo( mesh(), phi, runTime )
    
    from Foam.finiteVolume.cfdTools.general.include import setInitialDeltaT
    runTime = setInitialDeltaT( runTime, adjustTimeStep, maxCo, maxDeltaT, CoNum )
    
    from Foam.OpenFOAM import ext_Info, nl
    ext_Info() << "\nStarting time loop\n" << nl
    
    while runTime.run() :
        
        adjustTimeStep, maxCo, maxDeltaT, piso, nCorr, nNonOrthCorr, \
            momentumPredictor, transonic, nOuterCorr, correctPhi, checkMeshCourantNo = readControls( runTime, mesh() )
        
        maxAlphaCo, alphaCoNum, meanAlphaCoNum = alphaCourantNo( runTime, mesh(), alpha1, phi )
        
        CoNum, meanCoNum = CourantNo( mesh(), phi, runTime )
        
        from Foam import fvc
        # Make the fluxes absolute
        fvc.makeAbsolute( phi, U )
        
        runTime = setDeltaT( runTime, adjustTimeStep, maxCo, CoNum, maxAlphaCo, alphaCoNum, maxDeltaT )
        
        runTime.increment()
        ext_Info() << "Time = " << runTime.timeName() << nl << nl
        
        timeBeforeMeshUpdate = runTime.elapsedCpuTime()

        # Do any mesh changes
        mesh.update()
        
        if mesh.changing():
            ext_Info() << "Execution time for mesh.update() = " << runTime.elapsedCpuTime() - timeBeforeMeshUpdate << " s" << nl
            gh.ext_assign( g & mesh.C() )
            ghf.ext_assign( g & mesh.Cf() )
            pass
        
        if mesh.changing() and correctPhi:
            cumulativeContErr = fun_correctPhi( runTime, mesh(), phi, p, p_rgh, rho, U, cumulativeContErr, nNonOrthCorr, pRefCell, pRefValue)
            pass

        # Make the fluxes relative to the mesh motion
        fvc.makeRelative( phi, U )

        if mesh.changing() and checkMeshCourantNo:
           from Foam.dynamicFvMesh import meshCourantNo
           meshCoNum, meanMeshCoNum = meshCourantNo( runTime, mesh(), phi )
           pass

        twoPhaseProperties.correct()
     
        alphaEqnSubCycle( runTime, piso, mesh(), phi, alpha1, rho, rhoPhi, rho1, rho2, interface )
        
        UEqn = _UEqn( mesh(), alpha1, U, p, p_rgh, ghf, rho, rhoPhi, turbulence, g, twoPhaseProperties, interface, momentumPredictor )

        # --- PISO loop
        for corr in range( nCorr ):
            cumulativeContErr = _pEqn( runTime, mesh(), UEqn, U, p, p_rgh, gh, ghf, phi, alpha1, rho, g, \
                                       interface, corr, nCorr, nNonOrthCorr, pRefCell, pRefValue, cumulativeContErr )
            pass
        
        turbulence.correct()

        runTime.write()

        ext_Info() << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << \
              " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << nl
        
        pass

    ext_Info() << "End\n" <=", "010701" ):
   if __name__ == "__main__" :
      import sys, os
      argv = sys.argv
      os._exit( main_standalone( len( argv ), argv ) )
      pass
   pass
else:
   from Foam.OpenFOAM import ext_Info
   ext_Info()<< "\nTo use this solver, It is necessary to SWIG OpenFoam1.7.1 or higher \n "
   pass

#--------------------------------------------------------------------------------------

Output:

Wireframe View:

At 1 timestep

At 5 timestep

At 10 timestep

At 15 timestep

At 30 timestep



Volume view:

At time step 1

At time step 5

At time step 15

At time step 25

Two Liquid Mixing simulation using OpenFoam, PythonFlu, Mayavi2

May 18, 2011

To see Floating object simulation using openFoam and Pythonflu

https://dhastha.wordpress.com/2011/06/02/floating-object-simulation-using-openfoam-pythonflu/

To see Dambreak simulation using openFoam and PythonFlu

https://dhastha.wordpress.com/2011/04/17/dambreak-simulation-using-openfoam-pythonflu-mayavi2/

Short description:

I am going to demonstrate twoStream simulation which means that how two fluids are mixing using OpenFoam, PythonFlu, and Mayavi2. I am going to use OpenFoam and Pythonflu for pre-processing, simulation and finally Mayavi2 for post-processing.

Problem definition:

This problem is a mixing of two fluid in a two-dimensional geometry. The geometry has closed walls and two fluids are separated at initial state.We are examining the mixing of  flow at steady state. The geometry is illustrated below:

How to run:

To run this case you should install OpenFoam 1.7.1, Pythonflu and Mayavi2.

Download case file and solver from here

1. Extract the twoLiquidMixingFlux.tar in home directory

2. Open your terminal and type following commands

 cd twoLiquidMixingFlux
 blockMesh
 setFields

3. Run the solver

 python twoLiquidMixingFlux.py

The script as follows

#!/usr/bin/env python

#----------------------------------------------------------------------------
## Copyright (C) 2011 Dhasthagheer

#----------------------------------------------------------------------------
def _createFields( runTime, mesh, g ):
    from Foam.OpenFOAM import ext_Info, nl
    from Foam.OpenFOAM import IOdictionary, IOobject, word, fileName
    from Foam.finiteVolume import volScalarField
    ext_Info() << "Reading field p_rgh\n" << nl
    p_rgh = volScalarField( IOobject( word( "p_rgh" ),
                                      fileName( runTime.timeName() ),
                                      mesh,
                                      IOobject.MUST_READ,
                                      IOobject.AUTO_WRITE ),
                            mesh )

    ext_Info() << "Reading field alpha1\n" << nl
    alpha1 = volScalarField( IOobject( word( "alpha1" ),
                                       fileName( runTime.timeName() ),
                                       mesh,
                                       IOobject.MUST_READ,
                                       IOobject.AUTO_WRITE ),
                             mesh )

    ext_Info() << "Reading field U\n" << nl
    from Foam.finiteVolume import volVectorField
    U = volVectorField( IOobject( word( "U" ),
                                  fileName( runTime.timeName() ),
                                  mesh,
                                  IOobject.MUST_READ,
                                  IOobject.AUTO_WRITE ),
                        mesh )

    from Foam.finiteVolume.cfdTools.incompressible import createPhi
    phi = createPhi( runTime, mesh, U )

    ext_Info() << "Reading transportProperties\n" << nl
    from Foam.transportModels import twoPhaseMixture
    twoPhaseProperties = twoPhaseMixture( U, phi )

    rho1 = twoPhaseProperties.rho1()
    rho2 = twoPhaseProperties.rho2()

    from Foam.OpenFOAM import dimensionedScalar
    Dab = dimensionedScalar( twoPhaseProperties.lookup( word( "Dab" ) ) )

    # Read the reciprocal of the turbulent Schmidt number
    alphatab = dimensionedScalar( twoPhaseProperties.lookup( word( "alphatab" ) ) )

    # Need to store rho for ddt(rho, U)

    from Foam.OpenFOAM import scalar
    rho = volScalarField ( word( "rho" ), alpha1 * rho1 + ( scalar(1) - alpha1 ) * rho2 )
    rho.oldTime()

    # Mass flux
    # Initialisation does not matter because rhoPhi is reset after the
    # alpha1 solution before it is used in the U equation.
    from Foam.finiteVolume import surfaceScalarField
    rhoPhi = surfaceScalarField( IOobject( word( "rho*phi" ),
                                           fileName( runTime.timeName() ),
                                           mesh,
                                           IOobject.NO_READ,
                                           IOobject.NO_WRITE ),
                                 rho1 * phi )

    # Construct incompressible turbulence model
    from Foam import incompressible
    turbulence = incompressible.turbulenceModel.New( U, phi, twoPhaseProperties )

    ext_Info() << "Calculating field g.h\n" << nl
    gh = volScalarField ( word( "gh" ), g & mesh.C() )
    ghf = surfaceScalarField ( word( "ghf" ), g & mesh.Cf() )

    p = volScalarField( IOobject( word( "p" ),
                                  fileName( runTime.timeName() ),
                                  mesh,
                                  IOobject.NO_READ,
                                  IOobject.AUTO_WRITE ),
                        p_rgh + rho * gh )

    pRefCell = 0
    pRefValue = 0.0
    from Foam.finiteVolume import setRefCell, getRefCellValue
    pRefCell, pRefValue = setRefCell( p, p_rgh, mesh.solutionDict().subDict( word( "PIMPLE" ) ), pRefCell, pRefValue )

    if p_rgh.needReference():
        p.ext_assign( p + dimensionedScalar( word( "p" ),
                                             p.dimensions(),
                                             pRefValue - getRefCellValue( p, pRefCell ) ) )
        p_rgh.ext_assign( p - rho * gh )
        pass

    return p_rgh, alpha1, U, phi, twoPhaseProperties, rho1, rho2, Dab, alphatab, rho, rhoPhi, turbulence, gh, ghf, p, pRefCell, pRefValue

#--------------------------------------------------------------------------------------
def alphaEqn( mesh, phi, alpha1, alphatab, Dab, rhoPhi, rho, rho1, rho2, turbulence ):
    from Foam import fvm
    from Foam.OpenFOAM import word
    tmp = turbulence.ext_nut()
    alpha1Eqn = fvm.ddt( alpha1 ) + fvm.div( phi, alpha1 ) - fvm.laplacian( Dab + alphatab * tmp() , alpha1, word( "laplacian(Dab,alpha1)" ) )

    alpha1Eqn.solve()

    rhoPhi.ext_assign( alpha1Eqn.flux() * ( rho1 - rho2 ) + phi * rho2 )
    from Foam.OpenFOAM import scalar
    rho.ext_assign( alpha1 * rho1 + ( scalar( 1 ) - alpha1 ) * rho2 )

    from Foam.OpenFOAM import ext_Info, nl
    ext_Info() << "Phase 1 volume fraction = " << alpha1.weightedAverage( mesh.V() ).value()  \
               << "  Min(alpha1) = " << alpha1.ext_min().value() << "  Max(alpha1) = " << alpha1.ext_max().value() << nl

    pass

#----------------------------------------------------------------------------------------
def fun_UEqn( mesh, U, p_rgh, ghf, rho, rhoPhi, turbulence, twoPhaseProperties, momentumPredictor, finalIter ):
    from Foam.OpenFOAM import word
    from Foam.finiteVolume import surfaceScalarField
    from Foam import fvc
    muEff = surfaceScalarField( word( "muEff" ),
                                twoPhaseProperties.muf() + fvc.interpolate( rho * turbulence.ext_nut() ) )

    from Foam import fvm, fvc

    UEqn = fvm.ddt( rho, U ) + fvm.div( rhoPhi, U ) - fvm.laplacian( muEff, U ) - ( fvc.grad( U ) & fvc.grad( muEff ) )

    if finalIter:
        UEqn.relax( 1.0 )
        pass
    else:
        UEqn.relax()
        pass

    if momentumPredictor:
        from Foam.finiteVolume import solve
        solve( UEqn  == fvc.reconstruct( ( - ghf * fvc.snGrad( rho ) - fvc.snGrad( p_rgh ) ) * mesh.magSf() ),
               mesh.solver( U.select( finalIter ) ) )
        pass

    return UEqn

#--------------------------------------------------------------------------------------
def fun_pEqn( runTime, mesh, UEqn, U, p, p_rgh, gh, ghf, phi, rho, finalIter, corr, nCorr, nNonOrthCorr, pRefCell, pRefValue, cumulativeContErr ):
    rAU = 1.0/UEqn.A()

    from Foam import fvc
    rAUf = fvc.interpolate( rAU )

    U.ext_assign( rAU * UEqn.H() )

    from Foam.finiteVolume import surfaceScalarField
    from Foam.OpenFOAM import word
    phiU = surfaceScalarField( word( "phiU" ),
                               ( fvc.interpolate( U ) & mesh.Sf() ) + fvc.ddtPhiCorr( rAU, rho, U, phi ) )

    from Foam.finiteVolume import adjustPhi
    adjustPhi(phiU, U, p)

    phi.ext_assign( phiU - ghf * fvc.snGrad( rho ) * rAUf * mesh.magSf() )

    from Foam import fvm
    for nonOrth in range( nNonOrthCorr + 1 ):
        p_rghEqn = fvm.laplacian( rAUf, p_rgh ) == fvc.div( phi )
        p_rghEqn.setReference( pRefCell, pRefValue )

        p_rghEqn.solve( mesh.solver( p_rgh.select( finalIter and corr == nCorr-1 and nonOrth == nNonOrthCorr) ) )

        if nonOrth == nNonOrthCorr:
           phi.ext_assign( phi - p_rghEqn.flux() )
           pass
        pass

    U.ext_assign( U + rAU * fvc.reconstruct( ( phi - phiU ) / rAUf ) )
    U.correctBoundaryConditions()

    from Foam.finiteVolume.cfdTools.incompressible import continuityErrs
    cumulativeContErr = continuityErrs( mesh, phi, runTime, cumulativeContErr )

    p == p_rgh + rho * gh

    if p_rgh.needReference():
       from Foam.OpenFOAM import pRefValue
       p.ext_assign( p + dimensionedScalar( word( "p" ),
                                            p.dimensions(),
                                            pRefValue - getRefCellValue(p, pRefCell) ) )
       p_rgh.ext_assign( p - rho * gh )
       pass

    return cumulativeContErr

#--------------------------------------------------------------------------------------
def main_standalone( argc, argv ):

    from Foam.OpenFOAM.include import setRootCase
    args = setRootCase( argc, argv )

    from Foam.OpenFOAM.include import createTime
    runTime = createTime( args )

    from Foam.OpenFOAM.include import createMesh
    mesh = createMesh( runTime )

    from Foam.finiteVolume.cfdTools.general.include import readGravitationalAcceleration
    g = readGravitationalAcceleration( runTime, mesh )

    from Foam.finiteVolume.cfdTools.general.include import readPIMPLEControls
    pimple, nOuterCorr, nCorr, nNonOrthCorr, momentumPredictor, transonic = readPIMPLEControls( mesh )

    from Foam.finiteVolume.cfdTools.general.include import initContinuityErrs
    cumulativeContErr = initContinuityErrs()

    p_rgh, alpha1, U, phi, twoPhaseProperties, rho1, rho2, Dab, \
    alphatab, rho, rhoPhi, turbulence, gh, ghf, p, pRefCell, pRefValue = _createFields( runTime, mesh, g )

    from Foam.finiteVolume.cfdTools.general.include import readTimeControls
    adjustTimeStep, maxCo, maxDeltaT = readTimeControls( runTime )

    from Foam.finiteVolume.cfdTools.incompressible import CourantNo
    CoNum, meanCoNum = CourantNo( mesh, phi, runTime )

    from Foam.finiteVolume.cfdTools.general.include import setInitialDeltaT
    runTime = setInitialDeltaT( runTime, adjustTimeStep, maxCo, maxDeltaT, CoNum )

    from Foam.OpenFOAM import ext_Info, nl
    ext_Info() << "\nStarting time loop\n" << nl

    while runTime.run():
        pimple, nOuterCorr, nCorr, nNonOrthCorr, momentumPredictor, transonic = readPIMPLEControls( mesh )
        adjustTimeStep, maxCo, maxDeltaT = readTimeControls( runTime )
        CoNum, meanCoNum = CourantNo( mesh, phi, runTime )

        from Foam.finiteVolume.cfdTools.general.include import setDeltaT
        runTime = setDeltaT( runTime, adjustTimeStep, maxCo, maxDeltaT, CoNum )

        runTime.increment()
        ext_Info() << "Time = " << runTime.timeName() << nl << nl

        # --- Pressure-velocity PIMPLE corrector loop
        for oCorr in range( nOuterCorr ):
            finalIter = oCorr == nOuterCorr-1

            twoPhaseProperties.correct()

            alphaEqn( mesh, phi, alpha1, alphatab, Dab, rhoPhi, rho, rho1, rho2, turbulence )

            UEqn = fun_UEqn( mesh, U, p_rgh, ghf, rho, rhoPhi, turbulence, twoPhaseProperties, momentumPredictor, finalIter )

            # --- PISO loop
            for corr in range( nCorr ):
                cumulativeContErr = fun_pEqn( runTime, mesh, UEqn, U, p, p_rgh, gh, ghf, phi, rho, \
                                              finalIter, corr, nCorr, nNonOrthCorr, pRefCell, pRefValue, cumulativeContErr )
                pass

            turbulence.correct()
            pass

        runTime.write()
        ext_Info() << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << "  ClockTime = " << runTime.elapsedClockTime() << " s" << nl << nl

        pass

    ext_Info() << "End\n" << nl      import os     return os.EX_OK #-------------------------------------------------------------------------------------- from Foam import FOAM_REF_VERSION if FOAM_REF_VERSION( ">=", "010701" ):
   if __name__ == "__main__" :
      import sys, os
      argv = sys.argv
      os._exit( main_standalone( len( argv ), argv ) )
      pass
   pass
else:
   from Foam.OpenFOAM import ext_Info
   ext_Info()<< "\nTo use this solver, It is necessary to SWIG OpenFoam1.7.1 or higher \n "
   pass

#--------------------------------------------------------------------------------------

4. Copy and save following code as mayaviFoam.py in same directory

#!/usr/bin/env mayavi2

import os
import glob
from os.path import join, abspath, dirname
from enthought.mayavi.scripts import mayavi2
from enthought.mayavi.sources.vtk_file_reader import VTKFileReader
from enthought.mayavi.modules.surface import Surface

def setup_data(fname):
    """Given a VTK file name `fname`, this creates a mayavi2
    reader for it and adds it to the pipeline.  It returns the reader
    created.
    """
    r = VTKFileReader()
    r.initialize(fname)
    mayavi.add_source(r)
    return r

def surface():
    """Sets up the mayavi pipeline for the visualization.
    """
    # Create an outline for the data.
    s = Surface()
    mayavi.add_module(s)

def load_vtk_file():
    mayavi.new_scene()
    filepath = glob.glob("VTK/*_0.vtk")
    vtk_file_name = os.path.basename(filepath[0])
    data_dir = "VTK"
    fname = join(data_dir, vtk_file_name)
    r = setup_data(fname)
    surface()

@mayavi2.standalone
def main():
    if os.path.exists("VTK"):
        load_vtk_file()

    else:
        os.system('foamToVTK -ascii')
        load_vtk_file()

if __name__ == '__main__':
    main()

5. Run post processing script mayaviFoam.py

 python mayaviFoam.py

Output:

At 0 timestep
fs

At time step 50

A course in CFD with Opensource Softwares

May 4, 2011

Hi friends, I started a course in CFD with Opensource softwares. Here , I am going to share my day by day experience. I hope this will help those who all are interest in CFD.

Purpose of this course:

1. To give an introduction to Open source software for CFD

2. To give an introduction to OpenFOAM in order to ‘get started’

3. To introduce how to modify OpenFOAM for specific purposes

4. To increase the use of OpenFOAM in India

How to learn more after this course:

1. Learn by doing

2. Exchange knowledge with other OpenFOAM users at the Forum

3. For more advance courses, http://www.openfoam.org

4. Manuals and source codes are available in http://www.openfoam.org

5. The D’oxygen manual

6. The OpenFOAM wiki (www.openfoamwiki.net)

Syllabus for first and second days:

1. An introduction to how to install and use standard solvers, utilities and libraries of OpenFoam

2. An introduction to the organization of the code and cases

3. We will use the OpenFoam User Guide and Programmers guide as base, but I will add my personal experience

Syllabus for third and fourth days:

1. We will have a deeper look into the code and make our own solvers, utilities and libraries. For this we must know how to compile all part, of the code. In particulars we will have a look at turbulence models and boundary conditions.

2. We will use the OpenFoam user guide, and programming guide as reference material. It is good to invest in a book in C++, like ‘C++ Direkt’ by Jan Skansholm (student litteratur)

3. After these days you should be able to investigate the code and find out what it does. You should also be able to make simple modification of the code.

Syllabus for 5th day:

1. Programming in OpenFOAM

2. Examples of different ways of using OpenFOAM.

Learning Outcomes:

1. Learn how to download, install, compile and run OpenFOAM solvers and utilities

2. Learn how to implement solvers and utilities

3. Learn how to implement a turbulence model

4. Learn how to implement boundary condition

5. Learn the basics of C++ and object orientation

6. Learn how to do CFD  with OpenFOAM together with python, paraview, mayavi2

7. Learn basics of Linux, Doxygen, Compilation procedure, Debugging, Version Control System and VTK

8. Learn how to use OpenFoam by doing a project work

Reference:

* http://www.openfoam.org

* OpenFoam userguide

* OpenFoam programming guide

* OpenFoam Doxygen

* C++ direkt, Jan Skansholm, Studentlitteratur

* An introduction to Computational Fluid Dynamics, H K Versteeg & W Malalaskera

* Computational Methods for Fluid Dynamics , J.H. Ferziger & M. Peric

How to compile and Install OpenFOAM 1.6

April 30, 2011

Step1:

To get the sources you will need to do the following

cd ~/OpenFOAM/
git clone http://repo.or.cz/r/OpenFOAM-1.6.x.git/
cd OpenFOAM-1.6.x
git pull

Step 2:

Download Third party source from http://space.dl.sourceforge.net/project/foam/foam/1.6/ThirdParty-1.6.General.gtgz

cd ~/OpenFOAM/
tar -zxf ThirdParty-1.6.General.gtgz

Step 3:

In ~/OpenFOAM/OpenFOAM-1.6.x/etc/settings.sh change the line (line no 98)

compilerInstall=OpenFOAM

to

compilerInstall=System

Step 4:

Next source the OF setting by

. ~/OpenFOAM/OpenFOAM-1.6.x/etc/bashrc

Step 5:

Next step is to compile OF.

Set the number of processors to use, in this case 4

export WM_NCOMPPROCS=4

Compile the ThirdParty stuff

cd ~/OpenFOAM/
ln -s ThirdParty-1.6 ThirdParty-1.6.x
cd ThirdParty-1.6.x
./Allwmake

Assuming you have no errors, we can now compile OpenFOAM re-source the settings (to pick up newly created directories)

cd ~/OpenFOAM/OpenFOAM-1.6.x
. etc/bashrc

Then compile away

./Allwmake

Thats it. You have done.

CFD Analysis Process

April 27, 2011

CFD:

Computational Fluid Dynamics (CFD) provides a qualitative (and sometimes even quantitative) prediction of fluid flows by means of
• mathematical modeling (partial differential equations)
• numerical methods (discretization and solution techniques)
• software tools (solvers, pre- and postprocessing utilities)

CFD Analysis Process:

The general process for performing a CFD analysis is outlined below so as to provide a reference for understanding the various aspects of a CFD simulation. The process includes:

1. Forumulate the Flow Problem
2. Model the Geometry and Flow Domain
3. Establish the Boundary and Initial Conditions
4. Generate the Grid
5. Establish the Simulation Strategy
6. Establish the Input Parameters and Files
7. Perform the Simulation
8. Monitor the Simulation for Completion
9. Post-process the Simulation to get the Results
10. Make Comparisons of the Results
11. Repeat the Process to Examine Sensitivities
12. Document

In further detail, these steps include:

   1.      Formulate the Flow Problem

The first step of the analysis process is to formulate the flow problem by seeking answers to the following questions:
* what is the objective of the analysis?
* what is the easiest way to obtain those objective?
* what geometry should be included?
* what are the freestream and/or operating conditions?
* what dimensionality of the spatial model is required? (1D, quasi-1D, 2D, axisymmetric, 3D)
* what should the flow domain look like?
* what temporal modeling is appropriate? (steady or unsteady)
* what is the nature of the viscous flow? (inviscid, laminar, turbulent)
* how should the gas be modeled?

   2.      Model the Geometry and Flow Domain

The body about which flow is to be analyzed requires modeling. This generally involves modeling the geometry with a CAD software package. Approximations of the geometry and simplifications may be required to allow an analysis with reasonable effort. Concurrently, decisions are made as to the extent of the finite flow domain in which the flow is to be simulated. Portions of the boundary of the flow domain conicide with the surfaces of the body geometry. Other surfaces are free boundaries over which flow enters or leaves. The geometry and flow domain are modeled in such a manner as to provide input for the grid generation. Thus, the modeling often takes into account the structure and topology of the grid generation.

   3.      Establish the Boundary and Initial Conditions

Since a finite flow domain is specified, physical conditions are required on the boundaries of the flow domain. The simulation generally starts from an initial solution and uses an iterative method to reach a final flow field solution.

   4.      Generate the Grid

The flow domain is discretized into a grid. The grid generation involves defining the structure and topology and then generating a grid on that topology. Currently all cases involve multi-block, structured grids; however, the grid blocks may be abbuting, contiguous, non-contiguous, and overlapping. The grid should exhibit some minimal grid quality as defined by measures of orthogonality (especially at the boundaries), relative grid spacing (15% to 20% stretching is considered a maximum value), grid skewness, etc… Further the maximum spacings should be consistent with the desired resolution of important features. The resolution of boundary layers requires the grid to be clustered in the direction normal to the surface with the spacing of the first grid point off the wall to be well within the laminar sublayer of the boundary layer. For turbulent flows, the first point off the wall should exhibit a y+ value of less than 1.0.

   5.      Establish the Simulation Strategy

The strategy for performing the simulation involves determining such things as the use of space-marching or time-marching, the choice of turbulence or chemistry model, and the choice of algorithms.

   6.      Establish the Input Parameters and Files

A CFD codes generally requires that an input data file be created listing the values of the input parameters consisted with the desired strategy. Further the a grid file containing the grid and boundary condition information is generally required. The files for the grid and initial flow solution need to be generated.

   7.      Perform the Simulation

The simulation is performed with various possible with options for interactive or batch processing and distributed processing.

   8.      Monitor the Simulation for Completion

As the simulation proceeds, the solution is monitored to determine if a “converged” solution has been obtained, which is iterative convergence. Further discussion can be found on the page entitled Examining Iterative Convergence.

   9.      Post-Process the Simulation to get the Results

Post-Processing involves extracting the desired flow properties (thrust, lift, drag, etc…) from the computed flowfield.

  10.      Make Comparisons of the Results

The computed flow properties are then compared to results from analytic, computational, or experimental studies to establish the validity of the computed results.

  11.      Repeat the Process to Examine Sensitivities

The sensitivity of the computed results should be examined to understand the possible differences in the accuracy of results and / or performance of the computation with respect to such things as:
* dimensionality
* flow conditions
* initial conditions
* marching strategy
* algorithms
* grid topology and density
* turbulence model
* chemistry model
* flux model
* artificial viscosity
* boundary conditions
* computer system

Further information can be found on the pages entitled Verification Assessment and Validation Assessment.

  12.      Document

Documenting the findings of an analysis involves describing each of these steps in the process.