Eulerian-Lagrangian Formulation for Compressible Navier-Stokes Equations

The Eulerian-Lagrangian formulation of the (inviscid) Euler dynamics in terms of advected Weber-Clebsch potentials (Lamb, 1932), was extended by Constantin to cover the viscous Navier-Stokes dynamics (Constantin, 2001). Numerical studies (Ohkitani & Constantin, 2003), of this formulation of the Navier-Stokes equations concluded that the diffusive Lagrangian map becomes non-invertible under time evolution and requires resetting for its calculation. They proposed the observed sharp increase of the frequency of resettings as a new diagnostic of vortex reconnection. In previous work we were able (Cartes et al., 2007; 2009) to complement these results, using an approach that is based on a generalised set of equations of motion for the Weber-Clebsch potentials, that turned out to depend on a parameter τ, which has the unit of time for the Navier-Stokes case. Also to extend our formulation to magnetohydrodynamics, and thereby obtain a new diagnostic for magnetic reconnection. In this work we present a generalisation of the Weber-Clebsch variables in order to describe the compressible Navier-Stokes dynamics. Our main result is a good agreement between the dynamics for the velocity and density fields that come from the dynamics of Weber-Clebsch variables and direct numerical simulations of the compressible Navier-Stokes equations. We first present the inviscid Eulerian-Lagrangian theory, then Constantin’s extension to viscous fluids and derive our equations of motion for the Weber-Clebsch potentials that describe the compressible Navier-Stokes dynamics. Then, performing direct numerical simulations of the Taylor-Green vortex, we check that our formulation reproduces the compressible dynamics.


Introduction
The Eulerian-Lagrangian formulation of the (inviscid) Euler dynamics in terms of advected Weber-Clebsch potentials (Lamb, 1932), was extended by Constantin to cover the viscous Navier-Stokes dynamics (Constantin, 2001). Numerical studies (Ohkitani & Constantin, 2003), of this formulation of the Navier-Stokes equations concluded that the diffusive Lagrangian map becomes non-invertible under time evolution and requires resetting for its calculation. They proposed the observed sharp increase of the frequency of resettings as a new diagnostic of vortex reconnection. In previous work we were able (Cartes et al., 2007;2009) to complement these results, using an approach that is based on a generalised set of equations of motion for the Weber-Clebsch potentials, that turned out to depend on a parameter τ, which has the unit of time for the Navier-Stokes case. Also to extend our formulation to magnetohydrodynamics, and thereby obtain a new diagnostic for magnetic reconnection. In this work we present a generalisation of the Weber-Clebsch variables in order to describe the compressible Navier-Stokes dynamics. Our main result is a good agreement between the dynamics for the velocity and density fields that come from the dynamics of Weber-Clebsch variables and direct numerical simulations of the compressible Navier-Stokes equations. We first present the inviscid Eulerian-Lagrangian theory, then Constantin's extension to viscous fluids and derive our equations of motion for the Weber-Clebsch potentials that describe the compressible Navier-Stokes dynamics. Then, performing direct numerical simulations of the Taylor-Green vortex, we check that our formulation reproduces the compressible dynamics. Will-be-set-by- IN-TECH here p is the pressure field. Now the equations for evolution of the vorticity ω = ∇×u field D t ω = ω ·∇u , ( 2 ) where D t is the convective derivative A well known consequence of this equation is the preservation of vorticity lines (Helmholtz's theorem).
In other words the intersections of surfaces λ =const. and μ =const. are the vorticity lines. If vorticity lines follow Euler equations and are preserved, then the fields λ and μ follow the fluid. Clebsch variables can also be used to find a variational principle for Euler equations. We can write a Lagrangian density for Euler equations L = |u| 2 2 + λ∂ t μ , ( 8 ) and the variations of L in function of the fields λ, μ and φ give us the system equations δL δμ = −D t λ = 0( 9 ) δL δλ = D t μ = 0 δL δφ = ∇·u = 0.
we can obtain the evolution equation for u

