A high accuracy Leray-deconvolution model of turbulence and its limiting behavior. 7 0 William Layton∗ Roger Lewandowski† 0 2 n a J Abstract 2 1 In 1934 J. Leray proposed a regularization of the Navier-Stokes equations whose limits were weak solutions of the NSE. Recently, a modification of the 1 v Leray model, called the Leray-alpha model, has atracted study for turbulent 0 flow simulation. One common drawback of Leray type regularizations is their 4 low accuracy. Increasing the accuracy of a simulation based on a Leray regu- 0 1 larization requires cutting the averaging radius, i.e., remeshing and resolving 0 on finer meshes. This report analyzes a family of Leray type models of ar- 7 bitrarily high orders of accuracy for fixed averaging radius. We establish the 0 / basic theory of the entire family including limiting behavior as the averaging h p radius decreases to zero, (a simple extension of results known for the Leray - model). We also give a more technically interesting result on the limit as h t the order of the models increases with fixed averaging radius. Because of a this property, increasing accuracy of the model is potentially cheaper than m decreasing the averaging radius (or meshwidth) and high order models are : v doubly interesting. i X r MCS Classification : 76D05, 35Q30, 76F65, 76D03 a Key-words : Navier-Stokes equations, Large eddy simulation, Deconvolution models. 1 Introduction In 1934 J. Leray [Leray34a], [Leray34b] studied an interesting regularization of the Navier-Stokes equations (NSE). He proved that the regularized NSE has a unique, smooth, strong solution and that as the regularization length-scale δ 0, the → regularized system’s solution converges (modulo a subsequence) to a weak solution oftheNavier-Stokes equations. This modelhasrecently beenattracting new interest as continuum model upon which large eddy simulation can be based (see, e.g., the work of Geurts and Holm [GH03] and Titi and co-workers [CHOT05], [CTV05], ∗Department of Mathematics, University of Pittsburgh, Pittsburgh, PA 15260, USA; [email protected], www.math.pitt.edu/˜wjl, partially supported by NSF grant DMS-0508260. †IRMAR, UMR 6625, Universit´e Rennes 1, Campus Beaulieu, 35042 Rennes cedex FRANCE; [email protected],http://perso.univ-rennes1.fr/roger.lewandowski/ 1 [VTC05]). If w denotes a local spacial average of the velocity w associated with filter length-scale δ , the classical Leray model is given by ∂ w+w w ν w+ q = f, (1.1) t ·∇ − △ ∇ w = 0. ∇· Leray chose w = g ⋆w, where g is a Gaussian1 with averaging radius δ. In (1.1), ν δ δ is the kinematic viscosity, (w,q) denotes the model’s velocity and pressure, Ω is the flow domain and f is the body force, assumed herein to be smooth and divergence free. We take Ω = (0,2π)3, the initial condition w(0,x) = u (x) in Ω , 0 and impose periodic boundary conditions (with zero mean) on the solution (and all problem data) w(t,x+2πe ) = w(t,x) and w(t,x)dx = 0. j ZΩ The Leray model is easy to solve using standard numerical methods for the Navier-Stokes equations, e.g., [LMNR06], and, properly interpreted, requires no ex- tra or ad hoc boundary conditions in the non-periodic case. However, as an LES model it has three main shortcomings: The Leray model’s solution can suffer an accumulation of energy around the • cutoff frequency (Geurts and Holm [GH03]). The Gaussian filter is expensive to compute. • Themodel(1.1)haslowaccuracyonthesmoothflowcomponents,e.g., u NSE • || − w = O(δ2) , Section 4. || The accumulation of energy around the cutoff length-scale can be ameliorated by new ideas in eddy viscosity which focus its effects on the smallest resolved scales or by time relaxation models with similar motivations, [Guer], [SAK01a], [SAK01b], [SAK02]. The expense of computing the filtered velocity is reduced (as proposed by Geurts and Holm [GH03]) if the Gaussian filter is replaced by a differential filter (in- troduced into LES by Germano [Ger86]). The resulting combination of Leray model plus differential filter, called the Leray-alpha model, has attracted an explosion of recent interest because of its theoretical clarity and computational convenience. For recent work on the Leray-alpha model, see Geurts and Holm [GH03], [GH05], Guer- mond and Prudhomme [GP05], Cheskidov, Holm, Olson and Titi [CHOT05], Ilyin, Lunasin and Titi [ILT05], Chepyzhov, Titi and Vishik [CTV05], [VTC05] (among other works). Finally, there is the issue of the accuracy of (1.1) upon the large scales / smooth flow components which must be improved if (1.1) is to evolve from a descriptive regularization to a predictive model. 1 By other choices of convolution kernel, differential filters, sharp spectral cutoff and top-hat filters can be recovered. 2 In this paper we complement the above work on (1.1) by presenting the analysis of a new (and related) family of Leray-deconvolution models which have arbitrarily high orders of accuracy and include the Leray-alpha model as the zeroth order (N = 0) case. To present the Leray-deconvolution models, some preliminaries are needed. Although any reasonable filter can be used, for definiteness, we fix the averaging to be by differential filter. Thus, given 2π-periodic free divergence field w with zero mean value, its spacial average over O(δ) length-scales, denoted w, is the unique 2π-periodic solution of the Stokes problem2 (1.2) δ2 w+w+ π = w in R3, w = 0, w = 0. − △ ∇ ∇· ZΩ It can be shown thatπ is a constant in the equation above andtherefore the pressure term disappaers. Given a filter, the (ill-posed) deconvolution problem is then: given w, find a (useful and stable) approximation to w. Let D denote the map: w chosen approximation to w, introduced in section N → 2. The approximate de-convolution operators, N = 0,1,2, , we consider have the ··· asymptotic accuracy and stability properties (see Section 2) (1.3) φ = D φ+O(δ2N+2) , for smooth φ, N (1.4) ||DN||L(H0→H0) ≤ N +1 < ∞, where (1.5) H = v L2 (IR3), 2π periodic, v = 0, v = 0 . 0 loc ∈ − ∇· (cid:26) ZΩ (cid:27) The Leray-deconvolution model is then w = 0 and ∇· (1.6) ∂ w+(D w) w ν w+ q = f in R3 (0,T), w = 0, t N ·∇ − △ ∇ × ZΩ subject to initial and 2π periodic boundary conditions. Because of (1.3), the formal − accuracy of (1.4) on the smooth flow components as an approximation of the Navier- Stokes equations is O(δ2N+2),N = 0,1,2, , Section 4. When N = 0, D w = w so 0 ··· (1.6) reduces to the Leray-alpha regularization. The behavior of the model as N increases for δ fixed is beyond known Leray- type theories and relevant for practical computation. Indeed, decreasing δ for fixed N requires reducing the computational meshwidth (a process which increases the storage and computing time dramatically) and resolving (1.1) or (1.2). On the other hand, increasing N for δ fixed requires only the solution of one additional Poisson or Stokes-type problem per deconvolution step. The main theoretical contribution of this paper is, in Section 3, to resolve this limiting behavior of the model. We 2 This filter is very close to the Gaussian filter; it is the first sub-diagonal Pad´e approximation thereof (as in the rational model [GL00]); in the Camassa-Holm/alpha model [FHT01], (1.2) is often called a Helmholz filter. One possible extension of this differential filter to non-periodic boundary condition is given in Manica and Merdan [MM06]. 3 first prove existence and uniqueness of a smooth solution to (1.6), then we show that (modulo a subsequence), for fixed δ, as N the model solution w = w(N) → ∞ converges to a weak solution of the Navier-Stokes equations. To our knowledge, this is the first result on the limiting behavior of a family of turbulence models as the order of accuracy of a family of models on the large scales increases. The difference between the two limiting cases is, loosely speaking, that as δ 0 , w w in every → → reasonable sense and the deconvolution process inherits this: D w w for N N → fixed and as δ 0. However, since deconvolution is ill-posed and the deconvolution → operators D are only an asymptotic (in δ for very smooth solutions and for N N fixed) inverse, the limiting behavior as N is more delicate. → ∞ AnothernewfeaturesincludeanO(δ2N+2)errorboundfortheenergynormofthe model’s error (section 4): w(δ) u , provided the underlying solution of the NSE || − || Navier-Stokes equations is sufficiently smooth, and an estimate of the time averaged error < (u w) 2 >1 both for general weak solutions of the Navier-Stokes ||∇ NSE − ||H 2 equations and for weak solutions having the k−5 energy spectrum typical observed 3 in fully developed turbulent flows. (These first of the two estimates is connected to related results for approximate deconvolution models in [LL06a], [LL06b], [LL03] and [DE06] and the second is an extension of the authors work in [LL06b].) Of course, ultimately, practical computations require analytic guidance in bal- ancing the computational meshwidth with δ , N and other model and algorithmic parameters. 1.1 Notation and preliminaries We use to denote the L2(Ω) norm and associated operator norm. We always ||· || impose the zero mean condition φdx = 0 on φ = w,p,f and w . Recall that Ω 0 R (1.7) H = w L2 (IR3), 2π periodic, w = 0, w = 0 . 0 loc ∈ − ∇· (cid:26) ZΩ (cid:27) We also define the space function (1.8) H = w H1 (IR3), 2π periodic, w = 0, w = 0 . 1 loc ∈ − ∇· (cid:26) ZΩ (cid:27) We shall define the space H for every s in the same way. We can thus expand the s velocity in a Fourier series w(x,t) = w(k,t)e−ik.x, where k Z3 is the wave-number. ∈ k X The Fourier coefficients abre given by 1 w(k,t) = w(x,t)e−ik.xdx. (2π)3 ZΩ Magnitudes of k are definbed by k = k1 2 + k2 2 + k3 21, | | {| | | | | |} k = max k , k , k . ∞ 1 2 3 | | {| | | | | |} 4 2π The length-scale of the wave number k is defined by l = . Parseval’s equality k ∞ | | implies that the energy in the flow can be decomposed by wave number as follows. For w H , 0 ∈ 1 1 1 w(x,t) 2dx = w(k,t) 2 = (2π)3 2| | 2| | ZΩ k X b 1 = w(k,t) 2 ,where k Z3.  2| |  ∈ k k−1<|k|≤k X X Moreover, for s IR,  b  ∈ (1.9) H = v Hs (IR3), 2π periodic, v = 0, w(0) = 0, k 2s w(k,t) 2 < . s loc ∈ − ∇· | | | | ∞ ( ) k X We define the H norms by b b s (1.10) w 2 = k 2s w(k,t) 2, || ||s | | | | k X where of course w 2 = w 2. It can be showbn that when s is an integer, w 2 = || ||0 || || || ||s sw 2 (see [DG95]). ||∇ || 1.2 About the filter Let w H and w H be the unique solution to the Stokes problem. 0 1 ∈ ∈ (1.11) δ2 w+w+ π = w in R3, w = 0, w = 0. − △ ∇ ∇· ZΩ It is usual in deconvolution studies to denote the filtering operation by G so that w = Gw. Writing w(x,t) = w(k,t)e−ik.x, k X it is easily seen that π = 0 and b ∇ w(k,t) (1.12) w(x,t) = e−ik.x 1+δ2 k 2 k | | X b Then writing w = G(w), we see that in the corresponding spaces of the type H , s the transfer function of G, denoted by G is the function 1 G(k) = , b 1+δ2 k 2 | | and we also can write on Hs spacbes type (1.13) δ2 w+w = w in R3, w = 0, w = 0. − △ ∇· ZΩ Moreover, one notes that the transfer function depends only on the modulus of the wave vector k. Therefore, by noting k = k , we shall write in the following G(k) | | instead of G(k). b b 5 2 Approximate de-convolution operators The de-convolution problem is central in image processing, [BB98]. The basic prob- lem in approximate de-convolution is: given w solve approximately for w: (2.1) Gw = w. Exact de-convolution is typically ill-posed. We consider the van Cittert [BB98] approximate de-convolution algorithm and associated operators, introduced into LES modeling by Adams, Kleiser and Stolz, e.g., [AS01], [AS02], [SA99], [SAK01a], [SAK01b], [SAK02]. For each N = 0,1, , w = D w is computed as follows. N N ··· Definition 2.1 (van Cittert approximate de-convolution algorithm) Setw = 0 w, for n = 1,2, ,N 1, ··· − perform w = w + w Gw n+1 n n { − } Set w = D w N N By eliminating the intermediate steps, we find N (2.2) D w = (I G)nw. N − n=0 X The approximate deconvolution operator D is a bounded operator (Lemma 2.1) N which approximates the unbounded exact deconvolution operator to high asymp- totic accuracy O(δ2N+2) on subspaces of smooth functions (Lemma 2.2). We begin by summarizing from [DE06], [BIL06] a few important, known properties of the approximation deconvolution operator D . N Lemma 2.1 [Stability of approximate de-convolution] Let the averaging be defined by (1.2) . Then D is a self-adjoint, positive semi-definite operator on H and N 0 D = N +1. N || || Lemma 2.2 [Error in approximate de-convolution]Let G be given by the differential filter (1.2). For any w H , 0 ∈ (2.3) w D w = (I ( δ2 +1)−1)N+1w N − − − △ = ( 1)N+1δ2N+2 N+1( δ2 +1)−(N+1)w. − △ − △ 2.1 Some calculations Since we consider the periodic problem, it is insightful visualize the approximate de-convolution operators D in terms of the transfer function of the operator D . N N Since these are functions of δk , where k = k , and not δ or k, it is appropriate to | | record them re-scaled by d k δk. ← 6 1 Since G = ( δ2 +1)−1 , we find, after rescaling, G(k) = . With that, the − △ 1+k2 transfer function of the first three deconvolution operators are b D = 1, 0 1 2k2 +1 D = 2 = , and c1 − k2 +1 k2 +1 k2 k2 2 c D = 1+ + . 2 k2 +1 k2 +1 (cid:18) (cid:19) These three are plotted tocgether with the transfer function of exact de-convolution (k2 +1) (in bold). Figure 1: Exact and approximate (N =0,1,2)) de-convolution operators More generally, we find N 1 n k2 N+1 D (k) := = (1+k2) 1 . N 1+k2 − 1+k2 " # n=0(cid:18) (cid:19) (cid:18) (cid:19) X d The large scales are associated with the wave numbers near zero (i.e., k small). | | Thus, the fact that D is a very accurate solution of the de-convolution problem N for the large scales is reflected in the above graph in that the transfer functions have high order contact near k = 0. The key observation in studying asymptotics of the model as N is that (loosely speaking) the region of high accuracy grows → ∞ (slowly) as well as N increases. The regularizationin thenonlinearity involves a special combination of averaging and deconvolution that we shall denote by H . N Definition 2.2 The truncation operator H : H H is defined by N 0 0 → H w := D w = (D G)w. N N N ◦ 7 Proposition 2.1 For each N = 0,1, , the operator H is positive semi-definite N ··· on H . The operator H : H H is a compact operator. Moreover, H maps 0 N 0 0 N → continuously H onto H . Further, 0 2 H w = D w, N N ||HN||L(H0→H0) = 1 , for each N ≥ 0, ||HN||L(H1→H1) = 1 , for each N ≥ 0. Proof These properties are easily read off from the transfer function H (k) of N H which we give next. For example, compactness follows since H (k) 0 as N N | | → k . b | | → ∞ b Remark 2.1 Similarly, by the same proof, H maps H into itself (see Lemma 2.5 N s below), and is compact. Since H (k) 1, one has for each s 0 N ≤ ≥ (cid:12) (cid:12) (cid:12) (cid:12) (2.4) b H w w . (cid:12) N(cid:12) s s || || ≤ || || The Fourier coefficients/transfer function of the operator H are similarly easily N calculated to be (after rescaling by k δk) ← k2 N+1 (2.5) H (k) = 1 . N − 1+k2 (cid:18) (cid:19) b They are plotted below for a few values of N. Figure 2: Transfer function of HN, N =0,10,50) These plotsarerepresentative ofthebehavior ofthewholefamily. Examining the abovegraphs,weobservethatH (u)isveryclosetouforthelowfrequencies/largest N 8 solution scales and that H (u) attenuates small scales/high frequencies. The break- N point between the low frequencies and high frequencies is somewhat arbitrary. The following (from [LN06a]) is convenient for our purposes and fits our intuition of an approximate spectral cutoff operator. We take for k the frequency for which H c N most closely attains the value 1. 2 b Definition 2.3 (Cutoff-Frequency) The cutoff frequency of H is N 1 k := greatest integer H−1 . c N 2 (cid:18) (cid:18) (cid:19)(cid:19) b In other words, the frequency for which H most closely attains the value 1. N 2 From the above explicit formulas, itbis easy to verify that the cutoff frequency grows to infinity slowly as N for fixed δ and as δ 0 for fixed N. Other prop- → ∞ → erties (whose proofs are simple calculations) of the operator H ( ) follow similarly N · k2 N+1 easily from H (k) = 1 . N − 1+k2 (cid:18) (cid:19) b Lemma 2.3 For all N 0, H (k) has the following properties: N ≥ b 0 < H (k) 1, N ≤ H (k) 0 as k , N | | → b | | → ∞ H (0) = 1. N b Let kc be cutoff frequency, thenb k as N for fixed δ and as δ 0 for N fixed. c | | → ∞ → ∞ → For any fixed value of k (or bounded set of values of k ) H (k) 1 as N for δ fixed and as δ 0 for N fixed N → → ∞ → Proof. bThis follows from (2.5). To extract strong convergence of deconvolution, some restriction to the large scales is needed (as shown clearly in the above two figures). One way to do this is to consider the action of deconvolution operators on trigonometric polynomials. Definition 2.4 The space IP denotes the 2π-periodic, divergence-free, wit hzero n mean value, vector, trigonometric polynomials of degree n and Π : H IP n 0 n ≤ → is the associated orthogonal projection. The model error is driven by the error in approximate deconvolution. With this definition a critical convergence result on approximate deconvolution is possible. 9 Proposition 2.2 H is a compact, symmetric positive semi-definite operator. For N n fixed, H is strictly positive definite on IP and N n H I on IP in the strong operator norm on H . N n 0 → Given n there is an N (n) large enough such that 0 1 Π v 2 H2v , for all N N (n), || n || ≤ || N || ≥ 0 Further, 1 v H2v 2 v for all v P . || || ≤ || N || ≤ || || ∈ n Proof. This follows from (2.5). For example, Compactness follows since H (k) 0 as k . N → → ∞ Other properties of the operator H ( ) follow similarly easily from its transfer N · fbunction. Proposition 2.3 Let Π denote the orthogonal H projection into kc 0 span eik·x : k k . For all w H : c 0 { | | ≤ } ∈ (2.6) (H w,w) C Π w 2, N ≥ || kc || (w H w,w) C (I Π )w 2. − N ≥ || − kc || Proof The claims follow from the definition of the cutoff frequency, the explicit formula for the transfer function and direct calculation. The following properties are important to the analysis that follows. Lemma 2.4 Let s IR, w H . Then H (w) H . s N s+2 ∈ ∈ ∈ Proof. Since H (w) satisfies periodic boundary conditions, it is divergence free N and has zero mean value by construction. We check the regularity. It is easy checked by using a Taylor expansion that there exists a constant C = C(δ,N) such that C(δ,N) (2.7) k Z3, k = 0, 0 Hˆ (k) , with k = k. N ∀ ∈ 6 ≤ ≤ k2 Therefore, for w H , w(x) = w(k)e−ik·x ∈ s k∈Z (2.8) (H (k))2k2(s+P2) w(k) 2 C(δ,N)2k2s w(k) 2. N b | | ≤ | | Recall that here w(0) = 0. Therefore in (2.8), k = 0. Hence, by definition b 6 (2.9) b ||HN(w)||Hs+2 ≤ C(δ,N)||w||Hs, and the proof is finished. Remark 2.2 The constant C(δ,N) in the above proof blows up as δ 0 or as → N . → ∞ 10