Large-eddy estimate of the turbulent dissipation rate using PIV

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website. • The final author version and the galley proof are versions of the publication after peer review. • The final published version features the final layout of the paper including the volume, issue and page numbers.

many particle images. The particle motion and thus the velocity field are averaged over the extent of the interrogation window. For a measurement of small-scale turbulence quantities, the velocity field must be resolved accordingly. For example, the dissipation rate ǫ = ν i,j (∂u i /∂x j ) 2 , where ν is the kinematic viscosity, involves the sum of squared derivatives of the components u i of the velocity field. It requires the resolution of the velocity field down to the Kolmogorov scale η, where the velocity field is smooth and where derivatives can be estimated from finite differences. Normally, the linear dimension L of the interrogation window is much larger than η, so that the magnitude of derivatives is underestimated. The obvious cure is to make the interrogation window as small as η, at the expense of a very limited view of the velocity field (Tanaka and Eaton 2010). On the other hand, when the size of the interrogation window is well within the inertial range, a now popular approach involves the assumption that the statistics of the resolved velocity field ū is universal, with a universal relation between the measured derivatives of the spatially averaged field ū and the true dissipation rate, ǫ LE (Sheng et al. 2000), with S = ( i,j S ij S ij ) 1/2 , and S ij = 1 2 (∂ū i /∂x j + ∂ū j /∂x i ), so that in homogeneous turbulence �S 2 � = 1 2 i,j s i,j , with s i,j = �(∂ū i /∂x j ) 2 �. For interrogation windows with inertialrange dimensions, finite velocity differences u scale as �(�u) 2 � ∝ ǫ 2/3 L 2/3 , so that ǫ LE becomes independent of the window size. This approach is attractive as it provides a means of estimating the dissipation rate from PIV measurements that still provide information about the large-scale structure of the velocity field.
The commonly used value of the Smagorinsky constant is C Sm = 0.17 (Sheng et al. 2000;Lavoie et al. 2007), but (1) ǫ LE = 2 3/2 C 2 Sm L 2 S 3 , Abstract The result of a particle image velocimetry (PIV) measurement is a velocity field averaged over interrogation windows. This severely affects the measurement of small-scale turbulence quantities when the interrogation window size is much larger than the smallest lengthscale in turbulence, the Kolmogorov length. In particular, a direct measurement of the dissipation rate demands the measurement of gradients of the velocity field, which are now underestimated because the small-scale motion is not resolved. A popular procedure is to relate the statistical properties of the measured, but underresolved gradients to those of the true ones, invoking a large-eddy argument (Sheng et al. in Chem Eng Sci 55(20): [4423][4424][4425][4426][4427][4428][4429][4430][4431][4432][4433][4434]2000). We argue that the used proportionality constant, the Smagorinsky constant, should depend on the window overlap, on the used elements of the strain tensor, and on the way in which derivatives are approximated. Using an analytic description, PIV measurements of velocity fields from a kinematic simulation and experiments in a synthetic jetdriven turbulent flow with zero mean velocity, we propose new values for this constant.