Weber transformation
Let us note a i as the initial coordinate (at t = 0) of a fluid element and X i (a, t) its position at time t and note A i (x, t) the inverse application: At time t Eulerian coordinates are by definition the variables x i = X i (a, t) then the Lagrangian velocity of a fluid element isũ and its acceleration Newton equations for the fluid element are where the forces F i X (a, t) are given by and p (X (a, t)) is the pressure field in Eulerian coordinates. Therefore the movement equations for the fluid elements are For an incompressible fluid, the transformation matrix, between Lagrangian and Eulerian coordinates, verifies this value is fixed from the relation between the volume elements in the two coordinate systems. We also note that this transformation is always invertible. From Eq. (16) we perform a coordinate transformation for the derivatives of p using Eulerian-Lagrangian Formulation for Compressible Navier-Stokes Equations www.intechopen.com to obtain wherep (a, t) is the pressure field in Lagrangian coordinates. We multiply Eq. (19) with the inverse coordinate transformation ∂X i ∂a j in order to obtain which is the Lagrangian form for the dynamic equations. The left hand side of this equation can be written as and Eq. (20) becomes where the termq (a, t) is given bỹ Now let us integrate Eq. (22) over t, maintaining a i fixed whereφ is written asφ This equation system (24) is called Weber transformation (Lamb, 1932). Now we perform a coordinate transformation identifyingμ i (a, t)=a i and the initial velocityλ i (a, t)=ũ i 0 (a) we obtain the evolution equations for the fields in Lagrangian coordinates 112 Hydrodynamics -Optimizing Methods and Tools

www.intechopen.com
Eulerian-Lagrangian Formulation for Compressible Navier-Stokes Equations If we now go to the Eulerian coordinates, identifying μ i (x, t)=A i (x, t) and λ i (x, t)= u i 0 (μ(x, t)), we obtain the Weber-Clebsch transformation Using the convective derivative, the dynamic equations for the Clebsch variables Eq. (27) can be written in Eulerian coordinates as The Weber-Clebsch transformation Eq. (28) and its evolution laws Eq. (29) are very similar to Clebsch variables Eq. (4) and the system (9). An important difference is the number of potential pairs. If we use Clebsch variables Eq. (4) to represent the velocity field u we have the problem that u is restricted to fields with mean helicity of value zero (Grossmann, 1975). In fact, writing h in terms of Clebsch variables the term λ∇μ is perpendicular to ∇λ ×∇μ and then their scalar product is zero. For the other terms, we integrate by parts but in a periodic domain the first term in the right hand side is zero. We also know that ∇·ω = 0 and therefore we have and we finally get If we consider now two pairs of Clebsch variables, for each component of the velocity field, we have and we arrive to a system of equations of second degree in its unknowns, this system does not have an analytic solution and we don't have a systematic way to find λ i and μ i for an arbitrary velocity field u.
If we use now the same number of pairs as spatial variables (three in this case), we get the Weber-Clebsch transformation with this representation we can write an arbitrary velocity field defining, at t = 0: which is completely equivalent to the Weber transformation.

Constantin's formulation of Navier-Stokes equations
Here we will recall Constantin's extension for the Eulerian-Lagrangian formulation of Navier-Stokes equations. The departing point (Constantin, 2001), is the expression for the Eulerian velocity u = u 1 , u 2 , u 3 from the Weber-Clebsch transformation The fields in this equation admit the same interpretation as in the Weber transformation: λ m are the Lagrangian velocity components, μ m are the Lagrangian coordinates and φ fixes the incompressibility condition for the velocity field. In a way similar to the Weber transformation, we have the Lagrangian coordinates a i = μ i (x, t) and the Eulerian coordinates x i = X i (a, t).
If we now consider the first term of the right hand side in Eq. (39) as a coordinate transformation, it is possible to write their derivatives in Lagrangian coordinates, as in Eq.
114 Hydrodynamics -Optimizing Methods and Tools

