An engineer's pipe flow
From Gerris
Revision as of 08:34, 12 February 2008 Traumflug (Talk  contribs) (Finishing the textbook case.) ← Previous diff 
Revision as of 10:32, 12 February 2008 Traumflug (Talk  contribs) (small syntax error) Next diff → 

Line 273:  Line 273:  
== About mesh refinement ==  == About mesh refinement ==  
  The examples used on this page use a <code>[[GfsRefineSolid]] = 7</code> to adjust the mesh roughly to geometry. The drawback with this is, around simulation time t = 1.5 the simulated flow with Re = 2485 "explodes". Obviously, numerical problems are the culprit and can't handle the developing turbulence.  +  The examples used on this page use a <code>[[GfsRefineSolid]] 7</code> to adjust the mesh roughly to geometry. The drawback with this is, around simulation time t = 1.5 the simulated flow with Re = 2485 "explodes". Obviously, numerical problems are the culprit and can't handle the developing turbulence. 
[[Image:Doppelbogenwatervelocityexploded.pngthumbleft"Explosion" at about t = 1.5 (velocity)]]  [[Image:Doppelbogenwatervelocityexploded.pngthumbleft"Explosion" at about t = 1.5 (velocity)]]  
Line 280:  Line 280:  
<br clear=all>  <br clear=all>  
  To solve the issue, I reran the same simulation with <code>[[GfsRefineSolid]] = 8</code>. Now, the simulation easily passes the t = 1.5 mark. Drawback is, it does so after a CPU time of no less than 140,000 seconds (about 38 hours). This is an about 20fold CPU work increase for a single refinement level.  +  To solve the issue, I reran the same simulation with <code>[[GfsRefineSolid]] 8</code>. Now, the simulation easily passes the t = 1.5 mark. Drawback is, it does so after a CPU time of no less than 140,000 seconds (about 38 hours). This is an about 20fold CPU work increase for a single refinement level. 
[[Image:Doppelbogenwatervelocityrefined.pngthumbleftRefined simulation at t = 2.0 (velocity)]]  [[Image:Doppelbogenwatervelocityrefined.pngthumbleftRefined simulation at t = 2.0 (velocity)]] 
Revision as of 10:32, 12 February 2008
All these oceans with their tides and ships cuising in it are impressive, aren't they? Did you notice how stiff and probably refreshing the breeze around MS Tangaroa is?
Well, as an ordinary mechanical engineer I often have much less impressive things to work out. Like water taps, curved pipes or some sort of flow meter. These smallsized internal flows are, what this tutorial is about. I'll take an example on a simple pipe with two bows in it.
In this tutorial I'll explain how to get more or less arbitrary geometry of internal flows into the Gerris flow solver and how to compute the fluid's behavior there. I'll discuss on what's all this fuss with the Reynolds Number is and how it translates to simple engineering tasks. Finally, I'll do some basic measures in the results and will compare against what a textbook on hydraulics considers as the correct solution.
This tutorial was initially written by Markus Hitter at jumpING in February 2008
Computer Aided Design
Undoubtly, different engineers have different tastes about what a good software package for doing CAD is. So, I'll explain the steps common to all modelers, only.
The few steps mentioned here sound simple, but can be a lot of work, of course:
 Create a model of where in your design the fluid flows.
 Create a cube of arbitrary size to cover this fluid model.
 Substract the modeled fluid flow from this cube.
 Make sure the inlet and outlet of this flow reaches the cube's sides.
The result might look like this:
Additional requirements for your CAD work:
 All of the cube's edges should be of equal size (or it wouldn't be a cube).
 It's only one inlet or outlet per cube side allowed.
 Avoid trapped volumes as it's a waste of computing time.
If you have this, save the geometry in STL format. For the sample scripts in this tutorial I'll use doppelbogen.stl
as the file name.
From CAD geometry to Gerris geometry
This step is quick and pretty straightforward as it's mostly a single step to resize the geometry to fit Gerris' fluid domain and to convert from STL to GTS format. Gerris comes along with all the tools needed.
The only caveat is, you have to make the cube slightly larger than the fluid domain as many CAD applications write triangualted data pretty exact, but not exact enough to be unambiguous. Small errors at the 15th digit behind the decimal point are common.
The transform
tool has a handy normalize
option which does resizing and relocating in one step. All put together, a single shell line to do the conversion and transformation looks like this:
tr d '\r' < doppelbogen.stl  \ stl2gts  transform revert normalize  \ transform scale 1.0001 > doppelbogen.gts
The commands in detail:

