Thursday, October 31, 2013

Promo 3D print

In these days there is a promo on  Shapeways 3D print site.

"Limited Time Only: New customers Save $15 on orders $75+ with promo code ghost "

Maybe it's a good occasion to print some pitot batch :-)) , look at this site

Wednesday, October 30, 2013

GPS altitude introduction

The global positioning system  is not usually included in the basic topics about air data measurement, nonetheless this instrument is useful to some auxiliary operations that should be carried during normal instruments operation.
Common approach is to log the altitude as above ground level, the departure runway is set then as the reference zero altitude above ground level. Same approach is used with barometer altimeters.
Sometimes this procedure is carried out in an inaccurate way and both the GPS and the barometer are unable to provide mean sea level altitude, not to mention the heavy impact on barometer measures accuracy that can be caused by a coarse, or absent, local weather compensation.

GPS works natively with a ECEF reference frame, the chipset converts this coordinate system to those of a selected ellipsoid. As example I will use a Venus638flpx GPS chipset by Skytraq, get the datasheet here , at page three you can read that GPS datum is by default WGS-84 . This application note at page 44 , indicates that Venus chipset can manage many reference ellipsoids.

WGS-84 provides an approximate model of earth surface, that is modelled as an ellipsoid. Returning to altitude measurements they are refered to the geoid datum. This reference system define a zero elevation surface, the geoid, that would coincide with the mean ocean surface of the Earth if the oceans and atmosphere were in equilibrium. Consider the following figure that depict current altitude reference systems as introduced until now.
Figure F19.1 Altitude datum representation for a flying platform.
H GPS altitude, h altitude , N undulation, AGL above ground level altitude, MSL mean sea level altitude

In the figure H represents graphically the altitude measure provided by the example GPS unit, h is the altitude as globally defined. Deviation between H and h is N or undulation, such deviations are caused by the natural variations of gravity field that as shown \(\epsilon\) can also deviate the gravity vector. Note that geoid surface can either be over or under the WGS-84 surface.

By figure inspection mean sea level altitude can be calculated as \(h=H-N\) therefore there is a need for an accurate knowledge of the geoid altitude.

This information is available through gravity models as EGM 2008 edition by U.S. National Geospatial-Intelligence Agency, on their site there is freely downloadable software in theSoftware and Coefficients for WGS 84 Geoid Undulation Computations by Harmonic Synthesis” section of this page. This software can calculate for us the undulation N in meters.

Is to be noted that the value of altitude supplied by  GPS units can be of different types.
The altitude can be referred to selected ellipsoid, in our example WGS-84, this value don't need to be elaborated further. On the other side, by default on Venus units, altitude reading is compensated for the deviation from ellipsoid, this altitude reading is defined on the manual as MSL. Refer to this document page 8 for typical NMEA $GPGGA position fix message format.
The heigh above WSG-84 geoid H is needed so if our fix message is as follow
$GPGGA,utc,latitude,N,longitude,E,1,#sats,HDOP,57,M,40.3,M,,*crc
                                                                                   ^         ^
GPS altitude H is then 57,2+40,3=97,5 m
Using the coordinates of my RC Group airfield into the EGM08 model I get the following output “45.231386 9.882303 40.488”
First two rows are decimal latitude and longitude and the last is N in meters.
If my GPS altitude reading at the same coordinates is 97,5 m then the altitude, or at last elevation,
is 97,5-40,488=57,012 m MSL.
As alternative you can also use this online calculator, note that this calculator is implementing an older EGM96.
With the same coordinates the output is the following

<<Your Input Coordinates and GPS Height:

Latitude = 45.2313861111111° N = 45° 13' 52.99" N
Longitude = 9.88230277777778° E = 9° 52' 56.29" E
GPS ellipsoidal height = 97.500 (meters)
Geoid height = 40.137 (meters)
Orthometric height (height above mean sea level) =57.363 (meters)
(note: orthometric height = GPS ellipsoidal height - geoid height) >>
There is a difference of 0,351m due to different N value between the two calculations.