www.intechopen.com
Eulerian-Lagrangian Formulation for Compressible Navier-Stokes Equations In the same way it is possible to write the derivatives of the Eulerian coordinate in terms of Lagrangian coordinates We also have the relation for the commutators Using relations Eq. (40), Eq. (41) and Eq. (42) we can compute the commutators between ∂ x and ∂ a ∂ ∂a i , Introducing the displacement vector ℓ m = μ m − x m which relates the Eulerian position x to the original Lagrangian position μ, we can express the commutator Eq. (43) as The term C m,k;i is related to the Christoffel coefficients Γ m ij of the flat connection in R 3 by the formula We consider now the diffusive evolution of our fields, with that goal in mind we define the operator where ν is the viscosity and u is the Eulerian velocity. When the operator Eq. (46) is applied over a vector or a matrix each component is taken in an independent way. Constantin imposes that the coordinates μ i are advected and diffused so they follow We also need a coordinate transformation that can be invertible at any time t, that condition is always satisfied when the diffusion is zero (ν = 0) and the fluid is incompressible, because the fluid element volume is preserved by the coordinate transformation, and therefore

115
Eulerian-Lagrangian Formulation for Compressible Navier-Stokes Equations In order to get the evolution for the λ i fields we apply D t on Eq. (39), and using the relation we obtain We also have, from Navier-Stokes equations: We compute the term △u i with the transformation Eq. (39), using Eq. (47) and regrouping the terms we obtain Now we split the φ field to obtain the pressure equation where c is a constant.
To obtain the λ l dynamics we have to invert the transformation matrix in Eq. (53) if the determinant of the transformation matrix follows then it will be impossible to perform a coordinate transformation. Therefore, if the matrix is invertible, the λ l dynamics is written as We have to remark that the dynamics of u is completely described by Eq. (47), (57) and the incompressibility condition for u, thus Eq. (54) becomes an identity.

Generalisation of Constantin's formulation
We begin with the Weber-Clebsch transformation for the velocity field u 116 Hydrodynamics -Optimizing Methods and Tools

www.intechopen.com
Eulerian-Lagrangian Formulation for Compressible Navier-Stokes Equations and perform a variation on the Weber-Clebsch transformation Eq. (58) to obtain the relation (Cartes et al., 2007) here δ represents a spatial or temporal variation. In the system (59) it is already possible to see that we have three equations (δu) and six unknowns to find (δλ i and δμ i , δφ is fixed by the continuity equation).
In order to write the temporal evolution of u, in terms of Weber-Clebsch potentials, we use the convective derivative D t and the identity We compute now the convective derivative for u, which can be written as a function of the potentials For that purpose we write in gradient form and noting that Finally we regroup the gradients We must note that this expression is very similar to Eq. (59), the only difference is given by the term 1 2 u 2 , that comes from the commutator between the gradient and the convective derivative.

General formulation for the compressible Navier-Stokes equations
We now consider the compressible Navier-Stokes equations with a general forcing term f where w is the enthalpy and ρ the density field. For this work we will consider, for simplicity and without loss of generality, a barotropic fluid, then the relation for the enthalpy w is where Ma is the Mach number for a flow of density ρ 0 = 1a n dv e l o c i t yu ∼ 1. In this approximation we suppose the density field ρ is very near to the uniformity state and consequently the Mach number is small. The usual compressible Navier-Stokes equations are obtained when the forcing term f is the viscous dissipation The idea is to find the evolution equations, in the most general way, for the potentials Eq. (58), now we replace D t , in the equations of motion Eq. (65), by its expression Eq. (64) and we define To wit we made the separation in Eq. (64) and Eq. (65) between gradient and non-gradient terms here G is an arbitrary gauge function, which comes from the fact that the separation in gradient and non-gradient terms is not unique.
The equation system (70) has 3 linear equations and 6 unknowns L i , M i .I no r d e rt os o l v et h i s system with f we must remark that, when ν = 0, the fields λ i and μ i follow Euler dynamics If we are in the overdetermined case (more equations than unknowns), in general, equation has no solution. Then we consider only the under determined case (more unknowns than equations). In order to obtain evolution equations in the same way as (Constantin, 2001) we look for advection diffusion equations. With that goal in mind we introduce L i and M i ,definedby The terms L i and M i must verify where G is an arbitrary scalar function, linked to the old function G, by the relation The process used to obtain L i and M i consists in solving the linear system (73).

