Divergence-Free Wavelet Analysis of Turbulent Flows - Springer Link
Journal of Scientific Computing, Vol. 17, Nos. 1–4, December 2002 (© 2002)
Divergence-Free Wavelet Analysis of Turbulent Flows Cem M. Albukrek, 1 Karsten Urban, 2 Dietmar Rempfer, 3 and John L. Lumley 1 Received August 27, 2001; accepted (in revised form) November 3, 2001 In this paper we study the application of divergence-free wavelet bases for the analysis of incompressible turbulent flows and perform several experiments. In particular, we analyze various nominally incompressible fields and study the influence of compressible perturbations due to experimental and computational errors. In addition, we investigate the multiscale structure of modes obtained from the Proper Orthogonal Decomposition (POD) method. Finally, we study the divergence-free wavelet compression of turbulent flow data and present results on the energy recovery. Moreover, we utilize wavelet decompositions to investigate the regularity of turbulent flow fields in certain non-classical function spaces, namely Besov spaces. In our experiments, we have observed significantly higher Besov regularity than Sobolev regularity, which indicates the potential for adaptive numerical simulations. KEY WORDS: Wavelets; turbulence.
1. INTRODUCTION The study, analysis and simulation of turbulent flows is probably one of the most challenging problems in engineering, scientific computing and applied mathematics. Wavelets offer some potential for both the analysis and the simulation of turbulent flows. This is mainly due to a good localization in both frequency and physical domain which entails a good separation of features associated with different characteristic length scales. Localization in physical space can be enhanced by employing compactly supported wavelets which all subsequent investigations will be confined to. A second important feature of wavelets is their compression property. By 1
Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, New York 14853. E-mail: {cma13;jll4}@cornell.edu 2 RWTH Aachen, Institut für Geometrie und Praktische Mathematik, Templergraben 55, 52056 Aachen, Germany. E-mail: [email protected] 3 Illinois Institute of Technology, 3110 S. State St., Chicago, Illinois 60616. E-mail: RempferD @iit.edu 49 0885-7474/02/1200-0049/0 © 2002 Plenum Publishing Corporation
50
Albukrek et al.
this we mean that often only relatively few but judiciously chosen terms in an infinite wavelet expansion suffice to recover the function within a desired tolerance. This relies on the so called cancellation properties of wavelets combined with the fact that certain weighted sequence norms of the expansion coefficients are equivalent to the norms for certain smoothness spaces such as Sobolev or Besov spaces (such estimates will be called norm equivalences). Such highly economical representations are therefore particularly promising for representing functions with a wide range of characteristic length scales such as those representing turbulent flow fields. The primary objective of this paper is to demonstrate the usefulness of wavelet concepts and their key mechanisms for the analysis and compression of turbulent flows. In particular, wavelet concepts are used to analyze the Proper Orthogonal Decomposition (POD) method [21], which aims at approximating the infinite-dimensional nonlinear dynamical system describing turbulent flows by some low-dimensional model. A variety of approaches for this reduction can be found in the literature. In e.g., [2, 16–18, 20, 23], wavelets already have been used for the analysis and simulation of turbulent flows. In most of the cited papers, the data analysis has been performed by means of orthogonal wavelets for periodic vorticity fields. For the simulation of periodic flows mainly exponentially decaying vaguelettes have been employed so far, see e.g., [17]. In [18] interpolatory wavelets have been used for the simulation. In this paper we propose to utilize compactly supported divergencefree wavelets [19, 24, 25] for the analysis of velocity fields. Since we can therefore work with the velocity/pressure formulation of the Navier– Stokes equations these tools work naturally also for three-dimensional problems. Since the incompressibility constraint is satisfied by the basis functions, such systems are natural candidates for a quantitative analysis of the role of incompressibility and its perturbation. Specifically, we will be concerned with the following issues. First, divergence-free wavelets and a wavelet basis for an appropriate complement space give rise to a Helmholtz-type decomposition for spaces of vector fields [25]. Using the divergence-free wavelet transform, we can split flow data obtained from experiments as well as direct numerical simulations into its solenoidal and compressible parts. In particular, this tool offers insight into the stability of the POD method, which is a well-known technique for constructing a low-dimensional turbulence model. POD modes are given as the eigenfunctions of the autocorrelation tensor for an ensemble of flow realizations. They can be interpreted as the coherent structures in a turbulent flow. In order to produce an incompressible flow model based on the POD-Galerkin method, the corresponding data fields have to be divergence-free. However, it turns out that all the considered flows do in fact contain some compressible components. By comparing the POD modes obtained from the available data and those from its divergencefree part we study the quantitative influence of this violation on the incompressible flow model.
Divergence-Free Wavelet Analysis of Turbulent Flows
51
Second, the above mentioned ability of wavelet expansions to separate data into contributions from different characteristic length scales will be used to provide insight into the multiscale structure of various flows and, in particular, of the coherent structures. Finally, we exploit the compression properties of divergence-free wavelet bases. Here the term compression refers to the reduction of data sets and should not be confused with compressibility as a flow property. We compress several turbulent flow fields in the sense that we keep possibly few terms in their wavelet expansion while retaining most of the energy. Such a reduction without loosing essential information is closely related to the regularity of the underlying field in a certain scale of Besov spaces [14]. Due to the above mentioned norm equivalences we are able to quantify the regularity of given data sets in the relevant (Besov) scale and draw conclusions about their adaptive recovery, see [5]. This latter aspect is particularly interesting with regard to adaptive direct numerical simulations. This paper is organized as follows: In Section 2, we review the relevant wavelet theory that is needed in this paper and Section 3 is devoted to a brief review of the POD method. Our particular implementation is described in Section 4. Sections 5 and 6 contain our numerical results. We will frequently use the following notation throughout the remainder of this paper. A ’ B means that there exists constants C \ c > 0 such that cA [ B [ CA and the constants are independent of all the parameters A and B may depend on.
2. WAVELETS AND SIGNAL ANALYSIS In this section, we review the relevant properties of wavelet bases that will be needed in the sequel.
2.1. Biorthogonal Wavelets Let us denote by t, t˜ ¥ L2 (R) a pair of functions that are refinable, i.e., there exist some sequences a :=(ak )k ¥ Z , a˜ :=(a˜k )k ¥ Z ¥ a2 (Z) such that t(x) := C ak t(2x − k),
t˜(x) := C a˜k t˜(2x − k)
k¥Z
(1)
k¥Z
for x ¥ R. Moreover, t and t˜ are dual or biorthogonal, i.e., (t( · ), t˜( · − k))0, R =d0, k , k ¥ Z, where ( · , · )0, R refers to the usual L2 (R)-inner product. We will call t and t˜ dual scaling functions. Defining for j, k ¥ Z the scaled and translated version of a function f: R Q R by fj, k (x) :=2 j/2f(2 jx − k),
x¥R
(2)
52
Albukrek et al.
the pair of dual scaling functions generate a sequence of nested spaces ˜j :=closL (R) span{t˜j, k : k ¥ Z}, defined by Vj :=closL2 (R) span{tj, k : k ¥ Z}, V 2 which are referred to as biorthogonal multiresolution analysis, [22]. ˜ j are defined by Vj+1 =Vj À Wj and The complement spaces Wj , W ˜ ˜ ˜ ˜ j and V ˜j + Wj . Finally, g, g˜ ¥ L2 (R) Vj+1 =Vj À Wj , j ¥ Z, as well as Vj + W ˜ j= are called biorthogonal wavelets if Wj =closL2 (R) span{gj, k : k ¥ Z}, W closL2 (R) span{g˜j, k : k ¥ Z} and the systems {gj, k : j, k ¥ Z}, {g˜j, k : j, k ¥ Z} form Riesz bases for L2 (R). Biorthogonal wavelets are defined by the relations g(x) := C bk t(2x − k), k¥Z
g˜(x) := C b˜k t˜(2x − k)
(3)
k¥Z
where bk :=(−1) k a˜1 − k , b˜k :=(−1) k a1 − k , see (1). Biorthogonal wavelets can be constructed from different families of scaling functions. One prominent example is formed by using cardinal B-splines of order d with nodes at the integers as primal scaling functions, which will be denoted by d t. In [6], a whole family d, d˜ t˜ for d+d˜ even, of dual scaling functions for cardinal B-splines was constructed. With increasing d˜, the polynomial exactness of the system induced by d, d˜ t˜ grows at the expense of increasing size of their supports. The corresponding biorthogonal wavelets will be denoted by d, d˜ g and d, d˜ g˜, respectively. However, in many places we will drop these indices for notational convenience, except when their specific values matter. Nowadays biorthogonal wavelet bases are also available in L2 ([0, 1]), [10], in L2 (W) for a wide class of bounded domains W … R n and even manifolds, see e.g., [3, 4, 11–13]. We will denote scaling and wavelet systems in L2 (W) by Xj :={tj, k : k ¥ Ij }, Uj :={gj, k : k ¥ Jj }, respectively, where Ij and Jj are finite index sets. The spaces Vj and Wj are defined analogously. Finally, the full wavelet basis will be denoted by U :={gl : l ¥ J}, where l=(j, k) and J=1j Jj (and J=Z 2 if W=R). 2.2. Wavelet Transform Given a function vj+1 =; k ¥ Ij+1 cj+1, k tj+1, k ¥ Vj+1 with single scale coefficients cj+1, k , the refinement equations (1), (3) can be used to split vj+1 into a coarse part (average) vj ¥ Vj with coefficients cj, k and a detail part wj expanded in terms of the wavelets gj, k with coefficients dj, k . In fact, these new coefficients are obtained by the recursion cj, k = C m ¥ Ij+1
2 −1/2a˜m − 2k cj+1, m ,
dj, k = C
2 −1/2b˜m − 2k cj+1, m
(4)
m ¥ Ij+1
Repeating this process results in a transformation that is often called decomposition or discrete wavelet transform (DWT) algorithm. Conversely, given the set of coarse scaling coefficients, c0, k , k ¥ I0 , and detail coefficients, dj, k , k ¥ Jj , j \ 0, the refinement equations (1), (3) can be used to
Divergence-Free Wavelet Analysis of Turbulent Flows
53
recover the scaling function coefficients recursively on the highest level through cj+1, k = C 2 −1/2ak − 2m cj, m + C 2 −1/2bk − 2m dj, m , m ¥ Ij
k ¥ Ij+1
(5)
m ¥ Jj
This process reconstructs a signal given by its coarse scale parts and the corresponding details. It is therefore also called reconstruction, or inverse discrete wavelet transform (IDWT). 2.3. Divergence-Free Wavelets Now, we briefly review the construction of divergence-free wavelets as introduced in [19, 24, 25]. The key point is that even for 1D biorthogonal wavelet systems there is a simple mechanism that relates different wavelet bases by differentiation and integration. Then, by fixing one component of a vector wavelet, one can arrange the remaining ones in order to derive divergence-free wavelets. From now on, we concentrate on biorthogonal B-spline wavelets only, i.e., the scaling functions t=td are cardinal B-splines of order d. It turns out that they are quite suitable for the construction of divergence-free wavelets due to certain relations involving their derivatives and primitives, [19]. To be precise, the following relations hold: d
tŒ(x)=d − 1 t(x+m(d − 1)) − d − 1 t(x − m(d − 1))
(6)
for the scaling functions, where m(d) :=d mod 2 as well as d, d˜
gŒ(x)=4 d − 1, d˜+1 g(x)
(7)
for the wavelets. We will frequently use the following short hand notations for the involved functions: t (1) :=d t, and g (1) :=d, d˜ g, as well as t (0) :=d − 1 t, and g (0) :=d − 1, d˜+1 g. A formula corresponding to (7) also holds for wavelet systems on the interval, see [25]. In fact, viewing Uj as a column vector consisting of the (0) wavelets gj, k on [0, 1], we obtain dxd U (1) j =Dj U j , where Dj is invertible. For the scaling function systems Xj , we obtain a similar relation, namely d (1) (0) dx X j =Cj X j , where the matrices Cj have full rank (but may in general be rectangular). As a consequence, we obtain
1 dxd X 2=span(X d span 1 U 2 =span(U dx span
(1) j
(0) j
(1) j
(0) j
) (8) )
54
Albukrek et al.
where span always denotes the closure of the linear span. Now, we use a tensor product approach (here, for simplicity only in 3D) and obtain wavelet systems on [0, 1] 3 by Y (ij, 1e, i2 , i3 ) :=G j,(i1e)1 é G j,(i2e)2 é G j,(i3e)3
(9)
e ¥ E g :={0, 1} 3 0 {(0, 0, 0) T}, i1 , i2 , i3 ¥ {0, 1} and
˛ XU
G (ij, ne)n :=
(in ) j (in ) j
if en =0 if en =1
This results in a diagram depicted in Fig. 1. In particular, we obtain also wavelet bases in the spaces H(div; W) and H(curl; W) defined by H(div; W) :={f ¥ L2 (W) : div f ¥ L2 (W)} and similarly for H(curl; W). Moreover, for the space H(div; W) it is possible to identify and characterize the subspace V(div; W) :={f ¥ H(div; W) : div f=0} of divergence-free functions. This means that one can construct functions k df l such that df Y df :={k df l : l¥J }
(J df being an appropriate index set) is a wavelet basis for V(div; W) and that each vector field fj ¥ Wj :=span Y div j … H(div; W) can be decomposed into its divergence-free part and some complement in the following way, dc df dc fj =Q df j fj +Q j fj , where Q j , Q j denote the biorthogonal projection onto df dc dc W df j :=span Y j and W j :=span Y j , respectively. The latter is a stable (but dc not necessarily orthogonal) complement in the sense that Wj =W df j À Wj .
Fig. 1. Construction principle for wavelet bases in H(div; W) and H(curl; W), where S(Y) := span Y. The solid lines indicate that the images of the spaces on the left of the arrow under the corresponding differential operator coincide with the spaces generated by the functions on the right of the arrow. Dashed lines indicate that the functions form a wavelet system for the corresponding function space.
Divergence-Free Wavelet Analysis of Turbulent Flows
55
Let us finally mention that the corresponding multidimensional filters can easily be derived from those of the univariate functions (see e.g., [6]) by building tensor products as in (9) and performing the modifications described in (6) and (7). 2.4. Wavelet Thresholding A key property of wavelet bases is their ability to compress signals in the following sense. Once a wavelet representation of some data is given, one can obtain an economic representation by discarding all wavelet coefficients below some (possibly scale dependent) threshold or equivalently by keeping a certain number of (possibly weighted) largest coefficients. The main analytical tool for quantifying such concepts is a family of norm equivalences. In fact, under certain assumptions on the primal and dual functions, one can show, [9], that for some c, ˜c > 0 depending on the Sobolev regularity of t, t˜, one has
>C d g> l
l¥J
2
’ C 2 2s |l| |dl | 2,
l
s, W
s ¥ (−c˜, c)
(10)
l¥J
where l=(j, k), |l| :=j and H s(W) denotes the Sobolev space (of order s ¥ R, normed by || · ||s, W ). Similar relations hold for divergence-free functions in Sobolev spaces when the above described divergence-free wavelets are used. In fact, it has been shown in [19, 25] that
>
C dl k df l
l ¥ J df
>
2
’ C rs, |l| |dl | 2, s, W
l ¥ J df
˛ 1,1+2
rs, j =
2j+2
,
if s=0 if s=1 (11)
Once relations like (10) or (11) are available, the error caused by discarding all coefficients below some threshold can be estimated in any of the above norms. In the present context, we apply the following version of thresholding based on measuring accuracy in H s. Rather than discarding all coefficients below some threshold we keep as few coefficients as possible so that the norm of the compressed expansion is still o times the original norm, where o ¥ (0, 1) is a fixed preservation rate. Thus one looks for the smallest set L(o) … J df such that ; l ¥ L(o) rs, |l| |dl | 2 \ o ; l ¥ J df rs, |l| |dl | 2. In view of (11) this can in fact be done by first reordering the vector v :=(vl )l ¥ J df , vl :=`rs, |l| |dl | in decreasing order of the modulus of its components, then summing in that order the squares of its entries until the sum reaches o times the full sum on the right-hand side of (11) which we denote by ||v|| 2(s) for a moment. Clearly, if N(o) is the number of summands obtained in this way, the set L(o) consists of the N(o) largest weighted wavelet coefficients. The error by which the compressed representation ; l ¥ L(o) dl k df l deviates
56
Albukrek et al.
from the original field in the norm || · ||s, W , is proportional to `1 − o ||v||(s) and obviously also to the error of the best N(o)-term approximation, i.e., the best approximation that can be obtained by any linear combination of N(o) basis functions, [14]. 3. PROPER ORTHOGONAL DECOMPOSITION In general, turbulent shear flows contain discernible, self-sustaining components, called the coherent structures, which account for a major part of the turbulent kinetic energy. Hence, it is natural to look for a good approximation to these structures in order to use them as a basis for a flow simulation. One prominent example of such a basis can be constructed by the POD method. We briefly outline here the main ingredients that will be needed in the sequel. Let us consider a specific flow scenario in a domain W … R n for a time interval [0, T]. We assume that we are given samples of velocity fields uj, k :=u(xk , tj ) ¥ R n, which can be seen as samples on a grid Wh = ( j − 1) T {x1 ,..., xM } … W of ‘‘snapshots’’ of the flow at a time tj := N − 1 ¥ [0, T], j=1,..., N. In order to determine a basis which is in a certain sense optimal in representing the ensemble U=(uj, k )j=1,..., N; k=1,..., M , one can determine f=(f1 ,..., fM ) T ¥ R M × n by
3 T1 C 1 C u z 2 ;||z|| 4 N
f :=arg max z
M
2
j, k
j=1
k
2 2
(12)
k=1
The summation over the snapshots can be interpreted as a time averaging. Note that (12) can be rewritten as
3 zz Rzz 4 , T
f :=arg max z
T
1 N where R := C uj u Tj T j=1
and uj =(uj, 1 ,..., uj, M ) T. This implies that f is the eigenvector of R corresponding to the largest eigenvalue, which coincides with the square of the 1 largest singular value of `T U. The vector f determined by (12) is the best approximation to the ensemble by a single vector. One obtains the desired basis {f i : i=1,..., J [ N} by taking the eigenvectors of R corresponding to the J largest eigenvalues, i.e., the principle axes of the cloud of data points. These vectors are also called the POD modes. They can be interpreted as capturing the coherent structures of the flow. These (finitely many) modes are often subsequently used as (discretizations of ) test and trial functions in a spatial Galerkin discretization of the Navier–Stokes equations which results in a system of ordinary differential equations. Thus one obtains a very economical approximate solution that can be used to further study flow characteristics for the given environment for comparisons with other simulations.
Divergence-Free Wavelet Analysis of Turbulent Flows
57
4. IMPLEMENTATION AND VALIDATION Implementation. We implemented two-dimensional and three-dimensional versions of the divergence-free and the complement wavelet transform for realizing the compression and analysis purposes indicated above. For visualization purposes, only the two-dimensional results will be presented. The underlying code was written in Matlab. In addition, C++ routines from the MultiScale Library (MSLib), [1] were used. To study the possible range of data reduction (wavelet compression) it is important to test different filters that give rise to different orders of cancellation properties and ranges for the norm equivalences. Therefore the test data has to be brought into a filter independent format. We use the following approach. First, we compute the piecewise linear interpolant for the given discrete sample values. Of course, the resulting function has an infinite expansion in terms of the divergence-free wavelets. Therefore, in the next step, we compute the biorthogonal projection onto a multiresolution space of sufficiently high level exceeding the sample resolution. This level is chosen so that the maximum of the relative errors at the grid points stays below 0.2%. Validation. We have used two kinds of data for our tests. First, as a basis for comparison, we used data for a periodic flow in a box which has kindly been provided to us by K. Schneider, see [17]. We have also validated our code by comparing our results with the data reported in [17]. Then, we used several snapshots of a two-dimensional velocity field obtained from the Particle Imaging Velocimetry (PIV) measurements of the turbulent wake behind an airfoil with the parameters displayed in Table I. A slice of a snapshot is shown in Fig. 2. This data has been kindly made accessible to us by M. Glauser and C. S. Yao from NASA Langley Research Center. Of course, only for low Mach numbers, as in the present case, air flow can be assumed to be incompressible. Moreover, when taking two dimensional slices one cannot expect to obtain strictly divergence-free fields in the turbulent flow regime, hence this data is particularly suitable for the analysis of compressible components in a given flow.
Table I. Parameters of Our PIV Test Data Experimental setup and grid specifications Airfoil Wind tunnel Flow No. of grid points
Chord length=101.6 mm, thickness=12.7 mm, angle of attack=4.5° Height=108 mm, width=254 mm Free stream velocity=18.1 m/s, mach number=0.05, Reynolds number=135.000 x: 168; y: 94; grid spacing: 0.471 mm
58
Albukrek et al.
Fig. 2. Slice of the PIV data.
5. DATA ANALYSIS In this section, we summarize our numerical results obtained from the application of wavelet decompositions to several data.
5.1. Artificial Compressibility In order to investigate the data, we decompose flow fields in order to determine their incompressible and compressible parts. As the PIV data consists of 2D streamwise slices of the corresponding 3D incompressible turbulent flow field, it is not surprising that the 2D vector fields are not exactly divergence-free. For given flow data, we compute the piecewise linear interpolant fj for some high level j as described above. Then, we dc compute Q df j fj and Q j fj , namely the incompressible part and the complement. By using the norm equivalences (10), (11), we compute the corresponding norms. For 240 different realizations of the flow field, the incompressible part accounts for 93% of the total energy on average. The next issue is the distribution of this energy over the scales of the flow. In Table II, we display the L2 -norms of the wavelet projection of two different snapshots onto different levels. These norms can be interpreted as the energies contained on particular levels of the flow. It can be observed that the loss of energy due to removing the complement part occurs mostly on the finer scales of the flow, for which the cross-stream turbulent flow fluctuations are important. The divergence-free transform hence also offers
Divergence-Free Wavelet Analysis of Turbulent Flows
59
Table II. Energies of the Original Field and the Divergence-Free Part Snapshot 15
Snapshot 140
j
orig.
div-free
perc.
orig.
div-free
perc.
4 3 2 1 0
0.0070 0.0148 0.0238 0.0383 0.0372
0.0063 0.0133 0.0222 0.0363 0.0347
90.4% 89.6% 93.1% 94.6% 93.3%
0.0070 0.0134 0.0228 0.0343 0.0418
0.0063 0.0118 0.0212 0.0316 0.0400
90.4% 87.9% 92.9% 92.4% 95.8%
natural means to quantify the associated experimental measurement error over the scales of the flow. Next, we use the DWT in order to investigate the multiscale structure of the given data. As in image processing, we consider plots of wavelet coefficients over the different levels. This offers insight into the distribution of the flow over the scales. However, as shown in Fig. 3, right-hand side, the contributions on finer scales decay so rapidly that they are no longer visible in the plot. This indicates that the main part of the flow is concentrated on the coarser scales. To show the structure of the high scale contributions the coefficients are amplified to attain essentially the same order of magnitude as the coarse scale contributions. In this way an approximate Normalized Divergence free Wavelet Coefficients for Levels 1 to 8
Divergence free Wavelet Coefficients for Levels 3 to 8
10 50
20 100 30
150 40
50
200
60 250 50
100
150
200
250
10
20
30
40
50
60
Fig. 3. Divergence-free wavelet decomposition of the data f from [17], i.e., f= ; j, e, k dj, e, k k df j, e, k (see (9)). The sets of coefficients dj, e :=(dj, e, k )k are shown for j=−1,..., 4 (starting from j=−1 in the top left corner, larger rectangles indicate larger j). On each level the contributions for e=(0, 1) (lower left corner), e=(1, 1) (diagonal) and e=(1, 0) (upper right corner) are shown separately.
60
Albukrek et al.
self-similarity is visible, see Fig. 3, left. Roughly speaking this corresponds to a measurement in H 1 rather than in L2 . 5.2. Analysis of the POD Modes As already mentioned, the POD modes may be interpreted as the coherent structures, i.e., those structures which capture the ‘‘main features’’ of the flow. The divergence-free wavelet analysis offers a tool to gain insight into at least two interesting issues concerning POD: Influence of Artificial Compressibility on the POD Modes. In many cases, extensive analysis on a set of turbulent flow data is carried out before an appropriate turbulence model can be developed. The kind of flow configuration that we investigate here is usually modeled by the incompressible Navier–Stokes equations. The POD modes extracted from data are used as test and trial functions in a Galerkin scheme for constructing a low-order numerical approximation of the full set of the Navier–Stokes equations. As an outcome, this imposes the divergence-free requirement on the POD modes. For the POD functions to be incompressible, the underlying flow data also has to meet the divergence-free condition. However, as already mentioned, data obtained from experiments or direct numerical simulations may contain errors. These errors may cause the flow field under consideration to be no longer divergence-free. The use of POD modes extracted from such data as a divergence-free basis in an incompressible flow simulation constitutes a violation on the used model. The proposed Helmholtz decomposition based on divergence-free wavelets offers a tool for quantifying the associated error. As two-dimensional slices of the PIV data include a naturally built-in artificial compressibility, we chose to use this data for our analysis. To this end, we first created a new set of divergence-free fields which were obtained by the divergence-free projection of the data. Then we applied POD decomposition on both sets, namely the original and the divergencefree part. In the first step, we compare the turbulent kinetic energies associated with the decomposed fields using the POD coefficients of the two data sets. From the stability of the POD method, it is expected that the energies of the original and divergence-free field are asymptotically equal. Here, we are interested in quantitative results. We find exactly the same percentage of energy recovery (93%) as mentioned above when we compare the norms of the POD coefficients corresponding to the original and filtered (i.e., divergence-free) data sets. The top plot in Fig. 4 shows the turbulent kinetic energy captured by the POD modes of the original data and its divergence-free part in a cumulative sense. The lower plot shows the percentage of the total kinetic turbulent energy captured by each individual mode. In both plots, the lines for the original data and its divergence-free part are quite close indicating the high quantitative stability of the POD with respect to perturbations.
Divergence-Free Wavelet Analysis of Turbulent Flows
61
Table III. Distribution of L2 -Norm Energy over the Scales of the POD Eigenfunctions Corresponding to the Original PIV Data Set j
mode 1
mode 25
mode 50
mode 100
mode 200
mode 240
9 8 7 6 5 4 3 2 1 0
0.0628 0.0885 0.1703 0.5583 0.6961 0.1267 0.1058 0.0790 0.0461 0.0142
0.0956 0.1959 0.3172 0.6266 0.5527 0.3346 0.1661 0.0143 0.0044 0.0006
0.1186 0.2527 0.4929 0.6135 0.4307 0.3167 0.0889 0.0361 0.0064 0.0182
0.1596 0.3755 0.6116 0.5316 0.4292 0.1627 0.0924 0.0458 0.0038 0.0011
0.2439 0.5450 0.6323 0.4422 0.2710 0.1339 0.0839 0.0133 0.0043 0.0011
0.8838 0.4662 0.2340 0.1264 0.0637 0.0216 0.0090 0.0022 0.0004 0.0001
Multiscale Structure of POD Modes. There has been an ongoing debate about whether the coherent structures have multiscale structure or not. A wavelet-based multiscale analysis of the POD modes interpreted as the coherent structures can shed some light on this question. For this purpose, we have performed wavelet transforms on several POD modes of the PIV data. For each mode starting from the most to the least energetic one we have computed the energy distribution over the different scales of the flow. The results were compiled in Table III. In all columns the maximal value is printed in boldface. The results indicate that the energy for the first (i.e., most energetic) few modes is concentrated around the 5th level of the decomposition. For the less energetic modes, the energy is distributed more evenly between the finer scales of the flow. Assuming that the coherent structures are represented by the first 25 POD eigenfunctions which capture Turbulent Kinetic Energy Distribution over Number of Modes
Percentage of Total Energy
100 80 PIV Data Incompressible Part
60 40 20 0
20
40
60
80
100 120 140 160 Number of POD Eigenfunctions
180
200
220
240
Turbulent Kinetic Energy Distribution over Modes
Percentage of Total Energy
7 6 5 4 3 2 1 0
0
5
10
15
20 25 30 POD Eigenfunction
35
40
45
50
Fig. 4. Kinetic turbulent energy (top) and percentage of total energy (bottom) for the original PIV data (straight lines) and the divergence-free part (dashed lines).
62
Albukrek et al.
about 60% of the turbulent kinetic energy, it can be seen that about 6 scales contribute most to these structures. Hence we see that at least for our data, the coherent structures have multiscale structure. Let us mention that we obtain similar results also for the data in [17] as well as for the divergence-free part of the given data. A final remark can be made from the above study: Our analysis reveals that the deviation from incompressibility is essentially limited to the finer scales of the flow, while the coherent structures are associated with the coarser scales of the flow. Hence the compressibility error induced by the construction of a low dimensional, incompressible flow approximation based on POD is by far less significant due to the use of only the first few dominant modes that represent the coherent structures. 6. COMPLEXITY OF TURBULENT FLOWS In this section we use the strong analytical properties of wavelets in order to investigate the structure of a turbulent flow. We are interested in two issues. First we study how sparse wavelet expansions of typical flow instances are, i.e., how many degrees of freedom are needed to obtain a certain energy recovery. Second, we investigate how effective an adaptive wavelet method can be at best for simulating incompressible turbulent flows. 6.1. Data Compression
100
100
90
90
80
80
Percentage of Energy Recovery
Percentage of Energy Recovery
In order to investigate the quantitative aspects of the compression properties of divergence-free wavelets, we decompose several snapshots of the given flow ensemble. The compression rate is the relation between the number of coefficients kept and the amount of energy retained in the field. Figure 5 is a typical example of the compression performance of divergencefree wavelets on the snapshots of the PIV data. 99% of the energy in the L2 -norm can be recovered by dropping 97% of the wavelet coefficients. On the other hand, the same recovery in the H 1-norm requires keeping 10% of the coefficients. Let us mention that the compression rates for the
70 60 50 40 30 20
0
0
20 40 60 80 Compression Percentage
60 50 40 30 20
H1 norm L2 norm
10
70
10
100
0 500 1000 1500 2000 2500 3000 Number of non zero Coefficients (out of 31764)
Fig. 5. Compression percentage (left) and number of retained coefficients (right) over recovered energy in L2 (solid line) and H 1 (dashed line).
Divergence-Free Wavelet Analysis of Turbulent Flows
63
periodic data are even higher and are comparable to those reported in [17]. The results displayed in Figs. 3 and 5 and have been computed with the basis corresponding to d=d˜=3 but we observe similar curves also for other choices of these parameters. This is related to the smoothness of our data which will be investigated in the next subsection. 6.2. Potential for Adaptivity A central objective for further research is to use divergence-free wavelets for a fully adaptive numerical simulation of turbulent flows. The corresponding spatial Galerkin projection of the Navier–Stokes equation could be done with respect to trial spaces spanned by adaptively chosen divergence-free wavelet bases. Quite recently, the connection between adaptive methods and nonlinear approximation has been analyzed rigorously for a wide class of problems, [5, 7]. In particular, it was shown under which circumstances it can be proven that the use of adaptive methods offers an asymptotic gain of efficiency over uniform refinement strategies. This means that the rate of convergence (relating accuracy and degrees of freedom) is higher for the adaptive than for the non-adaptive method. This is closely related to the regularity of the flow fields in certain smoothness scales which will be explained next. In the context of incompressible flows smoothness is usually measured in the scale of Sobolev spaces H s, i.e., derivatives are measured in L2 . The Besov spaces B sq (Lp ) can be used to measure smoothness of order s ¥ R in Lp (where the additional parameter q ¥ (0, .) permits some fine-tuning). Besov spaces arise as interpolation spaces between Sobolev spaces, see e.g., [14]. The key here is that they can also be characterized by weighted sequence norms based on wavelet expansions similar to (10), [5, 14]. Now fixing s and decreasing p the spaces B sq (Lp ) become larger and admit an increasingly stronger singular behavior of its elements. The Sobolev embedding theorem tells us how far p can be decreased so that the spaces remain continuously embedded in L2 . The critical value p=y is given by the relation 1 s 1 = + y n 2
(13)
where n denotes the spatial dimension. Figure 6 shows a diagram of these spaces. The horizontal axis corresponds to 1/p, i.e., the measure of integrability and the vertical axis shows the smoothness, s. The vertical line through (12 , 0) corresponds to the Sobolev spaces where the smoothness is measured in L2 . The diagonal line shows the Besov spaces with the relation of s and y given by (13). Hence, the space B sy (Ly ) is (much) larger than H s when s, y are related by (13). Moreover, this gap increases when s grows and hence y decreases. The point is that Besov regularity according to the critical embedding line is (almost) characterized by nonlinear approximation and, in particular, by the best wavelet N-term approximation.
64
Albukrek et al.
Fig. 6. Diagram indicating smoothness spaces.
Roughly the error of best N-term approximation of u decays like N −s/n when u ¥ B sy (Ly ) for s and y related by (13). Thus high regularity in this scale permits high accuracy at the expense of relatively few degrees of freedom. Knowing about the Besov regularity of a flow would allow us to predict the size of the sets L(o) in the above thresholding scheme. The adaptive schemes in [5, 7] was shown to perform in a certain range asymptotically as well as the best N-term approximation. Thus, the more the norm of the approximated object in H s exceeds its norm in B sy (Ly ), the higher the efficiency gain offered by the adaptive method is. In particular, when the approximated function belongs to B sy (Ly ) but not to H s the gain is even reflected by an asymptotically higher convergence rate. Therefore we are interested in the norms of the flows in the spaces B sy (Ly ) for growing s. The fact that these norms can indeed be computed relies on the above norm equivalences. In fact, it is well-known that when s and y are related by (13) one has [14] (see also (10))
> C dg> l l
l¥J
y s
B y (Ly (W))
’ C |dl | y
(14)
l¥J
This relation is valid for a range of s that depends on the underlying wavelet basis, especially its regularity. In order to use (14) to measure the Besov regularity of flow data recall that the data is discrete and therefore has to be mapped first into a function corresponding to an element of a multiresolution space for a very high level. This function automatically belongs (for any set of coefficients) to a given Besov space as long as the underlying wavelets belong to this space and thus (14) holds for the corresponding (finite) array of wavelet coefficients. The same applies to the corresponding Sobolev norm given by (10). Thus as long as the wavelets have sufficiently high order, the two norms can be evaluated for the given flow data in the corresponding regularity range by computing the weighted discrete norms. Therefore we cannot draw any asymptotic conclusions but can only make quantitative comparisons between the two norms. We use a vector-valued wavelet decomposition with smooth, biorthogonal wavelets (d=6, d˜=8). The result shown in Fig. 7 demonstrates that Sobolev and Besov norms differ more and more for higher values of s. With increasing s, the Besov norms (although large) even seem to become
Divergence-Free Wavelet Analysis of Turbulent Flows
65
18
10
16
10
14
Weigted Coefficient Norm
10
12
10
10
10
8
10
6
10
Sobolev (Hs) Besov (Bs )
4
10
τ,τ
2
10
0
1
2
3
4
5
6
s
Fig. 7. Sobolev and Besov regularity of a PIV snapshot.
independent of s whereas the Sobolev norm increases rapidly. We obtained similar curves for all snapshots of the PIV data and other choices for d and d˜. This indicates a significant potential for wavelet compression of flow data that will be addressed in future research. As far as we know the above conclusions have no theoretical backup yet. To our knowledge, so far regularity results for the relevant Besov scales have been proven for the elliptic case [8], and for 1D conservation laws [15]. In the elliptic case, the Besov regularity was shown to be about twice the Sobolev regularity. In [15] it was shown that even though the solutions to 1D conservation laws may have no positive Sobolev regularity, the Besov regularity in the relevant scale (13) is arbitrarily high. In view of the results in [5] this implies that nonlinear or adaptive techniques may converge with arbitrarily high order. So far, no analogous results for the incompressible Navier–Stokes equations are known let alone the effect of turbulence. ACKNOWLEDGMENTS We are very grateful to Wolfgang Dahmen for many intensive and fruitful discussions. This work has been supported by U.S. National Science Foundation within the project Hierarchical modeling of incompressible turbulent flows using divergence-free vector wavelets, Grant No. CTS-9613948. REFERENCES 1. Barinka, A., Barsch, T., Urban, K., and Vorloeper, J. (2001). The Multilevel Library: Software Tools for Multiscale Methods and Wavelets, Version 2.1, RWTH Aachen, IGPM Preprint 205. 2. Berkooz, G., Elezgaray, J., and Holmes, P. J. (1992). Coherent structures in random media and wavelets. Physica D 61, 47–58.
66
Albukrek et al.
3. Canuto, C., Tabacco, A., and Urban, K. (1999). The wavelet element method I: Construction and analysis. Appl. Comp. Harm. Anal. 6, 1–52. 4. Cohen, A., and Masson, R. (1999). Wavelet methods for second order elliptic problems— preconditioning and adaptivity. SIAM J. Sci. Comp. 21(3), 1006–1026. 5. Cohen, A., Dahmen, W., and DeVore, R. A. (2001). Adaptive wavelet methods for elliptic operator equations—convergence rates. Math. Comp. 70, 27–75. 6. Cohen, A., Daubechies, I., and Feauveau, J.-C. (1992). Biorthogonal bases of compactly supported wavelets. Comm. Pure and Appl. Math. 45, 485–560. 7. Dahlke, S., Dahmen, W., and Urban, K. (2001). Adaptive Wavelet Methods for Saddle Point Problems—Optimal Convergence Rates, RWTH Aachen, IGPM Report 204. 8. Dahlke, S., and DeVore, R. A. (1997). Besov regularity for elliptic boundary value problems. Comm. Partial Diff. Eq. 22(1/2), 1–16. 9. Dahmen, W. (1997). Wavelet and multiscale methods for operator equations. Acta Numerica 6, 55–228. 10. Dahmen, W., Kunoth, A., and Urban, K. (1999). Biorthogonal spline-wavelets on the interval—stability and moment conditions. Appl. Comp. Harm. Anal. 6, 132–196. 11. Dahmen, W., and Schneider, R. (1999). Composite wavelet bases for operator equations. Math. Comput. 68, 1533–1567. 12. Dahmen, W., and Schneider, R. (1999). Wavelets on manifolds I: Construction and domain decomposition. SIAM J. Math. Anal. 31, 184–230. 13. Dahmen, W., and Stevenson, R. (1999). Element-by-element construction of wavelets satisfying stability and moment conditions. SIAM J. Numer. Anal. 37, 319–325. 14. DeVore, R. A. (1998). Nonlinear approximation. Acta Numerica 7, 51–150. 15. DeVore, R. A., and Lucier, B. (1996). On the size and smoothness of solutions to nonlinear hyperbolic conservation laws. SIAM J. Math. Anal. 27, 684–707. 16. Farge, M. (1992). Wavelet transforms and their applications to turbulence. Ann. Rev. Flu. Mech. 24, 395–457. 17. Farge, M., Schneider, K., and Kevlahan, N. (1999). Non-Gaussianity and coherent vortex simulation for two-dimensional turbulence using an adaptive orthogonal wavelet basis. Phys. Flui. 11(8), 2187–2201. 18. Griebel, M., and Koster, F. (2000). Adaptive wavelet solvers for the unsteady incompressible Navier–Stokes equations. In Malek, J., et al. (eds.), Advances in Mathematical Fluid Mechanics, Springer, pp. 67–118. 19. Lemarié-Rieusset, P. G. (1992). Analyses multi-résolutions non orthogonales, commutation entre projecteurs et derivation et ondelettes vecteurs a` divergence nulle (in french). Rev. Mat. Iberoam. 8, 221–236. 20. Lewalle, J. (1993). Wavelet transforms of the Navier–Stokes equations and the generalized dimensions of turbulence. Appl. Sci. Res. 51(1/2), 109–113. 21. Lumley, J. L. (1967). The structure of inhomogeneous turbulence. In Yaglom, A. M., and Tatarski, V. I. (eds.), Atmospheric Turbulence and Wave Propagation, Moscow, Nauka. 22. Mallat, S. (1989). A theory for multiresolution signal decomposition: The wavelet representation. IEEE Patt. Anal. Mach. Intell. 11, 674–693. 23. Meneveau, C. (1991). Analysis of turbulence in the orthonormal wavelet representation. J. Fluid Mech. 232, 469–520. 24. Urban, K. (1995). On divergence-free wavelets. Adv. Comput. Math. 4, 51–82. 25. Urban, K. (2001). Wavelet bases in H(div) and H(curl). Math. Comp. 70, 739–766.
Divergence-Free Wavelet Analysis of Turbulent Flows Cem M. Albukrek, 1 Karsten Urban, 2 Dietmar Rempfer, 3 and John L. Lumley 1 Received August 27, 2001; accepted (in revised form) November 3, 2001 In this paper we study the application of divergence-free wavelet bases for the analysis of incompressible turbulent flows and perform several experiments. In particular, we analyze various nominally incompressible fields and study the influence of compressible perturbations due to experimental and computational errors. In addition, we investigate the multiscale structure of modes obtained from the Proper Orthogonal Decomposition (POD) method. Finally, we study the divergence-free wavelet compression of turbulent flow data and present results on the energy recovery. Moreover, we utilize wavelet decompositions to investigate the regularity of turbulent flow fields in certain non-classical function spaces, namely Besov spaces. In our experiments, we have observed significantly higher Besov regularity than Sobolev regularity, which indicates the potential for adaptive numerical simulations. KEY WORDS: Wavelets; turbulence.
1. INTRODUCTION The study, analysis and simulation of turbulent flows is probably one of the most challenging problems in engineering, scientific computing and applied mathematics. Wavelets offer some potential for both the analysis and the simulation of turbulent flows. This is mainly due to a good localization in both frequency and physical domain which entails a good separation of features associated with different characteristic length scales. Localization in physical space can be enhanced by employing compactly supported wavelets which all subsequent investigations will be confined to. A second important feature of wavelets is their compression property. By 1
Sibley School of Mechanical and Aerospace Engineering, Cornell University, Ithaca, New York 14853. E-mail: {cma13;jll4}@cornell.edu 2 RWTH Aachen, Institut für Geometrie und Praktische Mathematik, Templergraben 55, 52056 Aachen, Germany. E-mail: [email protected] 3 Illinois Institute of Technology, 3110 S. State St., Chicago, Illinois 60616. E-mail: RempferD @iit.edu 49 0885-7474/02/1200-0049/0 © 2002 Plenum Publishing Corporation
50
Albukrek et al.
this we mean that often only relatively few but judiciously chosen terms in an infinite wavelet expansion suffice to recover the function within a desired tolerance. This relies on the so called cancellation properties of wavelets combined with the fact that certain weighted sequence norms of the expansion coefficients are equivalent to the norms for certain smoothness spaces such as Sobolev or Besov spaces (such estimates will be called norm equivalences). Such highly economical representations are therefore particularly promising for representing functions with a wide range of characteristic length scales such as those representing turbulent flow fields. The primary objective of this paper is to demonstrate the usefulness of wavelet concepts and their key mechanisms for the analysis and compression of turbulent flows. In particular, wavelet concepts are used to analyze the Proper Orthogonal Decomposition (POD) method [21], which aims at approximating the infinite-dimensional nonlinear dynamical system describing turbulent flows by some low-dimensional model. A variety of approaches for this reduction can be found in the literature. In e.g., [2, 16–18, 20, 23], wavelets already have been used for the analysis and simulation of turbulent flows. In most of the cited papers, the data analysis has been performed by means of orthogonal wavelets for periodic vorticity fields. For the simulation of periodic flows mainly exponentially decaying vaguelettes have been employed so far, see e.g., [17]. In [18] interpolatory wavelets have been used for the simulation. In this paper we propose to utilize compactly supported divergencefree wavelets [19, 24, 25] for the analysis of velocity fields. Since we can therefore work with the velocity/pressure formulation of the Navier– Stokes equations these tools work naturally also for three-dimensional problems. Since the incompressibility constraint is satisfied by the basis functions, such systems are natural candidates for a quantitative analysis of the role of incompressibility and its perturbation. Specifically, we will be concerned with the following issues. First, divergence-free wavelets and a wavelet basis for an appropriate complement space give rise to a Helmholtz-type decomposition for spaces of vector fields [25]. Using the divergence-free wavelet transform, we can split flow data obtained from experiments as well as direct numerical simulations into its solenoidal and compressible parts. In particular, this tool offers insight into the stability of the POD method, which is a well-known technique for constructing a low-dimensional turbulence model. POD modes are given as the eigenfunctions of the autocorrelation tensor for an ensemble of flow realizations. They can be interpreted as the coherent structures in a turbulent flow. In order to produce an incompressible flow model based on the POD-Galerkin method, the corresponding data fields have to be divergence-free. However, it turns out that all the considered flows do in fact contain some compressible components. By comparing the POD modes obtained from the available data and those from its divergencefree part we study the quantitative influence of this violation on the incompressible flow model.
Divergence-Free Wavelet Analysis of Turbulent Flows
51
Second, the above mentioned ability of wavelet expansions to separate data into contributions from different characteristic length scales will be used to provide insight into the multiscale structure of various flows and, in particular, of the coherent structures. Finally, we exploit the compression properties of divergence-free wavelet bases. Here the term compression refers to the reduction of data sets and should not be confused with compressibility as a flow property. We compress several turbulent flow fields in the sense that we keep possibly few terms in their wavelet expansion while retaining most of the energy. Such a reduction without loosing essential information is closely related to the regularity of the underlying field in a certain scale of Besov spaces [14]. Due to the above mentioned norm equivalences we are able to quantify the regularity of given data sets in the relevant (Besov) scale and draw conclusions about their adaptive recovery, see [5]. This latter aspect is particularly interesting with regard to adaptive direct numerical simulations. This paper is organized as follows: In Section 2, we review the relevant wavelet theory that is needed in this paper and Section 3 is devoted to a brief review of the POD method. Our particular implementation is described in Section 4. Sections 5 and 6 contain our numerical results. We will frequently use the following notation throughout the remainder of this paper. A ’ B means that there exists constants C \ c > 0 such that cA [ B [ CA and the constants are independent of all the parameters A and B may depend on.
2. WAVELETS AND SIGNAL ANALYSIS In this section, we review the relevant properties of wavelet bases that will be needed in the sequel.
2.1. Biorthogonal Wavelets Let us denote by t, t˜ ¥ L2 (R) a pair of functions that are refinable, i.e., there exist some sequences a :=(ak )k ¥ Z , a˜ :=(a˜k )k ¥ Z ¥ a2 (Z) such that t(x) := C ak t(2x − k),
t˜(x) := C a˜k t˜(2x − k)
k¥Z
(1)
k¥Z
for x ¥ R. Moreover, t and t˜ are dual or biorthogonal, i.e., (t( · ), t˜( · − k))0, R =d0, k , k ¥ Z, where ( · , · )0, R refers to the usual L2 (R)-inner product. We will call t and t˜ dual scaling functions. Defining for j, k ¥ Z the scaled and translated version of a function f: R Q R by fj, k (x) :=2 j/2f(2 jx − k),
x¥R
(2)
52
Albukrek et al.
the pair of dual scaling functions generate a sequence of nested spaces ˜j :=closL (R) span{t˜j, k : k ¥ Z}, defined by Vj :=closL2 (R) span{tj, k : k ¥ Z}, V 2 which are referred to as biorthogonal multiresolution analysis, [22]. ˜ j are defined by Vj+1 =Vj À Wj and The complement spaces Wj , W ˜ ˜ ˜ ˜ j and V ˜j + Wj . Finally, g, g˜ ¥ L2 (R) Vj+1 =Vj À Wj , j ¥ Z, as well as Vj + W ˜ j= are called biorthogonal wavelets if Wj =closL2 (R) span{gj, k : k ¥ Z}, W closL2 (R) span{g˜j, k : k ¥ Z} and the systems {gj, k : j, k ¥ Z}, {g˜j, k : j, k ¥ Z} form Riesz bases for L2 (R). Biorthogonal wavelets are defined by the relations g(x) := C bk t(2x − k), k¥Z
g˜(x) := C b˜k t˜(2x − k)
(3)
k¥Z
where bk :=(−1) k a˜1 − k , b˜k :=(−1) k a1 − k , see (1). Biorthogonal wavelets can be constructed from different families of scaling functions. One prominent example is formed by using cardinal B-splines of order d with nodes at the integers as primal scaling functions, which will be denoted by d t. In [6], a whole family d, d˜ t˜ for d+d˜ even, of dual scaling functions for cardinal B-splines was constructed. With increasing d˜, the polynomial exactness of the system induced by d, d˜ t˜ grows at the expense of increasing size of their supports. The corresponding biorthogonal wavelets will be denoted by d, d˜ g and d, d˜ g˜, respectively. However, in many places we will drop these indices for notational convenience, except when their specific values matter. Nowadays biorthogonal wavelet bases are also available in L2 ([0, 1]), [10], in L2 (W) for a wide class of bounded domains W … R n and even manifolds, see e.g., [3, 4, 11–13]. We will denote scaling and wavelet systems in L2 (W) by Xj :={tj, k : k ¥ Ij }, Uj :={gj, k : k ¥ Jj }, respectively, where Ij and Jj are finite index sets. The spaces Vj and Wj are defined analogously. Finally, the full wavelet basis will be denoted by U :={gl : l ¥ J}, where l=(j, k) and J=1j Jj (and J=Z 2 if W=R). 2.2. Wavelet Transform Given a function vj+1 =; k ¥ Ij+1 cj+1, k tj+1, k ¥ Vj+1 with single scale coefficients cj+1, k , the refinement equations (1), (3) can be used to split vj+1 into a coarse part (average) vj ¥ Vj with coefficients cj, k and a detail part wj expanded in terms of the wavelets gj, k with coefficients dj, k . In fact, these new coefficients are obtained by the recursion cj, k = C m ¥ Ij+1
2 −1/2a˜m − 2k cj+1, m ,
dj, k = C
2 −1/2b˜m − 2k cj+1, m
(4)
m ¥ Ij+1
Repeating this process results in a transformation that is often called decomposition or discrete wavelet transform (DWT) algorithm. Conversely, given the set of coarse scaling coefficients, c0, k , k ¥ I0 , and detail coefficients, dj, k , k ¥ Jj , j \ 0, the refinement equations (1), (3) can be used to
Divergence-Free Wavelet Analysis of Turbulent Flows
53
recover the scaling function coefficients recursively on the highest level through cj+1, k = C 2 −1/2ak − 2m cj, m + C 2 −1/2bk − 2m dj, m , m ¥ Ij
k ¥ Ij+1
(5)
m ¥ Jj
This process reconstructs a signal given by its coarse scale parts and the corresponding details. It is therefore also called reconstruction, or inverse discrete wavelet transform (IDWT). 2.3. Divergence-Free Wavelets Now, we briefly review the construction of divergence-free wavelets as introduced in [19, 24, 25]. The key point is that even for 1D biorthogonal wavelet systems there is a simple mechanism that relates different wavelet bases by differentiation and integration. Then, by fixing one component of a vector wavelet, one can arrange the remaining ones in order to derive divergence-free wavelets. From now on, we concentrate on biorthogonal B-spline wavelets only, i.e., the scaling functions t=td are cardinal B-splines of order d. It turns out that they are quite suitable for the construction of divergence-free wavelets due to certain relations involving their derivatives and primitives, [19]. To be precise, the following relations hold: d
tŒ(x)=d − 1 t(x+m(d − 1)) − d − 1 t(x − m(d − 1))
(6)
for the scaling functions, where m(d) :=d mod 2 as well as d, d˜
gŒ(x)=4 d − 1, d˜+1 g(x)
(7)
for the wavelets. We will frequently use the following short hand notations for the involved functions: t (1) :=d t, and g (1) :=d, d˜ g, as well as t (0) :=d − 1 t, and g (0) :=d − 1, d˜+1 g. A formula corresponding to (7) also holds for wavelet systems on the interval, see [25]. In fact, viewing Uj as a column vector consisting of the (0) wavelets gj, k on [0, 1], we obtain dxd U (1) j =Dj U j , where Dj is invertible. For the scaling function systems Xj , we obtain a similar relation, namely d (1) (0) dx X j =Cj X j , where the matrices Cj have full rank (but may in general be rectangular). As a consequence, we obtain
1 dxd X 2=span(X d span 1 U 2 =span(U dx span
(1) j
(0) j
(1) j
(0) j
) (8) )
54
Albukrek et al.
where span always denotes the closure of the linear span. Now, we use a tensor product approach (here, for simplicity only in 3D) and obtain wavelet systems on [0, 1] 3 by Y (ij, 1e, i2 , i3 ) :=G j,(i1e)1 é G j,(i2e)2 é G j,(i3e)3
(9)
e ¥ E g :={0, 1} 3 0 {(0, 0, 0) T}, i1 , i2 , i3 ¥ {0, 1} and
˛ XU
G (ij, ne)n :=
(in ) j (in ) j
if en =0 if en =1
This results in a diagram depicted in Fig. 1. In particular, we obtain also wavelet bases in the spaces H(div; W) and H(curl; W) defined by H(div; W) :={f ¥ L2 (W) : div f ¥ L2 (W)} and similarly for H(curl; W). Moreover, for the space H(div; W) it is possible to identify and characterize the subspace V(div; W) :={f ¥ H(div; W) : div f=0} of divergence-free functions. This means that one can construct functions k df l such that df Y df :={k df l : l¥J }
(J df being an appropriate index set) is a wavelet basis for V(div; W) and that each vector field fj ¥ Wj :=span Y div j … H(div; W) can be decomposed into its divergence-free part and some complement in the following way, dc df dc fj =Q df j fj +Q j fj , where Q j , Q j denote the biorthogonal projection onto df dc dc W df j :=span Y j and W j :=span Y j , respectively. The latter is a stable (but dc not necessarily orthogonal) complement in the sense that Wj =W df j À Wj .
Fig. 1. Construction principle for wavelet bases in H(div; W) and H(curl; W), where S(Y) := span Y. The solid lines indicate that the images of the spaces on the left of the arrow under the corresponding differential operator coincide with the spaces generated by the functions on the right of the arrow. Dashed lines indicate that the functions form a wavelet system for the corresponding function space.
Divergence-Free Wavelet Analysis of Turbulent Flows
55
Let us finally mention that the corresponding multidimensional filters can easily be derived from those of the univariate functions (see e.g., [6]) by building tensor products as in (9) and performing the modifications described in (6) and (7). 2.4. Wavelet Thresholding A key property of wavelet bases is their ability to compress signals in the following sense. Once a wavelet representation of some data is given, one can obtain an economic representation by discarding all wavelet coefficients below some (possibly scale dependent) threshold or equivalently by keeping a certain number of (possibly weighted) largest coefficients. The main analytical tool for quantifying such concepts is a family of norm equivalences. In fact, under certain assumptions on the primal and dual functions, one can show, [9], that for some c, ˜c > 0 depending on the Sobolev regularity of t, t˜, one has
>C d g> l
l¥J
2
’ C 2 2s |l| |dl | 2,
l
s, W
s ¥ (−c˜, c)
(10)
l¥J
where l=(j, k), |l| :=j and H s(W) denotes the Sobolev space (of order s ¥ R, normed by || · ||s, W ). Similar relations hold for divergence-free functions in Sobolev spaces when the above described divergence-free wavelets are used. In fact, it has been shown in [19, 25] that
>
C dl k df l
l ¥ J df
>
2
’ C rs, |l| |dl | 2, s, W
l ¥ J df
˛ 1,1+2
rs, j =
2j+2
,
if s=0 if s=1 (11)
Once relations like (10) or (11) are available, the error caused by discarding all coefficients below some threshold can be estimated in any of the above norms. In the present context, we apply the following version of thresholding based on measuring accuracy in H s. Rather than discarding all coefficients below some threshold we keep as few coefficients as possible so that the norm of the compressed expansion is still o times the original norm, where o ¥ (0, 1) is a fixed preservation rate. Thus one looks for the smallest set L(o) … J df such that ; l ¥ L(o) rs, |l| |dl | 2 \ o ; l ¥ J df rs, |l| |dl | 2. In view of (11) this can in fact be done by first reordering the vector v :=(vl )l ¥ J df , vl :=`rs, |l| |dl | in decreasing order of the modulus of its components, then summing in that order the squares of its entries until the sum reaches o times the full sum on the right-hand side of (11) which we denote by ||v|| 2(s) for a moment. Clearly, if N(o) is the number of summands obtained in this way, the set L(o) consists of the N(o) largest weighted wavelet coefficients. The error by which the compressed representation ; l ¥ L(o) dl k df l deviates
56
Albukrek et al.
from the original field in the norm || · ||s, W , is proportional to `1 − o ||v||(s) and obviously also to the error of the best N(o)-term approximation, i.e., the best approximation that can be obtained by any linear combination of N(o) basis functions, [14]. 3. PROPER ORTHOGONAL DECOMPOSITION In general, turbulent shear flows contain discernible, self-sustaining components, called the coherent structures, which account for a major part of the turbulent kinetic energy. Hence, it is natural to look for a good approximation to these structures in order to use them as a basis for a flow simulation. One prominent example of such a basis can be constructed by the POD method. We briefly outline here the main ingredients that will be needed in the sequel. Let us consider a specific flow scenario in a domain W … R n for a time interval [0, T]. We assume that we are given samples of velocity fields uj, k :=u(xk , tj ) ¥ R n, which can be seen as samples on a grid Wh = ( j − 1) T {x1 ,..., xM } … W of ‘‘snapshots’’ of the flow at a time tj := N − 1 ¥ [0, T], j=1,..., N. In order to determine a basis which is in a certain sense optimal in representing the ensemble U=(uj, k )j=1,..., N; k=1,..., M , one can determine f=(f1 ,..., fM ) T ¥ R M × n by
3 T1 C 1 C u z 2 ;||z|| 4 N
f :=arg max z
M
2
j, k
j=1
k
2 2
(12)
k=1
The summation over the snapshots can be interpreted as a time averaging. Note that (12) can be rewritten as
3 zz Rzz 4 , T
f :=arg max z
T
1 N where R := C uj u Tj T j=1
and uj =(uj, 1 ,..., uj, M ) T. This implies that f is the eigenvector of R corresponding to the largest eigenvalue, which coincides with the square of the 1 largest singular value of `T U. The vector f determined by (12) is the best approximation to the ensemble by a single vector. One obtains the desired basis {f i : i=1,..., J [ N} by taking the eigenvectors of R corresponding to the J largest eigenvalues, i.e., the principle axes of the cloud of data points. These vectors are also called the POD modes. They can be interpreted as capturing the coherent structures of the flow. These (finitely many) modes are often subsequently used as (discretizations of ) test and trial functions in a spatial Galerkin discretization of the Navier–Stokes equations which results in a system of ordinary differential equations. Thus one obtains a very economical approximate solution that can be used to further study flow characteristics for the given environment for comparisons with other simulations.
Divergence-Free Wavelet Analysis of Turbulent Flows
57
4. IMPLEMENTATION AND VALIDATION Implementation. We implemented two-dimensional and three-dimensional versions of the divergence-free and the complement wavelet transform for realizing the compression and analysis purposes indicated above. For visualization purposes, only the two-dimensional results will be presented. The underlying code was written in Matlab. In addition, C++ routines from the MultiScale Library (MSLib), [1] were used. To study the possible range of data reduction (wavelet compression) it is important to test different filters that give rise to different orders of cancellation properties and ranges for the norm equivalences. Therefore the test data has to be brought into a filter independent format. We use the following approach. First, we compute the piecewise linear interpolant for the given discrete sample values. Of course, the resulting function has an infinite expansion in terms of the divergence-free wavelets. Therefore, in the next step, we compute the biorthogonal projection onto a multiresolution space of sufficiently high level exceeding the sample resolution. This level is chosen so that the maximum of the relative errors at the grid points stays below 0.2%. Validation. We have used two kinds of data for our tests. First, as a basis for comparison, we used data for a periodic flow in a box which has kindly been provided to us by K. Schneider, see [17]. We have also validated our code by comparing our results with the data reported in [17]. Then, we used several snapshots of a two-dimensional velocity field obtained from the Particle Imaging Velocimetry (PIV) measurements of the turbulent wake behind an airfoil with the parameters displayed in Table I. A slice of a snapshot is shown in Fig. 2. This data has been kindly made accessible to us by M. Glauser and C. S. Yao from NASA Langley Research Center. Of course, only for low Mach numbers, as in the present case, air flow can be assumed to be incompressible. Moreover, when taking two dimensional slices one cannot expect to obtain strictly divergence-free fields in the turbulent flow regime, hence this data is particularly suitable for the analysis of compressible components in a given flow.
Table I. Parameters of Our PIV Test Data Experimental setup and grid specifications Airfoil Wind tunnel Flow No. of grid points
Chord length=101.6 mm, thickness=12.7 mm, angle of attack=4.5° Height=108 mm, width=254 mm Free stream velocity=18.1 m/s, mach number=0.05, Reynolds number=135.000 x: 168; y: 94; grid spacing: 0.471 mm
58
Albukrek et al.
Fig. 2. Slice of the PIV data.
5. DATA ANALYSIS In this section, we summarize our numerical results obtained from the application of wavelet decompositions to several data.
5.1. Artificial Compressibility In order to investigate the data, we decompose flow fields in order to determine their incompressible and compressible parts. As the PIV data consists of 2D streamwise slices of the corresponding 3D incompressible turbulent flow field, it is not surprising that the 2D vector fields are not exactly divergence-free. For given flow data, we compute the piecewise linear interpolant fj for some high level j as described above. Then, we dc compute Q df j fj and Q j fj , namely the incompressible part and the complement. By using the norm equivalences (10), (11), we compute the corresponding norms. For 240 different realizations of the flow field, the incompressible part accounts for 93% of the total energy on average. The next issue is the distribution of this energy over the scales of the flow. In Table II, we display the L2 -norms of the wavelet projection of two different snapshots onto different levels. These norms can be interpreted as the energies contained on particular levels of the flow. It can be observed that the loss of energy due to removing the complement part occurs mostly on the finer scales of the flow, for which the cross-stream turbulent flow fluctuations are important. The divergence-free transform hence also offers
Divergence-Free Wavelet Analysis of Turbulent Flows
59
Table II. Energies of the Original Field and the Divergence-Free Part Snapshot 15
Snapshot 140
j
orig.
div-free
perc.
orig.
div-free
perc.
4 3 2 1 0
0.0070 0.0148 0.0238 0.0383 0.0372
0.0063 0.0133 0.0222 0.0363 0.0347
90.4% 89.6% 93.1% 94.6% 93.3%
0.0070 0.0134 0.0228 0.0343 0.0418
0.0063 0.0118 0.0212 0.0316 0.0400
90.4% 87.9% 92.9% 92.4% 95.8%
natural means to quantify the associated experimental measurement error over the scales of the flow. Next, we use the DWT in order to investigate the multiscale structure of the given data. As in image processing, we consider plots of wavelet coefficients over the different levels. This offers insight into the distribution of the flow over the scales. However, as shown in Fig. 3, right-hand side, the contributions on finer scales decay so rapidly that they are no longer visible in the plot. This indicates that the main part of the flow is concentrated on the coarser scales. To show the structure of the high scale contributions the coefficients are amplified to attain essentially the same order of magnitude as the coarse scale contributions. In this way an approximate Normalized Divergence free Wavelet Coefficients for Levels 1 to 8
Divergence free Wavelet Coefficients for Levels 3 to 8
10 50
20 100 30
150 40
50
200
60 250 50
100
150
200
250
10
20
30
40
50
60
Fig. 3. Divergence-free wavelet decomposition of the data f from [17], i.e., f= ; j, e, k dj, e, k k df j, e, k (see (9)). The sets of coefficients dj, e :=(dj, e, k )k are shown for j=−1,..., 4 (starting from j=−1 in the top left corner, larger rectangles indicate larger j). On each level the contributions for e=(0, 1) (lower left corner), e=(1, 1) (diagonal) and e=(1, 0) (upper right corner) are shown separately.
60
Albukrek et al.
self-similarity is visible, see Fig. 3, left. Roughly speaking this corresponds to a measurement in H 1 rather than in L2 . 5.2. Analysis of the POD Modes As already mentioned, the POD modes may be interpreted as the coherent structures, i.e., those structures which capture the ‘‘main features’’ of the flow. The divergence-free wavelet analysis offers a tool to gain insight into at least two interesting issues concerning POD: Influence of Artificial Compressibility on the POD Modes. In many cases, extensive analysis on a set of turbulent flow data is carried out before an appropriate turbulence model can be developed. The kind of flow configuration that we investigate here is usually modeled by the incompressible Navier–Stokes equations. The POD modes extracted from data are used as test and trial functions in a Galerkin scheme for constructing a low-order numerical approximation of the full set of the Navier–Stokes equations. As an outcome, this imposes the divergence-free requirement on the POD modes. For the POD functions to be incompressible, the underlying flow data also has to meet the divergence-free condition. However, as already mentioned, data obtained from experiments or direct numerical simulations may contain errors. These errors may cause the flow field under consideration to be no longer divergence-free. The use of POD modes extracted from such data as a divergence-free basis in an incompressible flow simulation constitutes a violation on the used model. The proposed Helmholtz decomposition based on divergence-free wavelets offers a tool for quantifying the associated error. As two-dimensional slices of the PIV data include a naturally built-in artificial compressibility, we chose to use this data for our analysis. To this end, we first created a new set of divergence-free fields which were obtained by the divergence-free projection of the data. Then we applied POD decomposition on both sets, namely the original and the divergencefree part. In the first step, we compare the turbulent kinetic energies associated with the decomposed fields using the POD coefficients of the two data sets. From the stability of the POD method, it is expected that the energies of the original and divergence-free field are asymptotically equal. Here, we are interested in quantitative results. We find exactly the same percentage of energy recovery (93%) as mentioned above when we compare the norms of the POD coefficients corresponding to the original and filtered (i.e., divergence-free) data sets. The top plot in Fig. 4 shows the turbulent kinetic energy captured by the POD modes of the original data and its divergence-free part in a cumulative sense. The lower plot shows the percentage of the total kinetic turbulent energy captured by each individual mode. In both plots, the lines for the original data and its divergence-free part are quite close indicating the high quantitative stability of the POD with respect to perturbations.
Divergence-Free Wavelet Analysis of Turbulent Flows
61
Table III. Distribution of L2 -Norm Energy over the Scales of the POD Eigenfunctions Corresponding to the Original PIV Data Set j
mode 1
mode 25
mode 50
mode 100
mode 200
mode 240
9 8 7 6 5 4 3 2 1 0
0.0628 0.0885 0.1703 0.5583 0.6961 0.1267 0.1058 0.0790 0.0461 0.0142
0.0956 0.1959 0.3172 0.6266 0.5527 0.3346 0.1661 0.0143 0.0044 0.0006
0.1186 0.2527 0.4929 0.6135 0.4307 0.3167 0.0889 0.0361 0.0064 0.0182
0.1596 0.3755 0.6116 0.5316 0.4292 0.1627 0.0924 0.0458 0.0038 0.0011
0.2439 0.5450 0.6323 0.4422 0.2710 0.1339 0.0839 0.0133 0.0043 0.0011
0.8838 0.4662 0.2340 0.1264 0.0637 0.0216 0.0090 0.0022 0.0004 0.0001
Multiscale Structure of POD Modes. There has been an ongoing debate about whether the coherent structures have multiscale structure or not. A wavelet-based multiscale analysis of the POD modes interpreted as the coherent structures can shed some light on this question. For this purpose, we have performed wavelet transforms on several POD modes of the PIV data. For each mode starting from the most to the least energetic one we have computed the energy distribution over the different scales of the flow. The results were compiled in Table III. In all columns the maximal value is printed in boldface. The results indicate that the energy for the first (i.e., most energetic) few modes is concentrated around the 5th level of the decomposition. For the less energetic modes, the energy is distributed more evenly between the finer scales of the flow. Assuming that the coherent structures are represented by the first 25 POD eigenfunctions which capture Turbulent Kinetic Energy Distribution over Number of Modes
Percentage of Total Energy
100 80 PIV Data Incompressible Part
60 40 20 0
20
40
60
80
100 120 140 160 Number of POD Eigenfunctions
180
200
220
240
Turbulent Kinetic Energy Distribution over Modes
Percentage of Total Energy
7 6 5 4 3 2 1 0
0
5
10
15
20 25 30 POD Eigenfunction
35
40
45
50
Fig. 4. Kinetic turbulent energy (top) and percentage of total energy (bottom) for the original PIV data (straight lines) and the divergence-free part (dashed lines).
62
Albukrek et al.
about 60% of the turbulent kinetic energy, it can be seen that about 6 scales contribute most to these structures. Hence we see that at least for our data, the coherent structures have multiscale structure. Let us mention that we obtain similar results also for the data in [17] as well as for the divergence-free part of the given data. A final remark can be made from the above study: Our analysis reveals that the deviation from incompressibility is essentially limited to the finer scales of the flow, while the coherent structures are associated with the coarser scales of the flow. Hence the compressibility error induced by the construction of a low dimensional, incompressible flow approximation based on POD is by far less significant due to the use of only the first few dominant modes that represent the coherent structures. 6. COMPLEXITY OF TURBULENT FLOWS In this section we use the strong analytical properties of wavelets in order to investigate the structure of a turbulent flow. We are interested in two issues. First we study how sparse wavelet expansions of typical flow instances are, i.e., how many degrees of freedom are needed to obtain a certain energy recovery. Second, we investigate how effective an adaptive wavelet method can be at best for simulating incompressible turbulent flows. 6.1. Data Compression
100
100
90
90
80
80
Percentage of Energy Recovery
Percentage of Energy Recovery
In order to investigate the quantitative aspects of the compression properties of divergence-free wavelets, we decompose several snapshots of the given flow ensemble. The compression rate is the relation between the number of coefficients kept and the amount of energy retained in the field. Figure 5 is a typical example of the compression performance of divergencefree wavelets on the snapshots of the PIV data. 99% of the energy in the L2 -norm can be recovered by dropping 97% of the wavelet coefficients. On the other hand, the same recovery in the H 1-norm requires keeping 10% of the coefficients. Let us mention that the compression rates for the
70 60 50 40 30 20
0
0
20 40 60 80 Compression Percentage
60 50 40 30 20
H1 norm L2 norm
10
70
10
100
0 500 1000 1500 2000 2500 3000 Number of non zero Coefficients (out of 31764)
Fig. 5. Compression percentage (left) and number of retained coefficients (right) over recovered energy in L2 (solid line) and H 1 (dashed line).
Divergence-Free Wavelet Analysis of Turbulent Flows
63
periodic data are even higher and are comparable to those reported in [17]. The results displayed in Figs. 3 and 5 and have been computed with the basis corresponding to d=d˜=3 but we observe similar curves also for other choices of these parameters. This is related to the smoothness of our data which will be investigated in the next subsection. 6.2. Potential for Adaptivity A central objective for further research is to use divergence-free wavelets for a fully adaptive numerical simulation of turbulent flows. The corresponding spatial Galerkin projection of the Navier–Stokes equation could be done with respect to trial spaces spanned by adaptively chosen divergence-free wavelet bases. Quite recently, the connection between adaptive methods and nonlinear approximation has been analyzed rigorously for a wide class of problems, [5, 7]. In particular, it was shown under which circumstances it can be proven that the use of adaptive methods offers an asymptotic gain of efficiency over uniform refinement strategies. This means that the rate of convergence (relating accuracy and degrees of freedom) is higher for the adaptive than for the non-adaptive method. This is closely related to the regularity of the flow fields in certain smoothness scales which will be explained next. In the context of incompressible flows smoothness is usually measured in the scale of Sobolev spaces H s, i.e., derivatives are measured in L2 . The Besov spaces B sq (Lp ) can be used to measure smoothness of order s ¥ R in Lp (where the additional parameter q ¥ (0, .) permits some fine-tuning). Besov spaces arise as interpolation spaces between Sobolev spaces, see e.g., [14]. The key here is that they can also be characterized by weighted sequence norms based on wavelet expansions similar to (10), [5, 14]. Now fixing s and decreasing p the spaces B sq (Lp ) become larger and admit an increasingly stronger singular behavior of its elements. The Sobolev embedding theorem tells us how far p can be decreased so that the spaces remain continuously embedded in L2 . The critical value p=y is given by the relation 1 s 1 = + y n 2
(13)
where n denotes the spatial dimension. Figure 6 shows a diagram of these spaces. The horizontal axis corresponds to 1/p, i.e., the measure of integrability and the vertical axis shows the smoothness, s. The vertical line through (12 , 0) corresponds to the Sobolev spaces where the smoothness is measured in L2 . The diagonal line shows the Besov spaces with the relation of s and y given by (13). Hence, the space B sy (Ly ) is (much) larger than H s when s, y are related by (13). Moreover, this gap increases when s grows and hence y decreases. The point is that Besov regularity according to the critical embedding line is (almost) characterized by nonlinear approximation and, in particular, by the best wavelet N-term approximation.
64
Albukrek et al.
Fig. 6. Diagram indicating smoothness spaces.
Roughly the error of best N-term approximation of u decays like N −s/n when u ¥ B sy (Ly ) for s and y related by (13). Thus high regularity in this scale permits high accuracy at the expense of relatively few degrees of freedom. Knowing about the Besov regularity of a flow would allow us to predict the size of the sets L(o) in the above thresholding scheme. The adaptive schemes in [5, 7] was shown to perform in a certain range asymptotically as well as the best N-term approximation. Thus, the more the norm of the approximated object in H s exceeds its norm in B sy (Ly ), the higher the efficiency gain offered by the adaptive method is. In particular, when the approximated function belongs to B sy (Ly ) but not to H s the gain is even reflected by an asymptotically higher convergence rate. Therefore we are interested in the norms of the flows in the spaces B sy (Ly ) for growing s. The fact that these norms can indeed be computed relies on the above norm equivalences. In fact, it is well-known that when s and y are related by (13) one has [14] (see also (10))
> C dg> l l
l¥J
y s
B y (Ly (W))
’ C |dl | y
(14)
l¥J
This relation is valid for a range of s that depends on the underlying wavelet basis, especially its regularity. In order to use (14) to measure the Besov regularity of flow data recall that the data is discrete and therefore has to be mapped first into a function corresponding to an element of a multiresolution space for a very high level. This function automatically belongs (for any set of coefficients) to a given Besov space as long as the underlying wavelets belong to this space and thus (14) holds for the corresponding (finite) array of wavelet coefficients. The same applies to the corresponding Sobolev norm given by (10). Thus as long as the wavelets have sufficiently high order, the two norms can be evaluated for the given flow data in the corresponding regularity range by computing the weighted discrete norms. Therefore we cannot draw any asymptotic conclusions but can only make quantitative comparisons between the two norms. We use a vector-valued wavelet decomposition with smooth, biorthogonal wavelets (d=6, d˜=8). The result shown in Fig. 7 demonstrates that Sobolev and Besov norms differ more and more for higher values of s. With increasing s, the Besov norms (although large) even seem to become
Divergence-Free Wavelet Analysis of Turbulent Flows
65
18
10
16
10
14
Weigted Coefficient Norm
10
12
10
10
10
8
10
6
10
Sobolev (Hs) Besov (Bs )
4
10
τ,τ
2
10
0
1
2
3
4
5
6
s
Fig. 7. Sobolev and Besov regularity of a PIV snapshot.
independent of s whereas the Sobolev norm increases rapidly. We obtained similar curves for all snapshots of the PIV data and other choices for d and d˜. This indicates a significant potential for wavelet compression of flow data that will be addressed in future research. As far as we know the above conclusions have no theoretical backup yet. To our knowledge, so far regularity results for the relevant Besov scales have been proven for the elliptic case [8], and for 1D conservation laws [15]. In the elliptic case, the Besov regularity was shown to be about twice the Sobolev regularity. In [15] it was shown that even though the solutions to 1D conservation laws may have no positive Sobolev regularity, the Besov regularity in the relevant scale (13) is arbitrarily high. In view of the results in [5] this implies that nonlinear or adaptive techniques may converge with arbitrarily high order. So far, no analogous results for the incompressible Navier–Stokes equations are known let alone the effect of turbulence. ACKNOWLEDGMENTS We are very grateful to Wolfgang Dahmen for many intensive and fruitful discussions. This work has been supported by U.S. National Science Foundation within the project Hierarchical modeling of incompressible turbulent flows using divergence-free vector wavelets, Grant No. CTS-9613948. REFERENCES 1. Barinka, A., Barsch, T., Urban, K., and Vorloeper, J. (2001). The Multilevel Library: Software Tools for Multiscale Methods and Wavelets, Version 2.1, RWTH Aachen, IGPM Preprint 205. 2. Berkooz, G., Elezgaray, J., and Holmes, P. J. (1992). Coherent structures in random media and wavelets. Physica D 61, 47–58.
66
Albukrek et al.
3. Canuto, C., Tabacco, A., and Urban, K. (1999). The wavelet element method I: Construction and analysis. Appl. Comp. Harm. Anal. 6, 1–52. 4. Cohen, A., and Masson, R. (1999). Wavelet methods for second order elliptic problems— preconditioning and adaptivity. SIAM J. Sci. Comp. 21(3), 1006–1026. 5. Cohen, A., Dahmen, W., and DeVore, R. A. (2001). Adaptive wavelet methods for elliptic operator equations—convergence rates. Math. Comp. 70, 27–75. 6. Cohen, A., Daubechies, I., and Feauveau, J.-C. (1992). Biorthogonal bases of compactly supported wavelets. Comm. Pure and Appl. Math. 45, 485–560. 7. Dahlke, S., Dahmen, W., and Urban, K. (2001). Adaptive Wavelet Methods for Saddle Point Problems—Optimal Convergence Rates, RWTH Aachen, IGPM Report 204. 8. Dahlke, S., and DeVore, R. A. (1997). Besov regularity for elliptic boundary value problems. Comm. Partial Diff. Eq. 22(1/2), 1–16. 9. Dahmen, W. (1997). Wavelet and multiscale methods for operator equations. Acta Numerica 6, 55–228. 10. Dahmen, W., Kunoth, A., and Urban, K. (1999). Biorthogonal spline-wavelets on the interval—stability and moment conditions. Appl. Comp. Harm. Anal. 6, 132–196. 11. Dahmen, W., and Schneider, R. (1999). Composite wavelet bases for operator equations. Math. Comput. 68, 1533–1567. 12. Dahmen, W., and Schneider, R. (1999). Wavelets on manifolds I: Construction and domain decomposition. SIAM J. Math. Anal. 31, 184–230. 13. Dahmen, W., and Stevenson, R. (1999). Element-by-element construction of wavelets satisfying stability and moment conditions. SIAM J. Numer. Anal. 37, 319–325. 14. DeVore, R. A. (1998). Nonlinear approximation. Acta Numerica 7, 51–150. 15. DeVore, R. A., and Lucier, B. (1996). On the size and smoothness of solutions to nonlinear hyperbolic conservation laws. SIAM J. Math. Anal. 27, 684–707. 16. Farge, M. (1992). Wavelet transforms and their applications to turbulence. Ann. Rev. Flu. Mech. 24, 395–457. 17. Farge, M., Schneider, K., and Kevlahan, N. (1999). Non-Gaussianity and coherent vortex simulation for two-dimensional turbulence using an adaptive orthogonal wavelet basis. Phys. Flui. 11(8), 2187–2201. 18. Griebel, M., and Koster, F. (2000). Adaptive wavelet solvers for the unsteady incompressible Navier–Stokes equations. In Malek, J., et al. (eds.), Advances in Mathematical Fluid Mechanics, Springer, pp. 67–118. 19. Lemarié-Rieusset, P. G. (1992). Analyses multi-résolutions non orthogonales, commutation entre projecteurs et derivation et ondelettes vecteurs a` divergence nulle (in french). Rev. Mat. Iberoam. 8, 221–236. 20. Lewalle, J. (1993). Wavelet transforms of the Navier–Stokes equations and the generalized dimensions of turbulence. Appl. Sci. Res. 51(1/2), 109–113. 21. Lumley, J. L. (1967). The structure of inhomogeneous turbulence. In Yaglom, A. M., and Tatarski, V. I. (eds.), Atmospheric Turbulence and Wave Propagation, Moscow, Nauka. 22. Mallat, S. (1989). A theory for multiresolution signal decomposition: The wavelet representation. IEEE Patt. Anal. Mach. Intell. 11, 674–693. 23. Meneveau, C. (1991). Analysis of turbulence in the orthonormal wavelet representation. J. Fluid Mech. 232, 469–520. 24. Urban, K. (1995). On divergence-free wavelets. Adv. Comput. Math. 4, 51–82. 25. Urban, K. (2001). Wavelet bases in H(div) and H(curl). Math. Comp. 70, 739–766.