Overall measure accuracy is obtained by the sum of accuracy of ellipsoid model, geoid model and GPS measurement itself. Accuracy of WGS-84 and EGM08 are related to the geographic coordinates.

A strong conservative value for EGM08 accuracy is 0,7m at 99% of confidence. A full compilation of documents with accuracy considerations for different geographical areas can be found at this link, for example refer to table 4 page 143.
Returning to our example, the datasheet of Venus GPS page three indicates an accuracy of 2,5m CEP radius; that is an indication of 2D precision of the unit.
Using CEP and some assumptions on error distribution is possible to calculate spherical accuracy SEP. SEP is defined as the radius of sphere centered at the true position estimate in 3D with probability 50%.
Assuming a typical VDOP / HDOP 2, PDOP / HDOP = 2 lead to approximate SEP=5m, similar results here at page 4-5
Probability of 50% is a low value, since half of the measurements will lie outside this radius. The 99% Spherical Accuracy Standard is then calculated as \(SAS_{90}=2,2SEP=11m\)
In a worst case scenario and with a confidence level of 99% the uncertainty is \(u_{t}=11+0,7=11,7m\)
So with the cited GPS unit, with no particular measurement methods, the altitude should be reported as  \(57,0 (\pm 11,7)m\) with 99% confidence level; 94% of total altitude deviation comes from the GPS unit uncertainty.

Is good practice to check our instruments readings against well know referral points, they can be mainly of two kinds.
Landmarks altitude accuracy is variable, but usually they are easy to use. Most commons landmarks are airports or main city squares, is to note that every, even if really small, airports have a good indication of altitude since it is a useful information for aircraft altimeter calibration.
Geodetic networks can also be used for calibration, this kind of networks are used also for cadastral purposes and can reach a really good, let's say impressive, degree of accuracy. In Italy, for example, there is a network called IGM95,Italian counter part of EUREF89 or ETRF89,  by IGM , the approximative location of points can be found on the IGM website ; a non official compilation of IGM 95 points can be found here.
If you test your unit against a referral point you will have a clear feedback about your GPS units degree of accuracy. Note that the geoid deviation calculated internally by your unit is not necessarily equal to the here above calculated deviation, anyway the measurement will be off with both corrections. You should find instead that your measurements fall into the uncertainty interval calculated above.
If you plan to record flight altitudes at your flying site, for the best results, setup a referral point of  altitude. Is a good practice to record both, GPS altitude and GPS corrected altitude, in such a way you can correct your past measurements as soon as another more accurate referral point is available.
These are the premises needed to result interpretation of a next to come example on GPS measurements, of course I will use also a Venus GPS unit.

Monday, October 21, 2013

Barometric altitude primer part two


A precedent post introduced some math about atmospheric models, relevant issues have been pointed out and within this new post they will be solved.

Barometric altitude is wide used in general aviation aircraft, use of this kind of altimeter is standard and considered reliable, some of this reliability comes from the simplicity of the involved devices that in some cases are fully mechanical. In terms of relative accuracy the results are good, in the same airspace uncompensated altimeters readings from different aircrafts flying at the same altitude will be practically identicals, so it's a reliable instrument to separate airplanes flight paths.

Unfortunately some problem arise when it's time to land. An altimeter that follows model equations will report a wrong altitude when the plane is on the runway. This is due to the local weather conditions that lead to a different pressure at ground level.

For example consider a runway located at 0 m MSL(mean sea level), and with 101800Pa of reported sea level pressure, by the way that is a possible value for a good autumn day in Cremona. The ISA model sets pressure to 101325 at 0 m MSL, so our altimeter is reading 101800 Pa that corresponds to -39,5 m, that of course it's wrong.
To calculate this values launch barometricr2.sce and execute the following commands, that are included at the end of file.

<<-->function F=isaaltitude(x)
--> F=pressure_reading-p0isa*(1-0.0065*x/T0)^(g/Rair/0.0065);
-->endfunction
-->p0isa=101325
p0isa =
101325.
-->pressure_reading=101800
pressure_reading =
101800.
-->fsolve(0,isaaltitude)
ans =
- 39.465884
>>