Moore-Penrose solution
As the system (73) is under determined, we must impose additional restrictions to solve it.
The most straightforward way is to force the coefficients M i = 0 as in Constantin's formulation then we will have 3 equations for the 3 unknowns. Another, more general, method relies in the imposition of additional conditions on the solution's length.
For that purpose we use the Moore-Penrose algorithm (Ben-Israel & Greville, 1974;Moore, 1920;Penrose, 1955), which produces 3 additional conditions that allow us to solve this more general system (73).
For the under determined case, the Moore-Penrose general solution consists in finding the solution to the linear system (73) with the imposition that the norm is minimal. The constant τ is introduced here because λ i and μ i have different dimensions. In fact, the Weber-Clebsch transformation Eq. (58) means that the dimensions for λ i and μ i are where H represents the squared symmetric matrix and the arbitrary function G is given by in order to minimise the general norm Eq. (82) with the objective to achieve numerical stability in our simulations. Replacing these solutions L i and M i in Eq. (72) we arrive to the explicit evolution equations

Comparison of the invertibility conditions
Constantin's method will have problems when the determinant det(∇µ)=0 which is the case in a manifold of codimension 1. In three dimensional space the generic situation becomes that, for any point in the space (x 1 , x 2 , x 3 ), there is a time t * for which the determinant becomes zero.
In our more general formulation, with three equations and six unknowns, the inversibility of which corresponds to isolated points in a manifold of codimension 4 in space-time.
In consequence the condition det(∇µ)=0 will arrive more frequently because of its lower codimension, than the condition with a higher codimension, for det(H)andτ → 0isasingular limit.

120
Hydrodynamics -Optimizing Methods and Tools

Resettings
As we saw, when the determinant det(H) is zero the Weber-Clebsch potential evolution equations (83) are no longer defined. In order to avoid this situation, we follow (Ohkitani & Constantin, 2003) and we perform a resetting. More precisely, when the spatial minimum of the determinant where ǫ is a pre defined lower limit. We reset the fields in the following way where u i (t 0 ) are the components of the velocity field obtained from Eq. (58) in the instant t 0 , when the resetting is performed, and the fields λ i , μ i and φ i are generated in the same way as in section 4.2. It was already pointed (Cartes et al., 2007;Ohkitani & Constantin, 2003), that the vanishing of det(H) is related to intense particle diffusion that takes place near reconnection of vorticity lines in the case of incompressible fluids, that means the spatial position of the minima of det(H) are the places where the reconnections take place.

Numerical results
In this section we will show the results from numerical simulations of our formulation for compressible Navier-Stokes equations. We used pseudo-spectral methods because they are easy to implement and their high precision. The technical details of the implementation are described in section 6.

Taylor-Green flow
The Taylor-Green flow is an standard flow used in the study of turbulence (Taylor & Green, 1937). Its advantages are the existence of numerous studies, see for instance (Brachet et al., 1983) and references therein, which allow us to perform comparisons, at the same time we can economise memory and computation resources by using its symmetries (Cartes et al., 2007). The initial Taylor-Green condition is: As the length and the initial velocity are of order 1, the Reynolds number is defined as R = 1/ν.

Periodic field generation
Periodic fields are generated from the Weber-Clebsch representation Eq. (58) as we also impose that μ i p and the other fields λ i and φ in Eq. (58) are periodic. In order to generate an arbitrary velocity field u we can use We note that the non-periodic part of μ i in Eq. (88) is made in a way that μ i gradients are periodic.
The initial ρ is given by imposing the incompressibility condition over w at t = 0 and the relation Eq. (66).

