Optimizing Velocity Field Measurement with 3D-Printed Particles and MATLAB: A Cost-Effective System for Flow Visualization

Abstract

This article aims to highlight the importance of including quantitative measurements when conducting flow visualization tests, such as those performed in towing tanks, within fluid mechanics analysis. It investigates the possibility of measuring velocity fields with an economically accessible technique compared to other techniques that require large financial investments, such as traditional PIV. The development of a MATLAB R2024b code based on image recognition and the use of 3D-printed tracer particles is proposed. Code workflow and how to make a correct selection of the processing parameters and its activity are explained and demonstrated on artificial images, generated by a computer, as well as real images, obtained in a 2D-test in the tank, achieving an accuracy, in absolute values, of 95%. However, the proposed velocimetry system currently has one important limitation, the impossibility of distinguishing between particles in different planes, which limits the study to two-dimensional tests. Then, the opportunity to include this technique in the study of more complex tests requires further investigation.

1. Introduction

Experimental techniques are used in various disciplines for the validation of computational models, the study of complex phenomena, or problems that are difficult to predict. In the field of fluid mechanics, water tunnels or towing tanks provide valuable elements in fundamental studies and qualitative information by flow visualization [1,2,3,4]. However, in addition to qualitative data, it is important to be able to extract quantitative values and precise measurements of relevant magnitudes such as velocity or pressure to complement visual results, which is the aim of this paper.

Quantitative velocity measurements are a fundamental part of experimental fluid mechanics. Point measurements of velocity are possible thanks to techniques such as hot-wire anemometry or Doppler velocimetry, but undoubtedly the invention and development of techniques to measure complete scalar and vector fields instantaneously was a major achievement in this field. All the main techniques for flow visualization in water are examined in [5,6,7,8], with the latest giving particular attention to their potential for extracting quantitative information.

Most techniques are based on measuring the displacement of small markers in the fluid by comparing their positions at different instants of time. In these techniques, a region of the fluid is illuminated with successive laser sheets with known temporal separation. The light reflected by the particles is collected in a camera (or several), and the velocity is calculated by determining the particle-field displacements between two frames and the known time separation between two laser pulses. This causes these techniques to be grouped under the name Pulsed Laser Velocimetry (PLV) [9].

These can be grouped according to the type of marker used and its concentration in the fluid. For conventional Particle Image Velocimetry (PIV) the concentration is such that individual particles can be identified in an image, but not with sufficient certainty to track them independently between images [10]. Another technique called Particle Tracking Velocimetry (PTV) appears when this concentration is reduced, so it is possible to track an individual particle [11]. Finally, Laser Speckle Velocimetry (LSV) would represent the opposite extreme, where the density of the particles is so high that they cannot be distinguished individually from each other [12].

A step further, recently developed methods assume the conservation of brightness or intensity from frame to frame and use the temporal and spatial variation of intensity to infer motion, named Optical Tracking Velocimetry (OTV) [13,14]. However, these advanced techniques, based on laser technology, require significant financial investment and a high level of qualification and training of personnel in the use of the equipment.

The main novelty of this paper is to present an alternative method based on PTV theory to be used when faced with a lack of advanced equipment (laser beam, tracer particles, limited computational resources, etc.), providing a cost-efficiency, easy to use and accurate method to apply in generalized experiments for educational and investigation purposes.

To obtain that, a simplified particle tracking method using 3D-printed ABS particles and a user-friendly MATLAB interface is proposed. It enables accurate velocity field measurements in the university’s towing tank for the first time, demonstrating reliable performance with an accuracy of 86% and minimal computational and investment requirements.

In the next section, Section 2, the numerical values for the key parameters relevant to a particle tracking system are calculated for the case study. The workflow of the developed code is described in detail, along with the steps required for its validation and application in a towing tank experiment. The Section 3 provides a performance analysis of the proposed software, including the results of validation with a simulated experiment and the results obtained from the towing tank tests, while also discussing the main limitations of the system. Finally, the key findings and insights are summarized in the Section 4.

2. Materials and Methods

This chapter presents the proposed system for the determination of the velocity field in the towing tank, detailing its structure and operation, and the basis of the proposed experiment.

2.1. Experiment Characteristics Values

It is important to first understand the characteristics that such a system should have so that the results properly represent the flow motion.

First, we must consider how readily a particle can follow continuous fluid fluctuations. This characteristic can be measured using the Stokes number (Stk). It is defined as the relationship between the response time and the flow system time (Equation (1)) [15]. This number can also be understood as the ratio between the inertial and viscous forces that act on the particle.

(1)$Stk={\displaystyle \frac{t_0{u}_0}{l_0}}$

where $t_0$ is the relaxation time of the particle, $u_0$ is the velocity of the fluid away from an obstacle, and $l_0$ is the characteristic length of the obstacle or a characteristic dimension of the fluid such as the thickness of the boundary layer.

Particles with a low Stokes number follow the streamlines (perfect advection), whereas particles with a high Stokes number are dominated by their inertia and continue their original trajectory, potentially hitting an obstacle within the fluid stream [16]. If Stk < 0.1, the errors of follow-up are less than 1% [17].