To obtain a correct measurement it's possible to account for effective sea level pressure, try the following commands.

<<-->pressure_reading=101800;
-->p0isa=101800;
-->function F=isaaltitude(x)
--> F=pressure_reading-p0isa*(1-0.0065*x/T0)^(g/Rair/0.0065);
-->endfunction
-->fsolve(0,isaaltitude)
ans =
0.
>>

If our altimeter is setup with the former method and the airfield is located at 150 m of altitude, on the runway the altimeter reading will be 150 m.
For some flight activities, as soaring, is convenient to know altitude above ground level rather than respect mean sea level, pressure data is the same but the altitude information is used to zero the altimeter indication. The sea level pressure value is available trough control tower communication or from an automatic weather observation station. A classical instrument should be setup, with minor deviations, as per the following video.

V18.1 Kollsman altimeter setup procedure

Setup should be repeated on a regular basis during flight, or anyway when the pilot detects a change in weather condition. Every variation in outside, temperature pressure and humidity will impact on the altimeter reading, pilots developed some mnemonic rules to deal with this phenomenon, this rules are effectively commented in the following video.

V18.2 True altitude video

So if you are using a barometric altimeter to log your RC/FPV/other aircraft altitude you will need to take care of setup, pragmatically speaking if you flight in a stable weather day and for limited time intervals and you are happy with a non optimal degree of accuracy it is sufficient to setup your altimeter just once prior to take off. If you are willing to push to the edge of measurement state of the art than you need to send, on a regular basis, the airfield pressure, corrected to sea level, to the flying platform.

Above described procedures are quite standard, on the other hand in many digital instrumentation the accuracy is pushed further, with supplementary measurements of temperature and humidity and the use of more sophisticated models the overall uncertainty is reduced. I'll show how our simple model can take into account temperature deviation from ISA atmosphere. Here a link to a Scilab example of altitude calculation that compensates for ISA temperature deviation.

Air temperature is relevant because modifies the air density profile, following calculation treat two identical pressure readings of \(P_{1}=P_{2}=100129 Pa\) taken at two different temperatures \(T_{1}=14,35°C\) and \(T_{2}=34°C\). Pressure at sea level is the standard value of 101325 Pa in both cases.

Barometric altitude of P1 reading can be calculated as

<<-->exec('isatemp.sce', -1)
Input static pressure reading Pa : 100129
Pressure height 100.04 m >>

A check for standard thermal distribution

\(T_{sea}=T_{1}+altitude0,0065\) substituting our values leed to 15°C=14,35+0,65 that is exactly the temperature at sea level for ISA model. Air density \(\rho_{1}\) is then 1,213 \(\frac{kg}{m^3}\).

Now that the temperature profile have been checked it's safe to assume the ISA altitude of 100 m.

Focusing on \(P_{2}=P_{2}\) altitude

<<-->exec('isatemp.sce', -1)
Input static pressure reading Pa : 100129
Pressure height 100.04 m >>

The altitude is the same, it's common to say pressure altitude is the same

Check of termal distribution

\(T_{sea}=T_{2}+altitude0,0065\)

substituting our values lead to 34,65°C=34+0,65 that is 19,65°C higher than ISA expected temperature so we denote it as ISA +19,65°C or rounded to unit ISA+20°C
Test failed, the measurement of pressure is not carried out in ISA conditions, some form of compensation is needed.

Air is an ideal gas here hence \(P_{1}= \rho_{1}R_{air}T_{1}=P_{2}= \rho_{2}R_{air}T_{2}\) then \(\frac{T1}{T2}= \frac {\rho_{2}} {\rho_{1}}\) \(\rho_{2}=\rho_{1}\frac{T1}{T2}\) substituting \(\rho_{2}=\frac {1,213(273,15+15)} {(273,15+34)}=1,13 \frac{kg}{m^3}\)

Recalculate the altitude using the standard formula with our sea level temperature

