Here’s a quick tutorial on interFoam solver from the OpenFOAM toolbox. It demonstrates how much mess one makes when sneezing close to a bowl of soup.
The Case
Every time I sneeze during lunchtime my imagination runs wild with images of soup splattering all over the place. It wouldn’t be very polite to try it in real life so I decided to try sneezing virtually.
That means directing a stream of air of high velocity into a liquid at a standstill at a very close range and seeing what happens.
The hard part is finding out how fast a sneeze really is. There’s a bunch of references quoting a common urban myth that says it reaches speeds of over 160 km/h. But there’s also a single scientific paper stating that no more than 17 km/h can be squeezed out of a human. To be honest, I really don’t care how fast a sneeze is as long as I make a huge virtual mess out of mine. It turns out that in ideal conditions a 20 km/h sneeze wouldn’t make a lot so I happily used the myth numbers, 150 km/h more exactly.
The Geometry
I very carefully modeled a soup bowl from my china set with precise caliper measurements. FreeCAD was the program of choice. I created two parts, one for the bowl and one for the soup at a standstill. As the geometry is rather simple, I meshed both parts with Mesh workbench and exported both to STL (for more complicated geometry I would most probably use SALOME for meshing and face assignment). I’m used to ASCII STL format – I added a .ast when exporting and then renamed it back to stl. Also FreeCAD’s default export units are millimeters or whatever your settings are – you can change them or -in my case- use the surfaceFeatureConvert utility. I also named my solids by manually editing STL files; in the line that states ‘solid’ I added a name: solid soup. It doesn’t sound right, though.
The Mesh
Looking into makeSnappyMesh script reveals the meshing process:
Copy the saved geometry and scale it from millimeters:
surfaceConvert -scale 0.001 ../geometry/plate.stl constant/triSurface/plate.stl surfaceConvert -scale 0.001 ../geometry/soup.stl constant/triSurface/soup.stl
Surface features: there are no sharp features on the plate so this is actually not needed;
it’s included here anyway in case the geometry should be updated later
surfaceFeatureExtract
blockMesh contains the computational domain; in this very case, the plate is not entirely included.
When the wave returns it leaves it and that’s visible as a sharp end near the sneezer’s mouth.
blockMesh
snappyHexMesh is configured to make a lot of cells near plate walls. And to run in parallel. My i7 processor has 4 physical cores and 8 logical with hyperthreading. It turns out that using 4 cores is faster than using 8.
There’s a special refinement zone called inlet. It will refine a small rectangular region which will then be captured by topoSet and used as an inlet patch.
It makes sense to run checkMesh on a decomposed mesh and then reconstruct it. Instead of scrolling way back to checkMesh‘s output, save it to a file and display it at the end.
decomposePar -force mpirun -np 4 snappyHexMesh -parallel -overwrite mpirun -np 4 checkMesh -parallel > log.checkMesh reconstructParMesh -constant
topoSet will capture the refined inlet and create a cellSet. createPatch will – guess what – create a patch from that set.
topoSet createPatch -overwrite
Display the output and clean up
cat log.checkMesh rm log.checkMesh
Boundary Conditions
I had a lot of ideas how boundary conditions should really be stated but in the end it turns out that tried and tested ones from tutorials work best. Namely, tutorials/multiphase/interFoam/RAS/damBreak. There’s just one not-that-common-one, that is inlet velocity. It uses a tabulated value that changes over time. It’s also not perpendicular to inlet surface but it just means the numbers are a bit different. I simply made up sneeze durations and velocity profile. It doesn’t have to be realistic, it just has to make a big splash.
Initial conditions
There’s soup in the bowl. As stated earlier, I prepared geometry for the bowl and another one for the soup. setFields‘ surfaceToCell uses the soup STL surface to find all cells that contain the soup initially. It is also worth noting that it will replace the original alpha.water file with the new initial conditions so it’s wise to save your original file as alpha.water.orig.
Dynamic Mesh refinement
From OpenFOAM tutorials:
The two-phase algorithm in interFoam is based on the volume of fluid (VOF) method in which a specie transport equation is used to determine the relative volume fraction of the two phases, or phase fraction α, in each computational cell. Physical properties are calculated as weighted averages based on this fraction. The nature of the VOF method means that an interface between the species is not explicitly computed, but rather emerges as a property of the phase fraction field. Since the phase fraction can have any value between 0 and 1, the interface is never sharply defined, but occupies a volume around the region where a sharp interface should exist.
Interface in this case means water level or more precisely, soup level. Therefore, smaller cells more precisely define the location of the soup. That means a zillion cells in the mesh but only a small number of interesting ones – the ones that contain the interface or soup. It would be best to refine mesh dynamically based on how much soup each cell contains. Unfortunately, using dynamicRefineFvMesh in constant/dynamicMeshDict causes the solver to crash without stating what exactly went wrong. I’m still finding out how to solve this problem so any suggestions are most welcome.
Running
It will take a lot of time to solve the first half of a second since velocities are high but Courant number should be kept low. As the sneeze calms down all that’s left is splashing of soup which is a much slower process. The rest 2.5 seconds will take approximately the same computational time.
Postprocessing
Since simulation results are not very precise they should at least look nice. I imported an STL head model (that has no physical meaning and no mention in the simulation) and placed it approximately where the inlet should be. I used a plane source in ParaView to model a table and imported the same STL bowl I used in snappyHexMesh.
When displaying soup interface with coarser mesh (such as this one) the Contour filter makes horrible surfaces. They can be smoothed out by using the -guess what- Smooth filter. The larger the number of iterations, the smoother the soup. In my case it’s an asparagus soup so it should stay as smooth as possible.
Also, some additional light tweaking (View > Light Inspector) helps a lot. For the final animation I checked the Enable OSPRay for a more realistic rendering.
I saved every hundredth of a second of simulation time. That accounts for 300 frames which give 20 seconds of video at 15 frames per second. Not to say using compression for result files makes sense…
The Files
If you want to mess around with your interFoam, you should check the attachment. In it you’ll find the FreeCAD files, STL geometry, two OpenFOAM cases (one for meshing and one for simulation), scripts that run all required steps and ParaView state for the whole postprocessing bundle. I can’t attach 5 GB of zipped result files, though.
Have fun and watch your table manners.