When the Reynolds number associated with the particle is less than one, $t_0$ is calculated using Equation (2).

(2)$t_0={\displaystyle \frac{{\rho }_p{d}_p^2}{18\mu }}$

where ${\rho }_0$ is the density of the particles, $d_p$ is the diameter of the particles, and $\mu$ is the dynamic viscosity of the fluid.

Table 1 shows the calculation of the Stokes number for the experiment performed in this paper.

Second, the average number of particles per unit volume, C, determines the maximum number of potential points at which velocity can be measured, and therefore, it is related to the spatial resolution of the system [9]. If the particles are considered as randomly distributed points, the probability of finding k particles in volume V follows a Poisson distribution (Equation (3)).

(3)$Prob\left(k\mathrm{particles}\mathrm{in}V\right)={\displaystyle \frac{\left(C{V}^k\right)}{k!}}{e}^{-CV}$

where C is the average number of particles per unit volume and k is the number of particles to be found in the volume V.

Table 2 and Table 3 show the characteristic magnitudes and probability of the experiment carried out in this paper.

Finally, another characteristic aspect of these systems, also related to spatial resolution, is the sampling or imaging frequency. The rate at which sampling should be performed can be characterized by the mean distance between the nearest particles in the fluid volume and the Nyquist criterion. In communication theory, the Nyquist frequency is half the sampling frequency and indicates the maximum frequency that an analog signal can have before it can be properly converted into a digital signal [18].

In our study, to avoid an imprecise reconstruction of the velocity field, the distance traveled by each of the particles between two consecutive frames must be at least half the distance between two particles. As long as this is the case, the particle association between frames is guaranteed to be correct.

Table 4 shows the appropriate sampling interval for the experiment performed in this paper. Note that the high sampling rate and low particle speed make the exposition time blurry effect on particle images negligible.

Note that despite the limited values achieved compared to those of professional techniques, the results obtained with the proposed method are still successfully validated with potential flow simulation before its application to a real bidimensional experiment.

2.2. Code Development

The proposed particle tracking system is based on the development of a MATLAB code that uses image recognition to determine the trajectory and velocity of the 3D-printed tracer particles photographed. In this way, the aim is to include quantitative studies, in addition to qualitative ones, in the flow visualization tests that will be carried out in the tank. The code is structured in three interrelated parts with a user-friendly interface that allows the user to easily modify and select the processing parameters that best fit the study to be performed. Each of them is explained separately below.

The choice of MATLAB as the software to develop the tool is due to its flexibility, ease of use, universally free access to codes, and widespread documentation. MATLAB is a powerful tool to develop solutions from scratch and explore new algorithms; therefore, enabling a better understanding of the underlying mathematics and algorithms, making it ideal for research and education. Particle-tracking flow analysis enables the direct study of flow phenomena, as opposed to equation-based simulation methods. This makes it a complementary approach to more specialized software, such as Ansys and OpenFOAM, within broader fluid dynamics workflows.

The first part allows delimiting the region of interest for the study over an image and obtaining an approximate value of the radius of the particles that are in it, which will have to be adjusted later to facilitate the automatic detection of the same.

The second part sets the parameters to be used by the code in the image processing, which is crucial to ensure proper automatic detection of the particles in each image. To make it easier for the user to identify the most appropriate value of these parameters and how their variations and possible combinations affect image processing, a panel of buttons and text entries has been developed in the MATLAB App Designer application that allows them to be conveniently adjusted (Figure 1). The main parameters that affect the correct identification of particles and that can be modified are the image filter, radius range, and edge threshold.

• Image filter: Allows one to improve or modify the visualization of an image. The use of different types of filters is proposed. Gauss filter, sharpness filter or no filter.

• Radius range: Tells the program the approximate radius of the particles in the image to facilitate their detection. Thanks to the approximate value obtained in the first part of the code, there is an initial value available to begin adjusting it to the most appropriate one.

• Edge Threshold: The gradient value that a pixel must have to be considered part of the edge of a particle.

When these values are adjusted to achieve the lowest number of false positives and unidentified particles, the panel is closed, and the values are automatically saved to be used in the third part of the code. Figure 2a–c shows an example of sequential adjustment of the parameters for the correct detection of particles in the images.

Lastly, the third part automatically performs the analysis of the images. To conduct this, the code will use the values of the parameters selected in the previous part and different internal MATLAB functions. The process is as follows: the particles are identified with the internal function imfindcircles [19]. The algorithm uses the circular Hough transform (CHT). The main advantage of this tool is its robustness and, although it can adopt different approaches, its fundamental steps are as follows:

• Accumulator array computation: When a high-gradient pixel is detected, it is chosen as a candidate and allowed to cast “votes” in the accumulator array. A typical CHT pattern for casting votes is a circle of fixed radius around it (Figure 3a).

• Center estimation: The voting patterns of the candidate pixels belonging to the same circle tend to coincide in the center of the circle. By looking for this peak in the accumulator array, the center of the circle is determined (Figure 3b).