Introduction
Particle image velocimetry is an excellent and widely used tool to measure the velocity field of a flow. It is based on the correlation of interrogation windows (or volumes) containing it must be realized that its numerical value should depend on the degree of window overlap, on the way in which derivatives are approximated using finite differences and on which components of the strain tensor are used in the estimate of S. In two-dimensional PIV, we can measure s 1,1 , s 1,2 , s 2,1 , and s 2,2 , but the derivatives of the third velocity component are unavailable. By using isotropy and incompressibility, the missing contributions to S 2 can be related to the measured ones: Only using the diagonal ones we have which we will indicate by "1D", or using both diagonal and off-diagonal ones which will be denoted "2D". Isotropy only applies to the ensemble-averaged squared gradients (∂u i /∂x j ) 2 , so that only S 2 can be estimated, and not S 3 , and the assumption �S 3 � = �S 2 � 3/2 is unavoidable (Meneveau and Lund 1997).
While the isotropy relations apply to the true isotropic velocity field, they no longer hold for the quantities of the resolved field computed with finite differences, so that s 1,2 � = 2 s 1,1 . Therefore, Eqs. 2 and 3 do not give the same answer for S 2 , which must be compensated with a different constant C Sm . So far, this subtlety was missed in the applications of the large-eddy correction to measurement of the dissipation rate. In fact, the commonly used value of the Smagorinsky constant C Sm = 0.17 is based on a box filter in Fourier space, contrary to the real space box filter that is associated with PIV.
For turbulence that is not homogeneous, the large-eddy correction should apply locally, much as in numerical large-eddy simulations. For turbulence that is not isotropic, it is not possible to guess the missing contributions to S 2 from 2D PIV, and a measurement of all strain components, for example through tomographic PIV, is unavoidable.
The measured derivatives are affected by noise. Even for particle images with no noise, out-of-plane motion of particles results in noisy PIV vectors. Because the dissipation rate ǫ LE involves the sum of squared derivatives, noise leads to an over-estimate of ǫ LE . More averaging, for example by using larger interrogation window sizes L, or using a smaller window overlap when estimating derivatives, reduces the influence of noise, but at the expense of a larger systematic error. Therefore, a measurement of ǫ LE involves a trade-off between statistical errors due to noisy PIV vectors and systematic errors due to averaging and finite differences.
The purpose of this paper is to precisely quantify the systematic errors of large-eddy PIV. In fact, we will provide (2) �S 2 � = 15 4 (s 1,1 + s 2,2 ), values for the constant C Sm depending on the window overlap, on the way of estimating derivatives, and on the used gradients (Eqs. 2, 3). Another way of estimating ǫ is through measurement of the second-and third-order structure function, where for the longitudinal structure function G L 2 (r), the component of the velocity points in the same direction as the separation vector, i = j, while for the transverse i � = j (of course, G T 3 = 0). The inertial-range behavior of G L 2 , G L 2 (r) = C 2 ǫ 2/3 r 2/3 with C 2 = 2.12 (Yeung and Zhou 1997), or that of G L 3 (r) = −(4/5)r ǫ, can in principle be used to infer ǫ from measured structure functions. However, also these structure functions are affected by the intrinsic averaging of PIV, while at the same time, the inertial range is limited. We will also quantify the effects of this systematic error for measurements of ǫ through Ḡ L 2 . Starting from a model energy spectrum, Lavoie et al. (2007) have studied analytically the effect of averaging on the measured dissipation rate ǫ LE , on the statistical properties of the velocity field and on the measurement of decaying turbulence. This spectrum is characterized by a turbulent velocity and a dissipation rate ǫ 0 and the effects of averaging are quantified by the ratio ǫ LE /ǫ 0 . Because the spectrum is a second-order quantity, a caveat is that this analysis only works for second-order quantities, such as S 2 and the second-order structure function. Inspired by this work, we follow a similar approach, but in addition study the influence of the interrogation window overlap, the apparent small-scale anisotropy, the choice of the measured strain components and the way to approximate derivatives using finite differences.
The spectral transfer of the turbulent energy has recently been measured by Ni et al. (2014) using a filtered scale approach. For filter length-scales within the inertial range, this would provide the dissipation rate. An extensive study of several ways to measure the dissipation rate from PIV, among which the large-eddy method, and a method based on G L 2 , was done by Jong et al. (2009). They carefully distinguish between random and systematic errors and observe that methods that reduce random errors have large systematic errors and underestimate ǫ. Their preference is the estimate ǫ G 2 through G L 2 , which they deem less sensitive to the effects of averaging. Taking their work further, we will argue here that the finite-difference estimate of gradients and the influence of averaging on G 2 , which are considered systematic errors in Jong et al. (2009), can be evaluated a priori.
Tomographic PIV has the great advantage of providing complete information about the velocity field (Elsinga et al. 2006). The effect of resolution on an estimate of the dissipation using tomographic PIV was discussed by Worth et al. (2010). Approximate agreement with the spectral filtering approach of Lavoie et al. (2007) was found, but it should be realized that the latter pertains to a 2D slice of the velocity field. Tokgoz et al. (2012) studied the turbulent flow between the two rotating cylinders of a Taylor-Couette device and compared the dissipation rate using tomographic PIV to a direct measurement through the torque needed to rotate the cylinders. They also found a large effect of window size and window overlap, but large-eddy corrections were not considered. Throughout, we will use the Pao spectrum for the isotropic three-dimensional E(k) (Pope 2000), with the spectral constant C 3 = 1.62 (Pearson et al. 2002), the modulation of the large scales , and the modulation of the small scales with β = 5.2, c L = 6.78, c η = 0.4, and η the Kolmogorov scale. If we write the Reynolds number as Re = (5/3) 1/2 u 2 ǫ −2/3 0 η −2/3 , we realize that the spectrum depends on the large scales through u, and the small scales through ǫ 0 . Further, Summarizing, the energy spectrum is determined by the dissipation rate ǫ 0 and the turbulent velocity u (or, equivalently, by ǫ 0 and Re ).
From E(k), we can compute the effects of averaging, window overlap and finite differences on the estimate ǫ LE of the dissipation rate through the large-eddy approach and the estimate ǫ G 2 through the second-order structure function. So, given a value of ǫ 0 , the Reynolds number, and information about the interrogation windows, we can compute ǫ LS and ǫ G 2 , and precisely chart the systematic errors.
In addition to analytical expressions, we test our ideas using a kinematic simulation of PIV image pairs, followed by a PIV procedure to generate vector fields and an analysis of these vector fields. Kinematic simulation is the procedure to produce a three-dimensional velocity field starting from a model spectrum, for which we choose the same E(k) (Eq. 4). The generated vector field is random, isotropic and incompressible, but it lacks true dynamics, while its statistics is Gaussian. However, due to the random character of particle seeding, the PIV particle images suffer from the same particle loss as experimental ones, resulting in noisy vector fields. Finally, we will illustrate our findings by analyzing the turbulent flow field in an experiment where turbulence with near-zero mean velocity is generated using randomly driven synthetic jets; an experimental setup similar to the one described by Jong et al. (2009). 2 Using the spectrum to compute the effect of averaging and finite differences In homogeneous isotropic turbulence, the spectral tensor involving the Fourier transform u i (k) of the velocity component u i is from which the longitudinal and transverse one-dimensional spectra follow as respectively. In terms of the one-dimensional spectrum, We will first discuss the effect of averaging of the xcomponent of the velocity field over a two-dimensional interrogation window of size L, In spectral space, averaging amounts to multiplication of the Fourier transform u(k x , k y ) with the product of the sinc functions, H(k x ) H(k y ), with so that the spectral tensor of the resolved field becomes Analogously to Eq. 5, the one-dimensional spectrum is now Moving to polar coordinates and using k y = k r sin φ = (k 2 − k 2 x ) 1/2 cos φ, we can write the integral over k y and k z as From this spectrum, we can compute the second-order longitudinal structure functions of the averaged velocity field, while the transverse structure function involves the transverse spectrum (Eq. 6) Using these structure functions, we can assess the smallscale isotropy of the flow by computing the ratio where the derivative d/dx is again approximated by a finite difference. Since the combination of averaging and differencing is no longer isotropic for finite sized interrogation windows, this ratio will differ from unity. The small-scale motion of isotropic turbulent flows measured with PIV will therefore appear anisotropic.

