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
Dhastha
Trackbacks