• Radius estimation: This is usually conducted in a second independent step. The function returns a matrix that collects the x and y coordinates of the centers of the detected particles and a column vector with the radii of each one.

The next step is to link each particle in one image with the corresponding particle (itself) in the next image of the sequence. To conduct this, it is necessary to establish a criterion that tells the program which particle corresponds to which. The option chosen is a proximity criterion, expressed in Equation (4).

(4)$min(dist=\sqrt{{(x_2-x_1)}^2+{(y_2-y_1)}^2})$

$x_1$, $y_1$ and $x_2$, $y_2$ are the centers of the particles in the first and second images, respectively.

The next section of the code uses a for loop to draw the vectors that go from each center of an image to the nearest one in the next image. Thus, for each center and for each image, the trajectories followed by each of the particles along the sequence of images will be able to be drawn. A coherence analysis allows for the filtering of distance measurements with abnormal values due to possible errors in identifying one particle in one image but not in the following, based on a set threshold above which distance is not gathered. At this point, and applying Equation (5), it is easy to obtain the velocity of the particles at each point using input parameters such as the sampling interval and the calculated distances.

(5)$u(x,t)={\displaystyle \frac{\Delta (x,t)}{\Delta \left(t\right)}}$

In Equation (5) $\Delta (x,t)$ is the displacement of a marker, located in position x at instant t, in a small time interval $\Delta \left(t\right)$ separating two observations.

Finally, a series of images is proposed as the output of the code to allow easy analysis and interpretation of the results in a visual way. The first proposed image superimposes the vectors joining the corresponding centers of each pair of images. In this way, we have an approximate representation of the trajectory of the particles and, therefore, of the motion of the fluid. For the velocity field representation, the region of interest is discretized by applying a Delaunay triangulation [20] using vertices as the set of centers detected in the set of images that integrate the sequence, using an internal function Delaunay.

A Delaunay triangulation is a convex-connected lattice of triangles such that the circumscribed circumference of each triangle of the lattice does not contain any vertex of another triangle. This lattice has two interesting properties:

• The minimum angle within each triangle is maximized, that is, the triangles tend to be as equilateral as possible.

• Each point in the input set will have an edge joining it to its nearest point.

This type of mesh has been used to tessellate regions in numerous studies with different types of applications [21,22], including PTV techniques for 2D and 3D experiments [23,24,25].

The code continues by constructing a vector that stores the velocity value associated with each of the centers in the image sequence. Its values are normalized between 0 and 1 and a color from a previously selected range is associated with it. In this way, each vertex of the triangulation is assigned a color representative of the fluid velocity at that point. The color of the rest of the points within each cell of the mesh in the region of interest is deduced by interpolation using the internal function patch with the color associated with its vertexes. Note, the vertices corresponding to the last image of the sequence may be excluded from the triangulation, as they do not have a distance displaced (velocity value), and therefore, the color associated with them. The process is summarized in Figure 4.

In addition to the value of the distance between the nearest centers, the relative magnitude, defined in Equation (6), of the x component is stored for each of these vectors. Therefore, another vector field can be plotted that represents the contribution of the component x to the velocity at each point in the fluid domain, differentiating the velocity zones with a greater horizontal component from those with a greater vertical component.

(6)${\displaystyle \frac{\delta x}{\delta x+\delta y}}$

where $\delta x$ is the x component of the velocity vector and $\delta y$ is the y component of the same.

2.3. Validation with Synthetic Images

Another MATLAB code provided by the Universidad Politécnica de Madrid was used to recreate the theoretical movement of particles around two vortices with opposite directions of rotation from their potential equations. Together with velocity field results, images of the theoretical position of the particles at different time steps were extracted and introduced into the software developed to compare the results obtained.

The movement of particles around two vortices with opposite directions of rotation and their theoretical velocity field has been simulated in a computer using Equations (7) and (8) extracted from non-viscous incompressible stationary flow described with classical potential theory [26].

(7)$\Phi ={\displaystyle \frac{\Gamma }{2\pi }}\theta$

(8)$\Psi ={\displaystyle \frac{\Gamma }{2\pi }}ln\left(r\right)$

where $\Phi$ is the velocity potential, $\Psi$ is the stream function, $\Gamma$ is the vortex circulation and r is the distance from its center.

Images of the position of the particles in different time steps were extracted and introduced into the developed software, to extract the trajectory and velocity field with the proposed method and compare it with the theoretical one. The results are presented in Section 3.2.

2.4. Experiment in the Towing Tank