<<
->function F=isaaltitude(x)
--> F=pressure_reading-p0isa*(1-0.0065*x/T0)^(g/Rair/0.0065);
-->endfunction
-->T0=307.15
T0 =
307.15
-->pressure_reading=100129
pressure_reading =
100129.
-->p0isa=101325
p0isa =
101325.

-->fsolve(100,isaaltitude)
ans =
106.6349


So the pressure \(P_{2}\) do indicate a pressure altitude of 100 m identical to pressure altitude from \(P_{1}\) ), compensation for the different temperatures at sea level reveals that the second reading has been taken at 106,6 m MSL, to avoid notation confusion the second altitude is sometimes denoted as geometric altitude.
Between the non compensated reading and the compensated one there is a difference of 6,6 m that corresponds to a relative error of 6,2%.

Some different calculation methods have been considered, possible error sources have been introduced and basic numerical routines have been provided. The Scilab function fsolve have been used to shorten the solution but it is not usable directly on microcontrollers, in the next post will be presented a microcontroller like version of the barometric algorithm that use simpler form. Our analysis evidenced also the need for weather and altitude data at the airfield for an accurate altimeter operation.

Tuesday, October 15, 2013

Barometric altitude primer


This is the first post of a mini-series that treats altitude measurement, this series will cover from the basic theory to practical application and further to GPS sensor aiding.

The most used method for altitude measurement in airborne vehicles is by  mean of static pressure readings, this measure is referred as barometric altitude.

Working principle is simple and can be visualized recalling Stevin's law  for hydrostatic pressure calculation, as introductory example I will use the calculation of a bathtub water level.

Refer to the below figure for a section view of a bathtub. The Bathtub is filled with water at a temperature \(T_{W}\) of 40 degrees Celsius and current atmospheric pressure \(P_{atm}\) is 101325 Pa, our sensor is reading the pressure on the bathtub floor \(P_{F}\) that is 105000 Pa.



Figure 14.1 Section view of a bathtub

Water is considered still over all the tub volume, as the water pressure and themperature have a limited range of variation hence the density is considered constant and then \(\rho=f(T_{W})=992 kg/m^3\).


To calculate the pressure at the bathtub floor the definition of pressure is used,

Equation 14.1, definition of pressure

$$Pressure=\frac {Force over one square meter }{(m^2)}$$
Consider the force acting on the tub floor, the water layer close to the floor should withstand to the weight of all the water and air that lies on it.

The force on the floor ot \(FF\) is then is simply the force exerted by the weight of the water and air column. Pressure of air \(P_{atm}\) at the interface with water is known and equal to 101325 Pa hence using eq.14.1 \(P_{atm}=\frac {Force over one square meter }{(m^2)}\)


Regarding the water part of the column, given \(g\) as the gravity acceleration constant, \(h\) as the water column height and given \(S=one square meter\) as the column section hence
$$FF=P_{atm}S +(Weight of water for a column of one square meter of section S)g$$
$$FF=P_{atm}S + h\rho gS$$ 

therefore

$$PF=\frac {P_{atm} S+ h\rho gS}{S}= P_{atm} +\rho gh$$

So it is practical to write \(PF-P_{atm}=\rho gh\)
With the example data \(PF-P_{atm}=9929,81h=3675 \)
\(h=0,37 m\)
More appetible result is that the pressure \(P(z)\) at a determinated deep \(z\) is then
Equation 14.2
$$P(z)=P_{atm}+\rho gz $$
Calculate the atmospheric pressure at a given z altitude is exactly the same thing that calculate the pressure in a bathtub, it's necessary to calculate the force exerted by the over head air column weight.
To define how change the air density with altitude and have the chance to take this into account an explicit formulation for air density is needed.

To calculate the air density a widely used model are the model and the US 1962/72. For altitudes included in the troposphere, less than 11000m, they are equivalent, a huge number of altimeters use this models formulas to convert pressure readings in altitude indications.