Derivatives through finite differences
We take the velocity field ū(x, y) averaged over a twodimensional interrogation window and compute the central difference x as an approximation of the longitudinal derivative ∂u/∂x, with 0 < α < 1 the overlap factor, such that for α = 0.5 windows overlap 50 %, and 75 % for α = 0.25. In spectral space, the central difference amounts to multiplication by The derivative is sensitive to noise in PIV vector fields; with ǫ LE the sum of squared derivatives, noise will increase the estimated dissipation rate. The "least squares" approximation of the derivative is more effective in suppressing noise than the central difference; it is based on a least squares fit of a line to 4 points, which has the spectral transfer function The estimate of the transverse derivative ∂u/∂y is Now, the estimate of the average square has to be done simultaneously with the averaging over y in Eq. 7, with the result For isotropic and incompressible turbulence, �(∂u/∂y) 2 � = 2 (∂u/∂x) 2 , but this no longer holds for the discretized derivative of averaged velocity fields, so that for finite window sizes L, 2 y � = 2 2 x . The relation between 2 y and 2 x may depend both on the size of the interrogation window L, on the overlap factor α and on the approximation to the derivative.

Kinematic simulations
A kinematic simulation is based on generating random Fourier modes with a prescribed spectrum (Kraichnan 1970;Fung et al. 1992). Each realization of the simulated velocity field consists of a sum of Fourier components: in which v n and w n are spatial Fourier amplitudes, which depend on k n , the discrete wavenumber vector. Exactly how they depend on k n determines the spectrum and all other properties of the resulting flow field. The velocity field is made incompressible, ∇ · u = 0, by making sure that the vectors v n and w n are perpendicular to the wavenumber vector. The easiest way to do this is by defining two new vectors a n and b n and taking v n and w n as the cross-products of these vectors with the normalized vector k n = k n /k n , so: v n = a n ×k n and w n = b n ×k n . The directions of vectors (12) k n are chosen to be uniformly random, and the lengths are chosen according to a geometrical distribution, where k 0 and k N k are the smallest and largest wavenumbers present in the simulation, and N k is the total number of wavenumbers.
The velocity spectrum E(k) determines how the lengths of a n and b n depend on k n , The integral length-scale L int sets the length-scale for the kinematic simulations. For the energy spectrum, we use the same as that in the analytic computation (Eq. 4).
To generate PIV images, we randomly sprinkle particles in the generated velocity fields and advect each of them during 0.1 Kolmogorov times τ η = η 2/3 ǫ −2/3 0 , using a fourth-order Runge-Kutta procedure with an integration time step τ η /40. The image of a particle in a light sheet is calculated by means of a two-dimensional Gaussian pointspread function centered on the particle position, with the peak intensity dependent on the particle's position within the light sheet. The light sheet is homogeneous in the plane of the image and has a Gaussian intensity profile in the zdirection (normal to the plane), I(z) = I 0 exp(−(z/w) 2 ), where w is the width of the profile and I 0 is the intensity in the center. The intensity in each pixel of the eventual image is obtained by taking the sum of the point-spread functions of all particles integrated over each pixel. Much as in the experiments of Sect. 4.4, the digital images have a depth of 12 bits and dimension 1600 × 1200 pixels.

