Reconstruction of large-scale flow structures in a stirred tank from limited sensor data

We combine reduced order modeling and system identification to reconstruct the temporal evolution of large-scale vortical structures behind the blades of a Rushton impeller. We performed direct numerical simulations at Reynolds number 600 and employed proper orthogonal decomposition (POD) to extract the dominant modes and their temporal coefficients. We then applied the identification algorithm, N4SID, to construct an estimator that captures the relation between the velocity signals at sensor points (input) and the POD coefficients (output). We show that the first pair of modes can be very well reconstructed using the velocity time signal from even a single sensor point. A larger number of points improves accuracy and robustness and also leads to better reconstruction for the second pair of POD modes. Application of the estimator derived at Re = 600 to the flows at Re = 500 and 700 shows that it is robust with respect to changes in operating conditions.

To the best of our knowledge, only one study has applied POD to the three-dimensional (3D) flow field in a stirred tank. 14 The author used Large Eddy Simulation (LES) to simulate the flow agitated by a three-blade pitched turbine using a combined rotating/static mesh approach. In total, 1500 snapshots of the 3D velocity field were stored over four rotations. The POD was carried out separately for the inner (rotating) domain surrounding the impeller and the outer In the present study, we exploit this compact description of the flow, but we go a step further. In particular, we used system identification to construct a linear estimator of the temporal coefficients of the leading order modes using as inputs the velocity signals at one or more sensor points in the flow. The estimator has the form of a multiple-input multiple-output linear, time-invariant system. The system matrices and the order of the system are determined optimally using the N4SID algorithm. 15 The latter belongs to the family of subspace identification methods; the reader may refer to Reference 16 for a detailed description and to Reference 17 for an overview.
Although system identification has been applied to many different areas (see examples in Reference 16), its application to fluid flows is relatively recent. It was first applied by Guzm an-Iñigo et al. 18 to derive an estimator for the dynamics of perturbations in a laminar boundary layer flow. The perturbations were assumed to have very small amplitude, so their evolution was governed by the Navier-Stokes equations, linearized around the Blasius velocity profile. The assumption of small perturbations was relaxed in Reference 19, where the method was applied to reconstruct nonlinear flows exhibiting limit-cycle oscillations. In particular, the strongly periodic flow around a 2D aerofoil was considered. POD was used to decompose the perturbations around the time-average flow, and the temporal coefficients of the two leading modes were estimated very accurately (with fit values above 90%) using the signals of the two velocity components at a single probe point in the wake.
In both references cited earlier, the flow was 2D and laminar. In the present study, we apply system identification to the complex 3D flow inside an unbaffled stirred tank operating at the transitional Reynolds number of 600; this setup was also examined by Tamburini et al. 20 Although the majority of studies on stirred tanks have focused on baffled configurations due to their enhanced mixing properties, unbaffled tanks are widely used in industry due to the ease of cleaning and lack of dead zones behind the baffles when run at lower Reynolds numbers. We chose a Reynolds number of 600 for two reasons. First, flows in stirred tanks at low-to-moderate Reynolds number have important industrial applications, for example, in bioreactors where the shear rate has to be kept low to avoid cell damage. Second, at this Reynolds number, the transitional flow can be modeled accurately using Direct Numerical Simulaton (DNS). This is not critical for the application of the system identification algorithm; since we are interested in the dominant pair of modes, LES would have also been appropriate. This is however the first time the method is applied to complex transitional 3D flows, and we decided to employ an accurate simulation method. As will be seen later, this also allows us to investigate the asymptotic behavior of higher order POD modes.
The ability to estimate the temporal coefficients of the dominant POD modes from velocity signals at a very small number of sensor points opens new opportunities. For example, one can construct the estimator at nominal operating conditions and then apply it at offdesign conditions. In this case, the approach has predictive power, as only the input velocity signals are necessary in order to reconstruct the entire instantaneous velocity field of the large-scale structures at the new conditions. We explore the accuracy of the reconstructed flow at Re numbers 500 and 700.
A word of caution regarding the reconstructed velocity perturbations for transitional and/or turbulent flows from limited velocity measurements is warranted here. Such flows exhibit a wide range of spatial scales, the largest of which are encapsulated in the dominant POD modes. The reconstruction framework is therefore expected to work better for large scales that have a strong footprint at the probe locations. Smaller scales have a more narrow extent and are more difficult (or even impossible) to recover. For this reason, in this article, we aim at reconstructing the large-scale 3D vortical structures around the impeller blades.
We close this section by stressing that system identification can be used also for process optimization, such as mixing enhancement. For example, if we know that a control variable (such as injection velocity) affects the mixing performance in a stirred tank, then we can use this as an input variable and the POD coefficients of the scalar concentration as the output, and construct a linear, time-invariant system, which can be used for mixing optimization. It is very difficult to solve this optimization problem directly using the 3D unsteady velocity and scalar fields, but system identification offers a more tractable approach.
This article is organized as follows: the flow conditions, numerical methodology, and extraction of POD modes are detailed in Sections 2 and 3, respectively. In Section 4, we show that the form of the identified multi-input, multi-output system arises naturally from the flow equations and estimation theory. In the same section, we also present the steps for the application of N4SID algorithm. Results at design and off-design conditions are reported in Section 5, and we conclude in Section 6.