To calculte the air density two common models are the ISA model and the US 1962/72. For altitudes included in the troposphere, less than 11000m, they are equivalent, a huge number of altimeters use this models formulas to convert pressure readings in altitude indications.
Here below a short list of model main assumptions, reference altitude is the mean sea level altitude.
a- latitude 45°
b-No wind or other forms of eddies,
c-Thermal gradient constant across the trophosphere of 6,5 °C each 1000 meter

d-Dry air, perfect gas \(Cp/Cv=1,4\)

Molecular weight, from the standard air composition is 28,9644 kg/kmol

R=8,314462175 J/°K/mol
e-Standard sea level conditions, reference conditions
\(T_{0}=288,15°K\)
\(P_{0}=101325 Pa\)
All the altimeters installed on aircrafts have the possibility to compensate for local pressure conditions, many are capable to more sophisticated compensations for temperature and humidity.

The atmosphere is by far idealized so no local weather conditions are taken into account, all the coefficients are given to fit at best 45° latitude standard conditions, humidity impact directly into the altitude reading values. Practically speaking the barometric altitude reading can be amazingly off when the aircraft is flying through clouds in a really warm day, tens meters of error should not surprise us cause the instrument is operating with a wrong density profile.
There are techniques to compensate for a,d and e assumptions.
All the altimeters installed on aircrafts have the possibility to compensate for local pressure conditions, many are capable to more sophisticated compensations for temperature and humidity.

Leaving to this reference the derivation of barometric formula is possible to introduce a numeric Scilab example, in this example two different calculation methods are employed; Initially the barometric, or pressure, altitude is calculated using a user specified external pressure and formulas from the Goodrich 4081 Air data handbook equation 4.1 then the calculated altitude is used in the ISA/US1962/72 pressure formula. If the implementation of the direct formula is the same of the inverse formula then the input pressure must be equal to ISA pressure.
You can run the example, and examine the calculation details, here is the scilab barometricr1.sce file, and verify that results confirm the equivalence between the two methods. Just test what is the Scilab output with a 105000 Pa input.

The atmospheric model calculates the following data about air
-\(\rho=f(h)\)
-\(T=f(h)\)
-\(p=f(h)\)
As example calculation result have a look to the following figure that shows density vs altitude, you can visualize the graph setting graph=1 into the example Scilab file and executing it.

Figure 14.2 Atmosphere density vs altitue ISA or US1962/1972 model.

The model is completed with the viscosity calculation, carried out with the Sutherland's formula, as the other values you get the viscosity using the supplied example Scilab.

With this data about the atmosphere you have everything for calculate, in a practical closed form, all the performances parameters of a fixed wing aircraft. It is so easy to have an initial evaluation of the performance of a particular airplane configuration, take off distance, climbing characteristics and in general every performance parameter inside the flight envelope.

Atmospheric models are a coarse approximation of reality, anyway real time information on pressure, temperature and altitude can be used to adapt the model to current climate conditions.Into the next post a practical approach to solving the issue will be introduced.

To the next article 

Further reading
2007 David G.Hull, Fundamentals of Airplane Flight Mechanics, Springer Chapter 3 



Thursday, October 10, 2013

Festo BionicOpter

Forward backward and stationary fly.
Breathtaking realization of dragonfly by Festo

More media here

Have fun,
JLJ 

Wednesday, October 9, 2013

Pitot solid model

Preliminary design phase is exciting as challenging phase of every project, different solutions should be evaluated and eventually one should be selected.

So called virtual prototyping techs have been used for decades by now and are of proven efficiency in the preliminary design phases, virtual prototype for 8 mm pitot is shown on following video, just below another example of virtual prototype render.

Video P13.1 Animation of pitot model file

Figure P13.1 Render of flange virtual prototype
Referring to a medium DIY size project and setting the goal of global cost reduction it's quite common to observe that virtual prototyping is the right tool also for the home DIY maker projects.

Correlated effects of a massive use of virtual prototyping include the possibility to archive, search, modify and share your project at best.

The base technique that is used in virtual prototyping is solid modelling.

Solid models can speed up work and rise the overall accuracy level. A solid model can give precise information, among others, about weight, center of mass and inertia. When flying platforms are involved it's often necessary to have such information to feed simulators or as foundation for control related tasks or if you plan to use telemetry data for aerodynamics parameters identification.