Results
We will first discuss the results of the analytic model and those of the kinematic simulation. In both cases, the same three-dimensional spectrum E(k) is taken, with turbulence characteristics that are typical for the experiment, ǫ = 60 m 2 s −3 , η = 8.7 × 10 −5 m, and u = 1.9 m s −1 , so that Re = 460 and L int = 0.23 m.

Analytical
The influence of averaging on the spectrum and on the second-order structure function is shown in Fig. 1. The averaged spectrum shows the well-known sinc function modulation due to the spatial box averaging of the velocity field. Due to the relatively small Reynolds number, the second-order structure function Ḡ L 2 lacks a well-defined inertial range; it is even less pronounced in the structure function of the averaged velocity field. The model structure function never reaches to the theoretical prediction G L 2 (r) = C 2 ǫ 2/3 r 2/3 , with C 2 = 2.12 (Yeung and Zhou 1997). It illustrates that an error of about 30 % results if the dissipation rate is estimated from G L 2 , without allowing for the averaging of the velocity field, and without taking into account the finite Reynolds number. We will return to this issue in Sect. 4.2.
Starting from the large-eddy estimate Eq. 1, and the standard value of the constant C Sm = 0.17, we will now compute the influence of the window overlap and the method of estimating derivatives on ǫ LE . Practically, this was done by picking a value for ǫ 0 , and values for η and L, evaluating the integrals Eqs. 11, 13, and using Eqs. 1, 2 and 3 to compute ǫ LE . The result is shown in Fig. 2a, b for α = 0.5 and 0.25, respectively. Surprisingly, for the central  Fig. 1 a Influence of interrogation window averaging on the energy spectrum. Input parameters are u = 1.9 m/s, η = 8.7 × 10 −5 m, ǫ 0 = 60 m 2 s −3 , and L = 1.6 × 10 −3 m. b Influence of interrogation window averaging on the second-order structure function. The full line shows the structure function of the averaged velocity field, the dashed line that of the bare velocity field and the dash-dotted line the Kolmogorov prediction G 2 (r) = C 2 ǫ 2/3 r 2/3 , with C 2 = 2.12 differences approximation to the derivatives (Eq. 10) and half-overlapping windows (α = 0.5), ǫ LE ≈ ǫ 0 for window sizes L/η inside the inertial range. This result is surprising because the standard value of C Sm = 0.17 was computed assuming a box filter in Fourier space, and not in real space, while no allowance was made or window overlap or how to estimate derivatives.
Indeed, other values of α, the least squares approximation to the derivatives, and taking other components of the strain, result in large differences between the input and large-eddy dissipation rate. However, all results share the property that the correction is approximately independent of the interrogation window size for L inside the inertial range. This, of course, is the essence of the large-eddy method.
For window size L = 30 η the correction factors ǫ LE /ǫ 0 are listed in Table 1. After measuring ǫ LE using the standard value C Sm = 0.17, the results should be divided by the corresponding correction factor. It is perhaps more appropriate to define new effective values for the Smagorinsky constant C eff Sm , depending on the details of the PIV procedure; this has been done in the last column of Table 1.
When L is not inside the inertial range, the curves of Fig. 2 can still be used, albeit in an iterative fashion. When the Reynolds number is so small that an inertial range is completely absent, correction factors have to be recomputed using the procedures sketched in Sect. 2.1. However, in this case, the large-eddy method for measuring ǫ looses its charm.

