Closing the loop: nonlinear Taylor vortex flow through the lens of resolvent analysis
CClosing the loop: nonlinear Taylor vortex flowthrough the lens of resolvent analysis
Benedikt Barthel , Xiaojue Zhu , and Beverley J. McKeon Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, CA 91125, USA Center of Mathematical Sciences and Applications, and School of Engineering and Applied Sciences,Harvard University, Cambridge, MA 02138, USA Graduate Aerospace Laboratories, California Institute of Technology, Pasadena, CA 91125, USA
ABSTRACT
We present an optimization-based method to efficiently calculate accurate nonlinear models of Taylor vortex flow. We use theresolvent formulation of McKeon and Sharma (2010) to model these Taylor vortex solutions by treating the nonlinearity not asan inherent part of the governing equations but rather as a triadic constraint which must be satisfied by the model solution.We exploit the low rank linear dynamics of the system to calculate an efficient basis for our solution, the coefficients of whichare then calculated through an optimization problem where the cost function to be minimized is the triadic consistency ofthe solution with itself as well as with the input mean flow. Our approach constitutes, what is to the best of our knowledge,the first fully nonlinear and self-sustaining, resolvent-based model described in the literature. We compare our results todirect numerical simulation of Taylor Couette flow at up to five times the critical Reynolds number, and show that our modelaccurately captures the structure of the flow. Additionally, we find that as the Reynolds number increases the flow undergoesa fundamental transition from a classical weakly nonlinear regime, where the forcing cascade is strictly down scale, to a fullynonlinear regime characterized by the emergence of an inverse (up scale) forcing cascade. Triadic contributions from theinverse and traditional cascade destructively interfere implying that the accurate modeling of a certain Fourier mode requiresknowledge of its immediate harmonic and sub-harmonic. We show analytically that this finding is a direct consequence of thestructure of the quadratic nonlinearity of the governing equations formulated in Fourier space. Finally, we show that using ourmodel solution as an initial condition to a higher Reynolds number DNS significantly reduces the time to convergence.
Keywords: Taylor-Couette flow, low-dimensional models, nonlinear instability, transition to turbulence
Taylor-Couette flow (TCF), the flow between two concentric and independently rotating cylinders, is one of the canonicalproblems in fluid mechanics, and a paradigm for the study of linear stability, pattern formation, and rotationally driventurbulence. From the original investigations of Taylor (1923) to the pioneering experiments of Coles (1965), recent theoreticalanalyses (Gebhardt and Grossmann, 1993; Jones, 1981; Maretzke et al., 2014), high Reynolds number simulations (Ostillaet al., 2013; Ostilla-M´onico et al., 2014; Grossmann et al., 2016; Sacco et al., 2019), and experimental investigation (van Gilset al., 2011; Huisman et al., 2014; van Gils et al., 2012), TCF has remained a problem of interest for most of the last century.Perhaps the most well known characteristic of TCF is the incredibly rich array of stable flow states that exist over a range ofgeometries and relative rotation rates of the inner and outer cylinders (Coles, 1965; Andereck et al., 1986). We consider thecase of pure inner cylinder rotation, for which the problem is parameterized by the ratio of the inner and outer radii, η ≡ r i / r o ,and a single Reynolds number: R . In this case the laminar velocity profile becomes linearly unstable at a critical Reynoldsnumber: R c ∼ O ( ) . A centrifugal instability leads to the formation of a periodic array of toroidal vortex structures knownas Taylor vortices or Taylor rolls. For a given geometry, these steady, axisymmetric Taylor vortices are stable and exist forsome range of Reynolds number in what is known as the Taylor Vortex flow (TVF) regime. As the rotation rate of the innercylinder is increased further, the Taylor vortices experience a secondary instability giving rise to azimuthally traveling waveswhose phase speed is determined by the geometry but whose azimuthal periodicity is not unique (Coles, 1965). This regimeis known as Wavy Vortex flow (WVF) and is characterized by being time periodic in a stationary reference frame but steadyin a frame corotating with the traveling wave (Marcus, 1984). WVF is again stable for some range of Reynolds numbersbefore the traveling waves themselves become unstable, and a second temporal frequency arises causing the traveling wavesto become modulated in space and time in what is known as modulated wavy vortex flow (MWVF). As the driving of theinner cylinder is increased further still, the flow becomes disordered and begins to transition to turbulence. However, whilethe main sequence of transitions from laminar flow to the bifrucation to TVF at R c ≈ a r X i v : . [ phy s i c s . f l u - dyn ] F e b n to turbulence at R ≈ R ∼ O ( ) (Grossmann et al., 2016).In recent years there has been renewed interest (Dessup et al., 2018; Sacco et al., 2019) in the dynamics of TVF and WVFas a model system to study the self-sustaining process (SSP) proposed by Waleffe (1997). The SSP consists of streamwise rollswhich advect the mean shear giving rise to streaks of streamwise velocity, which become unstable to wave-like disturbances,which in turn nonlinearly interact to sustain the rolls (Waleffe, 1997; Hamilton et al., 1995). Dessup et al. (2018) performeddirect numerical simulations (DNS) to show that the mechanism of transition from TVF to WVF follows the same path asdescribed in the SSP: the traveling waves of WVF arise due to an instability of the streamwise velocity component (streaks) ofTVF with the cross-stream velocity (rolls) playing a negligible role in the instability mechanism. Sacco et al. (2019) extendedthis line of study to higher Reynolds numbers and showed that despite their origin as a centrifugal instability, turbulent Taylorvortices are preserved in the limit of vanishing curvature and are thus not dependent on a rotational effects, and rather aresustained through a nonlinear feedback loop between the rolls and streaks.The SSP is believed to be one of the building blocks of turbulence and thus the study of self sustaining solutions of theNavier-Stokes equations (NSE), known as exact coherent states (ECS) has been of great interest to researchers since theywere first discovered by Nagata (1990). The field has grown immensely since in the intervening years and we make no attemptto summarize it all here in the interests of brevity. Of primary interest is the observation that many of these solutions resemblestreamwise elongated vortices and streaks, and thus resemble structures observed in experiments and simulations of turbulentflows (Beaume et al., 2015). This has led to the idea that ECS make up the phase space skeleton of turbulence, and that theobservation of these structures indicates the turbulent trajectory passing by one of these ECS solutions.Due to the significant mathematical simplification, it is useful to model these elongated structures as being infinitelylong in their streamwise extent. Illingworth (2020) studied the linear amplification mechanism of such streamwise invariantstructures to identify the most amplified spanwise length scales in both channel and plane Couette flow (PCF), and found thatthe latter was far more efficient in amplifying these structures. However no structures actually observed in channel or planeCouette flow are truly invariant in the streamwise direction. Additionally, ECS are generally unstable; while the turbulenttrajectory may visit these states, they do not actually persist in nature. TVF is thus a valuable test case to study the nonlineardynamics sustaining ECS in general since it is in fact a stable solution observed in experiment and due to the cylindricalgeometry is exactly streamwise constant.We aim to help bridge the gap between linear stability theory, which accurately predicts the genesis of Taylor vortices, andhigher Reynolds numbers where Taylor vortices are sustained by fully nonlinear mechanisms as shown by Sacco et al. (2019).The resolvent formulation of the NSE, which interprets the nonlinear term in the NSE as a forcing to the linear dynamics,offers a natural path from linear to nonlinear analysis and as such will form the basis of our modeling efforts. Resolventanalysis has historically been successful in exploiting the linear dynamics to identify structures in turbulence, which is bynature a nonlinear phenomenon (McKeon and Sharma, 2010). However, recently attempts have been made to explicitlycharacterize and quantify the influence of nonlinear dynamics within the resolvent framework. For example, Moarref et al.(2014) and McMullen et al. (2020) used convex optimization to compute reduced order representations of turbulent statisticsand Rigas et al. (2021) studied solutions of the Harmonic-Balanced Navier-Stokes equations to identify optimal nonlinearmechanisms leading to boundary layer transition. Additionally Morra et al. (2021) and Nogueira et al. (2021) have directlycomputed the nonlinear forcing statistics for minimal channel and Couette flow respectively and analyzed the efficacy of lowrank resolvent reconstructions in capturing the relevant dynamics.In this work we use optimization-based methods not only to capture statistics, but to explicitly model the self-sustainingnonlinear system. To the best of these authors’ knowledge this work is the first example in the literature of “ closing theresolvent loop ” where the feedback loop formulation of the NSE introduced by McKeon and Sharma (2010) is used to generatesolutions to the fully nonlinear NSE. In §2 we outline the resolvent formulation of the NSE and our optimization-based model.In §3 we present the results of our model, compare them to DNS, and present a detailed analysis of the nonlinear dynamics.We discuss applications of our results to DNS of turbulent TCF in §4 and conclude with some final remarks in §5. We study the flow of an incompressible Newtonian fluid with kinematic viscosity ν between two concentric cylinders usingthe NSE in cylindrical coordinates, ∂ ˜ u ∂ t + ˜ u · ∇ ˜ u − R ∇ ˜ u − ∇ ˜ p = ∇ · ˜ u = N r N θ N z DOF ( N r × N θ × N z )
100 101 768 512 39,714,816200 101 768 512 39,714,816400 129 768 512 50,724,864650 129 768 512 50,724,8641000 193 1024 640 126,484,4802000 193 1024 640 126,484,480
Table 1.
Numerical details of DNSon the domain r ∈ [ r i , r o ] , θ ∈ [ , π ] , z ∈ ( − L z / , L z / ) . We will consider the case where the outer cylinder is held fixedwhile the inner cylinder rotates with a prescribed azimuthal speed U i . The equations are nondimensionalized using the gapwidth d ≡ r o − r i and the azimuthal velocity of the inner cylinder U i . The Reynolds number is defined as R ≡ U i d / ν . In thesenondimensional variables the limits of the radial domain are given as a function of the radius ratio η by r i = η / ( − η ) and r o = / ( − η ) . Throughout this work we fix η = . η was chosen to allow for comparison to past studies such asOstilla et al. (2013) and since it allows for a larger range of Reynolds numbers for which TVF is a stable solution of (1).We decompose the state [ ˜ u , ˜ p ] = [ ˜ u θ , ˜ u r , ˜ u z , p ] into a mean and fluctuating component, [ ˜ u ( r , θ , z , t ) , ˜ p ( r , θ , z , t )] = [ U , P ( r )] + [ u ( r , θ , z , t ) , p ( r , θ , z , t )] (3)with ( · ) ≡ lim T , L → ∞ π T L (cid:90) π (cid:90) L (cid:90) T ( · ) d θ dz dt , (4)which upon substitution into(1) and averaging over θ , z , and t results in the mean momentum equation U · ∇ U − R ∇ U − ∇ P = − u · ∇ u (5) ∇ · U = . (6)Subtracting (5) from the full NSE then results in a governing equation for the fluctuations, where we have grouped thoseterms which are nonlinear in the fluctuations on the right hand sides in anticipation of the following analysis. ∂ u ∂ t + U · ∇ u + u · ∇ U − R ∇ u − ∇ p = − (cid:16) u · ∇ u − u · ∇ u (cid:17) (7) ∇ · u = U | r i = ˆ e θ , U | r o = u | r i = u | r o = (9) In order to validate our model solution we perform DNS of TCF for a range of Reynolds number, 100 < R < η = .
714 and an aspect ratio L z / d =
12. However, our analysis is focused primarily on the cases R = , , and 400. The details of the numerical method can be found in Verzicco and Orlandi (1996); van der Poel et al. (2015); Zhuet al. (2018), and the details of the simulations performed in this work are summarized in table 1. Our model solution will be based on the resolvent formulation of McKeon and Sharma (2010), which assumes that the onedimensional mean velocity profile U ( r ) is known. This allows us to write (7) and (8) as a balance between the linear dynamicsand the nonlinear term which we group into a forcing term denoted by f . L u = N ( u , u ) ≡ f . (10) e Fourier transform both the velocity u ( r , z , θ , t ) and the nonlinear forcing f ( r , z , θ , t ) in time as well as the homogeneousspatial directions, θ and z . For a function q ( r , z , θ , t ) the Fourier transform is defined asˆ q ( r , k , n , ω ) ≡ (cid:90) ∞ − ∞ (cid:90) ∞∞ (cid:90) π q ( r , z , θ , t ) e − i ( k z z + n θ − ω t ) d θ dz dt . (11)This results in system of coupled ordinary differential equations (ODEs) for the Fourier modes ˆ u k ( r ) and ˆ f k ( r ) parameterizedby the wavenumber triplet k ≡ [ k z , n , ω ] . The domain is periodic in the azimuthal direction and we formally consider the case L z = ∞ , so we have n ∈ Z and k z , ω ∈ R . ( L k − i ω M ) ˆ u k = ˆ f k (12)The explicit expressions for the Fourier transformed linear operator L k and the weight matrix M are given in appendix A. Ifwe then invert the linear operator on the left hand side of (12), we arrive at the characteristic input-output representation ofthe NSE of McKeon and Sharma (2010) where the linear dynamics represent a transfer function from the nonlinearity, whichis interpreted as a forcing, and the velocity.ˆ u k = ( L k − i ω M ) − ˆ f k (13)A singular value decomposition (SVD) of this transfer function which we call the resolvent and denote by H provides anorthonormal basis for the velocity as well as the forcing, where superscript H denotes the conjugate transpose. H k ≡ ( L k − i ω M ) − = Ψ k Σ k Φ H k (14)The columns of Ψ k and Φ k are denoted by ψ k , j and φ k , j are referred to as the resolvent response and forcing modesrespectively and are ordered by their singular values σ k , j such that σ k , ≥ σ k , ≥ ... ≥ σ k , j ≥ σ k , j + . Thus ψ k , represents themost linearly amplified structure at that wavenumber. Note that these modes are vector fields over r which are orthonormalwith respect to an L inner product over the three velocity components (cid:104) a , b (cid:105) ≡ (cid:90) r o r i a ∗ m ( r ) b m ( r ) r dr (15)with associated norm (cid:107) a (cid:107) ≡ (cid:104) a , a (cid:105) / (16)where summation over m is implied, such that (cid:104) ψ k , i , ψ k , j (cid:105) = (cid:104) φ k , i , φ k , j (cid:105) = δ i j . (17)In this basis, each Fourier mode of the velocity and forcing may be written asˆ u k = ∞ ∑ j = σ k , j χ k , j ψ k , j (18)ˆ f k = ∞ ∑ j = χ k , j φ k , j (19)where χ k , j ≡ (cid:104) φ k , j , ˆ f k (cid:105) represents the projection of the (unknown) forcing onto the forcing modes. The various flow states observed in TCF (TVF, WVF, and MWVF) may be defined by their spatio-temporal symmetries(Rand, 1982). We define TVF, the focus of this study, as a solution to (1) which is steady, axisymmetric, and axially periodicwith fundamental wavenumber β z , meaning we restrict ourselves to wavenumber vectors of the form k = [ k β z , , ] where k ∈ Z . This fundamental wavenumber β z is related to the axial height of the Taylor vortices and is generally constrained bythe experimental apparatus or computational box since the domain must contain an integer number of vortices. The resolventformulation assumes an infinite axial domain so the choice of β z is not immediately obvious. However, we found that the esults shown in this work are robust to changes in β z as long as π / (cid:46) β z (cid:46) π /
3. Therefore, we choose the axial periodicityof our model to match that observed in our DNS, allowing for a direct comparison between our model and the DNS. Thespecific values of β z are listed in table 2. Given these symmetries our model solution will consist of an expansion in Fouriermodes u ( r , z ) = N k ∑ k = ˆ u k ( r ) e ik β z z (20)where each Fourier mode ˆ u k is itself an expansion in resolvent modes given by (18). We truncate the model at N k Fouriermodes each of which is expanded in N kSVD resolvent modes such that the final form of the TVF solution is given by u ( r , z ) = N k ∑ k = N kSVD ∑ j = σ k , j χ k , j ψ k , j ( r ) e ik β z z . (21) At any given wavenumber, the forcing ˆ f k is given by a convolution of the interactions of all triadically compatible velocitymodes.ˆ f k = ∑ k (cid:48) (cid:54) = − ∇ · ( ˆ u k (cid:48) ˆ u k − k (cid:48) ) (22)This means that the forcing at a given wave number k contains only interactions between Fourier modes whose wavenumbersadd up to k . Throughout this work we use the terminology “k = k + k ” to refer to a single (resonant) triad involving thenonlinear interaction between Fourier modes with wavenumbers k and k forcing the Fourier mode with wavenumber k .Equating the two expressions for the forcing mode given by (19) and (22) and substituting (18) for the velocity modesgives ∞ ∑ j = χ k , j φ k , j = ∑ k (cid:48) (cid:54) = ∞ ∑ p = ∞ ∑ q = − σ k (cid:48) , p χ k (cid:48) , p σ k − k (cid:48) , q χ k − k (cid:48) , q ∇ · (cid:16) ψ k (cid:48) , p ψ k − k (cid:48) , q (cid:17) . (23)Projecting both sides of (23) onto each forcing φ k , i and dropping the summation symbols for simplicity gives χ k , i = χ k (cid:48) , p χ k − k (cid:48) , q N kk (cid:48) , ipq (24)where the N kk (cid:48) , ipq are called the interaction coefficients and are given by N kk (cid:48) , ipq = − σ k (cid:48) , p σ k − k (cid:48) , q (cid:104) φ k , i , ∇ · (cid:16) ψ k (cid:48) , p ψ k − k (cid:48) , q (cid:17) (cid:105) (25)which, critically, can be computed solely from knowledge of the linear operator H .Nonlinear interactions between the velocity fluctuations also appear in the divergence of the Reynolds stress on the righthand side of (5). This term is referred to as the “mean forcing” and is given by the sum of nonlinear interactions of all the ˆ u k and their complex conjugates ˆ u − k , which can be directly interpreted as (22) evaluated at k = : − u · ∇ u = ∑ k (cid:48) (cid:54) = − ∇ · ( ˆ u k (cid:48) ˆ u − k (cid:48) ) ≡ ˆ f . (26)At this point we would like to reiterate that the mean velocity profile is assumed to be known a priori. Thus the left hand sideof (5), and therefore the Reynolds stress divergence on left hand side of (26), is also known.We have thus reduced the NSE (under the assumption of a known mean velocity) to the infinite system of coupledpolynomial equations (24) for the complex coefficients χ k , j with the auxiliary condition that (26) is satisfied. While derivingan exact (nontrivial) solution to (23) may be a daunting task, we will demonstrate that approximate solutions can be efficientlycomputed by minimizing the residuals associated with (24) and (26). .5 Optimization Problem We have recast the NSE in the language of resolvent analysis as χ k , j − χ k (cid:48) , m χ k − k (cid:48) , n N kk (cid:48) , jmn = , ∀ k , j (27) u · ∇ u − f , k (cid:48) , mn χ k (cid:48) , m χ ∗ k (cid:48) , n = , (28) f , k (cid:48) , mn ≡ σ k (cid:48) , m σ k (cid:48) , n ∇ · (cid:16) ˆ ψ k (cid:48) , m ˆ ψ ∗ k (cid:48) , n (cid:17) , (29)where we have expanded the velocity Fourier modes in their resolvent basis according to (18), and summation over m , n and k (cid:48) is implied. We truncate the expansion at some number of harmonics, N k , of the fundamental wavenumber, and ateach retained harmonic we truncate the singular mode expansion at N kSVD such that the total number of retained modes is N = ∑ N k k = N kSVD . We seek to minimize the residuals (in the sense of the L norm) associated with (27) and (28), and thusformulate the following optimization problem:min χ k , j g ( χ k , j ) = ag ( χ k , j ) + ( − a ) g triad ( χ k , j ) . (30)The first and second terms on the right hand side in (30) are defined as the mean constraint, g ( χ k , j ) ≡ (cid:107) u · ∇ u − f , k (cid:48) , mn χ k (cid:48) , m χ ∗ k (cid:48) , n (cid:107)(cid:107) u · ∇ u (cid:107) , (31)and the triadic constraint, g triad ( χ k , j ) ≡ | χ k , j − χ k (cid:48) , m χ k − k (cid:48) , n N kk (cid:48) , jmn | . (32)The former represents the residual in the mean momentum equation (5), while the latter represents the residual in the equationfor the fluctuations (7). The user defined weighting parameter a ∈ ( , ) determines the relative penalization of each of thesetwo constraints in the residual.At this point we would like to highlight several important aspects of problem (30). First, we reiterate that the left handterm in the mean constraint (31) is a known function since the mean velocity profile is assumed to be known a priori. Second,we emphasize that we have assumed no closure model and made no modeling assumptions regarding the form of the nonlinearforcing in the derivation of (30). Finally, while in general the amplitudes χ k , j ∈ C , for the special case of steady axisymmetricsolutions considered here it can be shown that χ k , j ∈ R ∀ j , k meaning the optimization need only be carried out over a realvalued domain.Finally, we note that while the reformulation of the NSE (7) in the resolvent framework (27) is reminiscent of a Galerkinmethod (GM) where the governing equations are projected onto some predetermined set of basis functions, the currentapproach is appreciably different. Since we consider a steady process we can not integrate the equations forward in time aswould be generally done in a GM. Furthermore, since the resolvent framework provides a basis for both the velocity and thenonlinearity, (27) is an exact representation of (7) whereas a GM represents the governing equations only in an integral sense. We solve the optimization problem (30) using a trust-region algorithm built into Matlab’s f minunc function for R = , , and 400. The gradient and Hessian of (30) may be derived explicitly and are input to the algorithm to improve accuracy. Theweighting parameter a in (30) is set to 0.01 which means that the triadic constraint is penalized 99 times more heavily than themean constraint. This reflects the observation that the triadic constraint which encodes the fully nonlinear governing equationfor the fluctuations is far more complex than the mean constraint which, given the fact that the mean profile is known, issimply a least squares fit to a curve.While this value of a was found to lead to the most consistent and accurate results, the results are qualitatively robust tochanges in a as long as 0 . (cid:46) a (cid:46) .
8. If a (cid:46) . a (cid:38) .
8, the mean constraint isweighted too heavily and the optimization tends to over fit to the input mean and converge to a local minimum which is notself-sustaining and inaccurately reconstructs the velocity field.
100 200 400 β z N k N kDNS ∀ k ( k ≤ ) , ( k > ) ( k ≤ ) , ( k > ) N
48 80 138 g ∗ . × − . × − . × − e dns . × − . × − . × − e . × − . × − . × − e nse . × − . × − . × − Table 2.
Reynolds numbers, fundamental wavenumber β z , truncation values, degrees of freedom: N , final residuals(normalized by initial residual): g ∗ , and error metrics for the three model solutions presented.The optimization also requires an initial guess. For this we solve the rank 1 formulation of (30) using just one wavenumber,the fundamental, and one resolvent mode, in which case the minimum can be found analytically, resulting in an amplitude χ rank ≈ .
13. We then initialize the full optimization such that χ = χ rank and the remaining χ j > are assigned randomvalues between -0.01 and 0.01.We assess the convergence and accuracy of our model solution using three metrics. First, we compute the final minimumresidual of the cost function in (30) denoted g ∗ . Second, we compute the error of the model solution compared to the temporalaverage of the DNS solution e dns ≡ N k ∑ k = (cid:107) ˆ u k − ˆ u k , dns (cid:107) . (33)The norm is defined in (16). Third we quantify the error of our model solution in solving the underlying governing equations(7), the details of which are discussed in §3.2.We note that the error metric comparing our model to the DNS should be viewed with some caution since the Fourierdecomposition of the DNS involves some inherent uncertainty. While our model is formulated in Fourier space, the DNS towhich we compare our model utilizes a finite difference method in the axial direction. This means that the five or six Taylorvortices in the computational domain are not necessarily exactly the same size and the fundamental wavenumber β z can onlybe defined in an average sense, β z = π n roll L z (34)where n roll is the number of Taylor vortices contained in the domain and L z is the axial domain size. We use this average β z inthe construction of our model. To compute the DNS Fourier modes used in (33) we extract a single Taylor vortex whose sizeis closest to the average and perform a Fourier decomposition on this reduced domain. The largest difference between thisbest fit wavenumber and the average we observed was 0 . R = N k = N kSVD =
22 for k ≤ N kSVD =
10 for k > g ∗ with increased N . A detailed discussion of the choice of these particular truncation values is presented in§3.6. These truncation values, the total degrees of freedom, the axial wavenumber β z , as well as the error metrics for all threeReynolds numbers are summarized in table 2. The final result of the model is shown in figures 1 and 2 where we compare the model result and the DNS. We plot themean-subtracted azimuthal velocity u θ , and the azimuthal vorticity ω θ = ∂ u r ∂ z − ∂ u z ∂ r for R = ,
200 and 400. The modelsolution is axisymmetric and steady by construction, and thus the radial and axial velocity are linked through continuity;no information is omitted by plotting ω θ . As a comparison we show the azimuthal average of the mean subtracted DNS,however at this Reynolds number the flow is axisymmetric and steady, so the average field shown is representative of the ow at any azimuthal location and at any instance in time. There is good agreement between the resolvent model (top rows)and the DNS (bottom rows). The model accurately captures the dominant structure of the flow including the strong plumesof azimuthal velocity. The azimuthal vorticity exhibits a checkerboard pattern of regions of roughly constant vorticity ofopposing signs. Regions of higher vorticity are concentrated near the walls, while the larger segments in the bulk of thedomain have comparatively lower levels of vorticity. These results are in agreement with the DNS of Sacco et al. (2019), whofound that as the Reynolds number increases this concentration of vorticity at the walls is enhanced and the bulk becomesincreasingly “empty” of vorticity.A more quantitative assessment of the model’s accuracy is shown in figure 3 where we compare the individual Fouriermodes of the model solution to the Fourier modes computed from the DNS. For clarity of presentation we focus on R = ( k ≤ ) , However an analysis of the accuracy of all retained Fourier modes for all Reynolds numbers ispresented in figure 4 in §3.2. Compared to the DNS, the model slightly over-predicts the amplitude of the radial velocity forthe fundamental Fourier mode, but the wall parallel components of the fundamental are captured almost exactly. Most strikingis the good agreement of the higher harmonics. The largest scale dominates the contribution to the Reynolds stress divergenceand is thus determined primarily by the mean constraint, which as mentioned previously is relatively “easy” to solve. Howeverthe smaller scales require accurately approximating the solutions to the nonlinear triadic constraint, a much less trivial task.Furthermore, small deviations between the Fourier modes of the model and DNS are not necessarily indicative of errors inour model. For the reasons discussed above, a Fourier decomposition of the DNS incurs some inherent uncertainty. A morerigorous assessment of how accurately our model solves the governing equations is presented in §3.2. Additionally, figure 3acompares the DNS mean velocity profile, used as an input to the model, to the mean velocity profile computed from Reynoldsstress divergence of the model itself.1 R ∇ ¯ U model = N k ∑ = − ∇ · ( ˆ u k ˆ u − k ) (35)The input and output mean velocity profile show very good agreement, with only some mild discrepancy at the edge ofthe inner boundary layer. The error in mean velocity may be computed as e ≡ (cid:107) ¯ U DNS − ¯ U model (cid:107) (36)which is associated with residual of the mean constraint (31) in (30). However, note that while (36) is written in terms of themean velocity, (31) is written in terms of the Reynolds stress divergence. The values of e for all three Reynolds numbers aretabulated in table 2. We generally do not use (36) as one of the measures of convergence since very few modes are requiredto accurately capture the mean, and thus e reaches a minimum long before the full nonlinear flow is converged. This isconsistent with past studies which have have shown that the mean velocity profile of various flows may be accurately modeledusing the Reynolds stress divergence of a single resolvent or eigen mode (Mantiˇc-Lugo et al., 2014, 2015; Rosenberg andMcKeon, 2019a).Overall, the success of the model in capturing this fully developed TVF indicates that, despite its fully nonlinear nature,the full solution remains relatively low dimensional. Note that our model constitutes approximately a million-fold reductionin degrees of freedom compared to the DNS. Nevertheless, given the relative simplicity of the flow, the model reduction is notas drastic as one might expect from an analysis purely of the energetic content of the flow. At R =
400 the velocity associatedwith the third harmonic ( k = ) is two orders of magnitude less than the fundamental, and yet nine wavenumbers must beretained in order to achieve the convergence shown here. The dynamic importance of these energetically weak harmonics isdiscussed in §3.2. We have shown that our model accurately captures the structure of the TVF observed in the DNS. Now we analyze theaccuracy of our model viewed from the perspective of a self sustaining process. In other words we assess how accurately ourmodel approximates a solution to the governing equations (7). In the resolvent framework, the nonlinear term is interpreted asa forcing to the linear dynamics. As such, a solution is self-sustaining if the sum of all triadic interactions at a particularwavenumber provide the correct forcing for the response at that wavenumber. This means that we must haveˆ u k = H k ∑ k (cid:48) (cid:54) = − ∇ · ( ˆ u k (cid:48) ˆ u k − k (cid:48) ) ∀ k . (37)Note that (37) is simply a restatement of (13) with the nonlinear forcing written explicitly in terms of ˆ u k and that our modelwill generally not satisfy (37) exactly. igure 1. Mean subtracted azimuthal velocity u θ computed from our model (top row) and DNS (bottom row) at (from leftto right) R = , igure 2. Azimuthal vorticity ω θ = ∂ u r ∂ z − ∂ u z ∂ r computed from our model (top row) and DNS (bottom row) at (from left toright) R = , Figure 3.
Model solution (lines) compared to the DNS (symbols) at R = U , computed fromReynolds stress divergence of model compared to input mean velocity from DNS ( a ). First five Fourier modes of modelsolution, ˆ u θ , ˆ u r , ˆ u z ( b-d ), k = k = k = k = ere we refer to the the direct result of our model, i.e. the quantity on the left hand side of (37), as the “primary” velocity,and we denote right hand side of (37), computed from that model solution, as the “forced” velocity. In figure 4 we plot theazimuthal component of both the primary and forced Fourier modes for all the wavenumbers and for all three Reynoldsnumbers. For all Reynolds numbers and wavenumbers agreement between the primary (open circles) and forced mode (lines)is very good indicating that our model is indeed a close approximation of a solution to (7). Figure 4 also shows the Fouriermodes computed from the DNS for comparison (open squares). We see that there is growing discrepancy between the modelresult and the DNS modes with increasing wavenumber. However note that the discrepancy between the model and the DNS,which is only present in the higher harmonics, is two orders of magnitude smaller than the amplitude of the fundamental.This discrepancy is due to the structure of the nonlinear forcing and is discussed further in §3.3. Additionally, we note that forboth the model and DNS the mode shapes of the higher harmonics, k (cid:38) k indicating some level of universality as the length scale decreases.We quantify the total error in nonlinear compatibility as e nse = N k ∑ k = (cid:107) ˆ u k − H k ∑ k (cid:48) (cid:54) = − ∇ · ( ˆ u k (cid:48) ˆ u k − k (cid:48) ) (cid:107) . (38)This may be thought of as the total residual associated with how accurately our model solution approximates a solution to (7),or in other words it represents the final residual associated with the triadic constraint (32). For all three Reynolds numbersconsidered the error is O ( − ) , with the exact values listed in table 2. Here we investigate which individual triadic interactions are most important in sustaining the flow and how these vary withReynolds number. As previously noted, for R = , k (cid:38)
5, do not contribute significantly to theenergy content of the flow. Nevertheless, they play a crucial role in the nonlinear forcing of the larger structures and arenecessary to achieve the convergence shown in figures 1 and 2. To identify the mechanics underlying the forcing structure wecompute the individual terms in the summation on the right hand side of (37). These individual contributions of the forcedvelocity, defined asˆ v k , k (cid:48) ≡ − H k ∇ · ( ˆ u k (cid:48) ˆ u k − k (cid:48) ) , (39)represent the contribution of each individual triadic interaction in (37) and are shown in figures 5 and 6 for R =
100 and R =
400 respectively. The individual contributions (39) are plotted in dashed colored lines and the full Fourier mode is plottedin solid black. By definition the individual contributions (dashed colors) sum to the full Fourier mode (solid black). To clarifythe following discussion we define a “(forward) forcing cascade” as the forcing of mode k by interactions involving strictlymodes k ≤ k and an “inverse forcing cascade” as the forcing of mode k by interactions involving k > k .At R = ≈ . R c , we observe a forcing mechanism reminiscent of a weakly nonlinear theory where the harmonicsare all driven exclusively by a forward forcing cascade. The k = k = k = k = k = k do not contribute to the forcing of modes with wavenumber k < k .The higher Reynolds number model solutions do not exhibit the same unidirectional forcing cascade observed close to thebifurcation from laminar flow. We plot the same decomposition of (39) for R = ≈ R c in figure 6. Due to the significantdisparity in amplitude between the individual triadic contributions and the composite mode we plot the two against separatey-axis. For the harmonics ( k > ) , the pair of contributions due to triadic interactions involving the fundamental ( k = ± k = ( k + ) − k = ( k − ) +
1, have large, almost equal amplitudes but are of opposite sign and almost cancel. The samephenomenon is observed for the triads involving the k = ± Γ k , k − k (cid:48) , k (cid:48) ≡ (cid:104) ˆ v k , k (cid:48) , ˆ u k (cid:105)(cid:104) ˆ u k , ˆ u k (cid:105) . (40)These projections are plotted in figures 7 and 8 for R = R =
400 respectively. Each sub-figure corresponds to oneFourier mode, k , and each tile represents the contribution to the k ’th Fourier mode from the triadic interaction between k and k . Note that the k ’th column includes contributions from ± k since ˆ u − k = ˆ u ∗ k and therefore both represent the samestructure in physical space. As expected from figure 6, we observe pairs of strong negative and positive correlations from thetwo triads involving the fundamental, with less pronounced, but still evident, pairing between triads involving k = ±
2. In rder to quantify the total contribution of these pairs of triads we sum the columns in figures 7, 8, and the equivalent casefor R = k , it is the two pairs of triads k = ( k + ) − k = ( k − ) + k = ( k + ) − k = ( k − ) + R =
200 are driven almost entirely by the pair of triads involving k = ± k = ± k ≥
4. At R =
400 the forcing is more distributed among the various triadic interactionsindicating a higher degree of nonlinearity. However, the triads involving k = ± k = ± k = ± k = ±
1. Nevertheless, it is the triads involving the fundamental which have the largest amplitude contributions anddisplay the largest degree of destructive interference.At this point we would like to revisit the discrepancy between the Fourier modes of the model solution and the DNSobserved in figure 4. Recall that due to the structure of the forcing, an accurate reconstruction of a particular Fourier modewith wavenumber k requires accurate knowledge of its harmonic, k +
1. Practically, the model must be truncated at somepoint and so there will always be a maximum wavenumber k m whose harmonic k m + u k m since ˆ u k m + is unavailable to participate in the inverse forcing cascade described above.This error will then “back propagate” through Fourier space until it is outweighed by the influence of the large scales and theconstraint imposed on those large scales by the input mean flow. A detailed analysis of how the Fourier space truncationaffects the convergence of the optimization is beyond the scope of this work, however it is interesting to note that, while thehigher harmonics of our model solution deviate slightly from the DNS, they remain nonlinearly compatible to a very goodapproximation. In other words, the primary and forced Fourier modes of the model solution in figure 4 agree very well asquantified by the small residuals as defined by (38) and listed in table 2. Many studies have approached the nonlinear modeling of TVF through weakly nonlinear (WNL) theory, where the generalpremise is that the structure of the largest scale is given by the critical eigenmode and that the higher harmonics are allderived from that fundamental mode (Stuart, 1960; Yahata, 1977; Jones, 1981; Gallaire et al., 2016). Despite being formallyvalid for only a small range of Reynolds numbers close to R c , the mathematical difficulties associated with the nonlinearityof the NSE often necessitate the use of WNL methods outside this domain of validity (Gallaire et al., 2016). Our resultsilluminate the physical mechanisms which lead to the eventual failure of WNL theory as the Reynolds number increases.WNL theory proceeds by expanding the solution in an asymptotic series about the bifurcation point such that the leadingorder solution u is the laminar base flow and the O ( ε ) solution u is given by the critical eigen mode. The O ( ε ) solution u as well as the mean flow correction is then found by solving the linear system forced by the nonlinear self interactionof u . The higher order terms may then be similarly computed sequentially by solving a forced linear system of the form L k u k = f ( u , u , ..., u k − ) . At the lowest Reynolds number considered here, R = u k depends only on interactions between larger scales. However, as discussed in§3.3, at higher Reynolds numbers the forcing is dominated by pairs of triads, one of which involves ˆ u k + , a mechanism whichis impossible in the WNL formulation.This means that near the bifurcation from the laminar state a model solution of the nonlinear flow may be truncated at thehighest wavenumber of interest since a given wavenumber depends only on its sub-harmonics. We define such a flow to be inthe “weakly nonlinear” (WNL) regime. As the Reynolds number increases the forcing cascade is no longer only from largeto small scales and an equally important inverse cascade mechanism emerges. In this case we define the flow to be in the “fully nonlinear” (FNL) regime. For the case of η = .
714 considered here, this transition occurs around some 100 < R < R ∼ O ( ) they persist in the limit of zerocurvature i.e. in the absence of centrifugal effects (Nagata, 1990; Sacco et al., 2019). At sufficiently high Reynolds numbers,they find that the temporal evolution of the r.m.s. velocity associated with the Taylor vortex and the mean shear are perfectlyout of phase and fluctuate with a common characteristic frequency. Their finding indicates a regenerative self sustainingprocess similar to the framework suggested by Waleffe (1997); Hamilton et al. (1995). Since we consider steady TVF it isdifficult to make a direct comparison between their results and ours. However, it is possible that the transition from weakly to igure 4. Azimuthal velocity component of the model solution’s primary Fourier mode (open circles) and forced Fouriermode (lines) as well as the Fourier modes from DNS (open squares) for R =
100 (red), R =
200 (blue), and R =
400 (black).Top row k = −
3, middle row k = −
6, bottom row: k = − igure 5. Azimuthal velocity component of the forced Fourier modes at R = v k , k (cid:48) ,are shown in dashed colors and the full Fourier mode, ˆ u k , θ , is plotted in solid black. The sum of the individual triadcomponents, (dashed lines) add up to the total forced mode (solid black). The sum of the individual triad components,(dashed lines) add up to the total forced mode (solid black). Top row k = −
2, bottom row: k = − As described in §3.2, we observe that in the FNL regime a crucial component of the forcing at a given wavenumber k is thedestructive interference of the two triads k = ( k − ) + k = ( k + ) −
1. This pair of triads lead to velocity contributionswith large amplitudes but with opposite sign. This means that an accurate reconstruction of a Fourier mode with wavenumber k requires knowledge of both its subharmonic k − k +
1. In this section we show that for streamwiseconstant and spanwise periodic solutions, as considered in this work, this large amplitude destructive interference is a directconsequence of the structure of the Fourier representation of the nonlinear term u · ∇ u . In cylindrical coordinates the nonlinearinteraction between two axisymmetric Fourier modes a = [ a r , a θ , a z ] and b = [ b r , b θ , b z ] with axial wavenumbers k a and k b isgiven by f a , b ≡ a · ∇ b + b · ∇ a (41)where the axial derivative in the gradient operator is replaced by multiplication by ik b and ik a respectively. For clarity ofexposition we limit the following analysis to the azimuthal component of the forcing and note that analogous arguments holdfor the remaining two components. The forcing at wavenumber k due to the interactions of k − f + ≡ ˆ u · ∇ ˆ u k − + ˆ u k − · ∇ ˆ u . (42)Using the continuity equation to eliminate the axial velocity, the azimuthal component takes the form f + θ = (cid:18) ˆ u , r ˆ u (cid:48) k − , θ + ˆ u k − , r ˆ u (cid:48) , θ + ˆ u , θ ˆ u k − , r r + ˆ u k − , θ ˆ u , r r (cid:19) − ( k − ) ( r ˆ u , r ) (cid:48) ˆ u k − , θ r − ( r ˆ u k − , r ) (cid:48) ˆ u , θ r ( k − ) . (43) igure 6. Azimuthal velocity component of the forced Fourier modes at R = v k , k (cid:48) ,are shown in dashed colors and the full Fourier mode, ˆ u k , θ , is plotted in solid black. The sum of the individual triadcomponents, (dashed lines) add up to the total forced mode (solid black). The predicted canceling azimuthal velocitycontributions, ˆ v ± k , cancel , θ , derived in §3.5 are plotted in red and green. Top row k = −
3, middle row k = −
6, bottom row: k = − igure 7. Projections of the velocity due to individual triadic interactions onto the full Fourier mode, Γ k , k , k , at R = k = −
2, bottom row: k = − Figure 8.
Projections of the velocity due to individual triadic interactions onto the full Fourier mode, Γ k , k , k , at R = k = −
3, middle row k = −
6, bottom row: k = − igure 9. Projections of the velocity due to individual triadic interactions onto the full Fourier mode summed over commonwavenumbers. R =
100 in red, R =
200 in blue squares, and R =
400 in black circles. Top row k = −
3, middle row k = −
6, bottom row: k = − k + − f − ≡ ˆ u − · ∇ ˆ u k + + ˆ u k + · ∇ ˆ u − , (44)with the azimuthal component taking the form f − θ = (cid:18) ˆ u , r ˆ u (cid:48) k + , θ + ˆ u k + , r ˆ u (cid:48) , θ + ˆ u , θ ˆ u k + , r r + ˆ u k + , θ ˆ u , r r (cid:19) + ( k + ) ( r ˆ u , r ) (cid:48) ˆ u k + , θ r + ( r ˆ u k + , r ) (cid:48) ˆ u , θ r ( k + ) . (45)Here the superscript (cid:48) denotes partial derivatives with respect to r . Additionally, note that for the streamwise constantfluctuations considered here ˆ u − = ˆ u ∗ = [ ˆ u , r , ˆ u , θ , − ˆ u , z ] .For some integer wavenumber k >
1, the Fourier modes associated with the nearest neighbor wavenumbers k ± u ( r ) k ± ≡ (cid:90) ∞ − ∞ u ( r , z ) e i β z ( k ± ) z dz = (cid:90) ∞ − ∞ u ( r , z ) e i β z k ( ± k ) z dz . (46)Since the destructive interference is most pronounced for small scales, we formally consider the case of k (cid:29)
1, for which wecan expand (46) in a Taylor series about k − = u k ± = (cid:90) ∞ − ∞ u (cid:16) e i β z kz ± i β z k − e i β z kz + O ( k − ) (cid:17) dz = ˆ u k + O ( k − ) . (47)This observation that for k (cid:29) u k and ˆ u k ± differ by a quantity which is O ( k − ) meaning that for large values of k the shapeof the Fourier modes does not change drastically with increasing k . Figure 4 shows that this is indeed the case. Substitutingˆ u k ± = ˆ u k + O ( k − ) into (43) and (45) we find at leading order f + θ = f θ , eq − k ( r ˆ u , r ) (cid:48) ˆ u k , θ r + O ( k − ) (48) − θ = f θ , eq + k ( r ˆ u , r ) (cid:48) ˆ u k , θ r + O ( k − ) (49)where f θ , eq is the same for both triads and is given by f θ , eq = (cid:20) ˆ u , r ˆ u (cid:48) k , θ + ˆ u k , r ˆ u (cid:48) , θ + ˆ u , θ ˆ u k , r r + ˆ u k , θ ˆ u , r r + ( r ˆ u , r ) (cid:48) ˆ u k , θ r (cid:21) . (50)The only remaining terms are equal in magnitude but of opposite sign. Furthermore, since both ˆ u , r and ˆ u k , θ are boundedand nonzero these terms will scale proportionally with k and therefore are expected to have large amplitudes since k (cid:29) ± k ( r ˆ u , r ) (cid:48) ˆ u k , θ r .Similar expressions can be derived for the other two coordinates such that the two vector forcing terms which are expectedto cancel are given byˆ f ± k , cancel ≡ ∓ k ( r ˆ u , r ) (cid:48) r ˆ u k ≈ ∓ k ( r ˆ u , r ) (cid:48) r ˆ u k ∓ ( k (cid:29) ) , (51)where the approximate equivalence on the right hand side is due to (46). This prediction may be tested by computing thecorresponding velocity contributions to the Fourier mode with wavenumber k , given byˆ v ± k , cancel = H k ˆ f ± k , cancel , ( k (cid:29) ) , (52)and comparing its shape to that of the total velocity contributions ˆ v k , ± defined in (39). Figure 6 shows the ˆ v ± k , cancel alongsidethe individual ˆ v k , k (cid:48) , and we see that the former is a quite accurate approximation of ˆ v k , ± for k >
3. Surprisingly, we seethat the prediction holds at least approximately even for k = , k (cid:29)
1. Thesefindings establish that ˆ f ± k , cancel is indeed responsible for the large amplitude destructive interference characteristic of the fullynonlinear regime.Inspection of the spectral dynamics of the flow corroborate this finding. If we assume that the Fourier modes ˆ u k obey apower law (cid:107) ˆ u k (cid:107) ∼ k − p , (53)then, from (51), the forcing component ˆ f ± k , cancel must obey the power law (cid:107) ˆ f ± k , cancel (cid:107) ∼ k − p . (54)Thus the flow will be in the WNL regime as long as p (cid:29) p (cid:46)
1. In figure 10 we plot the norm of the Fourier modes computed from the DNS data for a range of Reynolds numbers.For all cases the Fourier modes decay in Fourier space faster than k − which is depicted by the dashed black line. However,the decay rate at R =
100 is significantly faster than for the higher Reynolds number cases which seem to converge to a decayrate which is roughly independent of Reynolds number. The inset of figure 10 shows the exponent of the best fit power lawfor all Reynolds numbers. For R =
100 we fit the power law only to k ≤ < k ≤
10 the norm of the Fouriercomponents remains roughly constant. At R = k − and the transition to the fully nonlinear regime is associated with the decay rateapproaching k − . Here we address how the particular truncation values N kSVD were chosen, and how the number of retained wavenumbers andresolvent modes at each wavenumber affects the accuracy of the model. At R =
100 the flow is in the weakly nonlinear regimeand thus the flow may be arbitrarily truncated in Fourier space with out appreciably impacting the accuracy of the retainedharmonics. Additionally, in this case the optimal resolvent mode is a good approximation of the flow and thus retaining onlya single harmonic with N SVD = N SVD = igure 10. Norm of the Fourier modes computed from DNS at R = , , , , , ∼ k − . Inset shows exponent of best fit power law as in (53). Power law fit performed over the range 1 < k < R =
100 and 1 < k <
10 for all R > R = R = N k and N kSVD uniformly until the residual no longer decreased appreciably with added degrees of freedom. For R =
400 this“full” convergence was achieved with N k = N kSVD =
22. In figure 11 we plot the expansion coefficients σ k , j χ k , j in (18)and the χ k , j in (19). The σ k , j χ k , j and χ k , j represent the projection of the velocity and nonlinear forcing on to their respectiveresolvent basis ψ k , j and φ k , j respectively. We also plot the singular values σ k , j on the right y-axis. Figure 11 reveals thatwhile the velocity is represented by only a few response modes, the nonlinear forcing has a significant projection onto manysuboptimal modes.This finding is in agreement with Symon et al. (2021) and Morra et al. (2021) who show that the nonlinear forcing hassignificant projection onto the sub-optimal resolvent forcing modes for a variety of flows even if the resolvent operator is lowrank. The former considered both ECS as well as flow in a minimal channel while the latter focused entirely on turbulentchannel flow. In fact, as also observed by Morra et al. (2021), the projection onto the first two forcing modes, χ k , and χ k , , ismuch lower than the projection onto many of the suboptimal modes.We interpret this significant projection of the nonlinear forcing onto a wide range of suboptimal modes as beingnecessitated by the roll off in singular values and the structured nature of the nonlinear forcing. If the nonlinear forcing wasunstructured white noise, the solution would be given by the sum of resolvent response modes weighted by their singularvalues. However, for flows with structured nonlinear forcing, such as the TVF discussed here, the projection of the velocityfield onto the response basis, ψ k , j does not necessarily decrease monotonically with j as evidenced by figure 11. Therefore, ifwe consider a certain suboptimal ψ k , j which has a significant projection onto the velocity field, σ k , j χ k , j but is associated witha small singular value σ k , j the nonlinear forcing must by definition have a large projection onto the associated φ k , j since if σ k , j χ k , j is significant and σ k , j is small, then χ k , j must be large. Thus the lack of roll off in the χ k , j is a measure of how muchthe structure of the forcing differs from white noise.From a practical point of view we see that for 2 ≤ k ≤ ( k = ) and the higher harmonics, k > j (cid:38)
10. Neglecting these suboptimal modes for the higher harmonics does not affect theaccuracy of the solution and the associated 30% reduction in degrees of freedom results in a 90% reduction in computational igure 11.
Expansion coefficients of the velocity σ k , j χ k , j (blue circles) and nonlinear forcing χ k , j (green squares), andsingular values σ k , j (red triangles) for R = k = −
3, middlerow k = −
6, bottom row: k = − th order tensors in (30) scale as N . Neglecting these negligible sub-optimal modesin the higher harmonics we arrive at the final truncation values cited in §2, N kSVD = ∀ k ≤ , N kSVD = ∀ k >
4. With thisreduction in degrees of freedom the computational complexity has decreased to a point where the optimization may be carriedout cheaply on a personal computer. It is these results using these values of N kSVD that are presented in §3 and §3.2. However,if only the large scales are desired or lower levels of convergence are acceptable, the solution is robust to significantly moretruncation in both Fourier space and the SVD. It is well known, if not entirely understood, that Taylor vortices persist well into the turbulent regime (Grossmann et al., 2016).While the nature of the Taylor vortices does evolve with increasing Reynolds number as discussed in this work and Saccoet al. (2019), the general structure does not deviate significantly from the form at R =
400 shown in figures 1 and 2. Giventhe significant model reduction achieved by our model, we now investigate whether the large scale Taylor vortex structurecan be precomputed using our approach and then used to initialize a DNS at a higher Reynolds number to reduce the timeto converge to a statistically stationary state. Similar ideas have been investigated by Rosales and Meneveau (2006), whoinitialized DNS and LES of isotropic decaying turbulence with both standard Gaussian and more realistic non-Gaussianvector fields. They found that the latter, which displayed some of the physical features associated with turbulence, led toshorter transition times before realistic decay rates were observed.We performed two sets of DNS of Taylor Couette flow for a range of Reynolds numbers from 400 to 2000: the firstusing a random perturbation as an initial condition and one using the R =
400 model solution as an initial condition. Thesimulations were run until the torque at the inner and outer cylinder agree to within 1%. The simulation was then continuedfor an additional 200 non-dimensional time units at which point the simulation was deemed to be converged. We define thepercent reduction in time to convergence between the two cases P ≡ T − T m T ×
400 650 1000 1500 2000 P
82% 80% 75% 72% 65%
Table 3.
Percentage reduction in convergence time using model TVF solution as initial condition compared to randomperturbation as a function of Reynolds number. All cases use the R =
400 model result as an initial condition.where T and T m are the time required to reach convergence with the random and model initial conditions respectively. Table 3summarizes the savings for all the Reynolds numbers we considered. As expected the percentage of run time saved decreasesas the Reynolds number increases because the Taylor vortices change slightly and, more crucially, because the flow becomesmore three-dimensional and time dependent. However, it is remarkable that even at the highest Reynolds number, R = We have presented a fully nonlinear reduced order model of Taylor vortex flow for Reynolds numbers up to five timesgreater than the critical value. The resolvent formulation allows the governing equations for the fluctuations about a knownmean velocity to be transformed into a set of polynomial equations. We approximate the solution to these equations byminimizing their associated residual in conjunction with a constraint which ensures the model generates Reynolds stressescompatible with the input mean velocity profile. We are able to generate model solutions which solve the NSE to a very goodapproximation and replicate the flow field computed through DNS at a tiny fraction of the computational cost. We believe thisis the first explicit example of “closing the resolvent loop” published in the literature, although Rosenberg (2018) presenteda similar analysis applied to ECS in a channel in his doctoral thesis which inspired this work. We analyzed the nonlinearinteractions driving the flow for a range of Reynolds numbers and identified the transition from a weakly nonlinear regimeclose to the bifurcation from the laminar state where the structure of the flow is accurately modeled by the linear dynamicsand the forcing cascade is purely from large scales to small scales. At higher Reynolds numbers we define a fully nonlinearregime where an inverse forcing cascade from small to large scales emerges to counter the cascade from large to small scales.The dominant nonlinear interactions at a given wavenumber k involve the pair of triadic interactions k = ( k ± ) ∓
1. For thehighest Reynolds number case the pair of triads k = ( k ± ) ∓ We thank Kevin Rosenberg for many inspiring and helpful discussions and gratefully acknowledge the support of the U.S.Office of Naval Research under grants ONR N00014-17-1-2307 and N00014-17-1-3022.
The Navier-Stokes operator in cylindrical coordinates linearized about a one dimensional azimuthal mean flow U ( r ) andFourier transformed in θ , z , and t is given by L k = inUr + R (cid:16) r − ∇ (cid:17) R (cid:16) inr (cid:17) − Ur ∂∂ r (cid:16) ∂ U ∂ r + Ur (cid:17) − R (cid:16) inr (cid:17) inUr + R (cid:16) r − ∇ (cid:17) inr inUr − R ∇ ik r + ∂∂ r inr ik , (56)where the Laplacian operator is defined as ∇ = ∂ ∂ r + r ∂∂ r − (cid:18) k + n r (cid:19) . (57)The weight matrix, M , is defined as M ≡ . (58)For the streamwise constant modes considered here the continuity equation may be used to eliminate pressure and theaxial velocity to write the system in terms of the radial and azimuthal velocity, ˆ q = [ ˆ u r , ˆ u θ ] such that˜ L k ˆ q = ˜ M ˆ f (59)where ˆ f = [ ˆ f r , ˆ f θ , ˆ f z ] and the linear operators take the form˜ L k = (cid:34) − R ∇ k Ur r ∂∂ r (cid:0) rU (cid:1) − R ∇ (cid:35) , (60)˜ M ≡ (cid:20) − k z − ik z ∂∂ r (cid:21) . (61) EFERENCES
Andereck, C. D., Liu, S. S., and Swinney, H. L. (1986). Flow regimes in a circular Couette system with independently rotatingcylinders.
Journal of Fluid Mechanics , 164:155–183.Beaume, C., Chini, G. P., Julien, K., and Knobloch, E. (2015). Reduced description of exact coherent states in parallel shearflows.
Physical Review E , 91(4):043010.Coles, D. (1965). Transition in circular Couette flow.
Journal of Fluid Mechanics , 21(3):385–425.Dessup, T., Tuckerman, L. S., Wesfreid, J. E., Barkley, D., and Willis, A. P. (2018). Self-sustaining process in Taylor-Couetteflow.
Physical Review Fluids , 3:123902.Gallaire, F., Boujo, E., Mantic-Lugo, V., Arratia, C., Thiria, B., and Meliga, P. (2016). Pushing amplitude equations far fromthreshold: application to the supercritical Hopf bifurcation in the cylinder wake.
Fluid Dynamics Research , 48(6):061401.Gebhardt, T. and Grossmann, S. (1993). The Taylor-Couette eigenvalue problem with independently rotating cylinders.
Zeitschrift f¨ur Physik B Condensed Matter , 90(4):475–490.Grossmann, S., Lohse, D., and Sun, C. (2016). High–Reynolds Number Taylor-Couette Turbulence.
Annual Review of FluidMechanics , 48:53–80.Hamilton, J. M., Kim, J., and Waleffe, F. (1995). Regeneration mechanisms of near-wall turbulence structures.
Journal ofFluid Mechanics , 287:317–348.Huisman, S. G., van der Veen, R. C. A., Sun, C., and Lohse, D. (2014). Multiple states in highly turbulent Taylor-Couetteflow.
Nature Communications , 5:3820.Illingworth, S. J. (2020). Streamwise-constant large-scale structures in Couette and Poiseuille flows.
Journal of FluidMechanics , 889.Jones, C. A. (1981). Nonlinear Taylor vortices and their stability.
Journal of Fluid Mechanics , 102:249–261.Mantiˇc-Lugo, V., Arratia, C., and Gallaire, F. (2014). Self-Consistent Mean Flow Description of the Nonlinear Saturation ofthe Vortex Shedding in the Cylinder Wake.
Physical Review Letters , 113(8):084501.Mantiˇc-Lugo, V., Arratia, C., and Gallaire, F. (2015). A self-consistent model for the saturation dynamics of the vortexshedding around the mean flow in the unstable cylinder wake.
Physics of Fluids , 27(7):074103.Marcus, P. S. (1984). Simulation of Taylor-Couette flow. Part 2. Numerical results for wavy-vortex flow with one travellingwave.
Journal of Fluid Mechanics , 146:65–113.Maretzke, S., Hof, B., and Avila, M. (2014). Transient growth in linearly stable Taylor–Couette flows.
Journal of FluidMechanics , 742:254–290.McKeon, B. J. and Sharma, A. S. (2010). A critical-layer framework for turbulent pipe flow.
Journal of Fluid Mechanics ,658:336–382.McMullen, R. M., Rosenberg, K., and McKeon, B. J. (2020). Interaction of forced Orr-Sommerfeld and Squire modes in alow-order representation of turbulent channel flow.
Physical Review Fluids , 5(8):084607.Moarref, R., Jovanovi´c, M. R., Tropp, J. A., Sharma, A. S., and McKeon, B. J. (2014). A low-order decomposition ofturbulent channel flow via resolvent analysis and convex optimization.
Physics of Fluids , 26(5):051701.Morra, P., Nogueira, P. A. S., C., A. V. G., and Henningson, D. S. (2021). The colour of forcing statistics in resolvent analysesof turbulent channel flows.
Journal of Fluid Mechanics , 907:A24.Nagata, M. (1990). Three-dimensional finite-amplitude solutions in plane Couette flow: bifurcation from infinity.
Journal ofFluid Mechanics , 217:519–527.Nogueira, P. A. S., Morra, P., Martini, E., Cavalieri, A. V. G., and Henningson, D. S. (2021). Forcing statistics in resolventanalysis: application in minimal turbulent Couette flow.
Journal of Fluid Mechanics , 908.Ostilla, R., Stevens, R. J. A. M., Grossmann, S., Verzicco, R., and Lohse, D. (2013). Optimal Taylor–Couette flow: directnumerical simulations.
Journal of Fluid Mechanics , 719:14–46.Ostilla-M´onico, R., van der Poel, E. P., Verzicco, R., Grossmann, S., and Lohse, D. (2014). Exploring the phase diagram offully turbulent Taylor–Couette flow.
Journal of Fluid Mechanics , 761:1–26.Rand, D. (1982). Dynamics and symmetry. Predictions for modulated waves in rotating fluids.
Archive for Rational Mechanicsand Analysis , 79(1):1–37.Rigas, G., Sipp, D., and Colonius, T. (2021). Nonlinear input/output analysis: application to boundary layer transition.
Journal of Fluid Mechanics , 911.Rosales, C. and Meneveau, C. (2006). A minimal multiscale Lagrangian map approach to synthesize non-Gaussian turbulentvector fields.
Physics of Fluids , 18(7):075104.Rosenberg, K. (2018).
Resolvent-based modeling of flows in a channel . PhD thesis, California Institute of Technology.Rosenberg, K. and McKeon, B. J. (2019a). Computing exact coherent states in channels starting from the laminar profile: Aresolvent-based approach.
Physical Review E , 100(2):021101.Rosenberg, K. and McKeon, B. J. (2019b). Efficient representation of exact coherent states of the Navier–Stokes equations sing resolvent analysis.
Fluid Dynamics Research , 51(1):011401.Sacco, F., Verzicco, R., and Ostilla-M´onico, R. (2019). Dynamics and evolution of turbulent Taylor rolls.
Journal of FluidMechanics , 870:970–987.Stuart, J. T. (1960). On the non-linear mechanics of wave disturbances in stable and unstable parallel flows Part 1. The basicbehaviour in plane Poiseuille flow.
Journal of Fluid Mechanics , 9(3):353–370.Symon, S., Illingworth, S. J., and Marusic, I. (2021). Energy transfer in turbulent channel flows and implications for resolventmodelling.
Journal of Fluid Mechanics , 911.Taylor, G. I. (1923). VIII. Stability of a viscous liquid contained between two rotating cylinders.
Philosophical Transactions ofthe Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character , 223(605-615):289–343.van der Poel, E. P., Ostilla-M´onico, R., Donners, J., and Verzicco, R. (2015). A pencil distributed finite difference code forstrongly turbulent wall-bounded flows.
Computers & Fluids , 116:10–16.van Gils, D. P. M., Bruggert, G., Lathrop, D. P., Sun, C., and Lohse, D. (2011). The Twente turbulent Taylor-Couette (T3C)facility: strongly turbulent (multiphase) flow between two independently rotating cylinders.
The Review of ScientificInstruments , 82(2):025105.van Gils, D. P. M., Huisman, S. G., Grossmann, S., Sun, C., and Lohse, D. (2012). Optimal Taylor–Couette turbulence.
Journal of Fluid Mechanics , 706:118–149.Verzicco, R. and Orlandi, P. (1996). A Finite-Difference Scheme for Three-Dimensional Incompressible Flows in CylindricalCoordinates.
Journal of Computational Physics , 123(2):402–414.Waleffe, F. (1997). On a self-sustaining process in shear flows.
Physics of Fluids , 9(4):883–900.Yahata, H. (1977). Slowly-Varying Amplitude of the Taylor Vortices near the Instability Point. II: Mode-Coupling-TheoreticalApproach.
Progress of Theoretical Physics , 57(5):1490–1496.Zhu, X., Phillips, E., Spandan, V., Donners, J., Ruetsch, G., Romero, J., Ostilla-M´onico, R., Yang, Y., Lohse, D., Verzicco, R.,Fatica, M., and Stevens, R. J. A. M. (2018). AFiD-GPU: A versatile Navier–Stokes solver for wall-bounded turbulentflows on GPU clusters.
Computer Physics Communications , 229:199–210., 229:199–210.