If it's needed to know if you new instrument will withstand to a particular airspeed you can at least attain an envelope figure of structural limits by using your solid model and a free FEM code as Elmer.

The selection of static ports location problem can be effectively treated, with CFD is possible to evaluate different locations and check sensibility to angle of attack or sideslip.

Leaving apart the poetic of flight testing to build a physical prototype and test it, by far slower and expensive than virtual prototyping. As shown by those simple examples also in DIY cases virtual prototyping mean efficiency and cost reduction.

I heard that sometimes the creative process is hindered by to continuously have to draw every new sub-part so for the 8 mm pitot you can find the .stl 3D model file here.

This model is ready to be incorporated in your project, only material personalization is needed.

Wednesday, October 2, 2013

Pitot parts 3D print



Have a look to this link, I've used this service and it's fine.

The 3D models files are free online available here.

Have fun!

Flanged pitot assembly instructions

Figure P13.1 Flanged pitot

After some testing finally you can find a full descriptionof of assembly proccess  here, added new video at the end of the page. Just let me know if you encounter some problem.

Have fun,
JLJ

Tuesday, October 1, 2013

CFD Pipe


Computational fluid dynamics is a versatile and powerful tool, nevertheless it's scarcely used in DIY applications. Nowadays there are some free CFD software at home maker disposal, the calculation power of personal computers is enough to reasonably cope with CFD design tasks.

There are many CFD commercial packages available, as well as a lot of well documented simulation cases, just to stay in topic have a look to the pitot simulation here . CFD commercial codes, and other needed software to complete the toolchain is not in typical DIY maker budget range, hence the numeric example will be carried out using the free and well proven package Elmer CFD .
At a glance, returning to our typical basic air data problem, CFD can calculate the velocity, pressure, temperature and density field for our given geometry and so derived standard defined aerodynamics parameters as the lift and the drag. With an adequate calculation power it's possible to estimate position error for a particular aircraft and instrument.
For the sake of simplicity the problem solution with CFD techniques is split into six main subproblems that are defined here below; after general concepts a numeric example is inspected.

Goal of the example is to compare analytical results with CFD results.

All the files needed for the example Elmer CFD can be downloaded here.

CFD main phases

Preliminary study phase


It is necessary to study the problem analytically, in such a way the sought test cases are individuated. With the preliminary study done it is possible to individuate the assumptions and simplifications that can be acceptable. For example a simplification can regard a geometric complexity reduction or the decision to neglect thermal transfer equations, both the cases can yield to a dramatically reduction of modelling, meshing and simulation time. Study phase is fundamental to check simulation results.

Geometric modelling phase

Here a 3D or 2D model of the case to be simulated is draw, no matter what software or format you use as long as it is easy to pass the data to the mesher software.

Meshing

Mesher must add simulation related information to our model, the model should be divided in regions and surfaces, for the mesh to be usable all the boundaries relevant to the physical problem must be identified. This features will be used during the simulation setup phase to define fluids, solids and boundary conditions properties. Mesher must also divide the simulation physical domain into a discrete grid, in a CFD problem the physical variables are calculated, for example, on the vertex of the grid. Unfortunately the accuracy of the simulation result is dependent on the mesh geometry itself. Choice of the grid dimension is not trivial and once again a preliminary study can give some information, for example if our simulation aim to study a boundary layer to have an estimation of layer thickness can aid to choice a proper grid dimension.

Please have a look to this site  to have a good introduction to the meshing.

Simulation setup and run

At this point you define in fine detail the properties of all the fluids, solids and boundaries and define what kind of solver you want to use during the simulation, incompressible flow, compressible flow, steady solution, dynamic solution heat transfer and so on. You must also define initial conditions, these are the values that variables take at the beginning of simulation, you can calculate or neglect such a values only if you carried out a good preliminary study phase.