Using the second-order structure function
An alternative to the large-eddy method would be to use the second-order structure function. Jong et al. (2009) propose to take the value of the compensated G 2 (r)/r 2/3 in its maximum at, say, r = r * and use it to estimate ǫ. This estimate needs the value of the Kolmogorov constant C 2 whose value (and possible dependence on the Reynolds number) is well established (Yeung and Zhou 1997). When r * is inside the inertial range, they argue that the second-order structure function should not depend strongly on the window size since L ≪ r * . It is now straightforward to test this assumption through computation of Ḡ 2 of the resolved field. The result is shown in Fig. 3. At the Reynolds number considered here, it appears that the dependence on the window size is much stronger than that of the large-eddy method, so that the second-order structure function approach would necessitate an iterative procedure. Table 1 Correction factors ǫ LE /ǫ 0 and effective Smagorinsky constants for different combinations of parameters The large-eddy estimate ǫ LE has been computed with the (traditional) Smagorinsky constant C Sm = 0.17; in combination with the correction factor, this defines the effective Smagorinsky constant C eff Sm = 0.17/(ǫ LE /ǫ 0 ) 1/2 in the last column. The computation is for L = 30 η, but, as Fig. 2 illustrates, there is a slight dependence on the window size due to a Reynolds number effect   Fig. 2 Large-eddy corrected ǫ LE as a function of interrogation window size L for two different values of the window overlap, α = 0.5 (a) and α = 0.25 (b). We have used the standard value of the Smagorinsky constant C Sm = 0.17. The thick gray lines indicate ǫ LE = ǫ 0 , with ǫ 0 the input dissipation rate. The full lines use the longitudinal gradients s 11 , s 22 (Eq. 2), the dashed lines represent the two-dimensional approximation to the strain (Eq. 3). The lines indicated by "CD" are computed using the central difference approximation to the derivative (Eq. 10), those marked by "LS" use Eq. 12. Almost no correction is needed for α = 0.5 in combination with the central difference approximation. c Apparent small-scale anisotropy s 21 /s 11 , computed for α = 0.5. Parameters used for the model are u = 1.9 m/s, η = 8.7 × 10 −5 m, ǫ = 60 m 2 s −3 Page 7 of 9 89