Simulation results
The following simulations were made using a spatial resolution of 128 3 points, a Reynolds number of 200 and a Mach number of 0.3. We will compare the velocity field which comes from simulations made with the Weber-Clebsch potentials with the velocity field that comes from a direct numerical simulation of Navier-Stokes equations. Navier-Stokes equations are integrated using standard pseudo-spectral methods (Gottlieb & Orszag, 1977). The temporal scheme is Adams-Bashforth of order 2 (for details see section 6.3). In order to characterise and measure the precision of our algorithm for the Weber-Clebsch potentials, we compute the associated enstrophy which is defined as where E(k, t) is the energy spectrum and can be described from the velocity field in Fourier spaceû(k, t) as Then E(k, t) is obtained as the mean value over spherical shells with thickness △k = 1. This enstrophy is computed from the velocity field which comes from the Weber-Clebsch and direct Navier-Stokes simulations. Fig. (1) shows the temporal evolution of the enstrophy for different values of the parameter τ. We found good agreement between our formulation and the direct Navier-Stokes simulations. The spatial mean of the quantity ρ 2 /2, which represents the density field, can be seen in Fig.  (2). As our λ and μ fields evolved in time we had to reset them to be able to continue the simulation as the coordinate transformation becomes non-invertible. The temporal evolution of the interval between resettings is characterised by where t j is the resetting time, we fixed the value for the lower limit of det(H)a sǫ 2 = 0.01, is shown in Fig. (3). We can see that, for a given time, the interval is a growing function of τ. However the shape of △t is well preserved even when the range of τ goes through several orders of magnitude.

Conclusions and perspectives
We arrived to a good agreement between the derived generalised equations of motion for the Weber-Clebsch potentials that implying that the velocity field follows the compressible Navier-Stokes equations. These new equations were shown to depend on a parameter with the dimension of time, τ. Direct numerical simulations of the Taylor-Green vortex were performed in order to validate this new formulation. This Eulerian-Lagrangian formulation of compressible Navier-Stokes equations, allows us to study in detail the reconnection process, the turbulence generated by such process and the sound generated by those moving fluids using for example the two antiparallel vortex approach (Virk et al., 1995). This subject is known as aeroacoustics (Lighthill, 1952), which is relevant for aerodynamic noise production, and is a key issue in the design of air planes, turbines, etc.

Appendix -Numerical methods
The simulated equations are nonlinear partial differential equations solved by the pseudo-spectral methods. The flows in this work are periodic because we work in a periodic box.
A periodic field f verifies: f (x + L)= f (x) where L is the box periodicity length. In our simulations we choose L = 2π. In this representation a continuous function can be expressed by the infinite Fourier series Then we can define the scalar product by where the Fourier series coefficients are In a numerical simulation f is known by its values in a finite number of points over L with in a way that the distance between the points is The discretisation points are supposed to capture the shape of f (x).T h ex variable is defined in physical space and the points j = 0, 1, . . . , N − 1 are called collocation points. Then the Fourier coefficients Eq. (96) becomef with x n = 2πn N n = 0,1,...,N − 1 .
This is the discrete Fourier transformation (DFT). We projected f over a base formed by N sine and cosine functions. Then we can find an approximation for f , f N , by inverting Eq. (100) in the following way The points k = − N 2 ,..., N 2 − 1formthediscretespectral space. They characterise the functions from our projection base. We must note that the functions Eq. (100) and Eq. (102) are written in a more symmetric way than Eq. (94) and Eq. (96). A priori, to perform the summations Eq. (100) and Eq. (102) we must perform a number of operations of order O N 2 . This requirement was a handicap in the use of this method, until the introduction of Cooley and Turkey algorithm for the fast Fourier transformation (FFT) in 1965 (Cooley & Tukey, 1965). This FFT algorithm allows us to reduce the number of operations to O (N log 2 N). The use of spectral methods is justified by their convergence which is better than the convergence obtained with finite differences. That comes from the fact that, in a finite differences computation of order p, the approximation coefficients of a field f in a Taylor expansion of p + 1 points have an error of order O (△x p ). On the other hand with spectral is what we call an exponential convergence. From this point of view, the spectral methods gave us a considerable gain in memory use from a fixed precision.