To demonstrate the operation of the program with real images, a 2D test was carried out in the towing tank of the School of Aeronautics and Space Engineering (ETSIAE-UPM). The particles to be introduced into the fluid were printed with an Ultimaker S5 3D printer (UltiMaker, Madrid, Spain). ABS material was chosen because it has a density similar to that of water. This is a fundamental characteristic if we want to carry out future experiments in three-dimensional volumes. As stated in Equation (1), a lower diameter value would improve the advection of the particles; however, it introduces problems with respect to precision during the printing of spheres if the value is smaller than 5 mm. Therefore, the particle diameter used in this experiment is 5 mm (Figure 5). It was observed that if the particles were printed with a fill density of 100%, the particles sank into the water. Although for the proposed test with a two-dimensional sheet, this aspect is not decisive, it may be so for three-dimensional experiments. In this sense, a study of filling densities was carried out, and its results are presented later. The 2D test in this article was carried out with particles with an infill density of 20%. Although another density could have been selected, the important thing is to choose a common density for all particles to ensure similar behavior. The tank was filled with a sheet of water approximately 8 mm thick to ensure a two-dimensional movement of the particles (avoiding overlap in different planes) and a Sony (Sony DSC-QX100, Sony, Madrid, Spain) camera (Table 5) was arranged to cover a region of 0.25 m², illuminated with white light, into which 100 particles were poured. The simple experiment consisted of opening the drain and studying the movement. It follows a step-by-step summary:

• Particle printing.

• Fill the tank with a sheet of water.

• Distribute the particles in water evenly.

• Wait for 2 min for the water to settle.

• Place the camera in the right position.

• Take calibration images.

• Open the drain and start recording.

3. Results

This section details a performance analysis of the code when its parameters are modified as well as the accuracy achieved in the results with the analysis of synthetic images. Finally, the results of the experiment carried out in the Aerodynamics Laboratory tank are detailed.

3.1. Code Performance

Figure 6 shows the computation time for different sizes and numbers of particles, as well as for different numbers of images analyzed.

Figure 6a shows that the computation time remains constant with the particle size; therefore, this is not a critical parameter in the execution of the imfindcircles function as long as the range of radii is properly selected. From the error bars added, it can be observed that the variability in the result is quite high for the size of 5 mm, probably related to the fact that this value is the lower limit that allows setting the function for the detection of circles, that is, the function is being executed for its limiting value.

Figure 6b shows that the number of particles does have an important influence on the computation time, following a practically linear growth. The same happens for the number of images analyzed, Figure 6c.

The equation of the linear trend lines and, therefore, the slope value can be seen in Table 6, together with the statistic value $R^2$.

Figure 7 now analyzes how the selection of processing parameters affects the precision of the code.

A threshold of 0.1 can be used for particle identification in very uniform images generated by computer software. However, with real images, this value should be increased to avoid detection of false positives. A value of 0.3 is sufficient for images with a relatively uniform background such as that of the tank. At higher or lower values than the optimum (0.3), the number of false positives (blue) and undetected particles (orange) increases exponentially, Figure 7a. Increasing the range of radii will increase the number of false positives detected almost linearly for intervals of 1, 2, 3 and up to 4 values, Figure 7b, Table 6. It is, therefore, shown that to improve the accuracy of the code, in which images with a single particle size are used, setting a specific value for the size optimizes the identification of the particles. Figure 7c shows that the number of false positives remains low for radius values slightly larger or smaller than the particle size, allowing some flexibility in setting this value.

Finally, it is found that the uniformity of the tank background allows the processing of the images without the need to apply any type of filter to the image (Figure 7d). However, the sharp filter that improves the definition and highlights the edges of the particles gives the same good result. However, for other types of images, the application of one or more successive filters could improve the analysis.

3.2. Validation with Synthetic Images

The results obtained from the application of the code developed to images obtained with a computer-simulated experiment based on potential theory are discussed below.

Figure 8b,d,f represents the result obtained by processing 10 sequential images extracted from the simulated experiment, while Figure 8a,c,e represents the exact solution obtained in the simulation itself.

It can be seen that the trajectories and vector fields represented by both applications present important similarities. However, in the region closer to the center of the vortices, the proximity of the particles to each other gives rise to some observable association errors in the representation of the trajectories and in some points of the magnitude fields calculated.

To complete this presentation of results, Figure 9a,b includes the Delaunay triangulation used to obtain the velocity field and a superposition of the particle trajectories on this vector field, which allows us to check the accuracy of the estimation of the vortex centers by the proposed software.

Finally, in order to establish a quantitative comparison, the specific normalized value of the velocity corresponding to each particle in its last position has also been extracted in the simulation and with the code developed, Figure 10, and compared in Table 7.

The accuracy reaches 86.2% in terms of relative error, but the average error made to determine the speed is only 1.6%.

It is interesting to add the error distribution in the fluid domain, Figure 11. It can be seen that in the vicinity of high-gradient zones (center of the vortices), the error occurring is greater than in others, such as the central axis.

The origin of the discrepancy between the values obtained by the two codes may lie in the determination of the exact center of the particles in the different images (and, therefore, distances measured), influenced by variations in the brightness of the image or the shape and size of the pixels used which may affect algorithm CHT performance, Figure 12. However, these values demonstrate the validity of the developed program and the proposed experimental setup, and the number of particles used per surface region. Note that including more particles in the images analyzed would have improved the accuracy, but results already show high concordance.

3.3. Experiment in the Towing Tank

It is important to note that despite using ABS material if particles were printed with an infill density of 100%, they would sink in the water. Although for the proposed test, which is two-dimensional, this aspect is not decisive, it may be so for future three-dimensional experiments. In this sense, a study of the infill densities was carried out, as shown in Figure 13, in which the particle deposition time at the bottom of the tank was studied as a function of the infill density.