Kinematic simulations
All PIV image pairs are processed with the PivTec software (PivTec GMBH, Göttingen, Germany) with a different combination of interrogation window size (32 or 48 pixels), and window overlap factor α = 0.5 or α = 0.25. The resulting data are then used to compute dissipation rates four times, each time with a different combination of finitedifference scheme (central differences or least squares approach) and use or omission of off-diagonal derivatives (1D or 2D). In all cases, the Smagorinsky constant used is C Sm = 0.17. This yields the uncorrected dissipation rates, which are shown in Table 2. As can be seen they vary from approximately 30-90 m 2 s −3 . The corrected dissipation rates are shown in Table 2, below the uncorrected ones.
As can be seen, the corrected dissipation rates are all very close to the value used in the kinematic simulations: ǫ 0 = 60.1 m 2 s −3 . Furthermore, the values found using the central differences scheme are all slightly higher than those found using the least squares approach. This demonstrates that the least squares approach is indeed less prone to noise. Relatively large errors result in the case of large window overlaps (75 %, α = 0.25) when estimating gradients. The PIV vector fields of the kinetic simulations are noisy due to the finite number of randomly sprinkled particles and their out-of-plane motion. The combination of large window overlaps and noisy PIV vector fields is detrimental for the estimate of the dissipation rate.

Experiments
Approximately homogeneous isotropic turbulence with zero mean flow was created in a cubic box using an array of synthetic jets, in an arrangement inspired by the work of Hwang and Eaton (2004). The difference with their experiment is our usage of larger loudspeakers, which allowed us to reach larger Reynolds numbers.
The design of the turbulence chamber takes advantage of the property synthetic jets possess in which momentum transfer is possible without mass transfer occurring, in an average sense. The apparatus consists of a PVC cubic box with truncated corners, which were fitted with eight speaker-driven synthetic jets. Hwang and Eaton (2004) were able to reach a Re = 218; in our setup, we reach Re = 563. In this way, the inertial range is widened. The box, which is shown in Fig. 4, has a side length of 400 mm and speakers to 365 mm diameter (MTX Audio sub-woofer model RT15-04, Mitek Corporation, Phoenix, AZ, USA). Optical access was available on four of its six sides through Perspex windows.
The synthetic jets point toward the center of the box and are driven independently with random voltages generated by series of independent pseudo-random numbers. By employing a Gaussian filter, the driving voltages have Using the second-order structure function to estimate ǫ. a Compensated Ḡ L 2 (r)/r 2/3 with Ḡ L 2 of the averaged velocity field, computed using Eq. 8. The height h of the maximum (indicated by the gray line) provides the estimate ǫ G 2 = (h/C 2 ) 3/2 , with C 2 = 1.21 (Yeung and Zhou 1997). b ǫ G 2 /ǫ 0 as a function of the interrogation window size L/η. The used parameters are u = 1.9 m/s, η = 8.7 × 10 −5 m, ǫ 0 = 60 m 2 s −3 , and interrogation window size L = 1.6 × 10 −3 m Table 2 Dissipation rates computed from the kinematic simulations The uncorrected dissipation rates vary a lot, while most of the corrected dissipation rates are very close to the actual dissipation rate ǫ 0 = 60.1 m 2 s −3 1D/2D: omit or use off-diagonal components of the strain rate tensor. CD/LS: use central differences or least squares approach to approximate derivatives. The interrogation window size is L/η = 28, corresponding to 48 pixels with center frequency f 0 = 40 Hz and spectral width σ = 10 Hz. The eight instantaneous driving voltages add to 0, so that the pressure in the chamber stays approximately constant. Slight differences in speaker performance may result in a mean flow. To counteract this mean flow, speaker-specific coefficients were implemented in the driving algorithm, allowing fine tuning of the synthetic jets. A total of 1500 PIV image pairs were collected and processed using the PivTec software with a window size L = 1.3 × 10 −3 m and overlap factor α = 0.5. After correction with the appropriate C Sm we find for the central difference scheme ǫ LE = 61 m 2 s −3 and ǫ LE = 52 m 2 s −3 , for 1D and 2D gradients, respectively, while for the least squares approach ǫ LE = 54 m 2 s −3 and ǫ LE = 49 m 2 s −3 for 1D and 2D gradients, respectively. The least squares approach gives smaller values for ǫ after correction, which signifies the influence of noise on the measured vector fields. In addition, we find an apparent anisotropy s 2,1 /s 1,1 = 1.4, which should be compared to the value 1.6 of the analytical model. For ǫ = 54 m 2 s −3 we will now test the consistency between the measured and theoretical second-order structure function. The result is shown in Fig. 5a and indicates that the used value of ǫ might be a slight underestimate of the true one. Figure 5b shows the small-scale anisotropy ratio (Eq. 9). The finite-difference approximation of the derivatives in Eq. 9 results in an apparent anisotropy, as is demonstrated by comparing the experimental result to the one from the model. We approximated the derivative term in Eq. 9 with central differences, Figure 2 shows our main result and illustrates the advantage of the large-eddy PIV procedure to estimate the dissipation rate. At window sizes L/η 20, the correction factor or the effective Smagorinsky constant C eff Sm is almost independent of the window size. This, of course, is the essence of the large-eddy approach. For small Reynolds numbers, no sizable inertial range exists, and the result of large-eddy PIV now depends on the window size. This size dependence can be computed using the approach of this paper, but an iterative approach is needed for an estimate of the dissipation rate.

