My first EJS simulation: Stopwatch
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………
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.
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
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
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
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:
* 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
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
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.
mayaviFoam: A postprocessing script for OpenFOAM
mayaviFoam:
Hi, I have created a Postprocessing script called mayaviFoam. It is the a reader module to run with Mayavi2, an open-source, visualization application. It will convert OpenFOAM natural data format into VTK data format and the VTK data is loaded into Mayavi2 using VTKFileReader.
#!/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()
Save it as mayaviFoam.py
Procedure:
Open your terminal
1.
cd openfoam_tutorials/multiphase/interFoam/laminar/damBreak
2.
blockMesh
3.
setFields
4. Refere here
python dambreak.py
5.
python mayaviFoam.py
See the output like below
Dambreak simulation using OpenFOAM, Pythonflu, Mayavi2
In this tutorial am going to demonstrate dambreak simulation using Openfoam, pythonflu, mayavi2 in Ubuntu 10.10.
Install following required package
1. OpenFOAM (see the link: https://dhastha.wordpress.com/2011/03/28/openfoam-introduction/)
2. Mayavi2 (link: https://dhastha.wordpress.com/2010/12/18/mayavi2-introduction/)
3. Pythonflu (link: https://dhastha.wordpress.com/2011/04/17/pythonflu-intro/)
copy /opt/openfoam171/tutorials into specified folder.
Open your terminal and type
cd tutorials/multiphase/interFoam/laminar/damBreak
then type following commands
blockMesh setFields
Copy the following python code and run it as “python dambreak.py”
#!/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\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 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 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 )
pass
#--------------------------------------------------------------------------------------
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 )
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() ) + fvc.ddtPhiCorr( rAU, rho, U, phi ) )
from Foam.finiteVolume import adjustPhi
adjustPhi(phiU, U, p)
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 )
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
pass
#--------------------------------------------------------------------------------------
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 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 )
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() :
piso, nCorr, nNonOrthCorr, momentumPredictor, transonic, nOuterCorr = readPISOControls( mesh )
adjustTimeStep, maxCo, maxDeltaT = readTimeControls( runTime )
CoNum, meanCoNum = CourantNo( mesh, phi, runTime )
maxAlphaCo, alphaCoNum, meanAlphaCoNum = alphaCourantNo( runTime, mesh, alpha1, phi )
runTime = setDeltaT( runTime, adjustTimeStep, maxCo, CoNum, maxAlphaCo, alphaCoNum, maxDeltaT )
runTime.increment()
ext_Info() << "Time = " << runTime.timeName() << nl << nl
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 ):
_pEqn( runTime, mesh, UEqn, U, p, p_rgh, gh, ghf, phi, alpha1, rho, g, interface, corr, nCorr, nNonOrthCorr, pRefCell, pRefValue, cumulativeContErr )
pass
from Foam.finiteVolume.cfdTools.incompressible import continuityErrs
cumulativeContErr = continuityErrs( mesh, phi, runTime, cumulativeContErr )
turbulence.correct()
runTime.write()
ext_Info() << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << \
" ClockTime = " << runTime.elapsedClockTime() << " s" << nl << nl
pass
ext_Info() << "End\n" <=", "010700" ):
if __name__ == "__main__" :
argv = sys.argv
if len( argv ) > 1 and argv[ 1 ] == "-test":
argv = None
test_dir= os.path.join( os.environ[ "PYFOAM_TESTING_DIR" ],'cases', 'propogated', 'r1.7.0', 'multiphase','interFoam', 'laminar', 'damBreak' )
argv = [ __file__, "-case", test_dir ]
pass
os._exit( main_standalone( len( argv ), argv ) )
#--------------------------------------------------------------------------------------
then convert the openfoam data into VTK data by typing following command
foamToVTK
It will create a VTK files. then
cd VTK mayavi2 damBreak_0.vtk
It will load vtk files into mayavi2. Then add surface module
select Visualize->Modules->Surface
Use Mayavi object editor to change various timestep and see the instance changes
Followings are the images at different timestep
Pythonflu: Intro
Pythonflu:
Pythonflu is a python front-end to OpenFOAM(Open source CFD toolbox). It indent to define a new level of flexibility and user interaction scenario for numerical simulation software.
Installation from binaries
- Identify the following installation parameters :
- OpenFOAM version you have ( 1.7.1 or 1.6-ext )
- Linux distribution you use ( OpenSUSE 11.3, Ubuntu Maverick or Lucid )
- CPU type you posses ( i386 or amd64 )
- Make sure you have already installed corresponding version of OpenFOAM binary package (OpenCFD or Extended)
- Choose and download proper pythonFlu binary package from SourceForge
- Run the native package manager on the target package :
sudo rpm -i <URL to pythonFlu binary package>.rpm # for OpenSUSEsudo dpkg --install <path pythonFlu binary package>.deb # for Ubuntu
User Configuration
Before using pythonFlu do not forget to configure OpenFOAM itself
source /opt/openfoam171/etc/bashrc # for OpenCFD OpenFOAM-1.7.1source /usr/lib/OpenFOAM-1.6-ext/etc/bashrc # for Extended OpenFOAM-1.6
Getting Started
Create a working directory dedicated to the installed version of OpenFOAM:
mkdir -p ${FOAM_RUN}Copy the tutorial examples directory from the OpenFOAM distribution to your working directory :
cp -r ${FOAM_TUTORIALS} ${FOAM_RUN}Run the first example case of incompressible laminar flow in a cavity:
cd ${FOAM_RUN}/tutorials/incompressible/icoFoam/cavityblockMeshicoFlux
Dhastha
Scientific Visual Coder 