It is observed that a buoyancy of 100% of the particles is achieved with an infill density of 20% and 10%. Above 50%, they all settle to the bottom of the channel in a short time (3 s). On the other hand, for the fill densities between 50% and 20% the behavior is variable and with a longer deposition time. This would be the range within which the most suitable density for three-dimensional tests would be found. As mentioned above, floating particles with a fill density of 20% were selected for the proposed experiment. Other studies have also used floating particles for surface flow analysis using PTV [27].

Figure 14 represents the results obtained from the analysis of the video frames with the zenith view for the tank experiment.

Observing Figure 14a, it is easy to identify the region in which the drain is located, equivalent to a sink to which it would correspond with infinite velocity. One can also appreciate the exponential growth of the velocity of the particles as they approach the drain, with a significant increase in velocity in its vicinity (yellow zone).

However, it is also observed that the particles on the right (light blue) are attracted towards the drain at a higher velocity than their symmetrical counterparts on the left (dark blue). This is due to the slight inclination of the tank toward its end to favor emptying. Therefore, the movement of the particles on the right is favored by the action of gravity itself. The difference between one region and the other is 0.8 cm/s. Incorporating this system makes it possible for the first time to quantify these differences in the tank.

Referring to the representation of the velocity field with differential x and y components (Figure 14b), the radial characteristic of the flow around a sink can be appreciated, with regions of a significant horizontal component (bright yellow) in contrast to others with the most vertical component (dark red). On the far left, it is observed that the color pattern is very irregular. This is due to the low density of particles in the region, because of the low motion of the particles present and the absence of an incoming flow of them, and because of the nonaxial-symmetric behavior of the particles by the action of gravity, as already mentioned. The interpolation of data on larger triangles results in a lower accuracy of the results in that region of the fluid domain. The irregular polygonization observed in Figure 14b,d, therefore, identifies areas where the results are less reliable.

It should be noted that the characteristics of the flow tested, which involved the accumulation of a large number of particles without velocity on the drain grate, made it necessary to apply a mask over it in the images taken during pre-processing. A white rectangle was chosen to cover it, with a particle in its center to serve as the last point for the trajectory of the different particles that arrived at the drain. This avoids incorrect associations between them, although this point corresponds to a velocity value of 0 (blue) as opposed to the theoretically infinite value that is reached (yellow), Figure 14a.

The application of the code to a real experiment deserves some comments. Despite a precise selection of the processing parameters in the designed software, when analyzing a large number of images, it may appear that a particle identified in one image is not identified in the next. This would result in an incorrect association with an associated distance that may be very different from the real one if no other particle is found in the vicinity (the result would improve with increasing particle density). This would affect not only the representation of trajectories but also the velocity field, with velocity peaks in places where there are none. In the first step to introduce the filtering of these erroneous results, a threshold value was set for the measured distances (90 in this test) above which their value would be set to 0, thus solving the occurrences detected in the proposed example.

These occurrences were detected near the edges, far from the sink, where the velocities are low, which explains the gaps at the beginning of the trajectories of some of those identified in Figure 14c and dark areas in the same region in Figure 14b. Adjusting this threshold value will depend on the characteristics of the experiment in question; modifying the value 0 by another value that could be more accurate from the analysis of the velocity of the surrounding particles is proposed for future developments.

Accurate trajectories will be achieved as long as the Nyquist criterion defined in Table 4 is met. When two different particles are closer together than their displacements between two frames, the code may confuse the trajectories. This can be a problem in areas that tend to concentrate the particles, such as the center of the vortex.

As a sample, an illustration of the trajectory obtained for one of the particles in a previous analysis is included in Figure 15. Observing this type of phenomenon makes it possible to identify that the motion of these particles is below the Nyquist limit, and the sampling frequency should be increased, if possible, to solve it. This trajectory was obtained with a sampling frequency of 40 Hz, while those represented in Figure 14c were obtained at 1 Hz.

This experiment serves to demonstrate the capabilities of the developed code, as well as an example of its future applications, its potential, and its current limitations.

The major current limitation is the inability of the code to distinguish between particles in different planes. This limitation prevents its direct application to three-dimensional tests.

4. Conclusions

The aim of this article was to demonstrate the feasibility of developing a high-accuracy system for the representation of the velocity field with minimal investment in equipment. This system is intended for application in university towing tanks for both research and educational purposes.

This objective was achieved using a system based on the theory of Particle Tracking Velocimetry (PTV), consisting of tracking ABS particles of 5 mm diameter, 3D printed in-house, through image capture and simple, dynamic post-processing facilitated by a user-friendly code developed in MATLAB.

To evaluate the accuracy of the proposed system, we performed a simulation of particle movement and the representation of the velocity field using computer software based on the theory of potential flow. Synthetic images of the position of the particles at different instants were extracted and processed with the developed MATLAB code. The results were then compared with those of the simulation.

The developed code demonstrated a relative accuracy of 86% with these synthetic images. The computation time for sequences of 10 images did not exceed 10 seconds, and an appropriate selection of input parameters effectively prevented false positives and undetected particles. These results highlight the good performance of the proposed system.