Numeric calculation phase, usually stop when the variation of the solution across successive iterations is less than a threshold value specified during setup.
Data visualization and post processing
Here you can visualize, physical values as speed, pressure and density as vector and punctual values. If you want the value of some derived physical quantity you can extract data and compute your derived value, for example you can calculate the flowrate over a surface integrating the velocity associate to this surface.
Solution validation
After a result as been reach is common to run again the simulation using different meshes sizes and shapes. The results with different meshes should be the same, after that also a more selective constrain on convergence should be applied. It's to be noted that validation of the solution is reached if some analytical support is available or

Working example


Well, as CFD is a quite articulate topic, it's better to check general concepts introducing a simple example. The example consist in a pipe of 0,1 meter diameter and long 1 meter. A certain amount of air at 20°C and 101325 Pa of pressure is blowed into the pipe, with a constant speed over the inlet pipe section.
Preliminary study phase
As you know the interaction of air with the pipe walls will create a boundary layer, macroscopic effect is the velocity magnitude relation with the wall distance. As the flow pass through the pipe a variable speed profile along the pipe length is to expect, at least since the flow stream reach an Equilibrium.

Now that we have a well defined problem it's possible to use common available bibliography like this notes online.
Our CFD case will be identical to those depicted in reference figure 8-8, whereas a common problem was chose it's also available a closed form solution; a little investigation into the reference lead to focus on the equation 8-17 that state that fully developed velocity profile for a laminar flow into a circular pipe have a parabolic shape. In a fully developed flow region of our example pipe given pipe radius R=0.05 m, average velocity Vavg=0.1m/s and defined r as the position on radius the axial velocity u(r) can be writen as

Eq. P11.1
u(r) =2Vag(1-r^2/R)
substituting our parameter values 0.2(1-r^2/0.05^2) that yields for laminar flow so Re<2300 should be verified.
To calculate Reynolds number the below fluid properties at operating temperature and pressure are used
Name = "Air (room temperature)"
Reference Temperature = 20°C
v,Viscosity = 1.983e-5 Pas
Compressibility Model = Perfect Gas
Reference Pressure = 101325 Pa
rho, Density = 1.205 kg/m^3

Hence Re=rho*D*Vag/v=608, so it's safe to assume laminar flow condition.
Geometric modelling phase
Tridimensional model of the pipe is just a 1m cylinder with 0.1 m of diameter; you can use the really powerful and free of charge freecad. Generated model is then saved as .stp.format, like this.

After that is sufficient to  launch the ElmerGUI and just open the .stp file.

Meshing
Into ElmerGUI geometric viewer set the mesh max h to 0.1 and min to 0.005; after apply the mesh the cylinder looks like the figure here below.
Figure PF11.1 Cylinder, of course with opinable but easy to attain mesh

After meshing you will get three boundary surfaces and one volume as fluid, that in our case is air.
The three surface are exactly the inlet section, outlet secion and the wall.

Simulation setup and run

ElmerGUI can be used to introduce every parameter or setting that I've used, anyway for your convenience the configuration for the simulation phase is contained in this case.sif file, it's plain text so you can have a look into it.

Define the boundary conditions as per following table, an INLET with a constant velocity profile of magnitude 0,1m/s, the OULET with a constrain on the velocity component along y and z axis that is zero and finally on the cylinder surface a WALL BC. Pay attention to assing the condition to the correct surface number, coordinates as following 1 =x 2=y 3=z.

Boundary Condition 1
Target Boundaries(1) = 1
Name = "OUTLET"
Velocity 3 = 0
save scalars = logical true
Velocity 2 = 0
End

Boundary Condition 2
Target Boundaries(1) = 2
Name = "WALL"
Noslip wall BC = True
End

Boundary Condition 3
Target Boundaries(1) = 3
Name = "INLET"
Velocity 3 = 0
Velocity 1 = 0.1
Velocity 2 = 0
End


Define also the fluid properties as

Fluid properties, from case.sif
Name = "Air (room temperature)"
Reference Temperature = 20°C
Viscosity = 1.983e-5 mPa
Compressibility Model = Perfect Gas
Reference Pressure = 101325 Pa
Density = 1.205 kg/m^3