Pseudo-spectral methods
The pseudo-spectral methods are based in the computation of the approximation of a determined function interpolating over a collocation point set, that means the differential equation will be exactly solved over the collocation points.
We chose this method in order to compute the convolutions in physical space. As base we use the DFT from trigonometric functions which corresponds to the collocation points. In the case of non-linear terms we use two inverse FFT (O(N log 2 N) operations) in order to have these terms in physical space. Then we compute the product (N operations) and we perform one FFT (O(N log 2 N) operations) in order to go back into the spectral space. For example, if we compute we have The summation over n includes all the terms where p + q − k ≡ O [N]. Then the approximated coefficients for the term w (k) are composed by all the exact coefficients plus other terms for which the correspondent function e ikx n [N] can not be distinguished from the function e ikx n , this phenomena is called aliasing, the Fourier modes with higher wave numbers are taken for modes with lower wave numbers. Even with these defects the pseudo-spectral methods have interesting properties, for example when we have to deal with multiple dimensions. In fact the DFT in 3D can be factorised as e i k· x = e ik x x e ik y y e ik z z . Moreover a DFT in 3D can be computed as a succession of DFT in 1D, so for N 3 point we have to perform 3N 2 DFT in 1D over N points and then 3N 3 log 2 N operations, on the other hand if we compute the products as convolutions in spectral space we need to perform N 6 operations. This passage by physical space and the use of pseudo-spectral methods allow us to gain in computation time but generates the problem of aliasing.

Aliasing correction
The only way to correct the aliasing error over the convoluted term is eliminating the aliased terms which, in spectral space, belong to |k|≤N/2. Over a grid of N points, working modulo N,t h ev a l u e so fk will belong to the interval − N 2 , N 2 . The method which allows us to solve the problem implies the elimination of the part of the spectrum which lays outside the interval ] − k max , k max [, with the condition (for nonlinearities of order 2) 2k max − N < −k max ,thismeansk max < N 3 . W eremovethenallthe values from the spectrum whose wave numbers are bigger than N 3 and smaller than − N 3 .I n this way all the replicated values are fixed to zero in each time step and they are no longer a problem. The aliasing correction is very expensive, because we lose one third of, otherwise, useful modes, but these computations are completely equivalent to the Galerkin truncation.

Integration by parts and conservation of energy
Let us suppose for a moment that our product computation is not dealiased.W ec o n s i d e rt h e product of the quantity f by ∂ x g and perform an integration by parts in spectral space. We perform the summation in the interval ] − N/2, N/2[ T h et e r mw h i c hf o r b i d su st oi n t e g r a t eb yp a r t si na ne x a c tw a yi s(n − j) [N] .I f w e d o a dealiasing over the term the summation is now over the interval ] − k max , k max [ and Let us remember that the integration by parts is a necessary step in the computation of the amount of energy present in our system, and it must be preserved if there is not diffusion. This physical requirement is satisfied by the dealiased spectral methods.

Temporal scheme
We made a pseudo-spectral solver for the equations, with periodic boundary conditions and a FFT base which allows us to integrate the partial differential equations where we call L the linear operator in Fourier space and N the nonlinear term. The resolution method is a second order finite difference. The time step is made by the explicit Adam-Bashforth method a t+△t = 1 − ν k 2 2 △t a t + △t 3 2 N(a t ) − 1 2 N(a t−△t ) 1 + ν k 2 2 △t . (108)

www.intechopen.com
To begin with the temporal integration we did an Euler time step in the following way a △t = a 0 + △tN(a 0 ) 1 + νk 2 △t , where a 0 is the initial condition. The constant evolution of the calculation capacity of the modern computers implies in a permanent effort to adjust the existing numerical codes, or to create new codes following new points of view, aiming to adequately simulate fluid flows and the related transport of physical properties. Additionally, the continuous improving of laboratory devices and equipment, which allow to record and measure fluid flows with a higher degree of details, induces to elaborate specific experiments, in order to shed light in unsolved aspects of the phenomena related to these flows. This volume presents conclusions about different aspects of calculated and observed flows, discussing the tools used in the analyses. It contains eighteen chapters, organized in four sections: 1) Smoothed Spheres, 2) Models and Codes in Fluid Dynamics, 3) Complex Hydraulic Engineering Applications, 4) Hydrodynamics and Heat/Mass Transfer. The chapters present results directed to the optimization of the methods and tools of Hydrodynamics.