Usage examples¶
This section documents the example scripts shipped with the source code.
Lid-driven cavity flow¶
The openfoam_cavity_tutorial.py
script provides an example of low-level manipulation of OpenFOAM cases. In this
example we shall re-create the initial example from the OpenFOAM users’
guide. It’s worth reading over that section first before trying to follow the
transliteration below.
Firstly, we need to import some things from the firefish.case
module:
from firefish.case import Case, FileName, FileClass
The Case
class encapsulates an OpenFOAM case
directory. We don’t want to overwrite an existing case and so we write a little
convenience wrapper function:
def create_new_case(case_dir):
# Check that the specified case directory does not already exist
if os.path.exists(case_dir):
raise RuntimeError(
'Refusing to write to existing path: {}'.format(case_dir)
)
# Create the case
return Case(case_dir)
The mutable_data_file()
method will return a context
manager which can be used to manipulate an OpenFOAM data file. The data file is
created if it does not yet exist, its contents are parsed into a dictionary and
the dictionary is returned from the context manager. One the context is left the
dictionary is re-written to disk.
The upshot of this is that the programmer is insulated from manipulating OpenFOAM data files directly. Let’s write the controlDict file from the tutorial:
def write_initial_control_dict(case):
# Control dict from tutorial
control_dict = {
'application': 'icoFoam',
'startFrom': 'startTime',
'startTime': 0,
'stopAt': 'endTime',
'endTime': 0.5,
'deltaT': 0.005,
'writeControl': 'timeStep',
'writeInterval': 20,
'purgeWrite': 0,
'writeFormat': 'ascii',
'writePrecision': 6,
'writeCompression': 'off',
'timeFormat': 'general',
'timePrecision': 6,
'runTimeModifiable': True,
}
with case.mutable_data_file(FileName.CONTROL) as d:
d.update(control_dict)
Well-known file names are available through the FileName
class.
The blockMeshDict file is the next one to be created. This is an example of a relatively complex file. The complexity is somewhat hidden by the mapping to and from the Python domain but there is still some subtlety. Notice particularly the rather odd way in which the boundary dictionary is specified:
def write_block_mesh_dict(case):
block_mesh_dict = {
'convertToMeters': 0.1,
'vertices': [
[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0],
[0, 0, 0.1], [1, 0, 0.1], [1, 1, 0.1], [0, 1, 0.1],
],
'blocks': [
(
'hex', [0, 1, 2, 3, 4, 5, 6, 7], [20, 20, 1],
'simpleGrading', [1, 1, 1],
)
],
'edges': [],
# Note the odd way in which boundary is defined here as a
# list of tuples.
'boundary': [
('movingWall', {
'type': 'wall',
'faces': [ [3, 7, 6, 2] ],
}),
('fixedWalls', {
'type': 'wall',
'faces': [
[0, 4, 7, 3],
[2, 6, 5, 1],
[1, 5, 4, 0],
],
}),
('frontAndBack', {
'type': 'empty',
'faces': [
[0, 3, 2, 1],
[4, 5, 6, 7],
],
}),
],
'mergePatchPairs': [],
}
with case.mutable_data_file(FileName.BLOCK_MESH) as d:
d.update(block_mesh_dict)
At this point in the tutorial we’re ready to run the blockMesh command which is one function call away:
case.run_tool('blockMesh')
We’re close to being able to run the icoFoam utility. The transport properties need to be defined:
from firefish.case import Dimension
with case.mutable_data_file(FileName.TRANSPORT_PROPERTIES) as tp:
tp['nu'] = (Dimension(0, 2, -1, 0, 0, 0, 0), 0.01)
We also need to create the initial conditions. Notice that we have to specify a different class when creating them:
def write_initial_conditions(case):
# Create the p initial conditions
p_file = case.mutable_data_file(
'0/p', create_class=FileClass.SCALAR_FIELD_3D
)
with p_file as p:
p.update({
'dimensions': Dimension(0, 2, -2, 0, 0, 0, 0),
'internalField': ('uniform', 0),
'boundaryField': {
'movingWall': { 'type': 'zeroGradient' },
'fixedWalls': { 'type': 'zeroGradient' },
'frontAndBack': { 'type': 'empty' },
},
})
# Create the U initial conditions
U_file = case.mutable_data_file(
'0/U', create_class=FileClass.VECTOR_FIELD_3D
)
with U_file as U:
U.update({
'dimensions': Dimension(0, 1, -1, 0, 0, 0, 0),
'internalField': ('uniform', [0, 0, 0]),
'boundaryField': {
'movingWall': {
'type': 'fixedValue', 'value': ('uniform', [1, 0, 0])
},
'fixedWalls': {
'type': 'fixedValue', 'value': ('uniform', [0, 0, 0])
},
'frontAndBack': { 'type': 'empty' },
},
})
Before we can run icoFoam, we must create the mysterious fvSolution file:
def write_fv_solution(case):
fv_solution = {
'solvers': {
'p': {
'solver': 'PCG',
'preconditioner': 'DIC',
'tolerance': 1e-6,
'relTol': 0,
},
'U': {
'solver': 'smoothSolver',
'smoother': 'symGaussSeidel',
'tolerance': 1e-5,
'relTol': 0,
},
},
'PISO': {
'nCorrectors': 2,
'nNonOrthogonalCorrectors': 0,
'pRefCell': 0,
'pRefValue': 0,
}
}
with case.mutable_data_file(FileName.FV_SOLUTION) as d:
d.update(fv_solution)
And also the equally mysterious fvSchemes file:
def write_fv_schemes(case):
fv_schemes = {
'ddtSchemes': { 'default': 'Euler' },
'gradSchemes': { 'default': 'Gauss linear', 'grad(p)': 'Gauss linear' },
'divSchemes': { 'div(phi,U)': 'Gauss linear', 'default': 'none' },
'laplacianSchemes': { 'default': 'Gauss linear orthogonal' },
'interpolationSchemes': { 'default': 'linear' },
'snGradSchemes': { 'default': 'orthogonal' },
}
with case.mutable_data_file(FileName.FV_SCHEMES) as d:
d.update(fv_schemes)
The example script includes a main()
function which performs all of
these steps:
def main(case_dir='cavity'):
# Create a new case file, raising RuntimeError if the directory already
# exists.
case = create_new_case(case_dir)
# Add the information needed by blockMesh.
write_initial_control_dict(case)
write_block_mesh_dict(case)
# At this point there is enough to run blockMesh.
case.run_tool('blockMesh')
# Update the physical properties.
with case.mutable_data_file(FileName.TRANSPORT_PROPERTIES) as tp:
tp['nu'] = (Dimension(0, 2, -1, 0, 0, 0, 0), 0.01)
# Write the fvSolution and fvSchemes files.
write_fv_solution(case)
write_fv_schemes(case)
# Write the initial conditions for the p and U fields.
write_initial_conditions(case)
# Run the icoFoam application.
case.run_tool('icoFoam')
After the example script is run, paraFoam may be run to inspect the solution.
Supersonic flow over wedge¶
The openfoam_supersonic_wedge.py
script provides an example of setting up a compressible flow solver in OpenFoam.
As in openfoam_cavity_tutorial.py
we set up the OpenFOAM case directory using the firefish.case framework.
For flows with a Mach number above 0.3 compressible effects become non negligible. A compressible solver must therefore be used. In this case we use rhoCentralFoam. The control file must therefore be set accordingly:
def write_control_dict(case, n_iter):
"""Sets up the control dictionary.
In this example we use the rhoCentralFoam compressible solver"""
# Control dict from tutorial
control_dict = {
'application': 'rhoCentralFoam',
'startFrom': 'startTime',
'startTime': 0,
'stopAt': 'endTime',
'endTime': n_iter,
'deltaT': 0.001,
'writeControl': 'runTime',
'writeInterval': 1,
'purgeWrite': 0,
'writeFormat': 'ascii',
'writePrecision': 6,
'writeCompression': 'off',
'timeFormat': 'general',
'timePrecision': 6,
'runTimeModifiable': True,
'adjustTimeStep' : 'no',
'maxCo' : 1,
'maxDeltaT' : 1e-6,
}
with case.mutable_data_file(FileName.CONTROL) as d:
d.update(control_dict)
The mesh needs to be set up using the blockMeshDict. The mesh consits of three blocks in order to model the upstream and downstream portions as well as the wedge itself. The numbering order in which the vertices are set is very important! We declare a block via:
hex', [0, 7, 2, 1, 8, 15, 10, 9], [40, 40, 1],
'simpleGrading', [1, 1, 1]
The first line declares which vertices make up the corners of the block. This explanation best describes the order in which the vertices should be listed. The second part of the statement describes the cell density within the block in each of the three directions. The last part is used when we want the mesh density to vary within the block.
We must also set the thermodynamic properties of the gas. In this case the properties have been chosen so that the gas has a ratio of specific heats of 1.4 and that, if the temperature is 1K, then the speed of sound is 1m/s. As this is a commonly used fluid it can be done using the write_standard_thermophysical_properties function in the firefish.case module. We do this via:
write_standard_thermophysical_properties(case, StandardFluid.DIMENSIONLESS_AIR)
As this is a compressible flow we must also set the initial value of the temperatue field. Notice also that for the velocity we have set a slip boundary condition on the solid walls. This is because we are using an inviscid solver. When we move to a viscous solver we must set a no slip boundary condition on the solid walls.
def write_initial_conditions(case):
"""Sets the initial conditions"""
# Create the p initial conditions
p_file = case.mutable_data_file(
'0/p', create_class=FileClass.SCALAR_FIELD_3D
)
with p_file as p:
p.update({
'dimensions': Dimension(1, -1, -2, 0, 0, 0, 0),
'internalField': ('uniform', 1),
'boundaryField': {
'inlet' : {'type' : 'fixedValue', 'value' : 'uniform 1'},
'outlet': {'type': 'zeroGradient'},
'fixedWalls': {'type': 'zeroGradient'},
'frontAndBack': {'type': 'empty'},
},
})
# Create the U initial conditions
U_file = case.mutable_data_file(
'0/U', create_class=FileClass.VECTOR_FIELD_3D
)
with U_file as U:
U.update({
'dimensions': Dimension(0, 1, -1, 0, 0, 0, 0),
'internalField': ('uniform', [2, 0, 0]),
'boundaryField': {
'inlet' : {'type' : 'fixedValue',
'value' : ('uniform', [2, 0, 0])},
'outlet': {
'type': 'zeroGradient'
},
'fixedWalls': {
'type': 'slip'
},
'frontAndBack': {'type': 'empty'},
},
})
# Create the T initial conditions
T_file = case.mutable_data_file(
'0/T', create_class=FileClass.SCALAR_FIELD_3D
)
with T_file as T:
T.update({
'dimensions': Dimension(0, 0, 0, 1, 0, 0, 0),
'internalField': ('uniform', 1),
'boundaryField': {
'inlet' : {'type' : 'fixedValue', 'value' : ('uniform', 1)},
'outlet': {
'type': 'zeroGradient'
},
'fixedWalls': {
'type': 'zeroGradient'
},
'frontAndBack': {'type': 'empty'},
},
})
The example script includes a main()
function which performs all of
these steps. A boolean value can be passed to main()
in order to reduce the number of iterations
and so speed up automatic testing.
def main(case_dir='wedge', n_iter=10):
#Try to create new case directory
case = create_new_case(case_dir)
# Add the information needed by blockMesh.
write_control_dict(case, n_iter)
write_block_mesh_dict(case)
#we generate the mes1h
case.run_tool('blockMesh')
#we prepare the thermophysical and turbulence properties
write_standard_thermophysical_properties(case, StandardFluid.DIMENSIONLESS_AIR)
write_turbulence_properties(case)
#we write fvScheme and fvSolution
write_fv_schemes(case)
write_fv_solution(case)
write_initial_conditions(case)
case.run_tool('rhoCentralFoam')
After the example script is run, paraFoam may be run to inspect the solution.
Snappy Hex Mesh Example¶
The snappy_hex_example.py
script provides an example of running snappyHexMesh.
In order to be able to run snappyHexMesh we need to set up a control dictionary even though it plays no part in the actual mesh generation process. Likewise, in order to use paraFoam, we need to set up fvSchemes and fvSolution.
For snappyHexMesh to work we must have an underlying block mesh. This is generated in make_block_mesh and follows the same procedure as in previous examples.
For the actual mesh generation we first of all load a geometry using the new Geometry class:
rocket = Geometry(GeometryFormat.STL,'example.stl','example',case)
The idea behind this class is that, when we support more geometry file types in the future, it will abstract away the need to worry about wether something is an STL or OBJ etc.
We next scale and translate the rocket:
rocket.scale(0.5);
rocket.translate([0.5,2,2])
The Geometry class also contains mesh quality settings for this particular geometry.
Now that the geometry has been loaded we use it to initialise the SnappyHexMesh class:
snap = SnappyHexMesh(rocket,4,case)
This creates a new SnappyHexMesh class based on the example geometry and with a surface refinement level of 4. The SnappyHexMesh class automatically sets the mesh generation settings to a set of default values. These can be altered:
snap.snap=True
snap.snapTolerance = 8;
snap.locationToKeep = [0.0012,0.124,0.19]
snap.addLayers=False
Once the mesh generator is set up we can make the mesh via the call:
snap.generate_mesh()
Several things happen all at once here: * Surface features are extracted from the geometry (saving the STL in the trisurfaces directory if it has not already done so) * The mesh quality settings are written to a dictionary * The snappyHexMesh dictionary is written using the attributes of the SnappyHexMesh class * snappyHexMesh is run as a tool within the case directory
Once this has run the mesh can be viewed via paraFoam
Join STL Example¶
The join_stl_example.py
script
provides an example of combining multiple STL files into a single geometry and then
generating a mesh through snappyHexMesh. It is worth having a look at
snappy_hex_example.py
first in order to
get a more detailed overview on how SnappyHexMesh works.
It is now very straightforwards to generate a mesh made up from multiple .STL files.
Firstly one needs to make a list containing the paths of each STL file:
path_list = ['STLS/nosecone.stl', 'STLS/upperTube.stl', 'STLS/lowerTube.stl', 'STLS/finCan.stl', 'STLS/tail.stl']
One then needs to make a list of the human-readable names that correspond to each file:
part_list = ['nosecone', 'upperTube', 'lowerTube', 'finCan', 'tail']
When examining the mesh in paraFoam or when getting force outputs it will be these names that appear.
Once this has been done the Geometry classes can be loaded by a single call of:
parts = load_multiple_geometries(GeometryFormat.STL,path_list,part_list,case)
This produces a list of firefish.geometry.Geometry objects which can scaled, translated and rotated independently as required using the normal geometry functions.
We now use this list of Geometry objects to initialise SnappyHexMesh:
snap = SnappyHexMesh(parts,4,case)
As in snappy_hex_example.py
, we can now
alter the settings of Snappy Hex Mesh by altering the attributes of the SnappyHexMesh class. Once
we’ve updated these we generate the mesh via:
snap.generate_mesh()
This generates four different meshes: the blank block mesh, the castellated mesh, the snapped mesh and a mesh with layer addition. In order to use the final mesh as the starting point of our simulation we perform some trickery to delete the meshes we don’t want and move the mesh we do want into the constant folder
def getTrueMesh(case):
#the proper mesh is in the final time directory, delete the one in constant
os.chdir(case.root_dir_path)
call (["rm", "-r", "constant/polyMesh"])
call (["mv", "0.002/polyMesh", "constant/"])
call (["rm", "-r", "0.001/"])
call (["rm", "-r", "0.002/"])
os.chdir("../")
Kinematics Example¶
The kinematics_example.py
script
uses the new firefish.kinematics framework to implement the kinematics example in kinematics_basis.py
We firstly define a class that inherits from firefish.kinematics.KinematicBody. We want to model a cylindrical rocket whose principal moments of inertia vary as the rocket burns its motor. In order to vary the prinicpal moments of inertia automatically during the timestepping routine we must overload the update_moi() function.
class CylinderRocket(KinematicBody):
"""A rocket that is a cylinder with evenly distributed mass"""
def update_moi(self):
radius = 0.3
height = 2
Ixx = (self.mass/12.0)*(3*radius**2 + height**2)
Iyy = (self.mass/12.0)*(3*radius**2 + height**2)
Izz = self.mass*radius**2 / 2.0
self.MoI = [Ixx,Iyy,Izz]
In the main function we now undergo the time stepping routine. For each time step we must pass the forces acting on the rocket along with the torques to the routine. These forces and torques must be given in the body coordinate system. The example here burns the motor for fifty seconds and then lets it coast
def main():
"""Run through the simulation with a 50s motor burn"""
initialMass = 100
initialInertias = [0, 0, 0]
rocket = CylinderRocket(initialMass,initialInertias)
dt =1; duration = 100; gravity = 9.81;
simulation = KinematicSimulation(rocket,gravity,duration,dt)
while simulation.tIndex*dt <= duration:
thrust = 0
if simulation.tIndex*dt <= 50:
thrust = 2000.0
forces = [2.0, 0.0, thrust]
torques = [0.0, 0.0, 0.0]
mdot = 0.1
i = simulation.tIndex
simulation.time_step(forces,torques,mdot)
return simulation.posits[100,2] #z position