tr d '\r'
is put in to convert DOS line endings to Unix line endings, e.g. if your CAD runs on Windows. It doesn't hurt if the STL already has Unix line endings. If you ever happen to use a traditional Macintosh text editor, use atr '\r' '\n'
instead. 
stl2gts
does, as the name already hints, the conversion from STL format to GTS format. 
transform revert normalize
handily makes the geometry fit into the minimal Gerris fluid domain, a single GfsBox. For therevert
see next paragraph. 
transform scale 1.0001
oversizes the cube slightly (see above).
Depending on your CAD, the revert
option in one of the transform
commands might be needed to switch the triangles "inside out". The fastest way to find out wether this option is needed or not is to just proceed with it and to redo this step in case you run into calculation errors later. I'll refer to this when explaining how to calculate the flow.
As all commands read from standard input and write to standard output, there's no need for intermediate files.
The result of this step is a file with the name doppelbogen.gts
.
A simple calculation script
There are plenty of examples and descriptions on how to put a Gerris calculation script together, so I'll discuss only few topics here. First, have a look at what I usually use to get started:
1 0 GfsSimulation GfsBox GfsGEdge {} { GfsTime { end = 1 } GfsRefine 4 GfsRefineSolid 7 GfsSolid doppelbogen.gts # The Reynolds number is Re = L * U / Nu. GfsSourceViscosity 1e3 GfsOutputTime { istep = 1 } stdout GfsOutputScalarStats { step = 0.1 } stdout { v = Vorticity } GfsOutputScalarStats { step = 0.1 } stdout { v = Velocity } GfsOutputScalarStats { step = 0.1 } stdout { v = P } GfsOutputTiming { start = end } stdout GfsOutputSimulation { step = 0.1 } result%2.1f.gfs } GfsBox { left = GfsBoundaryInflowConstant 1 right = GfsBoundaryOutflow }
Besides a bunch of GfsOutput... stuff, the experienced Gerris user will notice the following:
 The simulation domain is set to a single cube. As you've seen in the chapter before, there's a quick way to make a cube designed in CAD to fit into this single GfsBox. The only drawback is, Gerris insists (for good reason) on calculating a single GfsBox on a single processor. The advantage is, the other processor of your dual core chip will still be available for your entertainmant.
 The fluid domain is set to have an inflow boundary condition with fluid velocity = 1 to the left and an outflow boundary condition to the right. These two are the cube sides where our modeled pipe starts and ends.
 The other four sides of the cube need no boundary conditions as they are fully covered by our solid.
 The initial grid is with a
GfsRefine 4
rarther coarse (remember how tiny our pipe is). To adapt solid geometry better, the grid is refined with aGfsRefineSolid 7
. As we expect most of the vortices to be near the solid walls, this is fine.  There is no
GfsAdaptVorticity
, as this is easily a source of trouble, at least for beginners.  With
GfsSourceViscosity
, the viscosity of the fluid is set to something. I'll be back on this when discussing physical properties.
So, fire up your favourite text editor and copy the script above into a file with the name calculation.gfs
. Save the file in the same directory as the doppelbogen.gts
file.
The result of this step is a file with the name calculation.gfs
.
Calculating and viewing results
Ready to enter the number crunching performance test? If so, fire up a Terminal, go to the directory containing the calculation.gfs
and doppelbogen.gts
file and enter:
gerris3D calculation.gfs
That's all!
As of this writing, the calculation for the fluid's time from t = 0 to t = 1 takes about 40 minutes of wall clock time on my fastest CPU. Luckily, we've set the calculation script to write out result files with GfsOutputResult a few times in between, so you can launch GfsView3D
and look what Gerris is doing (and wether it is doing what you intend it to do) a lot earlier.
As Gerris continues to calculate the flow, play around with GfsView, click the buttons available in the toolbar, look which mouse button does what in the main view. The timestep viewed is set by which of the results file you open, all other variants shown (pressure, velocity, turbulence, ...) can be changed by doubleclicking on one of the items in the panel to the left.
If you think your latest results file didn't change significantly from the one before, stop the calculation by hitting ControlC in the Terminal. If you'd ever change your mind later, you'd simply restart the simulation from where you left:
gerris3D results<latest number>.gfs
One error commonly hits beginners: Gerris computes the "wrong side" of the solid. You easily spot it in GfsView: the pipe where in the fluid should flow is black and everywhere a concretelike solid should be shown colorful fluids appear: Gerris is actually computing something like a thick wire in the wind.
Another possible result of this mismatch is Gerris refusing to start doing fluid calculations at all:
GfsERROR **: root cell is entirely outside of the fluid domain the solid surface orientation may be incorrect aborting...
The reason for this is, different STL modeling software has different conventions about where the "inside" of a triangle representing a solid surface fraction is. The simple solution is, go back to the step converting geometry from STL to GTS and remove the revert
option. Insert it if you didn't use it.
Once you found out wether your CAD package needs the option or not you can stick with your choice for all other STL geometry exported from there.
Is it water or is it air: physical properties
So, now you can almost feel how the fluid flows through your pipe. If you followed the Tips and Tricks section, you've perhaps already some movie showing the development of the flow over time.
But wait! What is this fluid? Is it air? Is it water? Is it pasty like honey? Sure, you know the size of your geometry in CAD, but didn't we scale all this by some magnitudes?
To get some light in this confusion, some very knowledgeable people invented a characteristic number to compare flows against each other: the Reynolds Number. The simple rule is:
 same Reynolds number = same fluid behaviour
So, all we have to do is to calculate the Reynolds Number for our real geometry and to make sure our model features the same number. The number, as given in literature, is:
 Reynolds number = fluid velocity * characteristic length / kinematic fluid viscosity
 Fluid velocity is a number you should already have for your real flow. In our example calculation, the model fluid velocity is set to 1. There's only few reason to choose something different from 1 for the model.
 Characteristic length is some length of the geometry. For pipes, typically the pipe's diameter is used. Other dimensions work just as fine. For the real flow, you can easily measure a characteristic length. For the modeled flow, it's the same length, multiplied with the scaling factor used when converting CAD geometry to Gerris geometry.
 Kinematic fluid viscosity should be known for the real flow as a property of the matter inside your pipe. In the Gerris model, we use it to adjust the Reynolds number to match the real flow.
All this put together, we have to calculate the model's fluid viscosity by the following formula:
 model viscosity = model velocity * real viscosity / real velocity * model length / real length
 = model velocity * real viscosity / real velocity * scaling factor
The same with numbers for our example geometry, water as a fluid and an real fluid velocity of 0.5 m/s:
 model viscosity = 1 m/s * 1.006*10^{6}m^{2}/s / 0.5 m/s * 1 m / 0.050 m
 = 4.024 * 10^{5}m^{2}/s
 Nota Bene: In literature you'll often see the term "dynamic viscosity", expressed in N*s/m^{2} or Pa*s. This is not to be confused with the "kinematic viscosity" we need here. However, both are related by density:
 kinematic viscosity = dynamic viscosity / density
Fine. Let's insert this number into our calculation script:
[...] GfsSourceViscosity 4.024e5 [...]
Rename the former results file for later access, then run the calculation again. Expect calculation time to almost double from before.
While the CPU converts electricity into crunched numbers, we can do a comparison of what we have done so far. The pipe's diameter, measured in the fluid domain size, is 0.1. With our intitial, more randomly choosen setting of viscosity we had a Reynolds number of
 Re = 1 * 0.1 / 1e3 = 100
After reviewing viscosity and choosing it to fit water at 0.5 m/s, we have:
 Re = 1 * 0.1 / 4.024e5 = 2485
Some time later, we can view the difference in pictures:
Range blue to red = 0.0 to 2.0.
Range blue to red = 0.0 to 200.0.
Range blue to red = 1.0 to 4.0.
Of course, each pair of the snapshots was taken with the same scaling. The differences are obvious.
Back to reality: measuring results
Often, having a few colorful pictures is nice, but concrete numbers are more meaningful. If you know where to measure from the beginning or if you want to manifest numbers over time in a results file, there's GfsOutputLocation
, for example:
GfsOutputLocation { istep = 1 } stdout 0.45 0.2 0.0
The drawback of this method is, you have to (re)start a calculation each time you change your mind. If you have a results file already, there's another, more interactive method to do some measurements:
 Open the results file of the timestep you want to measure in GfsView3D.
 Click "Linear" in the toolbar.
 Doubleclick the entry to the left which just appeared.
 In the "Linear" tab, select the type of measurement you want to make (pressure, velocity, ...).
 Change the colorscaling or leave it as is, it doesn't matter.
 In the "2D Plane" tab, select the Zcoordinate you want to measure (or adjust the plane to your needs).
 Controlclick into the fluid domain in the main view and read coordinates as well as measurement at the bottom of the window.
Taking concrete measurements is easy, isn't it? Did you wonder why I didn't use dimensions with results until now? Well, the numbers given are in Gerris' unity domain and before attaching a m/s
to the velocity number, we have to scale the number back to the reality's coordinate system.
Fortunately, scaling back is even easier than adjusting the Reynolds number shown above. If we still know the characteristic lengths and the reference fluid velocities (I called them just "velocity" above) of both, model and real flow, the following formulas (with values taken from a random place of our example) do it:
 real velocity = model velocity / model reference velocity * real reference velocity
 = 1.3 / 1 * 0.5 m/s = 0.65 m/s
as well as
 real vorticity = model vorticity / model reference velocity * real reference velocity * model characteristic length / real characteristic length
 = 25.2 / 1 * 0.5 m/s * 1 / 0.05 m = 252 s^{1}
For pressure, we additionally need the fluid's density:
 real pressure = model pressure * real density * (real reference velocity)^{2}
 = 0.74 * 1000kg/m^{3} * (0,5 m/s)^{2} = 185 N/m^{2} = 185 Pa
Not that surprising, time is scaled as well:
 real time = model time * real characteristic length / model characteristic length / real reference velocity
 = 1 * 0.05 m / 1 / 0.5 m/s = 0.1 s
The latter result means, the whole example simulation represents just the tenth of a second of the real world. Not much for over an hour of computing time, but try to think how long it would take to do such measurements in the laboratory ...
What the textbook says to such a case
When researching for this tutorial, I found a good scriptum from the Austrian Institute of Hydraulics and Rural Water Management, Vienna, in the Net. The following snippet discusses losses in bowed pipes:
Well, I hope your knowledge of german is sufficient to unterstand all the formulas ...
h_{r} is the "height pressure loss". The formula to convert it to a pressure difference is:
 Δp = ζ_{r} * ρ * v^{2} / 2
with:
 Δp = pressure difference
 ζ_{r} = drag coefficient (see table)
 ρ = density (~1000 kg/m^{3} here)
 v = mean fluid velocity (0.5 m/s here)
Looking at the textbook's snippet, our sample is most closely to a "Doppelkrümmer",so:
 ζ_{r} = 4.0 * (0,131 + 0,159 * (0.005 / 0.005)^{3.5}) * 90° / 90°
 = 4.0 * (0.131 + 0.159) = 1.16
 Remark: This is likely a border case for this formula.
With this ζ_{r}, we get a Δp of:
 Δp = 1.16 * 1000 kg/m^{3} * (0.5 m/s)^{2} / 2
 = 145 kg/m/s^{2} = 145 Pa
Discussion of the differences between Gerris and the textbook
About mesh refinement
The examples used on this page use a GfsRefineSolid 7
to adjust the mesh roughly to geometry. The drawback with this is, around simulation time t = 1.5 the simulated flow with Re = 2485 "explodes". Obviously, numerical problems are the culprit and can't handle the developing turbulence.
To solve the issue, I reran the same simulation with GfsRefineSolid 8
. Now, the simulation easily passes the t = 1.5 mark. Drawback is, it does so after a CPU time of no less than 140,000 seconds (about 38 hours). This is an about 20fold CPU work increase for a single refinement level.