Analysis of a Reduced-Order Approximate Deconvolution Model and its interpretation as a Navier-Stokes-Voigt regularization

We study mathematical and physical properties of a family of recently introduced, reduced-order approximate deconvolution models. We first show a connection between these models and the NS-Voigt model, and that NS-Voigt can be re-derived in the approximate deconvolution framework. We then study the energy balance and spectra of the model, and provide results of some turbulent flow computations that backs up the theory. Analysis of global attractors for the model is also provided, as is a detailed analysis of the Voigt model's treatment of pulsatile flow.


1.
Introduction. Over the past several decades, many large eddy simulation (LES) models have been introduced and tested, with varying degrees of success; see e.g. [8,41] for extensive overviews. Each of these models has its share of good and bad properties, and with each generation of models, the amount of "good" tends to increase while the amount of "bad" tends to decrease. Two attractive models of the current generation are approximate deconvolution models (ADMs) [1,44], and Voigt regularization models [12,33], both of which have recently seen significant interest from both mathematicians and engineers. These models are known to be well-posed [6,9,18,33], have attractive physical properties (such as energy and helicity conservation and cascades [33,35]), and have performed well in (at least some) numerical simulations [31,45], including when applied to the Euler equations for ideal fluids [12].
where D denotes the operator norm of D. These assumptions are true for most common types of approximate deconvolution defined with the Helmholtz filtering operation, including van Cittert [18], multiscale [19], Tikhonov, and Tikhonov-Lavrentiev [43], see also the review in [5]. While we will often be general in our use of D, in some circumstances such as computations and energy spectra calculations, it is necessary to select a specific operator. In these cases, we will use the most common type of approximate deconvolution, van Cittert, and we will denote it by D = D N , where N ∈ N denotes the order of deconvolution. This operator is defined in Section 2.
The RADM system (1.3)-(1.4) was derived from the fourth-order version of the ADM by making the approximation D ≈ F −1 in the viscous term, which reduced the fourth order system to second order. This system was proved to be well posed in [40] and of high order consistency with the spatially filtered NSE. In [24], with D chosen to be van Cittert approximate deconvolution, an efficient finite element implementation of the RADM performed very well on the Re τ = 180 turbulent channel flow benchmark test, as well as in a 2D benchmark flow problem with a contraction and two outlets.
Plan of the Paper. It is the purpose of this paper to study some fundamental mathematical and physical properties of the model (1.3)-(1.4). After providing mathematical preliminaries in Section 2, we show a strong connection between the RADM (1.3) and the Navier-Stokes Voigt (NS-Voigt) model (3.5) in Section 3. In fact, we show that the RADM can be considered as a NS-Voigt-deconvolution model. That is, if the regularization operator F −1 is replaced with D −1 F −1 , then the RADM is recovered after a change of variables. This is interesting and important since a significant amount of theory has been developed for the NS-Voigt model by Titi and coworkers in e.g. [28,33,36] as well as in [6]. We also show that NS-Voigt can be derived using the approximate deconvolution approach originally developed by Stolz and Adams [1,2,44,45]. In this respect we establish a strong tie between Voigt-type models and approximate deconvolution type models in general.
An energy analysis of the RADM is given in Section 4. This includes an energy balance and a Kraichnan-type analysis of the energy spectra of the model. Some numerical results are also given here, to back up the theory. Our results reveal that the RADM is more computable than the NSE and has spectral scaling similar to the Leray-α model [17], which matches that of true fluid flow on larger scales in the inertial range, but smaller scales are transferred more quickly to even smaller scales. Moreover, our analysis shows that increasing the deconvolution order N shortens the inertial range and creates less energy in each scale.
Section 5 considers long-time behavior of the model, and proves existence of smooth global attractors. Section 6 considers the pseudo-parabolic nature of the system, and presents a study of some analytical properties and the model's treatment of pulsatile flow, which is an interesting example because it allows for exact analytical solutions. Finally, conclusions are drawn in the final section, and prospective future work is discussed.
2. Preliminaries. We present here notation and mathematical preliminaries to make the analysis of later sections smoother. As is typical with theoretical analyses of fluid models, we restrict our attention herein to problems involving periodic flows on the box Ω = (0, 2πL) 3 . We stress that all machinery requires the space-periodic setting to have the commutation between the deconvolution operator and general differential operators. The application of ADM methods to motion of fluid bounded by rigid walls is problematic, and the theory is still missing. Some results in special geometries can be found in [4,45]. In this framework, NS-Voigt represents one of the few models which is naturally set in a bounded domain with Dirichlet boundary conditions. Furthermore, we consider the system (1.3)-(1.4) to be non-dimensionalized by introducing L, U , and L/U for the characteristic length, velocity, and time scales, respectively. Thus, the viscosity ν > 0 can be considered dimensionless and representing the inverse of the Reynolds number of the flow. To simplify the notation, in the theorems we will suppose without loss of generality that L = 1.
We characterize the divergence-free spaces by using Fourier Series on the 3D torus. We denote by (e 1 , e 2 , e 3 ) the orthonormal canonical basis of R 3 , and by x := (x 1 , x 2 , x 3 ) ∈ R 3 the standard point in R 3 . Let T be the torus defined by T := R 3 /(2πLZ) 3 . We use · and (·, ·) to denote the L 2 (T) norm and inner product, and we impose the zero mean-value condition on velocity, pressure, and external force. We will use the customary Lebesgue L p (T) and Sobolev W k,p (T) spaces, and for simplicity (if not explicitly needed) we do not distinguish between scalar and vector valued real functions.
We define, for any exponent s ≥ 0, where W s,2 (T) 3 := W s,2 (T) 3 and if 0 ≤ s < 1 the condition ∇ · w = 0 must be understood in a weak sense. For w ∈ V s , we can expand the velocity field with Fourier series as follows is the wave-number vector, and the Fourier coefficients are given by where, as above, w H 0 := w . The inner products associated with these norms are We can finally characterize V s ⊂ W s,2 (T) 3 as follows: Recall that due to the zero mean value, the Poincaré inequality holds in V 1 : There exists λ > 0 such that for every φ ∈ V 1 , This inequality will be used extensively in our analysis, and allows us to use the more convenient gradient norm in place of the W 1,2 (T) norm.
2.1. Approximate deconvolution. We assume the deconvolution operator D is self-adjoint and commutes with both ∆ and F . These assumptions are known to be satisfied for van Cittert approximate deconvolution [18,35], multiscale deconvolution [19], Tikhonov, and Tikhonov-Lavrentiev [43]. See also the reviews in [5,14] for more results about these operators. We denote the filtering operator by A := F −1 = (−α 2 ∆ + I), and we introduce now a family of abstract deconvolution operators {D N } N , which are of order zero (as differential or pseudo-differential operators), and such that We note that van Cittert approximate deconvolution, forms such a family, and the other types of deconvolution mentioned above can also be written to satisfy the above limiting property. Although we will use a general deconvolution operator D when possible, in some analysis (e.g. energy spectra analysis) and for computations, it is necessary to specify a specific operator. In these cases, we will use van Cittert approximate deconvolution, as it is one of the most popular in the ADM community. The basic properties we will use of the operator D = D N are summarized in the following lemma, see [9]. Lemma 2.1. For each N ∈ N the operator D N : V s → V s is self-adjoint, it commutes with differentiation, and the following properties hold true: If D N denotes the symbol of the operator D N , then This lemma shows that in the specific case of practical use, we can assume that the upper and lower bounds for D N are d 0 = 1 and d 1 = N + 1, which reflects in similar bounds for D N .