| FLOW CONDITIONS AND COMPUTATIONAL DETAILS
The flow configuration considered is that of an unbaffled top-covered stirred tank, agitated by a six-blade Rushton impeller. A schematic is shown in Figure 1. The fluid is assumed to be incompressible and the Reynolds number, defined in the usual way as Re = ρND 2 /μ, where ρ is the density, N the impeller speed, and D the impeller diameter, is 600. The same geometry has been simulated with DNS by Tamburini et al. 20 and investigated experimentally by Scargiali et al. 21 ; the physical dimensions are shown in Table 1. Simulations for two additional cases, with Re = 500 and 700, were conducted in order to assess the ability of the reduced order model that was identified for Re = 600 to reconstruct the flow structures at off-design conditions.
Due to the absence of baffles, the geometry has rotational symmetry and it is more convenient to express and solve the flow equations in a coordinate system that rotates with the blades. In this system, shown in Figure 1, the impeller is held at a fixed position and the governing equations take the form where u is the velocity vector, t the time, p the static pressure, ν the kinematic viscosity, x = (x, y, z) the position vector, and finally ω is the angular velocity vector of the inertial coordinates. In the present case, this vector has a nonzero component only in the axial z direction, that is, ω = (0, 0, ω), which indicates that the rotation takes place in the (x, y) plane and ω is the angular velocity of the impeller. The two terms within square brackets in the momentum Equation (1a) represent the centrifugal and Coriolis acceleration terms. A Dirichlet boundary condition, u w = Àω Â x, is imposed at the top, bottom, and peripheral walls. This is known as the "Frozen Rotor" approach 22 and is widely used when the geometry has rotational symmetry.
The above set of equations is solved using our in-house code "Pantarhei." This is a finite volume, unstructured grid solver, which has been used extensively in simulating transitional and turbulent flows in various configurations, including stirred vessels. [23][24][25][26][27][28] The convection terms are discretized in space using a second-order central scheme. The orthogonal viscous terms are treated implicitly, while the convection, nonorthogonal viscous, Coriolis, and centrifugal terms are all treated explicitly in time using quadratic extrapolation. A thirdorder accurate scheme is used to discretize the transient term. The pressure-velocity coupling is treated using the fractional step method.
The code is parallelized with the aid of PETSc 29 and HYPRE libraries. 30 The computational domain was discretized using 100 cells in the radial direction, 516 in the tangential direction, and 149 in the vertical (axial) direction; the total number of cells was $7.6 Â 10 6 . This is between the third and fourth finest grids considered in Reference 20, with the authors choosing the third finest to perform their simula- To check the grid independence of the results, we also performed simulations with a finer grid, consisting of 120, 636, and 192 cells in the radial, tangential, and vertical directions respectively, resulting in a total number of $19.6 Â 10 6 cells. In Figure 2, we compare the profiles of time-averaged radial, tangential, and axial velocities along a line bisecting two blades at a height of 0.067 m, a few mm above the impeller disk, for the two grids. As can be seen, the matching of the time-averaged profiles is very good.
F I G U R E 1 Sketch of the stirred tank and the impeller geometry. We also checked the resolution with respect to the Kolmogorov length scale, η = (ν 3 /ϵ) 1/4 , which is the smallest scale of the flow. 31 The local ratio of the grid spacing to η was estimated from V 1=3 cell =η cell , where V cell is the cell volume. For the lower grid resolution, this ratio was found to be less than 2.0 everywhere, while the average value in the domain was 0.49. These metrics satisfy fully the resolution requirements for DNS 31 and are comparable to those of other DNS studies in stirred tanks. 20,32 As a final check, we computed the power number, N p ¼ P ρN 3 D 3 , where P is the time-average of total energy dissipation (i.e., due to mean and fluctuating fields) integrated over the tank volume. The value was found to be 2.1, which is close to 2.2 reported in References 20,21. Slight differences in geometry (e.g., the blade thickness) are likely to account for the small discrepancy. All the abovementioned tests confirm that the lower resolution grid is deemed sufficient, and this grid was used to obtain all the results presented henceforth.