Assign this fluid to the “body” entity that represent simulation volume.

Initial state is set for the entire fluid volume, our initial state is an uniform axial flow field of 0.1 m/s magnitude at a pressure of 101325 Pa.

Set then the solver to search for a steady solution using only Navier-Stokes equation (no flux eq) You can note into sif file that 'flowsolver' solver is active, that solver integrate air velocity over the inlet boundary, just usefull to check for the value of flowrate. Note also that simulation is setup to write the final step data into a file .vtu file.

Running the solver after iterations it converges as per figure P11.2


Figure PF11.2 Simulation convergence history


After convergence, look for threshold into the setup section in case.sif file, have been reach it's possible to check for the flowrate value. In detail in the solver output windows at the last iteration should read

...

SaveScalars: -----------------------------------------

SaveScalars: Saving scalar values of various kinds

SaveScalars: Performed 1 reduction operations

SaveScalars: Found 3 values to save in total

SaveScalars: Showing computed results:

SaveScalars: 1: convective flux: flow solution over bc 1 9.420126351287E-004

...

ResultOutputSolver: Saving in unstructured VTK XML (.vtu) format

VtuOutputSolver: Saving results in VTK XML format

VtuOutputSolver: Saving to file: ./case0001.vtu


ResultOutputSolver: -------------------------------------


ElmerSolver: *** Elmer Solver: ALL DONE ***

ElmerSolver: The end


The simulation calculated flowrate is “flow over bc1”  ,  a value of 9.4e-4 kg/s

Data visualization and post processing

As you already seen some data can be directly read in the simulation output screen, this possibility is quite comfortable during long simulation runs.

Anyway the final simulation step variables will be stored in case.ep file and case000.vtu

Example results can be visualized with the elmergui integrated postprocessor VTK or, for example, with an external application that can handle .vtu files. Have a look to Paraview for a free well suited postprocess package that can handle .vtu files.

It is straightforward to setup tridimensional views and sections, see following figure for a view of velocity_x contour in a section view of pipe


PF11.3 Velocity contours of a section view, inlet at the left side, outlet on the right side


Note on the above figure that in the inlet section the velocity have the attended constant magnitude , as expected the velocity profile evolves and at the outlet seem to be fully developed.


PF11.4 Velocity plot with multiple slices with Paraview, inlet at the left side, outlet on the right side

Paraview can also handle multiple views and data from parallel processing runs. With this postprocessor is easy to retrieve simulation data and setup custom data filters, for example in many cases there in no need to configure elmer to write to file common simulation values; for example in the following figure the velocity profile is directly selected and plotted into Paraview.



PF11.5 Venturi velocity profiles at different location, Paraview

Solution validation


At glance our problem have been defined with the following parameters

Geometry of cylindric pipe

L=Pipe length 1m
D=Internal diameter of pipe 0.1m

Fluid properties, from case.sif

Name = "Air (room temperature)"
Reference Temperature = 20°C
v,Viscosity = 1.983e-5 Pas
Compressibility Model = Perfect Gas
Reference Pressure = 101325 Pa
rho, Density = 1.205 kg/m^3

Fluid axial velocity is Vag=0.1 m/s
Re=rho*D*Vag/v=608 Re<<2300,  laminar flow

At the input section mass flowrate can be writen as

Eq P12.2
qm=qv*rho=Vag*pi*D^2/4*rho=9.4 e-4 kg/s

Good news are that calculated flowrate is the same that previously simulated flowrate, starting from this result is necessary to have an evidence that simulated velocity profile is that predicted by eq P12.1. Simulation was setup on purpose to save velocity values along the inlet section radius, you find this data into g.dat file, into profile compare spreadsheet you find the figures. The validation idea is to plot simulated values against calculated values, you can see the results in the following figure.



F11.6 Velocity profile comparison, evident agreement between two methods

By inspection of figure it's safe to say that there is an evident match, it's far more evident in the following figure that plot simulated velocity deviation from the ideal case.



F11.7 Absolute deviation from ideal case



Have fun,

JLJ