Conclusion
So far, the large-eddy estimate ǫ LE only involved the second-order quantity S 2 . We have studied the influence of averaging, finite differences and window overlap on ǫ LE . The restriction to S 2 is inevitable because in planar PIV, not all gradients of the velocity field are accessible, and only averaged squares of the missing gradients can be guessed on basis of isotropy and incompressibility. Therefore, the consistency of the large-eddy method with the second-order structure function is not surprising.
x y Fig. 4 Model of turbulence chamber. Approximately homogeneous and isotropic turbulence is generated using synthetic jets generated by large loudspeakers. For clarity, only one speaker is shown. The inner side length of the chamber is 40 cm. The jet orifices have a diameter of 4 cm Comparison of the measured second-order structure function to the one computed from the averaged velocity field. a Full lines (almost coincident) are measured Ḡ x,x 2 and Ḡ y,y 2 , the dashed line is the model structure function. b Isotropy relation Eq. 9 using Ḡ x,x 2 and Ḡ y,y 2 , the dashed line is for the analytical model. The parameters used for the model are u = 1.9 m/s, η = 8.7 × 10 −5 m, ǫ = 54 m 2 s −3 Page 9 of 9 89 We conclude that large-eddy PIV so far has not dealt with the true Smagorinsky model, but with a surrogate one, based on the second invariant of the strain tensor and assuming that �S 3 � = �S 2 � 3/2 . Because energy transport is associated with the true S 3 , it is expected that a large-eddy PIV method for the measurement of ǫ based on the third invariant is superior. Such a new method should be based on the experimentally accessible components of the strain tensor.
A point of concern is the effect of intermittency, which leads to a different scaling behavior of the second-order structure function. Instead of G 2 (r) ∝ r ζ 2 , with ζ 2 = 2/3 it is found that ζ 2 = 0.70 (She and Leveque 1994). This only leads to a weak window size dependence of the skewness, �S 3 �/�S 2 � 3/2 ∝ L −0.05 , and thus to a weak dependence of ǫ LE on the window size L.
Our analytic derivation and kinematic simulation are based on a model spectrum E(k), especially on its behavior at large k. It allows us to calculate the effect of averages, window overlap and finite differences on ǫ LE and G 2 . These results should not depend strongly on the precise form of E(k) at kη > 1.
We conclude that the large-eddy estimate of the dissipation rate using PIV provides a good value for ǫ, but that the used Smagorinsky constant should depend on the measurement conditions.