To demonstrate its performance, a real two-dimensional test was performed on the university towing tank, resulting in successful results. The velocity field in the water tank was successfully represented and differences in the emptying velocity between the left and right regions of the drain, attributed to the inclination of the tank, were estimated at 0.8 cm/s. Hence, we can conclude:

• ABS particles with a diameter of 5 mm can be effectively tracked to represent flow motion within a towing tank.

• The post-processing of particle images can be performed using a simple and user-friendly MATLAB code requiring minimal inputs and offering acceptable computation times.

• The system demonstrated a relative accuracy of 86%, making it applicable to a wide range of 2D experiments.

However, some important limitations must be considered and addressed in future research, and potential solutions are proposed for each limitation.

One major limitation is the inability of the system to distinguish particles located in different planes, restricting its application to two-dimensional tests. The use of a laser sheet is a conventional solution frequently mentioned in the literature. Alternatively, future studies could employ simultaneous image capture from two different positions, correlating the results to achieve a three-dimensional representation. Additionally, using particles of different colors could help identify erroneous trajectories and improve associations between frames. For the calibration of a future 3D system, studying a well-determined flow over a cylinder is proposed as a reference case.

Furthermore, according to the Nyquist criterion, in regions where the particle spacing is less than 2.5 mm, correct associations are not guaranteed, potentially leading to velocity computation errors. Increasing the sampling frequency could improve this limitation. It is also recommended to develop new algorithms to improve the discrimination of erroneous particle pairings beyond simply setting a threshold.

These limitations, along with a Stokes number of approximately 0.1 and a Poisson probability of 0.19, indicate that the performance of the system remains below that of professional techniques. However, the use of alternative materials or industrial particles, instead of the 3D printed ones used in this study, is believed to have the potential to further improve accuracy and particle behavior in the fluid.

Despite these constraints, the primary advantages of the proposed system are its cost-effectiveness and ease of application, which is in agreement with the benefits of using water tunnels for aerodynamic investigations. In contrast, a traditional PIV application would require equipment investments ranging from €10,000 to €100,000 for a professional setup. Additionally, the higher qualifications required for laser handling and microscopic particle manipulation, as well as the significant investment in a towing tank capable of integrating such PIV systems (approximately €200,000), must be considered.

In conclusion, the proposed system successfully met the expected requirements, demonstrating its viability as a low-cost alternative to velocity field analysis in educational and research contexts.

Author Contributions

Conceptualization, J.J.A.-M., Á.A.R.-S., M.J.C.-C., R.B.-M. and A.L.-M.; methodology, J.J.A.-M., Á.A.R.-S., R.B.-M., J.C.G.-M., E.B.-B., A.L.-M. and D.V.; software, J.J.A.-M. and A.L.-M.; validation, J.J.A.-M., Á.A.R.-S., M.J.C.-C., R.B.-M., J.C.G.-M., E.B.-B. and A.L.-M.; formal analysis, J.J.A.-M., Á.A.R.-S., R.B.-M., J.C.G.-M., E.B.-B. and A.L.-M.; investigation, J.J.A.-M., Á.A.R.-S., M.J.C.-C., R.B.-M. and A.L.-M.; resources, J.J.A.-M., Á.A.R.-S., R.B.-M. and M.J.C.-C.; data curation, J.J.A.-M., A.L.-M. and D.V.; writing—original draft preparation, Á.A.R.-S., R.B.-M. and A.L.-M.; writing—review and editing, J.J.A.-M., Á.A.R.-S., M.J.C.-C., R.B.-M., J.C.G.-M. and E.B.-B.; visualization, J.J.A.-M., Á.A.R.-S., M.J.C.-C., R.B.-M. and A.L.-M.; supervision, J.J.A.-M., Á.A.R.-S., M.J.C.-C. and R.B.-M.; project administration, Á.A.R.-S., M.J.C.-C. and R.B.-M. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are contained within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

The following abbreviations are used in this manuscript:

LSV	Laser Speckle Velocimetry
PIV	Particle Image Velocimetry
PLV	Pulsed Laser Velocimetry
PTV	Particle Tracking Velocimetry

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Button panel for the selection of processing parameters.

Figure 2. Example of 15-pixel particle detection setting: (a) maximum radius 5, minimum radius 5, threshold 0. (b) maximum radius 5, minimum radius 5, threshold 0.3. (c) maximum radius 15, minimum radius 15, threshold 0.3.

Figure 3. (a) Candidate pixel placed in a real circle (solid circle) next to the classical voting pattern (dashed circle). (b) Different voting patterns coincide at the center of the real circle. Extracted from [19].

Figure 4. Process summary.

Figure 5. (a) 3D printing of particles. ABS material. (b) Zenith view from a section of the tank’s bottom including the drain and suspended particles in a laminar sheet of water.

Figure 6. Computing time as a function of: (a) particle size; (b) number of particles; (c) number of images. Reference experiment: 9 images with 36 particles of 8 pixels. The computation times in (a) correspond to the identification time on a single image.