| REDUCED ORDER MODEL OF THE FLUCTUATIONS AROUND THE TIME-AVERAGE FLOW
We use POD to construct a reduced order model of the velocity fluctuations around the mean flow. POD is a modal decomposition technique that yields a set of orthogonal modes that optimally represent the kinetic energy of the perturbation field. Below we provide only a brief summary; for more details refer to the book by Holmes et al. 33 and the review by Taira et al. 34 Using the Reynolds decomposition, the velocity vector u(  (1), are provided as follows: where a i (t) is the temporal coefficient, ϕ i (x) is the spatial distribution of mode i, and n is the order of the truncated expansion. We define the inner product between two fluctuating vectors u 0 i and u 0 j as where T denotes the transpose operation and V the tank volume. If the norm is equal to twice the instantaneous turbulent kinetic energy integrated in the domain. The modes ϕ i (x) are obtained by maximizing the average projection of u 0 (x, t) to ϕ i (x); the projection is defined using the inner product (Equation 4). It can be shown 33 that this optimization problem can be reduced to an eigenvalue problem, which we solve using the method of snapshots. More specifically, we store the velocity fluctuations around the time- and then stack K of these vectors column-by-column to form the It is shown in References 33,34 that the POD modes can be obtained from the singular value decomposition (SVD) of the matrix where V is a 3 m Â 3 m weighting matrix that accounts for the grid nonuniformity. For a second-order finite volume discretization, V is a diagonal matrix with the cell volume in the main diagonal. In Equation (7), Σ is a diagonal matrix containing the singular values σ i , while the left and right singular vectors are stored in matrices Φ and Ψ, respectively. The eigenvalues λ i of the POD problem are λ i ¼ σ 2 i =K and represent the energy content of mode i, the spatial modes ϕ i (x) are the columns of V À1/2 Φ, and the temporal coefficients a i (t j ) The process to extract the POD modes is the following. For all Reynolds numbers, the simulations first run for 40 impeller revolutions to allow the flow to reach a statistically steady state (10 probe points located in several areas of the tank were used to verify that this is the case). For the main Re = 600 case, the flow was simulated for another 120 rotations. During the simulation, we saved approximately 100 snapshots of the velocity field in the whole domain per rotation (i.e., one every 57 time steps). This provided a dataset of 12,000 snapshots from which the spatial distribution of the POD modes ϕ i (x) (i = 1…n) was computed using SLEPC, 35,36 a library that provides parallel subroutines for solving SVD problems efficiently.
Details about the temporal and spatial mode distributions are provided in Section 5.1.
From the above, it is clear that the same grid used for the flow simulation was also used for the POD mode extraction, that is, we did not down-sample. The size of the dataset is determined by two factors. First, it must be longer (by at least an order of magnitude) than the period of the flow features being examined. Second, there must be enough data to both train and validate the estimator. As such, a dataset of 20-40 rotations would have been sufficient for extracting the leading order modes and model training. We would then need additional 20-40 revolutions for model validation. In the present study, we used a conservatively long dataset (120 rotations) to ensure that we capture slowly varying features also (although these are not elaborated on in this study).
We now define the state vector of the temporal POD coefficients a t ð Þ ¼ a 1 t ð Þ,ÁÁÁ, a n t ð Þ ð Þ T and obtain its evolution equation by performing Galerkin projection, that is, by substituting the expansion (3) into the fluctuating equation set (2) and taking into account that the basis is divergence-free and orthonormal. The resulting system takes the form where matrix A arises from the projection of the linear terms on ϕ i while the forcing term F(t) arises from the nonlinear terms that appear in the right-hand side of (Equation 2a) as well as the fact that we retain only a finite number of POD modes (for more details on the derivation of Equation (8) see Reference 19). This equation represents a reduced-order model of the fluctuations. In the next section, we describe how to obtain an estimator that approximates the evolution of a(t) using information from a few sensor points in the flow. Note that we never compute the matrix A explicitly, all we need to know is that this matrix exists.

| Form of a linear dynamic estimator
We seek an estimate a e (t) of the state vector a(t) from a set of velocity Þ T at p sensor points located at Þwithin the stirred tank. From expansion (3), we can write s , where matrix S depends on the sensor locations and the velocity components being recorded, while g(t) is noise arising from the fact that only a finite number of POD modes are retained in the state vector a(t). From linear estimation theory, 37 it is known that an optimal estimation a e (t) can be obtained from where matrix ℒ , known as Kalman filter gain, is obtained from the solution of a Riccati equation. The forcing term in this equation involves the difference between the actual measurements, s(t), and their estimations, Sa e (t). We can rearrange Equation (9) to get The above equation represents a linear state-space system with input s(t) and output a e (t). Note that S can be computed from ϕ i , but the matrices A, ℒ are not explicitly known. Equation (10) motivates the more general form of the estimator where we have introduced the state vector X e (t). This is again a statespace system with input s(t) and output a e (t) and becomes identical to Note that the order of system (11), that is the size of the vector X e (t), is not necessarily equal to the size n of a e (t), which gives additional flexibility (in this case H is not a square matrix). Again matrices A, B, and H are not explicitly known and we seek to compute these matrices from the discrete values of the POD coefficients a(t j ) and the sensor data s(t j ) using system identification; this process is explained later.

| System identification from DNS data
Since the available data a(t j ) and s(t j ) from the DNS simulations are known at the discrete time instants t j , we write system (11) in discrete form as where the primes indicate evolution matrices from t j to t j + 1 , for example, A 0 = e Aδt , where δt = t j + 1 À t j . Subspace system identification methods extract the matrices A 0 , B 0 , and H 0 from the input, s(t j ), and output, a(t j ), data such as the estimation a e (t j ) optimally approximates a(t j ) over the length of the training dataset (see below). For system (12), this is known as deterministic identification. For stochastic systems, as in our case, noise terms w(t j ) and v(t j ) are included, that is, and stochastic system identification methods can also provide the covariance matrices, Q 0 , S 0 , and R 0 , defined from where E denotes the expectation operator and δ pq is the Kronecker delta. This equation indicates that w(t p ) and v(t p ) are assumed to be white noise sequences.
The identification process consists of three procedural steps, which are shown schematically in Figure 3. In the first step, the DNS simulations are performed to generate the input and output data, s(t j ) and a(t j ), respectively. In the second step, the linear model matrices are computed by maximizing the fit between the output of the system and the model for part of the original data, known as training dataset.
In the third step, the performance of the generated model is assessed on a validation dataset, which is different to the training dataset. Finally, the percentage fit between the modeled and the true output for the ith mode was evaluated from where jjÁjj is the 2-norm of a vector, the subscript e denotes the estimated output, and the overbar denotes time-average (here equal to 0). A large value of FIT i (max 100) indicates good matching between a i (t) and a i,e (t). to see that the slow decay of the energy content in Figure 4A follows a À11/9 power law for mode numbers between ≈130 and 600 (shown as a dashed green line). This power law for POD modes was first derived theoretically in Reference 40 from the well-known κ À5/3 law for the energy density with respect to wavenumber κ for homogeneous isotropic turbulence. This law has also been observed evidence of structure and thus are far from the previously mentioned À11/9 power law, which is valid for homogeneous isotropic turbulence (refer to Figure 4A). As the mode number increases, the power spectrum becomes increasingly more noisy, with smaller energy content and less clear dominant frequencies, as evident for Mode 61. For higher modes, such as 252, the whole spectrum is filled and there is no dominant frequency, suggesting that the mode represents homogeneous isotropic turbulence and hence obeys the À11/9 power law.

| Characterization of POD modes
We examined also the probability density function plots of a i (t) (not shown) and found that these morph smoothly to Gaussian distribution with increasing mode number, as expected.

| Reconstruction of the first and second modes from point data
In this section, we consider the quality of reconstruction of the largescale structures behind the impeller blades that correspond to Modes 1 and 2, that is, the output vector in Equation (12) which consists of the temporal coefficients a e (t) = [a 1,e (t), a 2,e (t)] T . We start by employing as input the velocity signals from a single probe point, so s (t) = [u 1 (t), v 1 (t), w 1 (t)] T , and proceed to explore the effect of including information from more points. More specifically, we consider between two to six points; in the last case, there is one point in the wake of each blade, that is, s(t) = [u 1 (t), v 1 (t), w 1 (t), …, u 6 (t), v 6 (t), w 6 (t)] T .

| Reconstruction using the velocity signals from a single sensor point
The location of the probe point is important for the quality of reconstruction. 19 In the present work, this point was chosen to be at the position of maximum turbulent kinetic energy and is indicated with a black dot (Point 1) in Figure 5. in Equation (15). The results are shown in Figure 11A for both the training and the validation datasets (in computing the fit for the latter, we have excluded the first 500 data points to allow the model to "spin-up" from the initial zero state, see Figure 10). It can be seen that, for the training dataset, the fit for model orders between 2 and 9 is consistently over 80%. However, the performance of the generated models when tested on the validation dataset shows significant variations. The optimal fit values, 77% and 76%, are obtained for model orders 6 and 8, respectively. MATLAB identifies the eighthorder model as optimal. This is determined by examining the variation  Figure 11A.
It can be seen that the singular values drop by two orders of magnitude between the eighth and ninth order, indicating that orders higher than 8 do not offer significant advantage.
In Figure 12, we plot the spatial distribution of the velocity RMS due to the first two modes, , computed using the true (left panel) and the estimated (right panel) temporal coefficients a 1 (t) and a 2 (t). It can be seen that although there are slight deviations between the true and estimated coefficients over time, as evidenced in Figure 10, the time-average distributions are almost indistinguishable.

| Effect of the number of sensor points on the quality of reconstruction
We now investigate the quality of reconstruction as a function of additional probe points. To this end, we exploit the 60 rotational symmetry of the system and consider 2, 3, 4, 5 and 6 points placed at the locations of maximum turbulent kinetic energy. The points are located at the same radial distance from the impeller axis, but the azimuthal angle of each sensor location increases by 60 in the anticlockwise direction with respect to the preceding location; all probe points are indicated by black dots in Figure 5. The fit for the training dataset is very high for all cases considered. For example, using six probes, the fit is almost 90% as shown in Figure 11B. For the verification dataset, the fit reaches almost 80% and remains almost unaffected by the model order. Note also that a high fit value can be obtained even with very small models that contain only two states. This is in contrast to the performance with a single probe, which requires a model of higher order, as shown in Figure 11B. In general, as the number of sensor points increases, the optimal model order decreases; at the same time, the fit increases and the estimators become more robust.

| Reconstruction of the third and fourth modes
The ability of the identification algorithm to estimate the third-and fourth-order modes was also examined. For the training step, the same velocity signals from the same six probe points as in the previous section were employed as inputs and the temporal coefficients a 3 (t) and a 4 (t) were used as outputs. For the training datasets, the model fits were good, with values ≈85%. However, the performance on the validation dataset was lower, with fits typically around 50%.
The frequency of oscillation of the temporal coefficients, equal to 3F N (refer to Figure 7), was predicted very well, with the loss of accuracy coming from the error in the prediction of the local amplitude. The results for the training and validation datasets for the optimal eighthorder model are shown in Figure 13.
The poor fit in the validation dataset is probably related to the fact that the selected sensor points do not correspond to the locations of the largest amplitude of Modes 3 and 4, as can be clearly seen from Figure 5C,D. Additionally, the energy content of these modes is substantially smaller (about 5 times) than that of the first pair, refer to Figure 4A (this can also be noticed by comparing the peaks at frequencies F N and 3F N in the velocity spectrum plot shown in Figure 9).
The velocity RMS due to Modes 3 and 4, ffiffiffiffiffiffiffiffiffiffiffiffiffi u0 3&4 2 q , computed using the true and the estimated temporal coefficients, a 3 (t) and a 4 (t), are shown in the left and right panels, respectively, of Figure 14. As can be seen, the RMS of the velocity field corresponding to the second mode pair is predicted quite accurately; the error in the maximum value is less than 10%.
In Figure 15, identical, but quantitatively the estimated peak value is about 55% of the DNS value. This is expected since only four modes were used for reconstruction. It is interesting to notice that although in the full tank these four modes account for ≈34% of the total kinetic energy ( Figure 4A), the peak value at the plane considered is captured with significantly better accuracy. We have tried to estimate higher modes using data from the same sensor points, but the reconstruction becomes more erratic and the fit values drop. This is due to the fact that higher order modes have smaller energy content spread over a broader range of frequencies (Figure 7), thereby representing random, rather that organized, motions. The contribution of the white noise

| Toward instantaneous 3D velocity field reconstruction
Combining the estimated coefficients a 1,e (t) and a 2,e (t) using one probe point, we now attempt to reconstruct the instantaneous 3D flow field in the whole vessel from u e x, t ð Þ≈ u x ð Þþa 1,e t ð Þϕ 1 x ð Þþa 2,e t ð Þϕ 2 x ð Þ. Note that the time-average together with these two modes account for more than 90% of the total kinetic energy of the flow, as shown in Figure 4B.
We first consider the 3D flow close to the impeller and compute the λ 2 value for the instantaneous true and estimated flow fields at 20 rotations into the validation dataset. Recall that λ 2 is the second largest eigenvalue of the S 2 + Ω 2 , where S and Ω are the strain and rotation rate tensors, respectively. 43 Isosurface plots of negative λ 2 (that visualize vortex structures) are shown in Figure 16. The flow close to the impeller is dominated by the large vortex structures captured by Modes 1 and 2 and, since the temporal coefficients are captured well (Figure 11A), the instantaneous flow is also rather well reconstructed. There are small-scale discrepancies close to the blade surface, but these are due to the higher modes.
Capturing the instantaneous flow away from the impeller is much more challenging. Contours of the true and estimated velocity magnitude on a vertical slice midway between the impeller blades at the same time instant are shown in Figure 17. Although at the height of the impeller the reconstruction is reasonably good, away from it the flow contains detailed features that are not reproduced. This makes physical sense. Recall that the Modes 1 and 2 are predominantly present around the impeller and have also a smooth but weak footprint in the rest of the vessel ( Figure 6A,B). For this reason, the estimated velocity field in Figure 17B looks smooth away from the impeller, while the real one contains random, small scale, features. Such features in the bulk of the vessel cannot be captured by sensor points located close to the impeller because they have a narrow spatial extent (coherence).

| Performance at off-design conditions
The robustness and accuracy of reconstruction was also tested at off-design conditions, specifically at Reynolds numbers 500 and 700, obtained by adjusting the rotational velocity of the blades.
The reconstruction process at the new conditions was the following: we performed DNS simulations for the aforementioned Reynolds numbers and recorded the values of velocity at the same monitoring points. We then computed the temporal coefficients a 1 (t) and a 2 (t) by projecting the fluctuating velocity data obtained by DNS to the POD modes computed at the reference conditions, that is, Re = 600. Finally, we employed the model described by (12), which was trained in the reference case, to reconstruct the temporal coefficients using the input data for the new Reynolds numbers. It is important to note that this approach has truly predictive power, because we have trained a model using reference data (at Re = 600) and then use it to estimate the flow at off-design conditions (Re = 500 and 700). The only input required is information from sensor points at the new operating conditions.
The results using a single probe point and a sixth-order model are shown in Figure 18. The fits for modes (

| CONCLUSIONS
In the present study, the POD technique was applied to reconstruct from limited data the large-scale vortical structures behind the blades of a Rushton impeller rotating inside an unbaffled stirred tank at the transitional Reynolds number of 600. About 34% of the turbulent kinetic energy was accounted for by two pairs of POD modes, with sharp spectral peaks at 1.5F N and 3.0F N , respectively. Higher order modes were found to contain less than 1% of the total energy and have broader spectra. The À11/9 power law, corresponding to homogeneous isotropic turbulence, was found to hold for modes with orders greater than ≈120.
The N4SID identification algorithm was then used to construct a data-driven estimator, that is, a linear, time-invariant system that reproduces the observed input-output relation between the velocity time signals at discrete sensor points and the temporal coefficients of each mode pair. Data from 1 to 6 sensor points were used. The points were placed at the location of maximum turbulent kinetic energy at the wake of each blade. For the first pair of modes, the quality of reconstruction was very good, even with a single probe point, but the order of the estimator was high, between 6 and 8. Increasing the number of probe points to 6 increased the fit and resulted in lower order and more robust estimators. Using a single input point to predict the second pair of POD modes was not successful, but a six-point input predicted well the frequency of the mode and yielded a satisfactory reproduction of the time trace of the temporal coefficients. The prediction of the RMS velocity for Modes 1 and 2 was excellent, while for Modes 3 and 4 very good (less than 10% error). The reconstruction of the entire 3D velocity field using the leading order modes was found to be reasonably good at the impeller height. However, further away from the impeller, where the footprint of the leading order modes was weaker, the detailed flow features were not reproduced.
The robustness of the estimator at off-design conditions, Re = 500 and 700, was then examined. The estimator constructed for the nominal operating condition, Re = 600, was employed to reconstruct the temporal coefficients at the new conditions. The time signals of the Modes 1 and 2 were again very well reproduced, with small performance degradation, even using data from a single probe point.
We have demonstrated that system identification is capable of reconstructing the instantaneous evolution of 3D large-scale structures inside a stirred vessel using information from very few data points. Many questions remain to be answered. What is the performance of the algorithm at higher Reynolds numbers? How can we represent stochastically the truncated modes? What is the quality of reconstruction of other fields, such as Reynolds stresses? How accurate is the estimation of scalar mixing time using the reconstructed velocity field? These are topics of future research.