Experimental Evaluation of a 3D-Printed Fluidic System for a Directional Anemometer

An evolution of a previously proposed anemometer capable of detecting both the magnitude and the direction of the wind on a plane is proposed. The device is based on a recently formalized principle, consisting of combining the differential pressures measured across distinct diameters of a cylinder to estimate the wind velocity and incidence angle. Differently from previous sensors based on the same principle, the proposed anemometers use 3D printing to fabricate the channel structure that calculates the pressure combination in the fluidic domain. Furthermore, commercial sensors with low power consumption are used to read the two pressures that result from the fluidic processing. The whole fabrication procedure requires inexpensive equipment and can be adopted by small enterprises or research laboratories. Two original channel structures, predicted by previous theoretical work but never experimentally validated, are proposed. The results of detailed experiments performed in a wind tunnel are reported.


Introduction
Anemometers are devices designed to detect the velocity vectors of gas streams. Traditional fields of application for anemometers are weather forecasting, assistance in navigation and airport operations. Recently, the importance of detecting the magnitude and, possibly, the direction of air streams has emerged in several innovative fields. Precision agriculture [1] requires the monitoring of several parameters, among which the wind velocity is recognized to be an important factor [2,3]. Even the annual production of honey was found to be affected by the wind intensity [4]. In the field of robotics, anemometers are being proposed as fundamental sensors to be combined with gas detectors in order to enable autonomous vehicles to locate gas sources [5][6][7]. UAVs (unmanned aerial vehicles) represent another emerging field that is stimulating the development of a new generation of anemometers. UAVs can be used, for example, to produce 3D maps of the wind distribution for studies of atmosphere phenomena [8,9] and to predict the evolution of forest fires [10]. Wind is also a primary source of disturbance on UAV flight parameters: the prediction of wind effects has been proposed as a way to improve trajectory and attitude control [11]. Although wind intensity and direction can be inferred by the knowledge of the UAV attitude and other flight parameters (indirect approach) [10,12], direct measurement of the wind velocity by an onboard anemometer may help to improve the robustness of the wind estimate and is hindered only by the limitations of the currently available sensors [13].
Sensors 2020, 20, 4094 3 of 21 prototype developed in [25], and following devices based on the same principle [28][29][30], very small pressure differences were acquired by monitoring the air flow rate induced by the pressure difference inside a capillary duct. This operation was accomplished using non-commercial MEMS flow sensors previously developed by our group [31]. While in [25,[28][29][30], the position of the point pairs at which the pressure was sampled and the weights that characterize the linear combination of differential pressures were determined on the basis of a numerical optimization procedure, a theoretical explanation of the approach was given only in 2018 [32], using spectral analysis. Different approaches that use several differential pressure samples measured across a cylinder but infer the wind direction and velocity from analytical approximations of the pressure around the stagnation point have been proposed [8,33]. The use of a sphere in place of the cylinder was proposed in [34].
In this work, we present a further development of the work done around the DP approach formalized in [32]. The differences with respect to the previous practical implementations of the approach [25,[28][29][30] are the following: (i) use of 3D printing technology for the fabrication of the fluidic components that operate the linear combination of pressures; (ii) use of commercial differential pressure sensors in order to facilitate the development of low-cost products; (iii) experimentation of different angle/weight combinations than in previous works. The use of 3D printing has been already proposed in the field of anemometry to fabricate Pitot tubes [35], Cobra probes [36], strain-based anemometers [37,38], hotwire flow sensors [39] and a capacitive pressure-based anemometer [40].
In this work, 3D printing (stereolithography) is used to fabricate complex fluidic structures and for the large diffusion of this technique that should favor reproducing the results presented here by other research groups. Figure 1a depicts the cross-section of a cylinder, taken by cutting the cylinder with a plane perpendicular to the cylinder axis. In the rest of this paper, we will use the term "diameter" to indicate a straight-line segment that passes through the center of the circle and whose endpoints lie on the circumference. For the sake of simplicity, depending on the context, we will use the word "diameter" to mean also the width of the circle, equal to twice the radius. In Figure 1a, a reference axis is shown; the diameter that lies along the reference axis will be indicated with d 0 (reference diameter). The cylinder is exposed to an air stream (wind) whose velocity vector u w is parallel to the plane of the cross-section and forms an angle θ with the reference axis. An even number of additional diameters are considered together with the reference diameter in such a way that diameter d k forms an angle φ k with the reference axis. The configuration is symmetric with respect to d 0 , so that for each diameter d k , we have also a symmetrical diameter dk at angle -φ k .

Principles of Operation
Sensors 2020, 20, x FOR PEER REVIEW 3 of 21 small pressure differences were acquired by monitoring the air flow rate induced by the pressure difference inside a capillary duct. This operation was accomplished using non-commercial MEMS flow sensors previously developed by our group [31]. While in [25,[28][29][30], the position of the point pairs at which the pressure was sampled and the weights that characterize the linear combination of differential pressures were determined on the basis of a numerical optimization procedure, a theoretical explanation of the approach was given only in 2018 [32], using spectral analysis. Different approaches that use several differential pressure samples measured across a cylinder but infer the wind direction and velocity from analytical approximations of the pressure around the stagnation point have been proposed [8,33]. The use of a sphere in place of the cylinder was proposed in [34].
In this work, we present a further development of the work done around the DP approach formalized in [32]. The differences with respect to the previous practical implementations of the approach [25,[28][29][30] are the following: (i) use of 3D printing technology for the fabrication of the fluidic components that operate the linear combination of pressures; (ii) use of commercial differential pressure sensors in order to facilitate the development of low-cost products; (iii) experimentation of different angle/weight combinations than in previous works. The use of 3D printing has been already proposed in the field of anemometry to fabricate Pitot tubes [35], Cobra probes [36], strain-based anemometers [37,38], hotwire flow sensors [39] and a capacitive pressurebased anemometer [40].
In this work, 3D printing (stereolithography) is used to fabricate complex fluidic structures and for the large diffusion of this technique that should favor reproducing the results presented here by other research groups. Figure 1a depicts the cross-section of a cylinder, taken by cutting the cylinder with a plane perpendicular to the cylinder axis. In the rest of this paper, we will use the term "diameter" to indicate a straight-line segment that passes through the center of the circle and whose endpoints lie on the circumference. For the sake of simplicity, depending on the context, we will use the word "diameter" to mean also the width of the circle, equal to twice the radius. In Figure 1a, a reference axis is shown; the diameter that lies along the reference axis will be indicated with d0 (reference diameter). The cylinder is exposed to an air stream (wind) whose velocity vector uw is parallel to the plane of the cross-section and forms an angle θ with the reference axis. An even number of additional diameters are considered together with the reference diameter in such a way that diameter dk forms an angle φk with the reference axis. The configuration is symmetric with respect to d0, so that for each diameter dk, we have also a symmetrical diameter d-k at angle -φk.  A given diameter configuration is identified by the number of diameters and by the angles that each diameter forms with the reference diameter (d 0 ). The number of diameters is expressed by an integer N, such that the diameter set is d −N , ... , d −1 , d 0 , d 1 , ... , d N , forming angles -φ N , .....−φ 1 , 0, φ 1 , ...., φ N . Then, the total number of diameters is 2N + 1, with a total of N independent angles. Let us define the diametric pressure as the pressure difference developed by the wind across a given diameter. A diametric pressure is represented by the symbol P D (φ,θ,u w ), where φ is the angle formed by the diameter with d 0 , while θ and u w are the wind direction (see Figure 1a) and magnitude. Following [32], it is possible to associate a set of N+1 weights w 0 ,..w N to a given diameter configuration, such as that of Figure 1a, and calculate the following linear combination of the diametric pressures:

Principles of Operation
where w −k = w k . In [32], it was shown that, with a proper choice of angles φ k and weights w k , it is possible to obtain a P C function whose dependence on the wind direction (θ) is a good approximation of a cosine function. In this way, we can write the following: where H(u w ) is a monotonic function of the wind magnitude u w only. A 2D anemometer can be simply obtained by using two diameter arrangements like the one in Figure 1a, mounted in such a way that their reference diameters (d 0 ) are orthogonal. Consider a reference system with the x-axis parallel to the central diameter of one arrangement (the "X section") and the y-axis parallel to diameter d 0 of the other arrangement (the "Y section"). Then, applying (1) to both sections, we can obtain two combined pressures P X and P Y , respectively, that show the following dependence on u w and θ: where θ is the angle formed by the wind direction with the x-axis. Clearly, θ and u w can be easily found from P X and P Y by the following: Due to the monotonic behavior of H(u w ), it is possible to calculate u w using a calibration table or a proper fitting function. In [32], a theoretical explanation of Equation (2) is given on the basis of spectral analysis. Different sets of angles and weights that provide optimal approximation of the cosine behavior are reported in [32]. Arrangements with a larger number of diameters (i.e., larger N) generally give a better approximation and a larger interval of wind velocities where the approximation holds. Two different approaches are proposed: uniform weights (UW), characterized by w o = w 1 = .... w N , and uniform diameter spacing (UDS), where the diameters are evenly spaced across the 0-90 • interval while the weights are not uniform. In the case of UW, it is necessary to find the optimal angles that, for a given number N, provide an optimum cosine approximation. In the case of UDS, the angles are fixed, and optimization is obtained by a proper set of weights. More details are given in [32].
A straightforward implementation of Equation (1) is placing one pressure sensor at the end of each diameter of the set. Diametric pressures, obtained by simple differences of two couples of pressure readings, can be combined by simple algebraic calculations. In this way, a large number of pressure sensors is required. For example, in the case of only three diameters per section (N = 1), the total number of sensors is 12. This number can be halved by using differential pressure sensors to acquire the diametric pressures directly, but the sensor count is still relatively high, especially for complex diameter arrangements (N>1). An important simplification was proposed in [25,[28][29][30], where combination Equation (1) was performed in the fluidic domain, using a structure similar to that depicted in Figure 1b. An exhaustive explanation of the mechanism is given in [32]. Briefly, two cavities, H 1 and H 2 in Figure 1b, are connected to the cylinder's outer surface by means of channels of small diameter. The points where the channels end at the outer surface are the extremes of the diameters chosen to implement the approach described above. For the sake of clarity, this correspondence is shown only for diameter d 1 in Figure 1b. Considering laminar flow inside the channels, it can be shown that the pressure inside the cavities is a weighted average of the pressure values picked up at the points where the channels end on the outer surface [32]. Finally, if we take the pressure difference between chambers H 1 and H 2 , we obtain a quantity equal to P C (θ,u w ) in Equation (1) with the following: By this approach, only two differential pressure sensors are required, greatly simplifying the device, reducing size and cost, and simplifying calibration. Different weights can be obtained by using channels of different diameter or length. We have chosen to leave the channel diameter constant and vary the length, since this option requires a coarser fabrication resolution.

Geometry of the Fluid-Dynamic Elements (Heads)
In all previous experimental works, the channel structure was fabricated using a computer-controlled milling machine. The drawback of this subtractive manufacturing approach is that complex internal cavities and channel configurations could only be obtained by assembling separate cylinder sections together. Even with this solution, limitations still apply when full 3D structures should be fabricated. Furthermore, alignment between the X and Y sections of the anemometer and leakage containment are critical operations.
Here, we have used a 3D printer to fabricate in a single piece both the X and Y sections of the cylinder. The object obtained in this way has been called the "head" of the cylinder. We have chosen to test two configurations that were proposed in [32] but have not yet been used in experimental works. The transverse and longitudinal sections of the structures are shown in Figure 2. Note that only the lower section is shown in the figure (X section), while the upper section of the head is identical but rotated by 90 • along the cylinder axis.
The first structure, shown in Figure 2a, implements a UDS configuration with N = 1 (three diameters). Different channel lengths have been obtained by the particular shape of the cavities visible in the figure. The angle between the central diameter (d 0 ) and the lateral diameters (d 1 and d −1 ) is 45 • in this structure. The proper weight ratio w 1 /w 0 is cos(45 • ), obtained by setting the L A0 /L A1 ratio equal to cos(45 • ). The second structure, shown in Figure 2b, is a UW configuration with N = 3 (seven diameters). The requirement of uniform weights translates into equal length for all channels. As a result, the cavities are more regular than in the previous structure. Another difference is that the channels in the UW structure have an inclination of 30 • with respect to the cylinder transverse plane. We made this choice to reduce the probability of water droplets reaching the internal cavities through the channels. We have designed two different UW structures, differing only in the length of the channels (2 mm and 4 mm). The structure in Figure 2b is the 4 mm structure. In the 2 mm structure, the cavities are proportionally larger, due to the shorter channels. The angles formed by the diameters (d 1 , d 2 , d 3 ) with the central diameter are reported in Table 1, where data from structures presented in previous works are also included for comparison. Channel inclination was not applied to the UDS configuration ( Figure 2a) because, due to the different length, the channels would have crossed the outer surface at different distances from the cylinder top. This, combined with the finite length of the cylinder, would have introduced possible additional causes of non-ideality.
transverse sections, are used to connect the cavities to the pressure sensors, placed on the bottom of the cylinder. The upper section (Y section), placed over the structures shown in Figure 2, is connected to its individual pressure sensor through channels (indicated with dashed lines) that pass through the thick membrane that separates the two cavities of the lower section, avoiding any interaction with the X section. As a result, two pairs of holes are present on the bottom face of the heads. These hole pairs will be indicated as X and Y ports, respectively, in the remainder of this paper.  Other dimensions indicated in Figure 2 are reported in Table 2 for the three structures examined in this paper. The H7L and H7S structures are different only in the lengths of the channels that connect the internal cavities to the outer surface. Note that the actual channel length for these structures is greater than LB due to the 30° inclination shown in Figure 2b, since LB is the length of the channel projection on the transverse cross-section. The three heads that appear in Table 2 are composed by the stacking of two orthogonal channel/cavity structures with the corresponding geometry indicated in the table.  Table 1. Angles and weights used in the devices proposed in this paper and in previous works.
The cylinder diameter (D cy in Figure 2) is set to 20 mm for all types of structures. Vertical channels, visible in the longitudinal sections of Figure 2 and represented as circular holes in the transverse sections, are used to connect the cavities to the pressure sensors, placed on the bottom of the cylinder. The upper section (Y section), placed over the structures shown in Figure 2, is connected to its individual pressure sensor through channels (indicated with dashed lines) that pass through the thick membrane that separates the two cavities of the lower section, avoiding any interaction with the X section. As a result, two pairs of holes are present on the bottom face of the heads. These hole pairs will be indicated as X and Y ports, respectively, in the remainder of this paper.
Other dimensions indicated in Figure 2 are reported in Table 2 for the three structures examined in this paper. The H7L and H7S structures are different only in the lengths of the channels that connect the internal cavities to the outer surface. Note that the actual channel length for these structures is greater than L B due to the 30 • inclination shown in Figure 2b, since L B is the length of the channel projection on the transverse cross-section. The three heads that appear in Table 2 are composed by the stacking of two orthogonal channel/cavity structures with the corresponding geometry indicated in the table. Table 2. Dimensions used for the three heads tested in this work.

Head
Type

Fabrication of the Fluid-Dynamic Elements (Heads)
The fabrication of the three heads has been performed by means of stereolithography, using the 3D printer "FORM 2" of Formlabs [41]. The material used was the "Resin Clear version 04" (methacrylate-based photopolymer) provided by Formlabs. The vertical resolution (single-layer thickness) of the printing process was 50 µm, while the minimum feature size on the X-Y plane was 150 µm. The total fabrication time was around 4 h for all heads. After laser polymerization, the heads were rinsed in isopropanol for 20 min to remove excess resin. This step, performed by the form-wash tool provided by Formlabs, was particularly important to prevent channel clogging from excess resin. The size of the resin bath allowed the fabrication of up to four heads in parallel.
Design of the heads was carried out using the open source program FreeCAD [42], producing output files in STL (STereoLithography) format, ready to be imported into the proprietary printer management software. Figure 3 shows photographs of the three different heads used in the experiments described in this paper. Figure  LA0 (mm)

. Fabrication of the Fluid-Dynamic Elements (Heads)
The fabrication of the three heads has been performed by means of stereolithography, using the 3D printer "FORM 2" of Formlabs [41]. The material used was the "Resin Clear version 04" (methacrylate-based photopolymer) provided by Formlabs. The vertical resolution (single-layer thickness) of the printing process was 50 μm, while the minimum feature size on the X-Y plane was 150 μm. The total fabrication time was around 4 hours for all heads. After laser polymerization, the heads were rinsed in isopropanol for 20 min to remove excess resin. This step, performed by the formwash tool provided by Formlabs, was particularly important to prevent channel clogging from excess resin. The size of the resin bath allowed the fabrication of up to four heads in parallel.
Design of the heads was carried out using the open source program FreeCAD [42], producing output files in STL (STereoLithography) format, ready to be imported into the proprietary printer management software. Figure 3 shows photographs of the three different heads used in the experiments described in this paper.  Figure 4 illustrates the structure of the devices used in the experiments described in this paper. On the left, the internal structure of the heads is shown through a perspective view directly extracted from the FreeCAD design. The particular head shown in the figure is the H7L head. The X and Y sections and the pressure ports that allow access to the internal cavities of the individual sections are clearly visible. The heads are combined with other, simpler pieces to form the cylinder structure, as shown on the right. In particular, a top cylindrical piece, indicated with "upper extension" in the figure, has been placed on top of the heads to prolong the cylinder structure, reducing the effect of discontinuities at the top border of the heads. To further reduce discontinuity effects, the upper extension ended with a chamfer. A hollow cylinder, indicated with "lower extension" in Figure 4, was connected to the bottom face of the heads to continue the cylinder structure and cover the connection to the four pressure ports. Stacking of the lower extension, head, and upper extension forms a continuous cylinder structure. The upper and lower extensions were secured to the head by using adhesive tape in order to facilitate replacement of the heads. Connection to the pressure ports was accomplished by inserting stainless steel needles into the holes in the bottom faces of the heads,  Figure 4 illustrates the structure of the devices used in the experiments described in this paper. On the left, the internal structure of the heads is shown through a perspective view directly extracted from the FreeCAD design. The particular head shown in the figure is the H7L head. The X and Y sections and the pressure ports that allow access to the internal cavities of the individual sections are clearly visible. The heads are combined with other, simpler pieces to form the cylinder structure, as shown on the right. In particular, a top cylindrical piece, indicated with "upper extension" in the figure, has been placed on top of the heads to prolong the cylinder structure, reducing the effect of discontinuities at the top border of the heads. To further reduce discontinuity effects, the upper extension ended with a chamfer. A hollow cylinder, indicated with "lower extension" in Figure 4, was connected to the bottom face of the heads to continue the cylinder structure and cover the connection to the four pressure ports. Stacking of the lower extension, head, and upper extension forms a continuous cylinder structure. The upper and lower extensions were secured to the head by using adhesive tape in order to facilitate replacement of the heads. Connection to the pressure ports was accomplished by inserting stainless steel needles into the holes in the bottom faces of the heads, as visible in Figure 3b. The needles were secured to the heads by means of cyanoacrylate glue, carefully deposited around the needle's outer surface before insertion into the pressure ports and around the junction perimeter.

Architecture of the Prototypes
around the junction perimeter.
Connection to the pressure sensors was accomplished by means of silicone tubes. Two different diameters were required to fit the head pressure ports to the input ports of the pressure sensors, as schematically shown in the figure. The lower extension was split into two parts in order to simplify the manual connection of the head pressure ports to the silicone pipes. The cylinder was joined to a circular plate providing support to the pressure sensors. Differential pressure sensor model SDP 810 by Sensirion [43] was used to detect PX end PY (see Equation (3)). The full-scale pressure of sensors SDP 810 is ± 500 Pa, while the resolution is around 0.015 Pa (calculated considering the 16-bit output data resolution). Each SDP 810 sensor requires around 12.5 mW from a 3.3 V power supply. Note that the estimated PX and PY peak value for a wind velocity of 1 m/s is in the order of 1 Pa [32], thus pressure sensors with resolutions well below 1 Pa are required in order to enable operation at low wind velocities. Connection to the pressure sensors was accomplished by means of silicone tubes. Two different diameters were required to fit the head pressure ports to the input ports of the pressure sensors, as schematically shown in the figure. The lower extension was split into two parts in order to simplify the manual connection of the head pressure ports to the silicone pipes.
The cylinder was joined to a circular plate providing support to the pressure sensors. Differential pressure sensor model SDP 810 by Sensirion [43] was used to detect P X end P Y (see Equation (3)). The full-scale pressure of sensors SDP 810 is ±500 Pa, while the resolution is around 0.015 Pa (calculated considering the 16-bit output data resolution). Each SDP 810 sensor requires around 12.5 mW from a 3.3 V power supply. Note that the estimated P X and P Y peak value for a wind velocity of 1 m/s is in the order of 1 Pa [32], thus pressure sensors with resolutions well below 1 Pa are required in order to enable operation at low wind velocities.
The pressure sensors are read with a general purpose interface based on the MSP430i2041 microcontroller (Texas Instruments, Dallas, TX, USA). In particular, an I2C hub was required to independently communicate with the two sensors, which could not be connected to the same I2C bus since they are shipped with the same I2C address. The I2C hub was implemented as a distinct printed circuit board (PCB) module, based on the PCA9518 integrated circuit of Texas Instruments. Connection from the microcontroller unit (MCU) to a personal computer was implemented using the RN42 Bluetooth module by Microchip. The three components of the electronic interface shown in Figure 4 were designed by the authors and built as three separate printed circuit boards (PCBs).
A photograph of the whole device is shown in Figure 5, where the various components represented in Figure 4 are clearly indicated. The pressure sensors are read with a general purpose interface based on the MSP430i2041 microcontroller (Texas Instruments, Dallas, Texas, USA). In particular, an I2C hub was required to independently communicate with the two sensors, which could not be connected to the same I2C bus since they are shipped with the same I2C address. The I2C hub was implemented as a distinct printed circuit board (PCB) module, based on the PCA9518 integrated circuit of Texas Instruments. Connection from the microcontroller unit (MCU) to a personal computer was implemented using the RN42 Bluetooth module by Microchip. The three components of the electronic interface shown in Figure 4 were designed by the authors and built as three separate printed circuit boards (PCBs).
A photograph of the whole device is shown in Figure 5, where the various components represented in Figure 4 are clearly indicated.

Experimemental Setup and Data Processing
The measurements were performed using a small wind tunnel, consisting of a 1.2 m long circular pipe with 10 cm diameter. An electronically controlled fan allowed the air speed inside the tunnel to be set from 0.1 to around 8 m/s. Precise measurement of the air velocity was accomplished by means of a reference hot-wire anemometer (Model Testo AG 425, Testo SpA, Milan, Italy). Rotary components of the air flow due to the effect of the fan blades were reduced by filling a whole section of the pipe following the fan with an array of pipes of smaller diameter (2 cm), for a total length of 10 cm. The anemometers under test were inserted into the wind tunnel through a hole at nearly 20 cm from the exhaust end of the tunnel. The anemometers were secured to a sample holder that allowed complete 360-degree rotation around the cylinder axis, with a resolution of 1 degree. If not otherwise specified, the cylinder axis was perpendicular to the wind tunnel axis (i.e., to the wind velocity). Inclination tests were performed by introducing a tilt angle with respect to the perpendicular condition. A maximum inclination angle of 30 degrees was allowed by the sample-holder. In all the experiments, the sensitive sections of the cylinder were placed in the center of the wind tunnel, where the velocity is the highest. The reference hot-wire anemometer was placed as closely as possible to the height of the cylinder active sections. Displacements of ± 1 cm around the centerline did not cause significant differences in the responses of both the reference anemometer and devices under test. In the inclination tests, the height of the anemometer sample-holder was varied to keep the active section on the wind-tunnel centerline.

Experimemental Setup and Data Processing
The measurements were performed using a small wind tunnel, consisting of a 1.2 m long circular pipe with 10 cm diameter. An electronically controlled fan allowed the air speed inside the tunnel to be set from 0.1 to around 8 m/s. Precise measurement of the air velocity was accomplished by means of a reference hot-wire anemometer (Model Testo AG 425, Testo SpA, Milan, Italy). Rotary components of the air flow due to the effect of the fan blades were reduced by filling a whole section of the pipe following the fan with an array of pipes of smaller diameter (2 cm), for a total length of 10 cm. The anemometers under test were inserted into the wind tunnel through a hole at nearly 20 cm from the exhaust end of the tunnel. The anemometers were secured to a sample holder that allowed complete 360-degree rotation around the cylinder axis, with a resolution of 1 degree. If not otherwise specified, the cylinder axis was perpendicular to the wind tunnel axis (i.e., to the wind velocity). Inclination tests were performed by introducing a tilt angle with respect to the perpendicular condition. A maximum inclination angle of 30 degrees was allowed by the sample-holder. In all the experiments, the sensitive sections of the cylinder were placed in the center of the wind tunnel, where the velocity is the highest. The reference hot-wire anemometer was placed as closely as possible to the height of the cylinder active sections. Displacements of ± 1 cm around the centerline did not cause significant differences in the responses of both the reference anemometer and devices under test. In the inclination tests, the height of the anemometer sample-holder was varied to keep the active section on the wind-tunnel centerline.
Pressure data acquisition was controlled by a program running on a personal computer that communicates with the anemometer through the Bluetooth link. For each measurement, a total of 16 acquisitions is averaged across a time span of nearly 2 s. By this procedure, the total noise (peak-to-peak) superimposed on the differential pressure measurements was found to be 0.02 Pa in a condition of zero wind velocity. Raw data were processed by means of scripts written with the Python language. The data were preliminarily corrected in order to take into account the pressure drop occurring along the connecting pipe. It is worth mentioning that, in conformity with other products having such a small detection limit, the SDP 810 sensors are actually thermal flow meters that infer the differential data from the mass flow rate through a capillary pipe. This is analogous to the approach that we proposed in [25,[28][29][30], where we used highly miniaturized research-grade MEMS flow sensors [31]. The fact that a flow exists in the pipes that connect the head to the pressure sensors results in a significant difference between the actual differential pressure at the head ports and the sensor reading. A procedure to take into account this error is described in an application note available on the Sensirion website. Due to the irregular characteristics of the connecting pipes and the presence of bends, we preferred to directly calibrate the combination of SDP 810 sensors and connecting pipes. To this aim, we used an air flow line equipped with a precision mass flow controller to establish a constant differential pressure across two sections of a relatively large pipe (0.8 cm diameter), separated by a narrow obstruction. The differential pressure sensor to be calibrated (either the P X or P Y sensor) was connected across these two sections, through the actual pipe configuration that connects it to the pressure ports of the anemometer heads, including the stainless steel needles. The pressure across the same points is measured by means of a reference pressure sensor (another SPD 810 device) which, differently from the sensor to be calibrated, is connected by means of large and short pipes so that it can be considered free from pressure loss. Measurements performed at different pressures were collected into look-up tables, one for each channel, which were interpolated using spline functions. The offset of both sensors (of the order of 0.02 Pa) was acquired and subtracted from the measurements. The values of P X and P Y , acquired and corrected by this procedure, will be indicated with P XC and P YC , respectively, in the rest of this paper.

Results of Measurements
The three head configurations listed in Table 2 have been experimentally characterized and compared over a wind velocity range spanning from 0.9 to 8 m/s. The heads were combined with other components to form the device structure shown in Figures 4 and 5. All additional parts were re-used, so that switching from one configuration to another was accomplished by simply changing the head. This operation was simplified by the use of adhesive tape for the connection of the parts that form the cylinder. For the sake of simplicity, detailed data resulting from the experiments are shown only for heads H3L and H7L. Data obtained with head H7S are shown in the final plot regarding measurement errors, which better summarizes the wind sensor performance and allows comparison between the different configurations. Figure 6a,b show the dependence of P XC and P YC on the wind direction θ (see Figure 1 for the meaning of angle θ) for the H3L head. The experimental data, represented by square and circular symbols, refer to two different wind velocities, namely 1.7 ( Figure 6a where ρ air and µ air are the air density and dynamic viscosity, respectively, while D cy is the cylinder diameter, equal to the head diameter given in Table 2. The Reynolds number is presented to facilitate comparison with the theoretical results described in [32]. Solid lines represent the best sinusoidal fits. Note that P XC and P YC follow a sinusoidal behavior with sufficient accuracy at the smaller velocity, while significant deviations are visible at the higher velocity setting. Deviations are larger around the peaks of the respective sinusoids. Plots in Figure 6a,b is well representative of the behavior of the H3L head: approximation of the sinusoidal dependence is acceptable at low wind velocities, with errors that progressively grow as the velocity is increased.
The same experiments performed with the H7L head produced the results shown in Figure 7a,b. The main difference that can be observed is the much better approximation of the sinusoidal behavior. Even at the higher velocity, deviation from the sinusoidal fits are modest. Note that the peak magnitude of the pressure is approximately the same for the two heads under test. The H7S head behaves much like the H7L configuration, with slightly larger deviations from the sinusoid. Sensors 2020, 20, x FOR PEER REVIEW 11 of 21 Plots in Figure 6a,b is well representative of the behavior of the H3L head: approximation of the sinusoidal dependence is acceptable at low wind velocities, with errors that progressively grow as the velocity is increased.
The same experiments performed with the H7L head produced the results shown in Figure 7a,b. The main difference that can be observed is the much better approximation of the sinusoidal behavior. Even at the higher velocity, deviation from the sinusoidal fits are modest. Note that the peak magnitude of the pressure is approximately the same for the two heads under test. The H7S head behaves much like the H7L configuration, with slightly larger deviations from the sinusoid.    Figure 6a,b is well representative of the behavior of the H3L head: approximation of the sinusoidal dependence is acceptable at low wind velocities, with errors that progressively grow as the velocity is increased.

Plots in
The same experiments performed with the H7L head produced the results shown in Figure 7a,b. The main difference that can be observed is the much better approximation of the sinusoidal behavior. Even at the higher velocity, deviation from the sinusoidal fits are modest. Note that the peak magnitude of the pressure is approximately the same for the two heads under test. The H7S head behaves much like the H7L configuration, with slightly larger deviations from the sinusoid.  In order to understand how the observed discrepancies affect the sensor accuracy, we have used Equation (4) with P X = P XC and P Y = P YC to estimate the angle θ (measured angle) and the quantity H(u w ), which we expect to be independent of θ. The capability of the H3L head to detect the wind direction is represented in Figure 8a,b, where the results of experiments performed over the whole explored wind range are shown. Figure 8a presents the measured angle as a function of the actual wind direction, read on the sample holder goniometer. A substantial agreement between the actual and measured quantity can be observed. The measurement errors are better highlighted by Figure 8b, which shows the difference between the actual and measured angle across the whole 360-degree range of wind directions. An oscillatory trend for the error can be observed, except for the 0.9 m/s data, which, due to the much smaller signal level, are likely to suffer from additional error from noise and other inaccuracies of the pressure sensors. The angular error varies from roughly −11 to +7 degrees. This uncertainty band reduces significantly if the 0.9 m/s data are excluded.
where PX and PY are the pressures that would be given by an ideal sinusoidal behavior (solid lines in Figures 6 and 7), Δθ is the angular error, and ΔPY, ΔPX are the pressure errors with respect to the ideal (sinusoidal) dependence. The consequence of Equation (7) is that the error on one of the two pressures is weighted by the value of the other pressure. As a result, errors located around the peaks give a negligible contribution to the angular error, since the other pressure is nearly zero.  Angle estimation data for the H7L head are shown in Figure 9a,b. Comparison with the H3L head ( Figure 8) does not show an important improvement in terms of angular accuracy, in spite of the much better fit of the sinusoidal behavior illustrated by Figures 6 and 7. This can be explained by observing that the larger discrepancies visible in Figure 6 for the H3L heads are mostly located where the pressure value is close to the peak. With simple arguments of error propagation theory, we can find that the error on the estimated angle is given by the following: where P X and P Y are the pressures that would be given by an ideal sinusoidal behavior (solid lines in Figures 6 and 7), ∆θ is the angular error, and ∆P Y , ∆P X are the pressure errors with respect to the ideal (sinusoidal) dependence. The consequence of Equation (7) is that the error on one of the two pressures is weighted by the value of the other pressure. As a result, errors located around the peaks give a negligible contribution to the angular error, since the other pressure is nearly zero. The situation is different when the estimate of quantity H(u w ) is considered. In the case in which the sinusoidal laws represented by Equation (3) are strictly respected, the value of H(u w ), estimated from P XC and P YC through Equation (4), would be independent of the wind direction θ. This would allow the use of H(u w ) as an indication of the wind velocity magnitude, which could be estimated from H(u w ) with simple expressions, as suggested in [32], or using an empirical calibration curve. The results of the estimation of H(u w ) are shown in Figure 10a,b for the H3L and H7L head, respectively. Sensors 2020, 20, x FOR PEER REVIEW 13 of 21 The situation is different when the estimate of quantity H(uw) is considered. In the case in which the sinusoidal laws represented by Equation (3) are strictly respected, the value of H(uw), estimated from PXC and PYC through Equation (4), would be independent of the wind direction θ. This would allow the use of H(uw) as an indication of the wind velocity magnitude, which could be estimated from H(uw) with simple expressions, as suggested in [32], or using an empirical calibration curve. The results of the estimation of H(uw) are shown in Figure 10a,b for the H3L and H7L head, respectively.   The situation is different when the estimate of quantity H(uw) is considered. In the case in which the sinusoidal laws represented by Equation (3) are strictly respected, the value of H(uw), estimated from PXC and PYC through Equation (4), would be independent of the wind direction θ. This would allow the use of H(uw) as an indication of the wind velocity magnitude, which could be estimated from H(uw) with simple expressions, as suggested in [32], or using an empirical calibration curve. The results of the estimation of H(uw) are shown in Figure 10a,b for the H3L and H7L head, respectively.  The much larger dependence on the wind direction exhibited by the H3L head is clearly visible. Obviously, this is a consequence of the larger discrepancies exhibited by P XC and P YC for the H3L head. Differently from the case of the angle estimation, the expression of H(u w ) is more sensitive to errors located in the proximity of peaks.
It is interesting to observe how the quantity H(u w ) depends on the wind velocity magnitude. To eliminate the effect of the oscillations visible in Figure 10, we have averaged H(u w ) across the full 360-degree angle range. The result is shown in Figure 11 for both the H3L and H7L head. A straight line fit for the H7L data is indicated. that can be used to calculate the wind velocity from H(uw) is as follows: where A = 1.38 and β = 0.511 with uw in m/s and H in Pa. The accuracy that it is possible to obtain from the various heads tested in this paper is summarized in Figure 12, which shows the angular and magnitude error as a function of the magnitude of the wind velocity for heads H3L, H7L, and H7S. Each error point is calculated by taking the worst (largest) error measured by sweeping the wind direction across the full 360-degree range. As a measure of the magnitude error, we considered the maximum relative deviation of H(uw) from the average value, corresponding to the peak value of the oscillations visible in Figure 10. The excellent approximation of the fit indicates that H(u w ) follows a power law dependence, as predicted in [32]. The fit has not been drawn for the H3L data since they do not appear to be aligned along a straight line. Considering the power law estimated from the fit of Figure 11, the expression that can be used to calculate the wind velocity from H(u w ) is as follows: where A = 1.38 and β = 0.511 with u w in m/s and H in Pa. The accuracy that it is possible to obtain from the various heads tested in this paper is summarized in Figure 12, which shows the angular and magnitude error as a function of the magnitude of the wind velocity for heads H3L, H7L, and H7S. Each error point is calculated by taking the worst (largest) error measured by sweeping the wind direction across the full 360-degree range. As a measure of the magnitude error, we considered the maximum relative deviation of H(u w ) from the average value, corresponding to the peak value of the oscillations visible in Figure 10.  For velocities smaller than those represented in Figures 11 and 12, the pressure differences to be measured become so small that the errors from the sensors cannot be neglected. This is clearly proven by Figure 13a, where pressures P XC and P YC are shown as the function of the wind direction for the H7L head and a wind velocity of 0.3 m/s.  For velocities smaller than those represented in Figures 11 and 12, the pressure differences to be measured become so small that the errors from the sensors cannot be neglected. This is clearly proven by Figure 13a, where pressures PXC and PYC are shown as the function of the wind direction for the H7L head and a wind velocity of 0.3 m/s. The peak value of the pressure (0.05 Pa) at this velocity is less than three times the measured peak-to-peak noise (0.02 Pa). Figure 13b shows the angular and magnitude error plots for the H7L heads, where two points at 0.5 and 0.3 m/s have been added with respect to the same plots in Figure 12. The steep error increase for velocities smaller than 0.9 m/s is due to the dominance of the The peak value of the pressure (0.05 Pa) at this velocity is less than three times the measured peak-to-peak noise (0.02 Pa). Figure 13b shows the angular and magnitude error plots for the H7L heads, where two points at 0.5 and 0.3 m/s have been added with respect to the same plots in Figure 12. The steep error increase for velocities smaller than 0.9 m/s is due to the dominance of the sensor noise. From Equation (7), it is possible to find the peak-to-peak value θ n-pp of the noise superimposed on the angle estimation, considering that the equivalent noise pressures of sensors X and Y can be considered two independent random processes with same peak-to-peak value P n-pp : where θ n-pp is in radians and H is given by Equation (4). At small wind velocities, H decreases according to a square law, while P n-pp is clearly unaffected. Then, the fluctuation of the angle estimation grows larger and larger as u W is decreased and, eventually, becomes the dominant contribution to the error. The same consideration can be easily extended to the contribution of P n-pp to the estimation of H. Finally, we have investigated the effect of the inclination of the cylinder axis with respect to the wind velocity vector. The results of experiments performed with the H7L head and a velocity of 2.4 m/s are shown in Figure 14. In particular, Figure 14a shows pressures P XC and P YC as a function of the wind direction in the case of a cylinder axis perpendicular to the plane of the wind velocity (α = 0) and for an inclination of 30 degrees (the tip of the cylinder was moved downstream). The attenuation of the pressure curves is close to cos(30 • ), which suggests that the situation is equivalent to considering only the velocity component on the plane perpendicular to the cylinder axis. Estimation of the wind direction performed with the inclined cylinder, shown in Figure 14b, is marked by only a slightly larger error than in the case of zero inclination.
(α = 0) and for an inclination of 30 degrees (the tip of the cylinder was moved downstream). The attenuation of the pressure curves is close to cos(30°), which suggests that the situation is equivalent to considering only the velocity component on the plane perpendicular to the cylinder axis. Estimation of the wind direction performed with the inclined cylinder, shown in Figure 14b, is marked by only a slightly larger error than in the case of zero inclination.

Discussion
The three heads tested in this work implement two different configurations formalized in [32]. In particular, the H3L head follows the UDS approach, with N = 1 (three diameters) per section, while the H7L and H7S represent different cases of the UW configuration with N = 3 (seven diameters). The difference between H7S and H7L is the length of the channels that connect the internal cavities to the cylinder outer surface. The errors produced by the proposed anemometers and summarized in Figure 12 are quite larger than the prediction that can be found in [32] for the same channel configurations and the same Reynolds number range. This result is in line with previous experimental tests [25,[28][29][30], based on the channel configurations summarized in Table 1. The main

Discussion
The three heads tested in this work implement two different configurations formalized in [32]. In particular, the H3L head follows the UDS approach, with N = 1 (three diameters) per section, while the H7L and H7S represent different cases of the UW configuration with N = 3 (seven diameters). The difference between H7S and H7L is the length of the channels that connect the internal cavities to the cylinder outer surface. The errors produced by the proposed anemometers and summarized in Figure 12 are quite larger than the prediction that can be found in [32] for the same channel configurations and the same Reynolds number range. This result is in line with previous experimental tests [25,[28][29][30], based on the channel configurations summarized in Table 1. The main reason, as already pointed out in [32], is the short length of the channels, which is in contrast with our hypothesis of a fully developed laminar flow. The possibility of calculating the linear combination of pressures represented by Equation (1) requires a linear dependence between the pressure drop across the channels and the mass flow rate through them. This is guaranteed only if the flow is laminar along the full length of the channels. Even if the Reynolds number of the flow inside the channels is small enough for a laminar flow to exist, a region of perturbed flow is present at the point where the channels end at the cylinder's outer surface. The depth of this region can be several times the channel radius and grow longer at higher velocities. Comparing the H7L and H7S heads, which differ only in the channel length, it is clear that the longer channels of the H7L configuration produce smaller errors across the whole velocity range.
As far as the H3L head is concerned, deviations from the ideal behavior predicted in [32] for the same configuration are even larger and show up clearly in the poor approximation of the sinusoidal behavior visible in Figure 6. It should be observed that, in the UDS configurations, the approximation of the sinusoidal behavior relies on precise ratios between the weights, which, in the proposed implementation, would require a precise dependence between the fluidic conductance and the channel lengths. In turn, this needs a fully developed laminar flow inside the whole length of the channels. In addition, the H3L configuration is based on a smaller number of diameters with respect to heads H7L and H7S, and this contributes to increasing the error. Figure 12 confirms that the large deviations from the sinusoidal behavior have little impact on the angular accuracy, which is similar for the H3L and H7L heads, while it clearly impacts the independence of the H(u w ) quantity from the angle (Figure 12a). In this respect, the H3L head performances are very similar to those of the H7S head, which has much shorter channels. The practical consequence is that the H3L and H7S heads will provide an estimate of the wind magnitude that depends more on the angle than the H7L head. The latter presents also the most regular dependence of the H magnitude on the wind magnitude, allowing simple determination of the latter.
It is useful to compare the anemometers presented in this paper with earlier devices based on the same principle, fabricated with a fully custom procedure. The maximum angular error and relative magnitude deviation are reported in Table 3 for the H7L and H3L heads and for three different channel configurations proposed in previous works (lines 3-5 of the table). A direct comparison between fabrication approaches is not possible, due to the different channel configurations. However, three considerations can be drawn: (i) UDS approaches, which require the implementation of different weights in the pressure combination and provide worse accuracies, independently of the order (N) and of the fabrication approach; (ii) the performances obtained with the H7L head are only slightly worse than those obtained with similar UW configurations and custom fabrication technology; (iii) in all cases, the errors are significantly larger than the prediction shown in [32]. The results obtained with the H7L heads indicate that combination of 3D printing and off-the-shelf pressure sensors is a viable approach to fabricating low-cost directional anemometers for applications that do not require high accuracies. On the other hand, the fact that, both in this work and in previous studies, a large discrepancy with respect to the theory presented in [32] was found indicates that further investigation is required to improve this class of devices. For example, longer channels could be designed, exploiting the potentiality of 3D printing, which here allowed the fabrication of inclined channels (H7L and H7S heads), which is not possible with the milling machine used in our previous works [25,[28][29][30].
Preliminary investigation performed by means of 2D FEA (finite element analysis) has been performed on the H3L head in order to check how the hypotheses of (i) uniform pressure gradient along the channels and (ii) uniform pressure distribution inside H 1 and H 2 cavities are satisfied. Simulations have been executed using the COMSOL Multiphysics platform with the SST (shear stress transport) turbulence model. Figure 15a shows the pressure distribution for a velocity of 0.9 m/s and a wind direction θ = 0. The pressure inside the central channel (indicated by an arrow in Figure 15a) is shown in Figure 15b as a function of the distance from the channel endpoint at the cylinder's outer surface. The same plots for a velocity of 5.3 m/s are shown in Figure 15c,d, respectively. These simulations suggest that the pressure gradient inside the channel is not constant and that the pressure in the cavities is not uniform. These non-idealities are much more important at higher wind velocities. In particular, Figure 15d prove that most of the pressure drop occurs in the first quarter of the channel. The results shown in Figure 15a,d should be regarded as a qualitative demonstration that the pressure irregularity inside the channel structure can be the cause of the observed discrepancy with respect of the performances predicted in [32] and, in particular, of the unsatisfactory approximation of the sinusoidal dependence visible in Figure 6b. Quantitative predictions would require a full 3D model.
The results shown in Figure 15a,d should be regarded as a qualitative demonstration that the pressure irregularity inside the channel structure can be the cause of the observed discrepancy with respect of the performances predicted in [32] and, in particular, of the unsatisfactory approximation of the sinusoidal dependence visible in Figure 6b. Quantitative predictions would require a full 3D model. The resolution of the differential pressure sensors determines the lowest applicable wind velocity: Equation (9) indicates that the peak pressure difference (equal to H) should be much larger than the equivalent noise pressure of the sensors to maintain an acceptable angular resolution. A velocity of 0.3 m/s can be considered as a lower bound for the operating range of devices described The resolution of the differential pressure sensors determines the lowest applicable wind velocity: Equation (9) indicates that the peak pressure difference (equal to H) should be much larger than the equivalent noise pressure of the sensors to maintain an acceptable angular resolution. A velocity of 0.3 m/s can be considered as a lower bound for the operating range of devices described in this paper. Operation at even smaller velocities can be achieved only using pressure sensors with a lower equivalent noise.
The effect of the cylinder inclination appears only as a uniform attenuation of the pressure differences. A proportionality to the cosine of the inclination angle was found, as in [8]. The capability of detecting the wind direction is not impaired, at least up to a cylinder inclination of 30 degrees, which was the maximum allowed by the set-up.
Finally, it is important to point out that the approach presented in this work is highly scalable, opening the way to the fabrication of directional wind sensors with even smaller dimensions. Uniform scaling of the heads would result in a proportional variation in the Reynolds number, tied to the cylinder diameter. This should not be a problem, since the detection principle was demonstrated to be applicable over a wide range of Reynolds numbers. The magnitude of the differential pressures does not change significantly with the cylinder diameter, being always of the order of the dynamic pressure, equal to ρ(u w ) 2 /2. The only critical issue lies in the reduction in the channel diameter, which would increase the pressure loss when, as in the case of this paper, the pressure sensors are actually flow sensors. Sensors characterized by smaller leakage for the same differential pressure can be selected to mitigate this possible cause of signal loss. The selection of smaller pressure sensors can be also adopted to further reduce the overall dimensions of the device.