Known theory for RADM solutions.
We define now the notion of weak solution for the RADM.
, and if the following energy identity holds true: We have the following existence and uniqueness theorem, which is based on a standard Galerkin approach, together with the energy estimate (2.1) below, see [ The proof of this result is standard and follows the same guidelines of the results in [6,12], with the key point being the energy identity (2.1) below. From [40], with the assumptions on the data, we have that the solution v ∈ L ∞ (0, T ; V 1 ), and so all manipulations in this proof are justified. Multiplying (1.3) by Dv and integrating over the physical domain gives α 2 (∇v t , ∇Dv) + (v t , Dv) + (Dv · ∇Dv, Dv) − (q, ∇ · Dv) + ν ∇Dv 2 = (f, Dv), and then using a vector identity on the first two terms, and that ∇ · Dv = 0, both the pressure and nonlinear terms vanish, leaving us with Integrating in time over [0, T ] provides (2.1).
This equality shows that the Galerkin approximate solution v m = 1≤|k|≤m v k e ik·x is bounded uniformly in L ∞ (0, T ; V 1 ). Then, by making use of the comparison argument one can easily prove that In particular, this proves that one can extract from the Galerkin approximate functions a converging sub-sequence, and moreover shows also that the solution v ∈ C([0, T ]; V 1 ), which is a uniqueness class. The regularity of the time derivative allows us to show that an energy equality (instead of an inequality as for the NSE) holds true.
In the case of the problem without viscosity, i.e. the RADM-Euler equations by adapting techniques from [6,12] and by using the estimates on the operator D from Lemma 2.1, it is easy to infer the following result (the same applies also to the Euler-Voigt-deconvolution model, as discussed in the next section) with ∇ · f = 0. Then, there exists a unique solution v of the RADM Euler equa- 3. Remarks on the derivation of the model and a connection to NS-Voigt.
Interestingly, despite its derivation as an altered/reduced ADM, the proof for wellposedness of the system in [40] uses ideas of the well-posedness proof for NS-Voigt from [6]. This suggests a connection between the two models, and we show now that a NS-Voigt deconvolution model can be derived from (1.
and setting w = Dv yields This system (3.1)-(3.2) is precisely the NS-Voigt regularization except for the D −1 term included in the regularization, and thus we could refer to (3.1)-(3.2) as the NS-Voigt-deconvolution regularization. Since the usual regularization in NS-Voigt is given by just F −1 , and since D ≈ F −1 , this "new" regularization D −1 F −1 will provide a higher order of consistency for (3.1)-(3.2) to the NSE (which are recovered if no regularization is applied).
In this section we enlarge the derivation of the Voigt models as given in [40], and in particular we highlight the fact that Voigt models can be seen as ADMs. The Euler-Voigt equations read as and have been introduced in the LES community by Cao, Lunasin, and Titi [12]. They started from the inviscid "simplified Bardina", which is the same model studied in Layton and Lewandowski [34] (zeroth order Stolz and Adams) in the case ν = 0. The exact words are "...Therefore, we propose the inviscid simplified Bardina model as regularization of the 3D Euler equations that could be implemented in numerical computations of three dimensional inviscid flows." They then re-introduced the viscosity, obtaining the NS-Voigt equations (1.3)-(1.4). Later, in Larios and Titi [33], the model is justified as follows ". . . it was noted in Cao Lunasin and Titi [12] that formally setting ν = 0 in simplified Bardina amounts to simply adding the term −α 2 ∆u t to the left-hand side of the Euler equations yielding what they call the Euler-Voigt equations. Remarkably, if one reintroduces the viscous term −ν∆u the resulting equations happen to coincide with equations governing certain viscoelastic fluids known as Kelvin-Voigt fluids, which were first introduced and studied by A.P. Oskolkov [37]. " We show now that the NS-Voigt model can also be derived with usual modeling techniques and in the unified framework of deconvolution models. In addition to the motivation of the order reduction explained in [24], we follow now a pure approximate deconvolution derivation.
In order to the derive the model, apply the filter F =:

Rewriting the equation only in terms of the filtered variable
which is the starting point in the Adams and Stolz [1,44] approach.
Introduce now a family of abstract deconvolution operators {D N } N which are of order zero (explicit examples are those recalled in Section 2.1, see especially properties from Lemma 2.1) such that We then insert the deconvolution into the equations, which allows us to write Taking the zeroth order approximation, we get (in the new approximate variables Finally, by applying the operator A to both sides, we arrive at the NS-Voigt equations Following the same procedure, but using a higher accuracy operator D N (as in (3.1)) we obtain the NS-Voigt deconvolution model After applying A, we recover exactly the RADM model (1.3) studied in [24,40].

Energy analysis.
This section studies the treatment of energy of the RADM system. We begin with the energy balance, which is useful both in obtaining a priori estimates for existence theorems and as a starting point to study the RADM energy cascade. Then, a Kraichnan-type analysis of the RADM energy cascade is given, followed by some numerical experiments to test the energy cascade in forced, homogeneous, and isotropic turbulent flows.

Energy balance.
Energy is a quantity of fundamental importance in fluid flow modeling, both in terms of the physical relevance of its solutions and its potential for success in simulations. For a model to be successful in accurate coarse mesh simulations, it must dissipate sufficient energy to remain stable. We show here that the deconvolution in the RADM acts as an energy sponge, which seems to explain a smoothing effect from increasing deconvolution in the RADM that was observed in [24] with van Cittert approximate deconvolution, and in [26] with multiscale deconvolution.
We present now some integral identities and estimates, which -at the level of Faedo-Galerkin approximations-are the core to obtain existence of weak solutions. At the level of weak solutions (which are smooth enough) they are still valid and show the balance of the main flow quantities.
Remark 4.1. Equation (2.1) shows that the RADM conserves a "model energy" The conservation of energy is one of the outstanding open problems for the Navier-Stokes equations. The analysis of the conservation of the "model energy" for LES models started with the analysis in [4,9], and several ADM models share this property.
Remark 4.2. We note that since D ≥ 1, the balance (2.1) suggests that deconvolution in the RADM increases the effective viscosity. This is consistent with numerical results in [26,24], which reported a smoothing effect by increasing the deconvolution order.

4.2.
Energy spectra. This section presents an analytic and computational study of the RADM energy spectra. A widely accepted theory of turbulent flows starting with intuitions of L.F. Richardson and A.N. Kolmogorov is that energy is input at large scales, transferred from large to small scales through the inertial range, and finally dissipated at small scales by viscosity [8,21,22]. It is further believed that viscosity has no effect except on very small scales, on which it acts to decay exponentially fast.
To fully resolve a turbulent flow, one must resolve the NSE on this entire range of active scales, which is far beyond the limit of modern computational tools for high Reynolds number problems. Computing the NSE without full resolution is widely known to be unstable. For a model to be computable in practice, its range of active scales needs to be significantly smaller than for the NSE, and thus studying its energy spectra can help make such a determination.
Since the RADM takes a form similar (in some sense) to that of Leray-α, NS-α, and the simplified Bardina models, we adapt physics-based energy spectra analysis performed in [12,17,20] for these respective models. We begin with a decomposition of the velocity into its Fourier modes, which gives the energy balance 1 2 where D N (k) is the Fourier symbol of the deconvolution operator D N , and Note that T k − T 2k represents the net amount of energy transferred into wave numbers between [k, 2k). Time averaging the energy balance (4.2) gives the energy transfer equation Since the RADM model energy is defined by (4.1), as this is the conserved quantity in the absence of viscosity and external forces, we take D = D N and define the model energy of an eddy of size k −1 by Combining this with the energy transfer equation (4.3) gives the relation with the last relation coming from Lemma 2.1, since D N (k) ∼ (N +1) for sufficiently large k. If the wave number k belongs to the inertial range, then T k ≈ T 2k and there is no leakage of energy through dissipation. Then, we have that k is in the inertial range if This shows that increasing the deconvolution order N has the effect of shortening the inertial range of the RADM. This is consistent with the energy balance given above, which suggests that increasing N reduces the total energy and increases the effective viscosity.
To determine the kinetic energy distribution in the inertial range, we begin by defining the average velocity of a size k −1 eddy, and relating it to model energy via The total energy dissipation rate is determined by the energy balance to be M := ν ∇Dv 2 , and from the Kraichnan energy cascade theory gives that the corresponding turnover time for eddies of this size is Solving for E M (k) provides the relation As is the case for related models, such as those of α-type or the simplified Bardina model, the kinetic energy spectrum can be divided into 2 parts. If k α > O (1), then These scalings suggest that for larger scales in the inertial range, a k −5/3 roll-off of energy is expected, but for smaller wave numbers, the slope increases to k −3 . This implies that significantly less energy will be held in higher wave numbers, suggesting the model is indeed more easily computable than the NSE, which has a k −5/3 roll-off of energy through its entire inertial range.

4.3.
Numerical testing of energy spectra. To test some of these scaling results, a direct numerical simulation (DNS) of the non-dimensionalized RADM is performed for forced turbulence, which is one of the most extensively simulated turbulent flows (see e.g., [11,15,25,29,30]). In particular, we study the energy spectrum of the model for three-dimensional homogeneous and isotropic turbulent flow, on a periodic cubic box with a side length 2π, with Reynolds number Re = 200. We use a pseudo-spectral method for the spatial discretization of the model with full de-aliasing, and a second-order Adams-Bashforth scheme for the time stepping. In addition, we employ a forcing method used by Chen et al. [15,16] and She et al. [42]. The numerical forcing of a turbulent flow is the artificial addition of energy at the large scales in a numerical simulation. Statistical equilibrium is achieved by the balance between the input of kinetic energy through the forcing and its output through the viscous dissipation. Van Cittert deconvolution is used with order N = 0, 1, 2, and we will specify N in the plots.
We first begin by testing the scaling law in the previous section, which predicts that the kinetic energy scales with k −5/3 at the beginning of the inertial range, and transitions to k −3 near the end of the inertial range. A plot of the resulting energy spectra from a computation with N = 2 and α = 1/32 at the resolution of 128 3 is shown in Figure 1. Along with the spectrum, on the log-log plot are also lines with slopes -5/3 and -3, and we observe that the spectrum appears to be in agreement with these slopes near the beginning and end of the inertial range, respectively.
We also investigate the sensitivity of N on the energy spectra of the RADM. Here, we compute with 128 3 resolution, α = 1/16 and N = 0, 1, 2. Results are shown in Figure 2, and we observe (as our analysis predicts) the RADM energy curves shorten as N increases; that is, as N is increased, there is less (or equal) energy in each wave number across the entire spectrum. This shows there is less total energy in the N = 2 system than N = 1, and N = 1 compared to N = 0, which is in agreement with the above analysis.  Lastly, we use this experiment setup to test the accuracy of the RADM energy spectra against a DNS of the NSE. In Figure 3, we show the resulting energy spectra for the RADM with various values of α for N = 2 at the resolutions of 63 3 and 128 3 along with the DNS energy spectrum obtained from the NSE at 256 3 resolution. The energy spectra are provided for three cases of α = 1/8, 1/16, and 1/32. As α decreases, the energy spectra of the RADM become closer to the highly resolved DNS energy spectrum.

5.
Existence of global attractors for the RADM. The long-time behavior of evolution equations can be characterized by their attractors, which in many important cases are finite dimensional. In particular, the notion of global attractor (cf. [46]) seems the most relevant in the context of the NSE. At present it is wellknown that such an attractor exists for the 2D NSE, while the lack of well-posedness in the 3D case is reflected in partial results concerning weak or trajectory attractors. In the case of LES models, the solution has naturally better regularity properties (even if a certain hyperbolic behavior is expected), hence it is natural to expect existence of such an object. In particular, for the NS-Voigt model the analysis of these questions has been addressed in [27,28], where existence and regularity of the global attractor and also determining modes are studied. Since the NS-Voigt model is a peculiar system (pseudo-parabolic, see also results in Section 6.2) with also non-viscous features, the proof of existence of the global attractor uses techniques typical of hyperbolic systems and the semigroup associated with the solution results to be only asymptotically compact. In particular see the questions related with the regularity of the solution in the space variables recalled in Sec 6.1.
In this section we follow very closely the approach in [28] and we just highlight the changes needed to adapt the proof to the new RADM system. The main difference with respect to the latter reference is that, due to the presence of approximate deconvolution operators, we need to deal with the space-periodic case.
For simplicity of analysis and notation, following [28], in this section we consider the model (1.3)-(1.4) in the equivalent form v t + α 2 Av t + P (Dv · ∇Dv) + νADv = P f, where P is the Helmholtz-Leray orthogonal projection in L 2 (T) onto V 0 , and A := −P ∆ is the Stokes operator with periodic boundary conditions (not to be confused with the Helmholtz operator of the Sections 2 3, still denoted for historical reasons with the same symbol). The existence of a global attractor will follow from a wellknown result, for which we refer for instance to [46].
where Z(t) is compact in V 1 for each 0 < t 0 ≤ t. Assume also there exists a continuous k : [t 0 , ∞) × R + → R + such that for every r > 0, k(t, r) → 0 as t → ∞ and that The main result of this section is the following.
Theorem 5.2. Let v 0 ∈ V 1 and assume the forcing is only spatially dependent, the continuous semigroup generated by the problem (5.1) (the existence of which is guaranteed by Theorem 2.1 of [40]). Then S(t) has a global attractor.
Proof. As usual we begin by showing the existence of an absorbing ball. Then, we show that the assumptions of Theorem 5.1 are satisfied, to conclude that a global attractor exists. We proceed formally, but calculations again can be completely justified by a Faedo-Galerkin approximation. We multiply (5.1) by Dv and proceed as above for the energy inequality analysis, along with Cauchy-Schwarz, and Young inequalities on the forcing term to obtain the bound d dt Using that D is positive and self-adjoint, along with assumptions on its boundedness, gives and then applying Poincaré inequality provides the bound Setting c 0 = νd0 2 min{α −2 , λ}, we get d dt and thus from the Gronwall inequality, lim sup

LUIGI C. BERSELLI, TAE-YEON KIM AND LEO G. REBHOLZ
Again applying the lower boundedness assumption on D, this reduces to which implies that S(t) : V 1 → V 1 , t ∈ R + has as absorbing ball B 1 := B(0, M 1 ) ⊂ V 1 , and we have the uniform estimate for t large enough. We now show the assumptions of Theorem 5.1 are satisfied, following analysis similar to [28]. Observe that S(t) can be represented as where Y (t) is the semigroup generated by the linear (which is linear since also D is a linear operator when the van Cittert deconvolution is used) problem y t + α 2 Ay t + νADy = 0, (5.2) with v solution of (5.1). We take the L 2 -inner product of (5.2) with y to find d dt y(t) 2 + α 2 ∇y(t) 2 + 2ν ∇D 1/2 y 2 = 0, and then using that D is positive and bounded from below, along with the Poincaré inequality, we obtain as in the previous calculations, d dt y(t) 2 + α 2 ∇y(t) 2 + c 0 λ y(t) 2 + α 2 ∇y(t) 2 ≤ 0, which immediately implies Therefore the semigroup Y (t) is exponentially contractive. By using the properties of D, and standard inequalities in Sobolev spaces (cf. [28,Eq. (3.11)]) it follows that P (Dv · ∇Dv) V −1/2 = sup for all positive T , and with bounds independent of T , but depending essentially on D .
Then we take (5.4)-(5.5) and we test it by A 1/2 z (procedure which can be again justified by approximation with Faedo-Galerkin method). We observe that z V s = (A s z, z) = A s/2 z 2 , by definition. Hence, by noting that P f − P (Dv · ∇Dv) ∈ L ∞ (0, T ; V −1/2 ) we get 1 2 Using previous estimates we obtain d dt and the bound on v V 1 ≤ M 1 implies that z is bounded in V 3/2 . This means that Z(t) maps V 1 into V 3/2 , and thanks to the compact embedding V 3/2 ⊂ V 1 , Z(t) is compact. Thus the assumptions of Theorem 5.1 are satisfied, and we conclude that S(t) is asymptotically compact. By [46], we conclude that S(t) has a global attractor in V 1 ; the rest of the proof follows as in [28]. With the same arguments one can also show that if the force is smoother, say in P f ∈ V 1 , then the attractor is a compact set in V 2 .
6. Pipe flows. In this section we perform some analytical computations with very classical tools, in order to focus on some special analytical features of the solutions of Voigt models and also to better understand the possible role of the RADM (and thus the NS-Voigt system for D 0 = I) in modeling of visco-elastic flows. We recall that prior to turbulence modeling, NS-Voigt equations appeared as a linear model for certain non-Newtonian fluids, see a series of papers by Oskolkov [37,38] In order to understand, at least qualitatively, the effect of the (non-viscous) damping term −α 2 ∆u t on the behavior of solutions, we restrict to the very simple setting of: a) fully developed flows, that is the velocity is directed only along the axis of the pipe; and b) we restrict to the NS-Voigt case; that is, assume D = D 0 = I. The deconvolution order N > 0 could be used, but details would change for each specific deconvolution operator, and also the analysis will become (even more) technical and also boundary conditions for filtering at the wall would need to be specified. Under the two above assumptions we develop analytical solutions, and even if we are aware that these solutions are not turbulent and describe only the linear behavior, they constitute a good source of data for debugging complex codes, as only very few analytical solutions are known to exist. Moreover, analytical solutions can be used to gain some physical insight into the model's solutions. On the other hand, we recall that a remarkable feature of the NS-Voigt equations is that this is one of the few LES models which can be fully treated with Dirichlet boundary conditions and with given external force, see [12,32]. Since the NS-Voigt equations can be treated in the presence of solid boundaries with Dirichlet boundary conditions, it is natural to investigate some class of exact solutions for NS-Voigt, as in the case of the NSE, to directly compare the results.
6.1. Remarks on the regularity of the solution. In this section we make some remarks on the solution of the NS Voigt equations, but similar reasoning can be transferred to the RADM (at least in the space periodic-setting, while here we have a problem with solid boundaries and zero Dirichlet boundary conditions). In particular, as remarked in Larios [32] the NS-Voigt have a pseudo-parabolic behavior, as the abstract equations considered in Carroll and Showalter [13]. One of the most interesting point is that the solutions of the NS-Voigt preserve the same regularity as the initial datum. Contrary to the NSE, for the solution of the NS-Voigt it is not possible to show an improvement of the regularity, especially for what concerns the space-variables. We recall that for the NSE (or more generally parabolic equations), as soon as a weak solution is strong, then it becomes smooth (even C ∞ if the force and the boundary of the domain are so) for all positive times. This has particular consequences on the fine properties of weak solution of the NSE, when obtained as limits of solutions to the NS-Voigt equations as α → 0 + , see [10].
Here, we show with a particular simple example that the nonlinear NS-Voigt cannot produce an improvement of the regularity, by making some very explicit computations on the spectrum of the Laplace operator. This can be also used to shed some light in the different behavior of the NS-Voigt with respect to the NSE, concerning some very special properties.
Let us consider the NS-Voigt system in a special physical situation. We have a channel with cross section E and with axis directed along the x 3 -direction. We consider an incompressible fluid, modeled by NS-Voigt equations in a semi-infinite straight pipe P = E × R + ⊂ R 3 . In a reference frame with z (we call z := x 3 for historical reasons) directed along the pipe axis and x := (x 1 , x 2 ) belonging to an orthogonal plane, the NS-Voigt equations read as where v(t, x, z) and p(t, x, z) respectively denote velocity and pressure, and ν > 0 represents kinematic viscosity.
Remark 6.1. In this subsection and later on also in Section 6.2 we use the symbol v(t, x, z) = (u(t, x, z), v(t, x, z), w(t, x, z)) to remember that the velocity is a vector field and that we will look for special solutions, where the axial velocity w(t, x, z) is the only non-zero one.
By following classical calculations dating back to Hagen and Poiseuille, we look for fully-developed solutions (also named Poiseuille-type solutions): v(t, x, z) = (0, 0, w(t, x)) and where p 0 (t) denotes an arbitrary function of time. The Poiseuille-type ansatz implies that the convective term cancels out identically, and this motivates the statement that the results will concern only the linear behavior of solutions. Moreover, it is easy to deduce from the equations that the pressure is p(t, z) = −λ(t) z. Finally, the dependence of w on the space variables x 1 and x 2 allows us to consider a problem reduced to the cross-section E: In the so-called "direct problem," without external force, a given pressure-drop is assigned and the problem is the following: Given λ(t), find w(t, x) such that x ∈ E, t ∈ R + , w(t, x) = 0 x ∈ ∂E, t ∈ R + , w(0, x) = w 0 (x) x ∈ E, where ∆ x denotes the Laplacian with respect only to the variables x 1 and x 2 . In this particular setting the nonlinear term is identically vanishing, but we still have a boundary value problem. We will start by making some explicit calculations in the case λ(t) is extremely smooth, say λ ≡ 1 and we assign an initial datum w 0 (x) ∈ H 1 0 (E). This corresponds to the classical pressure drop problem with a constant gradient superimposed to sustain the flow, but as we will see this assumption is completely inessential. We have the following theorem.
The lack of regularization in the NS-Voigt equation is due to the special role of the inviscid regularization and it has some consequences when the Voigt model is used for theoretical purposes, see for instance [10].
Proof of Thm. 6.1. First, we note that the existence and uniqueness of solution with initial datum w 0 and external force f = (0, 0, 0) follows as in [12]. Here we make a more explicit construction in order to highlight the regularity properties. The proof is based on an explicit analysis of the frequency budget, working with a spectral basis to construct the solution. (In the space-periodic case this can be done even more explicit with a Fourier series expansion, but here we deal with the boundary value problem).
Let {φ m } m∈N denote the eigenfunctions of the Laplace operator with vanishing Dirichlet boundary conditions on E, and by standard spectral theory the scalar functions {φ m } are smooth, orthonormal with respect to the L 2 (E) scalar product, and the eigenvalues {λ m } are an increasing sequence of strictly positive real numbers. Hence, we construct the solution as w(t, x) = m≥1 c m (t)φ m (x) and it is immediate to check that we have to solve the following ordinary differential equation, for all m ≥ 1 Observe that, from the regularity hypotheses on the initial datum, we have It is straightforward to check that 1≤m≤N c m (t)φ m (x) → w(t, x) as N → +∞, and from the explicit expression of the solution we can see that even if β m would be zero, then w(t, x) H s = ∞, for all s > 1, and for all t ≥ 0, due to the fact that the exponential e − λmν 1+α 2 λm t does not produce any smoothing (as if α = 0), since The situation will be different for what concerns time regularity, since any timedifferentiation will produce the zeroth order term λmν 1+α 2 λm which will not affect the convergence of the series describing the solution.
6.2. Pulsatile flows. In this section we consider pulsatile flows, which are driven by a time-periodic force (the pressure gradient). This is a remarkable case motivated also by recent studies in [3,7] for heart-beat driven, human physiological flows, including blood and cerebrospinal fluid.
In the so-called 'direct problem,' a pulsatile pressure-drop is assigned in problem (6.2). Together with this problem, in many cases of physiological interest (when the measure of the pressure cannot be done directly with sufficient accuracy by phase-contrast Magnetic Resonance Imaging or Doppler ultrasonography) it will be also meaningful to study a problem with assigned flux: E w(t, x) dx = Φ(t), for some smooth-enough given scalar function Φ(t). The Poiseuille-type ansatz will lead to the problem: Given Φ(t), find (w(t, x), λ(t)) such that which is an inverse problem and it is linked to one of the nowadays classical Leray problems.
The study of this problem when E is a circle, λ(t) is time-periodic, and α = 0 dates back to Richardson (1926) and Sexl (1930), while a big interest for applications to bio-fluids raised especially after the work of the physiologist J. R. Womersley [47]. In particular, he quantitatively clarified why for the NSE, if the pulsation is fast enough, then the qualitative behavior of solution can change drastically. Instead of a parabolic type profile (as in the stationary case) a new phenomenon called "Annular effect" appears: If the frequency of pulsation of λ is large enough (with respect to the other relevant parameters of the problem), then the maximum of the velocity is not located along the axis of the channel, but near the wall. This is particularly of interest in blood modeling, where -at first approximationwe have a pulsatile flow (blood driven by heart beat) in a straight channel (larger arteries). This motivated the analysis of Womersley, who explained by analytical formulae the experimental observation that in larger arteries of the rabbit and the dog there is also a reversal of the flow. This suggests that pulsatile flows, even in the linear range, can be different and much more complicated than the steady ones. The analysis of Womersley was based on a single mode pressure-drop λ(t) = e int and a circular cross-section E = {x 2 1 + x 2 2 < R}. In cylindrical coordinates (r, θ) and after the separation of variables w(t, x) = e i nt W (r), the system (6.2) becomes: find W (r) such that i n W (r) − ν W (r) + W (r) r = 1 for 0 < r < R, The above equation for a cylindrical symmetric solution can be solved by means of an expression involving the zeroth order Bessel function J 0 (r). The most widely known contribution 1 of Womersley is that of determining the non-dimensional quantity (named later Womersley number) which determines the qualitative behavior of solutions, where R is the radius of the channel, ω the frequency of the pressure drop, and ν the kinematic viscosity. The Womersley number is the ratio of the transient or oscillatory inertia force to the shear force. When Wo ≤ 1, the frequency of pulsations is sufficiently low that a parabolic velocity profile has time to develop during each cycle, and the flow will be very nearly in phase with the pressure gradient. In this case using instantaneous pressure to modulate a parabolic profile will give a reasonable solution. When Wo ≥ 10 the frequency of pulsations is sufficiently large that the velocity profile is relatively flat or plug-like, and the mean flow lags the pressure gradient by about 90 degrees.
In addition to this peculiar effect, the Womersley solution has become a paradigm in the analysis of biological fluids (see Fung [23] and Quarteroni and Formaggia [39]) for its simplicity, being one of the very few exact unsteady solutions to the NSE, to debug complicated CFD flow, and also to provide an improved source of initial/boundary data, as well as a benchmark solution for pulsatile flows.
We now perform an analysis of the NS-Voigt equations for this problem. Assume that Φ(t) = e i ωt , and by the ansatz w(t, x) = e i ωt W (r) in the case of (1.3) we arrive at the following boundary value problem for a single ODE: Find W (r) such that i ω W (r) − α 2 W (r) + W (r) r − ν W (r) + W (r) r = 1 for 0 < r < R, The above equation can be again solved exactly by means of an expression involving the zeroth order Bessel function J 0 (r). By explicit calculations, we get which reduces to the classical solution obtained by Womersley in the case α = 0. Hence, we get that the natural non-dimensional quantity is the following (which we call α-Womersley number) The α-Womersley number is always smaller than the Womersley number of the flow with the same viscosity, radius of the pipe, and frequency. In fact and the difference can be significant if α is not properly chosen. This gives the hint that a not finely tuned choice of α can completely change the behavior of solutions, since the flow can pass from one regime to another, if Wo is large and α is not small enough.
Here, we want to give a more explicit and real solution by considering the 2D case, that is, the cross section E is the 1D interval E =:] − R, R[, and showing the effect of a pulsatile force cos(ωt), which is not a peculiar effect of the space dimension. In particular, we will see that the qualitative behavior is not better in a 2D pipe, with respect to the 3D one. We consider then the following toy problem x ∈] − R, R[, t ∈ R + , w(t, x) = 0 x = ±R, t ∈ R + , which corresponds to the pulsatile flow in a rectangular 2D channel. In this case the exact solution (which is real as can be checked by direct inspection) has the following explicit expression: w(t, x) = sin(ω t) ω If we set, in the case α = 0 (i.e. the NSE), ν = 1, R = 1, and ω = 144, we find Wo=12, and the solution at t = 0 shows the annular effect and reversal of the flow at distance about R/2, as shown in Figure 4.  As is evidenced in Figure 6.2, this puts strong limitations on the choice of α, and moreover such a small α moves the NS-Voigt systems closer and closer to the NSE. Moreover, in LES, the parameter α has the role of the smallest resolved length scale, and when looking at the system in this sense, the limitations for this model are those of a good resolution of the flow, otherwise transient features are lost. This could be used to speculate that the space-filter alone cannot be very good when studying time-evolution problems with transient or periodic effects, or at the very least the smallest resolved length scale should not be constant across the domain. This idea is mentioned in [24], and indeed better results for turbulent channel flow were found when α was chosen of the order of the local mesh width, which was much finer in the boundary layer.
7. Conclusions. We have presented a detailed mathematical analysis of several important aspects of the RADM. In particular, we presented a detailed analysis of the model's treatment of energy, which showed that deconvolution acts to drain energy from the system, and likewise affects the energy spectrum. We derived a scaling for energy in its Fourier modes, which shows that increasing deconvolution in the model shortens the inertial range, and that the inertial range is split into two parts: The larger scales of the inertial range have a k −5/3 scaling with energy, while the lower numbers have a faster k −3 scaling. Interestingly, this is the same spectral scaling as both the Leray-α and simplified Bardina models. We presented energy spectra from computations of forced, isotropic, homogeneous turbulence with varying N and α, which verified the predicted scalings, and showed good agreement on large scales with DNS of the NSE. In addition to the RADM energy study, we also proved existence of smooth global attractors, and provided an analytical study of pulsatile flows that shows how small α may be requested to study certain laminar flows.
There are several important directions for future work with the RADM. First and foremost, more benchmark testing needs done. In particular, since the model performed well with Re τ = 180 turbulent channel flow, future testing on the benchmarks of Re τ =395 and 590 are natural next steps. Turbulent flow around obstacles is another important test problem, as is applying the model to specific application problems such as flow through an aorta. The model could also be studied with different types of filtering and deconvolution.