Figure 7. Number of errors (false positives, blue; undetected particles, orange) as a function of varying the value of different processing parameters: (a) edge threshold, (b) range of radii, (c) search size and (d) filter applied to the image. Optimal parameters for the example: threshold 0.3, range 0, size 15, no filter.

Figure 8. Comparison of results obtained from computer simulation and processing with the proposed code of synthetic images extracted: (a,b) trajectories, (c,d) velocity field (total magnitude) and (e,f) velocity field (distinction between x and y components).

Figure 9. (a) Delaunay triangulation used to discretize the fluid domain. (b) Checking the estimation of vortex centers.

Figure 10. Comparison of the velocity values obtained for each particle at their last position in the computer simulation (a) and after analyzing images 9 and 10 with the proposed code (b).

Figure 11. Distribution of the absolute error occurred in the processing of synthetic images.

Figure 12. Sources of error in particle center determination.

Figure 13. Particle deposition time for different infill densities. ABS material.

Figure 14. Results obtained after image processing of the tank emptying experiment: (a) velocity field (total magnitude), (b) velocity field (distinction between x and y components), (c) trajectories and (d) Delaunay triangulation.

Figure 15. Error in tracking the trajectory of a particle due to its proximity to another particle.

Table 1

Characteristics magnitudes of the particles used in the test.

Density (kg/m3)	Diameter (m)	Water Dyn. Visc. (kg/(ms))	Relaxation Time (s)	Charact. Velocity (m/s)	Charact. Length (m)	Stokes Number
1000	0.005	0.001	1.39	0.01	0.1	0.139

Table 2

Characteristic magnitudes of the test.

Particle Num.	Study Region Area (m2)	Area Unit (m2)	Num. Area Units
100	0.25	6.25 × ${10}^{-4}$	400

Table 3

Probability of finding a particle in a unit area during the test.

C 1	V 1	k	Prob(k)
0.25	1	1	0.19

¹ Since the test proposed is 2D, units of volume should be replaced by units of surface area.

Table 4

Nyquist criterion for the test performed.

Sampling Rate (Hz)	Time (s)	Average Speed (m/s)	Average Distance (m)	Limit (m)
8	0.125	0.01	1.25 × ${10}^{-3}$	2.5 × ${10}^{-3}$

Table 5

Photo-taking device characteristics.

Brand	Model	Sensor (Number of Pixels)	Application
Sony	DSC-QX100	20.2 MP	Zenithal

Table 6

Trend lines and $R^2$ values.

Chart	Trend Line	${\mathit{R}}^2$
Figure 6a	y = −0.0004x + 0.0903	0.65
Figure 6b	y = 0.075x + 1.855	0.995
Figure 6c	y = 0.4542x + 0.3087	0.999
Figure 7b	y = 10.9x − 2	0.95931

Table 7

Errors made with the proposed code.

Absolute Error	Mean Absolute Error	Relative Error	Mean Relative Error
4.15%	1.6%	13.8%	5%

References

1. Lorincz, D.J. A Water Tunnel Flow Visualization Study of the F-15; NASA STI/Recon Technical Report NASA-CR-144878 NASA, Dryden Flight Research Center: Edwards, CA, USA, 1978.

2. Pullin, D.I.; Perry, A.E. Some flow visualization experiments on the starting vortex. J. Fluid Mech.; 1980; 97, 239. [DOI: https://dx.doi.org/10.1017/S0022112080002546]

3. Erm, L.P.; Ol, M.V. An Assessment of the Usefulness of Water Tunnels for Aerodynamic Investigations; Defence Science and Technology Organisation: Melbourne, VIC, Australia, 2013.

4. Rodríguez-Sevillano, A.A.; Casati-Calzada, M.J.; Bardera-Mora, R.; Ballesteros-Grande, L.; Martínez-García-Rodrigo, L.; López-Cuervo-Alcaraz, A.; Fernández-Antón, J.; Matías-García, J.C.; Barroso-Barderas, E. Exploring the Effectiveness of Visualization Techniques for NACA Symmetric Airfoils at Extremely Low Reynolds Numbers. Fluids; 2023; 8, 207. [DOI: https://dx.doi.org/10.3390/fluids8070207]

5. Von Funck, W.; Weinkauf, T.; Theisel, H.; Seidel, H.-P. Smoke Surfaces: An Interactive Flow Visualization Technique Inspired by Real-World Flow Experiments. IEEE Trans. Vis. Comput. Graph.; 2008; 14, pp. 1396-1403. [DOI: https://dx.doi.org/10.1109/TVCG.2008.163] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18988989]

6. Willert, C.; Waltenspül, S.; Schanz, D.; Kompenhans, J.; Schulze, M. Prandtl’s flow visualization film C1 revisited. Proceedings of the 13th International Symposium on Particle Image Velocimetry; Munich, Germany, 22–24 July 2019.

7. Werle, H. Hydrodynamic Flow Visualization. Annu. Rev. Fluid Mech.; 1973; 5, pp. 361-386. [DOI: https://dx.doi.org/10.1146/annurev.fl.05.010173.002045]

8. Clayton, B.; Massey, B. Flow visualization in water: A review of techniques. J. Sci. Instruments; 2002; 44, 2. [DOI: https://dx.doi.org/10.1088/0950-7671/44/1/302]

9. Adrian, R.J. Particle-Imaging Techniques for Experimental Fluid Mechanics. Annu. Rev. Fluid Mech.; 1991; 23, pp. 261-304. [DOI: https://dx.doi.org/10.1146/annurev.fl.23.010191.001401]

10. Brossard, C.; Monnier, J.C.; Barricau, P.; Vandernoot, F.X.; Sant, Y.L.; Champagnat, F.; Besnerais, G.L. Principles and Applications of Particle Image Velocimetry. Aerospace Lab.; 2009; 1, pp. 1-11.

11. Sokoray-Varga, B.; Józsa, J. Particle tracking velocimetry (PTV) and its application to analyse free surface flows in laboratory scale models. Period. Polytech. Civ. Eng.; 2008; 52, 63. [DOI: https://dx.doi.org/10.3311/pp.ci.2008-2.02]

12. Kowalczyk-Hernandez, M. Laser speckle velocimetry. Proc. SPIE; 1996; 2729, pp. 139-145. [DOI: https://dx.doi.org/10.1117/12.233000]

13. Tauro, F.; Tosi, F.; Mattoccia, S.; Toth, E.; Piscopia, R.; Grimaldi, S. Optical Tracking Velocimetry (OTV): Leveraging Optical Flow and Trajectory-Based Filtering for Surface Streamflow Observations. Remote Sens.; 2018; 10, 2010. [DOI: https://dx.doi.org/10.3390/rs10122010]

14. Schmidt, B.E.; Sutton, J.A. High-resolution velocimetry from tracer particle fields using a wavelet-based optical flow method. Exp. Fluids; 2019; 60, 37. [DOI: https://dx.doi.org/10.1007/s00348-019-2685-6]

15. Crowe, C. Multiphase Flow Handbook; Mechanical and Aerospace Engineering Series; CRC Press: Boca Raton, FL, USA, 2005.

16. University of Cambridge. 2004–2024. Available online: https://www.doitpoms.ac.uk/tlplib/powder/stokes_number.php (accessed on 8 October 2024).

17. Tropea, C.; Yarin, A.; Foss, J. Springer Handbook of Experimental Fluid Mechanics; Springer: Berlin/Heidelberg, Germany, 2007; [DOI: https://dx.doi.org/10.1007/978-3-540-30299-5]

18. Peak Studios. 2024; Available online: https://www.peak-studios.de/es/nyquist-frequenz/ (accessed on 8 October 2024).

19. The MathWorks, Inc. 1994–2024. Available online: https://es.mathworks.com/help/images/ref/imfindcircles.html (accessed on 8 October 2024).

20. Delaunay, B. Sur la sphère vide. Bull. L’Académie Sci. L’URSS. Cl. Sci. Mathématiques Na; 1934; 1934, pp. 793-800.

21. Weatherill, N.P. Delaunay Triangulation in Computational Fluid Dynamics. Comput. Math. Appl.; 1992; 24, pp. 129-150. [DOI: https://dx.doi.org/10.1016/0898-1221(92)90045-J]

22. Cho, H.G.; Kim, S.S.; Kim, Y.S.; Cho, M.G. A Geometric Compression Algorithm For Massive Terrain Data Using Delaunay Triangulation. 1998. Proceedings of the WSCG ’99: The 7th International Conference in Central Europe on Computer Graphics and Visualization and Interactive Digital Media ’99; Plzen, Czech Republic, 8–12 February 1999; Volume 1.

23. Ezra, E.; Keinan, E.; Liberzon, A.; Nahmias, Y. Development of Three-Dimensional Streamline Image Velocimetry Using Superimposed Delaunay Triangulation and Geometrical Fitting. J. Fluids Eng. Trans. ASME; 2016; 138, 011205. [DOI: https://dx.doi.org/10.1115/1.4031611]

24. Azadi, R.; Wong, J.; Nobes, D.S. Determination of fluid flow adjacent to a gas/liquid interface using particle tracking velocimetry (PTV) and a high-quality tessellation approach. Exp. Fluids; 2021; 62, 48. [DOI: https://dx.doi.org/10.1007/s00348-020-03103-5]

25. Chen, K.; Anthony, S.M.; Granick, S. Extending Particle Tracking Capability with Delaunay Triangulation. Langmuir; 2014; 30, pp. 4760-4766. [DOI: https://dx.doi.org/10.1021/la500323r] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/24734998]

26. Meseguer Ruiz, J.; Sanz Andrés, P. Aerodinámica Básica; 2nd ed. Garceta Grupo Editorial: Madrid, Spain, 2010.

27. Tang, H.-w.; Chen, C.; Chen, H.; Huang, J.-t. An improved PTV system for large-scale physical river model. J. Hydrodyn. Ser. B; 2008; 20, pp. 669-678. [DOI: https://dx.doi.org/10.1016/S1001-6058(09)60001-9]