Renato N Elias
email: renato@nacad.ufrj.br

Milton A Gonçalves
email: milton@nacad.ufrj.br

G A Coutinho

Paulo T T Esperança

D Martins

S Ferreira

Alvaro L G A Coutinho
email: alvaro@nacad.ufrj.br

Marcos A D Martins
email: marcos.martins@petrobras.com.br

Marcos D A S Ferreira

PROGRESSIVE WAVE SIMULATION USING STABILIZED EDGE-BASED FINITE ELEMENT METHODS

établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Flows involving waves occur in several problems in hydrodynamics. Sloshing of liquids in tanks, wave breaking in ships, offshore platforms motions and green water on decks are important examples of this class of problems. The main computational challenge when solving such highly nonlinear problem is determining the evolution of the interface location. There are a large number of numerical methods devoted to the computation of free-surface problems. These methods are frequently classified as interface tracking and interface capturing methods. Interface tracking methods are based on a Lagrangian framework where the moving interface or boundary is explicitly tracked by the computational grid or by the particles of meshless methods which must be deformed or moved in order to follow the fluid flow. The deforming-spatialdomain/stabilized space time finite element formulation (DSD/SST) proposed by Tezduyar [START_REF] Tezduyar | Finite Element Stabilization Parameters Computed from Element Matrices and Vectors[END_REF] and Tezduyar et al [START_REF] Tezduyar | Incompressible Flow Computations with Stabilized Bilinear and Linear Equal-Order-Interpolation Velocity-Pressure Elements[END_REF][START_REF] Tezduyar | Finite Element Methods for Flow Problems with Moving Boundaries and Interfaces[END_REF] is a mesh based example of interface tracking method. Particle methods, such as those of Koshizuka et al [START_REF] Koshizuka | A particle method for incompressible viscous flow with fluid fragmentation[END_REF] and Violeau and Issa [START_REF] Violeau | Numerical modelling of complex turbulent free-surface flows with the SPH method: an overview[END_REF] are examples of smoothed particle hydrodynamics (SPH) methods to the simulation of free-surface problems. However, these methods still present a high computational cost since they need to compute the interaction between the particles using search algorithms. As a compromise between the advantages offered by mesh based and meshless methods, Del Pin et al present in [START_REF] Del Pin | The ALE/Lagrangian particle finite element method: a new approach to computation of free-surface flows and fluid-object interactions[END_REF] the particle finite element method (PFEM) applied to free-surface flows. In this method the critical parts of the continuum are discretized with particles while the remaining parts are treated by a Lagrangian finite element formulation. Another technique mixing Lagrangian and Eulerian flavors was proposed in [START_REF] Takizawa | Computation of free-surface flows and fluid-object interactions with the CIP method based on adaptive meshless Soroban grids[END_REF] by Takizawa et al. In this work the authors enhanced the Constrained Interpolation Profile method (CIP) for solving hyperbolic equations with a meshless Soroban grid. The resulting formulation was used to treat fluidobject and fluid-structure interaction in the presence of freesurfaces.

As a cost effective alternative to interface tracking methods, interface capturing methods have emerged. Interface capturing methods are Eulerian in their concept, thus they rely on a unique and fixed computational grid to capture the interface evolution. In this class of methods the interface is represented by a scalar function which marks the regions filled with the fluids involved. In other words, the interface position is implicitly captured in a scalar marking function value and the interface evolution is determined by the additional cost of solving an advection equation for the marker. As opposed to interface tracking methods, interface capturing methods require little effort to represent all complicated features of moving interfaces. Additionally, the parallel implementation and postprocessing of interface capturing methods are straightforward. The main drawback of interface capturing methods is the need to average the fluid properties at the interface cells (elements) due to the discontinuity of the Eulerian representation of the interface. Moreover, the accuracy and computational cost of interface capturing methods are typically associated to grid resolution, properties of the marking function chosen to represent the interface and numerical methods for solving the fluid flow and marking function advection. The well known volume-of-fluid scheme (VOF), firstly proposed by Hirt and Nichols [START_REF] Hirt | Volume of Fluid (VOF) Methods for the Dynamics of Free Boundaries[END_REF] for Cartesian grids, is an interface capturing technique which employs a step function ranging from 0 to 1 to represent the fraction of fluid within the grid cells. In this sense, the interface is implicitly represented by the partially filled cells. The main issues associated to VOF methods include the difficulty in advecting a discontinuous step function and the accurate modeling of surface tension effects. Level set methods [START_REF] Sethian | Level Set Methods and Fast Marching Methods: Evolving Interfaces in Computational Geometry[END_REF] implement free-surface flows in a different manner than VOF by changing the marking function employed to represent the interface. Therefore, the fluids are associated to the range of the distance function signs while the interface is implicitly represented by the zero level set. However, the level set method suffers when the distance function looses its properties and must be rebuilt. In fact, the success of level set method lies in its ability of building and keeping a signed distance function without losing its properties. The enhanced-discretization interface-capturing method, firstly proposed by Tezduyar in [START_REF] Tezduyar | Enhanced-Discretization Interface-Capturing Technique[END_REF] and the work of Lohner et al presented in [START_REF] Lohner | On the Simulation of Flows with Violent Free-Surface Motion[END_REF], are both examples of unstructured grid formulations based on the finite element method to solve free-surface flows using a VOF marking function.

In this work we use our VOF edge-based stabilized finite element solver [START_REF] Elias | Stabilized Edge-Based Finite Element Simulation of Free-Surface Flows[END_REF] and [START_REF] Elias | Computational Techniques for Stabilized Edge-Based Finite Element Simulation of Free-Surface Flows[END_REF], to deal with complex wave phenomena. The main characteristics of our solver are: Streamline-Upwind/Petrov-Galerkin (SUPG), Pressure-Stabilizing/Petrov-Galerkin (PSPG) and Least Squares Incompressibility Constraint (LSIC) stabilized finite element formulation; implicit time marching scheme with adaptive time stepping control; advanced Inexact Newton solvers; edge-based data structures to save memory and improve performance; support to message passing and shared memory parallel programming models; and large eddy simulation extensions using a classical Smagorinsky model.

GOVERNING EQUATIONS

Incompressible Fluid Flow Let

sd n Ω ⊂ » be the spatial domain, where n sd is the number of space dimensions. Let Γ denote the boundary of Ω . We consider the following velocity-pressure formulation of the Navier-Stokes equations governing the incompressible flow of two immiscible fluids within an Arbitrary Lagrangian-Eulerian frame [START_REF] Donea | Finite Element Method for Flow Problems[END_REF]:

( ) on 0, mesh f t t ρ   ∂     + - ⋅ - -⋅ = Ω ×         ∂ u u u u f 0 ∇ ∇ σ ∇ ∇ σ ∇ ∇ σ ∇ ∇ σ (1)
0 on 0, f t   ⋅ = Ω ×   u ∇ ∇ ∇ ∇ (2)
where ρ and u are the density and velocity, mesh u is the mesh velocity, f is the body force vector carrying the gravity acceleration effect and σ σ σ σ is the stress tensor given as ( )

, p p = -+ u I T σ σ σ σ ( 3
)
where p is the pressure, I is the identity tensor, T is the deviatoric stress tensor

2 ( ) µ = T u ε ε ε ε (4) and ( ) u ε ε ε ε
is the strain rate tensor defined as

( ) 1 ( ) ( ) 2 T = + u u u ε ∇ ∇ ε ∇ ∇ ε ∇ ∇ ε ∇ ∇ (5)
In the present work a large eddy simulation (LES) approach to turbulence is considered by the use of a classic Smagorinsky turbulence model [START_REF] Smagorinsky | General Circulation Experiments with Primitive Equations, part i: the Basic Experiment[END_REF]. In this model, the viscosity µ is augmented by a subgrid-scale viscosity SGS µ proportional to the norm of the local strain rate tensor and to a filter width h defined here as the cubic root of the element volume, ( )

2 SGS ( ): ( ) S C h µ ρ = u u 2ε ε 2ε ε 2ε ε 2ε ε (6)
where S C is the Smagorinsky constant, taken as 0.1.

The essential and natural boundary conditions associated with equations ( 1) and ( 2) can be imposed at different portions of the boundary Γ and represented by,

on g Γ u g = = = = ( 7
)
on h ⋅ Γ n h σ = σ = σ = σ = (8)
where g Γ and h Γ are complementary subsets of Γ .

Let us assume that we have some suitably defined finitedimensional trial solution and test function spaces for velocity and pressure, h

S u , h V u , h p S and h h p p V S =
. The finite element formulation of equations ( 1) and ( 2) using SUPG and PSPG stabilizations for incompressible fluid flows can be written (see Tezduyar and Osawa [START_REF] Tezduyar | Finite Element Stabilization Parameters Computed from Element Matrices and Vectors[END_REF]) as follows:

Find h h S ∈ u u and h h p p S ∈ such that h h V ∀ ∈ u w and h h p q V ∀ ∈ : ( ) ( ) ( ) ( ) ( ) ( ) 1 1 : , , 1 el e el
e h h h h h h h h h n h h SUPG PSPG e h h h h e n h LSIC e h h mesh h h mesh h h mesh d t p d d q d q p d t Ω Ω Γ Ω = Ω = Ω ∂ ⋅ + ⋅ - Ω ∂ Ω - ⋅ Γ + ⋅ Ω + ⋅ + ∂ ⋅ + ⋅ -⋅ - Ω ∂ + ⋅ - - -             +                           ∫ ∫ ∫ ∫ ∑ ∫ ∑ ∫ u w u u f w u w h u w u u u f w u u u u u ρ τ ρ τ ρ ρ τ ρ ρ ∇ ∇ ∇ ∇ ε σ ∇ ε σ ∇ ε σ ∇ ε σ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ σ ∇ ∇ σ ∇ ∇ σ ∇ ∇ σ ∇ ∇ ∇ ∇ 0 h e ⋅ Ω = u ∇ ∇ ∇ ∇ (9)
In the above equation the first four integrals on the left hand side represent terms that appear in the Galerkin formulation of the problem (1)-( 8), while the remaining integral expressions represent the additional terms which arise in the stabilized finite element formulation. Note that the stabilization terms are evaluated as the sum of element-wise integral expressions, where n el is the number of elements in the mesh. The first summation corresponds to the SUPG term and the second to the PSPG term. We have evaluated the SUPG and PSPG stabilization parameters according to Tezduyar et al. [START_REF] Tezduyar | Incompressible Flow Computations with Stabilized Bilinear and Linear Equal-Order-Interpolation Velocity-Pressure Elements[END_REF], as follows:

1 2 2 2 2 2 4 9 h h mesh SUPG PSPG h h ν τ τ -     -           = = +                        u u ( 10
)
here u h is the local velocity vector, and ν is the kinematic viscosity.

In Equation ( 9), the last summation is the Least Squares Incompressibility Constraint (LSIC) term [START_REF] Tezduyar | Finite Element Methods for Flow Problems with Moving Boundaries and Interfaces[END_REF], added to prevent spurious oscillations at high Reynolds number flows. The LSIC stabilization parameter is,

2 h h mesh LSIC h τ - = u u (11)
The discretization of Equation ( 9) leads us to a nonlinear system of equations to be solved at each time step.

Interface Capturing

In volume-of-fluid (VOF) method [START_REF] Hirt | Volume of Fluid (VOF) Methods for the Dynamics of Free Boundaries[END_REF], a scalar marking function can be employed to capture the position of the interface between the fluids by simply using the fluids fraction relationship.

The volume-of-fluid can be stated as: assuming the value 1 in regions filled with fluid A, e.g., water, and the value 0 in regions filled with fluid B, e.g., air, the position of the fluid interface will be defined by the isovalue contour

c ( , )= t φ φ x , where c 0,1 φ   ∈   . The value c =0.5
φ is usually assumed.

Finally, the function ( ) φ x is driven by a velocity field u satisfying the following transport equation, given in conservative form as,

( ) (
)

0 mesh t φ φ ∂ + ⋅ - = ∂ u u ∇ ∇ ∇ ∇ (12)
In the volume-of-fluid (VOF) formulation the fluid density and viscosity, employed in the fluid flow solution, are interpolated across the interface as follows:

( , ) 1 ( , )

B A t t ρ φ ρ φ ρ   = + -   x x (13)
( , ) 1 ( , )

B A t t µ φ µ φ µ   = + -   x x (14)
where subscripts A and B denote the values corresponding to each fluid.

The finite element formulation of Equation ( 12) can be written as follows: Find

h h S φ φ ∈ , such that, h h w V φ ∀ ∈ , ( ) ( ) ( ) (
)

1
Ω Ω = Ω = Ω ∂ + - ⋅ Ω ∂ ⋅ Ω ∂ + - + ⋅ Ω ∂ ⋅ + ⋅ + ⋅ Ω =              +                    ∫ ∫ ∑ ∫ ∑ ∫ u u u u u u u φ φ φ φ τ φ φ δ φ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ ∇ (15)
where h S φ and h V φ are standard test and weight finite element spaces. The first two integrals represent the Galerkin formulation of Equation ( 12), while the first element-wise summation represents the SUPG and the second summation term is the nonlinear discontinuity-capturing term, useful when sharp gradients and/or boundary layers are present in directions other than the streamlines [START_REF] Galeão | A Consistent Approximate Upwind Petrov-Galerkin Method for Convection-Dominated Problems[END_REF]. The evaluation of SUPG τ and δ stabilization terms follow the definitions described in [START_REF] Brooks | Streamlineupwind/Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equation[END_REF] and [START_REF] Galeão | A Consistent Approximate Upwind Petrov-Galerkin Method for Convection-Dominated Problems[END_REF] respectively. The discretization of Equation ( 15) leads us to a non-linear ordinary differential equation system, due to the discontinuity-capturing term. In this work, we have adopted the YZ β discontinuity capturing term as proposed in [18] where δ is computed as:

1 2 2 3 1 1 1 ( ) ( ) 2 h h e e j j h R x β β φ δ φ φ φ φ - - - =       ∂       =          ∂         ∑ (16)
where ( ) e R φ is the element residual of Equation ( 15)

( ) ( )
( ) h e h h h h mesh R t φ φ φ ∂ = + ⋅ - ∂ u u ∇ ∇ ∇ ∇ (17)
Note that if 1 β = and the reference value 1 φ = , the YZ β discontinuity capturing term renders to the Consistent Approximated Upwind (CAU) method [START_REF] Galeão | A Consistent Approximate Upwind Petrov-Galerkin Method for Convection-Dominated Problems[END_REF]. In order to avoid non physical results, values lying outside the range [0,1] were truncated with the following function:

min max , 0 ,1 h h φ φ     =     (18)

SOLUTION PROCEDURE

The computational solution kernels consist of predictor multicorrector time integration schemes as described in [START_REF] Brooks | Streamlineupwind/Petrov-Galerkin Formulations for Convection Dominated Flows with Particular Emphasis on the Incompressible Navier-Stokes Equation[END_REF], [START_REF] Franca | Stabilized Finite Element Methods: ii. The incompressible Navier-Stokes equations[END_REF] for both incompressible fluid flow and interface transport equations. The generalized trapezoidal rule is employed in the time discretization. The nonlinearities due to the convective term on the Navier-Stokes equation are treated by an Inexact Newton-GMRES algorithm as described in Elias et al [START_REF] Elias | Parallel Edge-Based Solution of Viscoplastic Flows with the SUPG/PSPG Formulation[END_REF]. In this solution algorithm, at the beginning of the nonlinear iterations in each time step, the algorithm computes large linear tolerances, producing fast nonlinear steps. As the iterations progress towards the solution, the inexact nonlinear method adapts the GMRES tolerances to reach the desired accuracy. A nodal-block diagonal preconditioner is employed for the flow and a simple diagonal preconditioner for the marker. Moreover, for both the fluid flow and the marking function advection we use an adaptive time stepping procedure based on a Proportional-Integral-Derivative (PID) controller (see [START_REF] Valli | Control Strategies for Timestep Selection in FE Simulation of Incompressible Flows and Coupled Reaction-Convection-Diffusion Processes[END_REF]) for further details). Most of the computational effort spent in this solution procedure is due to the matrix-vector products within the GMRES driver for both flow and marker. To improve the computational efficiency with respect to standard element-byelement and sparse matrix-vector storage schemes, we adopt an edge-based data structure in order to minimize indirect memory addressing, diminish floating point operation counts (flops) and memory requirements, as described in Elias et al [START_REF] Elias | Inexact Newton-Type Methods for the Solution of Steady Incompressible Viscoplastic Flows with the SUPG/PSPG Finite Element Formulation[END_REF] and Coutinho et al [START_REF] Coutinho | Performance Comparison of Data-Reordering Algorithms for Sparse Matrix-Vector Multiplication in Edge-Based Unstructured Grid Computations[END_REF] for both the Navier-Stokes equations and the marking function advection. Further computational gains are obtained from data preprocessing performed by the EdgePack library -a package to improve cache reutilization based on reordering and grouping techniques [START_REF] Martins | Edgepack: A Parallel Vertex and Node Reordering Package for Optimizing Edge-Based Computations in Unstructured Grids[END_REF]. All computational kernels are parallelized covering most of distributed and shared memory systems as well as vectorized and pipelined processors [START_REF] Elias | Parallel Edge-Based Inexact Newton Solution of Steady Incompressible 3D Navier-Stokes Equations[END_REF].

Enforcing Mass Conservation

The most challenging feature for a good interface capturing method resides in its ability to preserve the mass of the species involved. According to Lohner [START_REF] Lohner | On the Simulation of Flows with Violent Free-Surface Motion[END_REF], in VOF methods the mass loss can be associated to reasons such as: interface smearing due to numerical diffusion of the step function, inexact divergence free velocity field and boundary conditions, undulations in the solution of the marking function advection. Level set methods suffer when the marking function loose its signed distance property. These problems have been reported by many authors and are still the subject of several researches. In this work we have followed the procedure proposed in [START_REF] Lohner | On the Simulation of Flows with Violent Free-Surface Motion[END_REF] and [START_REF] Elias | Stabilized Edge-Based Finite Element Simulation of Free-Surface Flows[END_REF] to overcome mass losses. In this procedure the mass lost/gained are found by comparing the expected value, composed by the initial mass plus the inlet and outlet fluxes, at the end of each time step. Therefore, the values to be added or removed are made proportional to the absolute value of the normal velocity of the interface given by ( )

h h mesh n u φ φ ∇ - ⋅ = ∇ u u ( 19
)
the amount computed from Equation ( 19) guarantees that the mass correction will act mainly in regions where the interface is moving faster, while keeping the stationary regions untouched. Therefore, the portion of mass correction corresponding to each element is and projected onto the global nodes by L 2 projection. In [START_REF] Elias | Stabilized Edge-Based Finite Element Simulation of Free-Surface Flows[END_REF], Elias and Coutinho evaluate different test problems with the help of the mass preservation algorithm showing its efficiency.

In order to guarantee the effectiveness of the mass correction method we must have an accurate procedure for computing the volume of the fluid phases. In this work we have adopted a method, described in [START_REF] Elias | Stabilized Edge-Based Finite Element Simulation of Free-Surface Flows[END_REF], which makes use of the interface explicitly determined for those elements lying on the interface. To introduce the method, firstly we need to identify the elements as filled, partially filled or empty. This can be done by computing an element key according to the nodal value of the marking function φ. Thus, if φ ≥ 0.5 we assign a value 1, otherwise, it is marked as 0. Finally, the element key is found by summing its nodal values. Note that the empty elements can be skipped from the computation since the void region will be the difference between the domain and the filled region. Moreover, for elements partially filled, the filling volume can be found from the sub-tetrahedra and wedges formed by the interface plane cutting the original element.

TEST PROBLEMS

Progressive Wave

This case represents the simulation of generation and propagation of a progressive regular wave, from a piston-type wave maker. The computational domain is presented in Figure 1. It has 120 meters length, 20 meters height and 0.5 meters width. Figure 1 also presents the initial condition and initial mesh configuration.

Figure 1. Initial condition and mesh -67,320 nodes and 254,052 tetrahedral elements.

The movement is prescribed by a periodic function at the left face, simulating a piston-type wave maker. To avoid mesh entangling we introduce a buffer region where the periodic function describing the nodes displacement suffers a damping effect, given by,

cos( ) Ae t - = + ω ϕ 0 x x ( 20
)
where A is the piston displacement amplitude, x and 0

x are position coordinates, ω is the piston angular frequency, t refers to the time and ϕ is the phase angle. In this problem, the piston moves periodically, with amplitude 0.40 m from its original position and the oscillation period is 4 s. All nodes, from the piston up to 20 m, also oscillate gradually according to Equation [START_REF] Elias | Parallel Edge-Based Solution of Viscoplastic Flows with the SUPG/PSPG Formulation[END_REF].

The time step was fixed in 0.005 seconds and the simulation time is 40 seconds. Linear GMRES tolerance was set to 10 -3 and 10 -5 for the Navier-Stokes equations and scalar marking function equation, respectively. The Navier-Stokes Inexact-Newton flow solver converges in 3 iterations while the scalar marking function equation converges in 6 iterations. Free surface elevation at t = 0s, t = 8s, t = 16s, t = 24s, t = 32s and t = 40s.

According to linear wave theory, the wave celerity can be computed as the ratio between its length and period [START_REF] Journée | Offshore Hydromechanics[END_REF]. The analytical and numerical results for the wave celerity are 6.17 Three positions were chosen to track pressure and velocity fields. Figure 4 shows the monitored points. The piston is in the origin of coordinate system. The points are at x coordinate 60.0 meters and y coordinates 2.5, 5.0 e 7.5 meters.

The analytical solution was calculated from equations given by the potential formulation for regular progressive wave propagation, known as linear theory [START_REF] Journée | Offshore Hydromechanics[END_REF]. The velocity components and pressure equations are shown below:

( ) ( ) ( ) cosh cos sinh a k h y u kx t ky ζ ω ω +     = - (21) ( ) ( ) ( )
sinh sin sinh a k h y v kx t ky ζ ω ω +     = - (22)
ρ ρ ζ ω +     = - + - (23)
In equations 21, 22 e 23, u and v are the velocity components, p is pressure and x and y are position coordinates. The wave amplitude and frequency are ζa e ω, respectively. Water depth is h and k is the wave number. The wave number is calculated from dispersion relationship, given by equation 24 .

2 tanh( )

k kh g ω = ( 24
)
where g denotes gravity.

Figure 5 shows the pressure field for all monitored points. The numerical simulation is compared to the linear theory. The transient phase of the pressure time history is due to the initial condition. In all three points, after the initial transient, a good agreement between the numerical and analytical pressure solutions is observed.

Mean pressure is presented in table 1 and the dynamic pressure amplitudes are presented in table 2. The values were taken in the last 5 minutes of simulation. Note that, there is a good agreement of the frequency, but the amplitudes are not so close. On the other hand, it is important to observe that linear theory can not cope with viscosity effects and non-linear free surface effects.

Progressive wave interaction with a fixed cylinder

This case refers to the interaction of a progressive wave with a fixed circular cylinder placed on the center of the domain, as shown in Figure 12. The ground based circular section structure has a depth of 15 m (5m in the air and 10 m under the water) and its radius is 5m. All dimensions of progressive wave original case domain were kept, except the width. This dimension is taken as 50 m, being a good choice to diminish the effects of lateral walls on the wave-cylinder interaction. The piston, represented by the face at y = 0 m, moves cyclically in x-direction with a 0.40 m amplitude and period of 4 seconds.

The fluid properties are: ρ water = 1000.0 kg/m 3 , µ water = 0.001 kg/(m s), ρ air = 1.0 kg/m 3 , µ air = 0.001 kg/(m s). The air viscosity is considered equal to the water in order to avoid unstable regions in the air phase.

Figure 12. Progressive wave over a cylinder. Free surface initial condition and domain.

Figure 13 shows the computational domain and the partitions for parallel processing. The computation was performed in parallel (MPI) in a SGI Altix XE 1200, employing 8 of its 88 Xeon cores (2.66 GHz). The maximum inexact-Newton linear tolerance was set to 10 -3 while the linear tolerance for the marking function equation was 10 -5 . The nonlinear loops were stopped after the relative residual as well as the relative step increment decreased 3 orders of magnitude. In addition, a fixed time step of 0.005 was employed. For 50 seconds of simulation it was spent about 70 hours of CPU time. Figure 15 presents the kinetic energy of the fluid. Since the fluid is initially at rest, the kinetic energy is introduced by the piston. The oscillatory movement of the piston, the wavecylinder and wave-opposite face interactions are captured in the kinetic energy oscillatory behavior. Figure 16 shows the water mass preservation with respect to the initial water mass, and reveals how efficiently the mass correction algorithm works. Note that most of the disturbances occur at the initial stage, up to 1 second, which coincides with the moment that the piston pushes the fluid, which is initially at rest. It is important to remind that the mass correction algorithm is based on the velocity field normal to the free surface, making the correction effects more pronounced when the interface is moving faster. The maximum computed water loss was around 0.005% which is one order of magnitude smaller than the adopted nonlinear stopping criteria. It means that material losses are negligible. The history of the horizontal hydrodynamic force exerted by the fluid over the cylinder is illustrated by Figure 17. The amplitude of the linear force evaluated using the MacCamy and Fuchs analytical formulation is equal to 612,299 N [START_REF] Sarpkaya | Mechanics of Wave Forces on Offshore Structures[END_REF] and by the WAMIT® code is 591,257 N, [START_REF] Wamit | Wamit 6.1 User Manual[END_REF]. These results, based on linear diffraction theory, are close to the numerical results obtained by the program. The causes of the changes in the time history of the horizontal force (between 37s and 45s for example) could at first be imputed to the wall reflections at the end of the channel, but deserves more investigation

CONCLUSIONS

In this work we have successfully extended our parallel incompressible flow solver to problems involving wave interactions with moving physical boundaries. This goal was reached by introducing Arbitrary Lagrangian-Eulerian formulation into the stabilized finite element formulation. We have validated the present formulation simulating a progressive wave generated from a piston-like wave maker and its interaction with a vertical circular cylinder. The main conclusions are that some important parameters in progressive wave behavior, as celerity, wavelength and height, are in good agreement with existing linear theoretical values. In the problem of the progressive wave interaction with a fixed circular cylinder, we have computed the horizontal force over the cylinder and the kinetic energy of the fluid, which is introduced by the piston movement. We noticed the linear force given by the MacCamy and Fuchs analytical formulation and the WAMIT® code are close to the numerical results obtained in the simulation. Furthermore, we verified in this case that the maximum loss computed is around 0.005%, which evidences the efficiency of global mass conservation procedure.

Figure 2 .

 2 Figure 2. Generation and propagation of a progressive wave.Free surface elevation at t = 0s, t = 8s, t = 16s, t = 24s, t = 32s and t = 40s.

  m/s and 6.42 m/s, which are in good agreement. Beyond the wave celerity, other comparison parameters were verified. The wave length and height provided by the linear wave generator theory ([START_REF] Hughes | Physical Models and Laboratory Techniques in Coatal Engineering[END_REF] and[START_REF] Dean | Water Wave Mechanics for Engineers and Scientists[END_REF]) are 24.67 m and 1.49 m, respectively. They compare well with the values obtained from the simulation, 24.70 m and 1.58 m.

Figure 3 .

 3 Figure 3. Progressive wave displacement.

Figure 4 .

 4 Figure 4. Monitored points. x = 60m and y = 2.5m, 5.0m, 7.5m.

Figure 5 .

 5 Figure 5. Pressure on monitored points.

y = 2

 2

Figure 6 .Figure 7 .Figure 8 .Figure 9 .Figure 10 .Figure 11 .

 67891011 Figure 6. Velocity. Component u. x = 60.0 and y = 2.5 m.

  Figure13shows the computational domain and the partitions for parallel processing. The computation was performed in parallel (MPI) in a SGI Altix XE 1200, employing 8 of its 88 Xeon cores (2.66 GHz). The maximum inexact-Newton linear tolerance was set to 10 -3 while the linear tolerance for the marking function equation was 10 -5 . The nonlinear loops were stopped after the relative residual as well as the relative step increment decreased 3 orders of magnitude. In addition, a fixed time step of 0.005 was employed. For 50 seconds of simulation it was spent about 70 hours of CPU time. Figure14presents five snapshots representing the time solutions at 10, 20, 30, 40 and 50 seconds.

Figure 13 .Figure 14 .

 1314 Figure 13. Progressive wave over a cylinder. Computational domain: Tetrahedral finite elements unstructured mesh with 195,426 nodes and 1,140,159 elements and partitioning for 8 processors.

Figure 15 .Figure 16 .

 1516 Figure 15. Kinetic energy of progressive wave.

Figure 17 .

 17 Figure 17. Force exerted by the progressive wave over the cylinder.

Table 1 .

 1 Mean pressure.

y =

2.5 m y = 5.0 m y = 7.5 m Linear amplitude of dynamic pressure (Pa)

			1371,88	2185,01	3914,47
	EdgeCFD dynamic pressure (Pa) amplitude	of	1953,80	2941,15	4831,95

Table 2 .

 2 Amplitude of dynamic pressure.

	Velocity field for each one of three monitored points is shown
	in Figures 6 -11. As it happens in pressure solution, the
	simulation needs an adjustment time to converge to steady-state
	solution.

ACKNOWLEDGMENTS

The authors would like to thank the financial support of the PETROBRAS Research Center; the MCT/CNPq Brazilian Council for Scientific Research and MCT/FINEP CT-AQUAVIÁRIO-2006. The Center for Parallel Computations (NACAD) and the Ocean Basin (LabOceano) at the Federal University of Rio de Janeiro provided the computational resources for this research.


E Tassi

Hamiltonian reduced fluid model for plasmas with temperature and heat flux anisotropies

We derive a Hamiltonian fluid model for strongly magnetized plasmas describing the evolution, for an arbitrary number of species, of density, velocity and electromagnetic fluctuations as well as of the temperature and heat flux fluctuations associated with motions parallel and perpendicular to the direction of a background magnetic field. The model is derived as a reduction of the infinite hierarchy of equations obtained by taking moments, with respect to Hermite-Laguerre polynomials in velocitymagnetic moment coordinates, of a Hamiltonian drift-kinetic system. A closure relation directly coupling the heat flux fluctuations in the directions parallel and perpendicular to the background magnetic field is shown to provide a fluid reduction which preserves the Hamiltonian character of the parent drift-kinetic model. An alternative set of dynamical variables is found, in terms of which the Poisson bracket of the fluid model takes a simple direct sum structure and permits an easy identification of the Casimir invariants. Such invariants, in the limit of translational symmetry with respect to the direction of the background magnetic field, turn out to be associated with Lagrangian invariants of the fluid model. It is pointed out that the coupling between parallel and perpendicular heat flux evolutions introduced by the closure, is necessary in order to guarantee the existence of a Hamiltonian structure with a Poisson bracket obtained as extension of a Lie-Poisson bracket.

I. INTRODUCTION

Kinetic theories provide an accurate description of continuous media such as gases and plasmas. Such models consist of evolution equations for distribution functions defined on particle phase space, possibly coupled self-consistently to equations governing the evolution of fields affecting the particle dynamics. Paradigmatic examples of kinetic theories are the Vlasov-Maxwell system, for plasmas in the presence of electromagnetic fields and the Vlasov-Poisson system, for plasmas in electrostatic fields, or for particle distributions subject to their own gravitational potential. A less accurate, but computationally less demanding, level of description of continuous media is provided by fluid models, which, in the Eulerian representation, consist of evolution equations for fluid dynamical variables depending only on spatial coordinates, such as particle densities, momenta, temperatures, etc..coupled with governing equations for the fields. The fluid dynamical variables usually correspond to moments of the kinetic distribution functions, obtained by integrating over velocity space the distribution functions, multiplied times non-negative integer powers of the velocity coordinates. Such powers determine the order of the moments.

Evolution equations for the fluid moments can then be derived directly from the evolution equations for the distribution functions. Because the evolution equation for each moment typically depends also on moments of higher order, this leads to an infinite hierarchy of coupled evolution equations for the fluid variables. In order to get a closed system for a finite number of fluid variables a reduction procedure, usually referred to as to closure, is required. Imposing a closure amounts to assume that a moment of a given order is determined by lower-order moments. Such lower-order moments become the dynamical variables of the fluid model which is then no longer equivalent to the original kinetic description but evolves only a finite number of dynamical variables depending only on the spatial coordinates.

Clearly, the choice of the closure relating a given moment to lower order moments, is crucial in determining the characteristics of the fluid model.

In the plasma physics literature several examples of closure relations exist, which satisfy different criteria. For instance in Refs. [1,2] the closure is such to guarantee energy conservation for the resulting fluid model. Fluid models retaining kinetic effects such as Landau damping have been presented in Refs. [3][4][5][6][7]. Non-dissipative closures preserving time-reversal symmetry of the original kinetic equation have been proposed in Ref. [8]. A fluid closure based on entropy production arguments has been proposed in Ref. [9].

In spite of a considerable effort devoted to identifying fluid closures, little is known about the consequences of such closures on some structure of the model equations. In particular, in many cases it is not known whether, in the limit where the parent kinetic model has a Hamiltonian structure (which is the case for instance for Vlasov-Maxwell and Vlasov-Poisson systems [10,11]), the closure is such to preserve the Hamiltonian structure in the fluid model as well.

This should be the case unless one voluntarily intends to introduce dissipation in the fluid model. In the mathematical physics literature, results related to this problem date back at least to Refs. [12,13], where the connection between the Poisson bracket for the Vlasov equation and the Poisson brackets in terms of the moments had been elucidated.

The connection between the Hamiltonian dynamics of the Vlasov moments and integrability has been discussed in Ref. [14]. The Hamiltonian structure of the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy of moment equations, on the other hand, had been presented in Ref. [15]. In Ref. [16] Hamiltonian structures for finite fluid reductions originating from the Vlasov-Poisson system are given. Some recent results in this context have also been obtained with regard to the identification of Hamiltonian fluid closures for three [17] and four-field [18] models obtained from the Vlasov-Poisson system, and two-field models obtained from drift-kinetic [19,20] and gyrokinetic systems [21]. A closure relation, preserving a Hamiltonian structure in a fluid model derived from a drift-kinetic system and retaining an arbitrary number of moments, was presented in Ref. [22].

The results presented in this paper belong to this line of research, devoted to identifying closures that preserve the Hamiltonian structure of a parent Hamiltonian kinetic theory. In particular, in the present work, we consider, as parent kinetic theory, a system of drift-kinetic equations coupled with electromagnetic fields through a quasi-neutrality condition and Ampère's law. Drift-kinetic equations (see, e.g. Ref. [START_REF] Hazeltine | Plasma Confinement[END_REF]) are of interest for low-frequency phenomena in plasmas where an intense component of the magnetic field is present. Indeed, under these assumptions, the fast gyration motion of particles around the intense component of the magnetic field can be averaged out because it takes place on too short time scales, when compared to the frequencies of interest. As a result, when compared to the six-dimensional phase space of the Vlasov distribution function, the drift-kinetic equation possesses the advantage of being defined on a reduced five-dimensional phase space, where the gyration angle coordinate has been ignored because irrelevant for low-frequency phenomena. The Hamiltonian drift-kinetic model adopted in this paper can be derived from the gyrokinetic model treated in Ref. [2] upon applying a number of simplifying assumptions consisting of neglecting finite Larmor radius effects, background inhomogeneities and dissipative terms.

The purpose of this paper is to show that a closure relation exists so that it is possible to derive, from such Hamiltonian drift-kinetic system, a Hamiltonian fluid model evolving six field variables for each particle species.

More precisely, indicating with N the number of particle species present in the plasma, we can construct a 6N -field Hamiltonian fluid model evolving, for each species and self-consistently with the electromagnetic fields, the density, the canonical fluid momentum, temperature and heat flux fluctuations originated from both the parallel and perpendicular particle motions. Such field variables are proportional to moments of a generalized distribution function taken with respect to Hermite and Laguerre polynomials of normalized velocity along the strong component of the magnetic field and of magnetic moment, respectively.

With respect to recent related results, such as those of Refs. [19][20][21][START_REF] Grasso | [END_REF], this model indeed accounts additionally also for moments taken with respect to the particle magnetic moment, thus yielding a fluid model that can describe also phenomena related to temperature and heat flux anisotropies.

The existence of the Hamiltonian structure in the resulting fluid model, besides assuring the absence of spurious dissipative effects and an unambiguous identification of the conserved total energy, can provide also information about further conservation laws of the system. Indeed, as is the case for the Hamiltonian structure of continuous media in the Eulerian description, also such fluid model turns out to possess a Poisson bracket of noncanonical type (see, e.g.

Refs. [25,[START_REF] Marsden | Introduction to Mechanics and Symmetry[END_REF]). Poisson brackets of this type are characterized by the presence of Casimir invariants, which constrain the dynamics of the system. Thanks to the Hamiltonian structure, we also provide Casimir invariants of the model, which become particularly relevant in the two-dimensional limit, when the system is translationally invariant with respect to the direction of the intense component of the magnetic field. In this limit, infinite families of Casimir invariants emerge. The identification of the Casimir invariants is made possible thanks to a change of variables in terms of which we are able to cast the Poisson bracket of the model in its simplest form, unveiling its direct sum structure. This alternative set of variables, which are referred to as to normal fields, turns out to be related to the existence of 6N Lagrangian invariants for the model in its two-dimensional limit.

The paper is organized as follows. In Sec. II the parent drift-kinetic model is introduced and its Hamiltonian structure is presented. In Sec. III we introduce the hierarchy of fluid equations, in terms of moments with respect to Hermite-Laguerre polynomials, and present the closure relation. The Hamiltonian structure of the resulting fluid model is presented in Sec. IV, whereas in Sec. V we introduce the normal field variables and present the Casimir invariants. We conclude in Sec. VI and provide the A where the eigenvalues of the matrix W (1) s are explicitly determined.

II. HAMILTONIAN PARENT DRIFT-KINETIC MODEL

In a Cartesian coordinate system (x, y, z) we consider the following system of drift-kinetic evolution equations :

∂g s ∂t + c B φ - v c A, g s + v ∂ ∂z g s + q s F s T s φ - v c A = 0, (1)
complemented by the quasi-neutrality relation

n 0 M c 2 B 2 ∆ ⊥ φ = - s q s dW s g s (2)
and by Ampère's law

c 4π ∆ ⊥ A - s q 2 s n 0 m s c A = - s q s dW s vg s . (3)
The index s = 1, • • • , N labels the N particle species present in the plasma and, in order to simplify the notation, the sum over s always means that the index s goes from 1 to N . The dynamical field variable g s is defined as

g s (x, y, z, v, µ s , t) = fs (x, y, z, v, µ s , t) + q s T s v c F s (v, µ s )A(x, y, z, t), ( 4
)
where fs is the perturbation of the distribution function for the s-th species with respect to the equilibrium Maxwellian

F s (v, µ s ) = n 0 m s 2πT s 3/2 exp - v 2 2v 2 ts - µ s B m s v 2 ts , (5)
with v indicating the velocity coordinate along z, which is the direction of the magnetic guide field, with µ s = m s (v 2 x + v 2 y )/(2B) indicating the magnetic moment of s-th particle species and with the constants q s , m s , T s , v ts , c and n 0 indicating the particle charge and mass, the equilibrium temperature for the s-th species, the thermal speed v ts = T s /m s , the speed of light and the equilibrium density, respectively. The magnetic field is given by B(x, y, z, t) = ∇A(x, y, z, t) × ẑ + B ẑ, with constant B, whereas φ = φ(x, y, z, t) indicates the electrostatic potential.

In Eq. ( 2) the constant M is defined as M = s m s . We indicated with the symbol ∆ ⊥ the perpendicular Laplacian

operator ∆ ⊥ = ∂ xx + ∂ yy , whereas the bracket [ , ] is defined as [f, g] = ∂ x f ∂ y g -∂ y f ∂ x g
, for two functions f and g.

Finally, the volume element dW s is defined as dW s = 2π(B/m s )dµ s dv. The spatial coordinates x, y and z vary over the 3-dimensional torus T 3 , whereas v ∈ R, µ s ≥ 0 for all s and t ≥ 0. The perturbations of the distribution functions fs are assumed to decay to zero sufficiently rapidly as v → ±∞ so that boundary terms vanish when integrating by parts with respect to v. The electrostatic potential φ is also assumed to have zero mean value, so that the Laplacian operator can be inverted in Eq. ( 2).

The model consisting of Eqs. ( 1), ( 2) and ( 3) can be derived from the gyrokinetic model described for instance in

Ref. [2], by taking the limit of vanishing Larmor radius, and in the absence of background magnetic inhomogeneities and dissipation. The model assumes the common δf approximation which is valid for plasmas close to equilibrium, i.e. fs F s and with weak variations along the guide field, i.e. ∂/∂z ∂/∂x and ∂/∂z ∂/∂y. With regard to Eq. ( 2) we observe that, with respect to the quasi-neutrality relation adopted in Ref. [2], the finite Larmor radius contribution depending on the perpendicular temperature has been neglected, whereas the leading order contribution polarization term has been retained. A concrete example where this assumption is justified corresponds to the case where two species are present in the plasma, with electrons retaining temperature and heat flux anisotropies, and ions following an isothermal description.

The parent drift-kinetic model ( 1), ( 2) and (3) possesses a Hamiltonian structure which is a straighforward extension of the Hamiltonian structure of the model treated in Ref. [START_REF] De Blank | Proceedings of the 28th Conference on Contr. Fusion and Plasma Phys[END_REF]. Such Hamiltonian structure consists of the Hamiltonian functional

H(g 1 , • • • , g N ) = s 1 2 d 3 xdW s T s g 2 s F s + q s g s φ - v c A , (6)
and of the noncanonical Poisson bracket

{F, G} = s {F, G} s = - s d 3 xdW s c q s B g s [F gs , G gs ] + v F s T s F gs ∂G gs ∂z , (7)
where the subscripts on the functionals F and G indicates functional derivatives. Indeed, the generic expression for a Hamiltonian system describing the dynamics of N field variables χ 1 , • • • , χ N is given by

∂χ i ∂t = {χ i , H}, i = 1, • • • , N (8)
where

H = H(χ 1 , • • • , χ N ) is the Hamiltonian functional and { , } a Poisson bracket. One obtains that, with N = N , by substituting χ 1 , • • • , χ N with g 1 , • • • , g N , H
with the Hamiltonian functional (6) and using as Poisson bracket the one indicated in Eq. ( 7), yields namely the drift-kinetic evolution equations (1).

We remark that, whereas the Poisson bracket (7) consists of the direct sum of Poisson brackets each of which is associated with a given species, the coupling between the different species in the dynamics comes from the last two terms in the Hamiltonian ( 6), containing φ and A, which involve all species according to ( 2) and ( 3), respectively. The first term in the Hamiltonian ( 6), on the other hand, comes from expanding the free energy functional of the original full (i.e. without the δf approximation) drift-kinetic model up to second order, about the equilibrium Maxwellian (5).

III. FLUID EQUATIONS

From the drift-kinetic equations (1) a hierarchy of evolution equations for dynamical variables depending only on the spatial coordinates x, y, z, as well as on time, can be obtained. As dynamical variables for a given species s, we choose the moments g mns defined as

g mns (x, y, z, t) = 1 n 0 √ m! dW s H m v v ts L n µ s B m s v 2 ts g s (x, y, z, v, µ s , t), (9)
where H m and L n are the Hermite and Laguerre polynomials of order m and n, respectively, with m and n non-negative integers. We recall that such polynomials can be defined in the following way

H m (x) = (-1) m e x 2 /2 d m dx m e -x 2 /2 , L n (x) = e x n! d n dx n x n e -x . (10)
The first Hermite and Laguerre polynomials, which correspond in particular to those of interest for the present analysis, are then given by

H 0 v v ts = 1, H 1 v v ts = v v ts , H 2 v v ts = v 2 v 2 ts -1, H 3 v v ts = v 3 v 3 ts -3 v v ts , L 0 µ s B m s v 2 ts = 1, L 1 µ s B m s v 2 ts = 1 - µ s B m s v 2 ts . (11)
The choice of the Hermite and Laguerre polynomials as "weight" functions, in Eq. ( 9), for the parallel velocity and magnetic moments, respectively, turns out to be practical due to the orthogonality relation

dW s H m v v ts H i v v ts L n µ s B m s v 2 ts L j µ s B m s v 2 ts F s = n 0 m!δ mi δ nj , (12)
which involves the equilibrium Maxwellian F s . The dynamical variable g s is decomposed over the basis of Hermite and Laguerre polynomials in the following way

g s (x, y, z, v, µ s , t) = +∞ m,n=0 1 √ m! H m v v ts L n µ s B m s v 2 ts F s (v, µ s )g mns (x, y, z, t), (13)
from which the expression (9) follows by making use of the relation (12).

We consider here the evolution of the moments g 00s , g 10s , g 20s , g 30s , g 01s and g 11s . Such moments are indeed also proportional to the physically most relevant dynamical variables typically adopted in anisotropic reduced fluid models for plasmas, in the presence of a strong guide field. The correspondence between the moments and the physical variables is given by :

g 00s = n s n 0 , g 10s = u s v ts + q s A m s v ts c , g 20s = 1 √ 2 T s T s , g 30s = 2 3 1 n 0 q s T s v ts , g 01s = - T ⊥s T s , g 11s = - 1 n 0 q ⊥s T s v ts , (14)
where, for each species s, n s = dW s fs are the density fluctuations, u s = (1/n 0 ) dW s v fs the velocity fluctuations,

T s = (1/n 0 ) dW s (m s v 2 -T s
) fs and T ⊥s = (1/n 0 ) dW s (µ s B -T s ) fs the parallel and perpendicular temperature fluctuations, respectively, whereas q s = (1/2) dW s (m s v 2 -3T s )v fs and q ⊥s = dW s (µ s B -T s )v fs are the parallel and perpendicular heat flux fluctuations, respectively.

By multiplying both sides of Eq. ( 1) for H m (v/v ts )L n (µ s B/(m s v 2 ts )) and integrating over dW s one obtains, for the six moments indicated in Eq. ( 14), the following set of evolution equations:

∂g 00s ∂t = - c B [φ, g 00s ] + v ts B [A, g 10s ] -v ts ∂g 10s ∂z + v 2 ts q s cT s ∂A ∂z , (15)
∂g 10s ∂t = - c B [φ, g 10s ] + √ 2 v ts B [A, g 20s ] + v ts B [A, g 00s ] - √ 2v ts ∂g 20s ∂z -v ts ∂g 002 ∂z -v ts q s T s ∂φ ∂z , (16)
∂g 20s ∂t = - c B [φ, g 20s ] + √ 3 v ts B [A, g 30s ] + √ 2 v ts B [A, g 10s ] - √ 3v ts ∂g 30s ∂z - √ 2v ts ∂g 10s ∂z + √ 2v 2 ts q s cT s ∂A ∂z , (17)
∂g 30s ∂t = - c B [φ, g 30s ] + 2 v ts B [A, g 40s ] + √ 3 v ts B [A, g 20s ] -2v ts ∂g 40s ∂z - √ 3v ts ∂g 20s ∂z , ( 18
)
∂g 01s ∂t = - c B [φ, g 01s ] + v ts B [A, g 11s ] -v ts ∂g 11s ∂z , ( 19
)
∂g 11s ∂t = - c B [φ, g 11s ] + √ 2 v ts B [A, g 21s ] + v ts B [A, g 01s ] - √ 2v ts ∂g 21s ∂z -v ts ∂g 01s ∂z , (20)
where use of the recursive relation

H m+1 (x) = xH m (x) -mH m-1 (x) has been made.
From Eqs. ( 18) and ( 20) it emerges that the system is not closed, because such two equations depend on g 40s and g 21s , which are undetermined. This leads to a closure problem, amounting to express g 40s and g 21s in terms of the six lower order moments. Recalling that we require to preserve, in the fluid model, the Hamiltonian character of the parent drift-kinetic model, we intend to show that closing the system ( 15)-( 20) by imposing the relations

g 40s = α s g 11s , g 21s = √ 2α s g 30s , (21)
with non-vanishing constant α s , leads to a fluid model with a Hamiltonian structure.

IV. HAMILTONIAN ANISOTROPIC FLUID MODEL

In order to simplify the notation we rename the moment variables in the following way:

g 0s = g 00s , g 1s = g 10s , g 2s = g 20s , g 3s = g 30s , g 4s = g 01s , g 5s = g 11s .
Applying the closure relation (21) we obtain that the anisotropic fluid model reads

∂g 0s ∂t = = - c B [φ, g 0s ] + v ts B [A, g 1s ] -v ts ∂g 1s ∂z + v 2 ts q s cT s ∂A ∂z , (22)
∂g 1s ∂t = - c B [φ, g 1s ] + √ 2 v ts B [A, g 2s ] + v ts B [A, g 0s ] - √ 2v ts ∂g 2s ∂z -v ts ∂g 0s ∂z -v ts q s T s ∂φ ∂z , (23)
∂g 2s ∂t = - c B [φ, g 2s ] + √ 3 v ts B [A, g 3s ] + √ 2 v ts B [A, g 1s ] - √ 3v ts ∂g 3s ∂z - √ 2v ts ∂g 1s ∂z + √ 2v 2 ts q s cT s ∂A ∂z , (24)
∂g 3s ∂t = - c B [φ, g 3s ] + 2α s v ts B [A, g 5s ] + √ 3 v ts B [A, g 2s ] -2α s v ts ∂g 5s ∂z - √ 3v ts ∂g 2s ∂z , ( 25
)
∂g 4s ∂t = - c B [φ, g 4s ] + v ts B [A, g 5s ] -v ts ∂g 5s ∂z , ( 26
)
∂g 5s ∂t = - c B [φ, g 5s ] + 2α s v ts B [A, g 3s ] + v ts B [A, g 4s ] -2α s v ts ∂g 3s ∂z -v ts ∂g 4s ∂z . (27)
In order to show that the system ( 22)-( 27) possesses a Hamiltonian structure, we need to show that it can be cast in the form (8). A natural choice for the 6N field variables is clearly given by χ is = g is , with i = 0, • • • , 5 and s labelling the N particle species.

A natural candidate as Hamiltonian functional is obtained directly from the expression of the Hamiltonian of the parent model ( 6) upon replacing g s with the following truncated version of the series (13):

g s (x, y, z, v, µ s , t) = F s (v, µ s ) g 0s (x, y, z, t)H 0 v v ts L 0 µ s B m s v 2 ts +g 1s (x, y, z, t)H 1 v v ts L 0 µ s B m s v 2 ts + g 2s √ 2 (x, y, z, t)H 2 v v ts L 0 µ s B m s v 2 ts + g 3s √ 6 (x, y, z, t)H 3 v v ts L 0 µ s B m s v 2 ts + g 4s (x, y, z, t)H 0 v v ts L 1 µ s B m s v 2 ts +g 5s (x, y, z, t)H 1 v v ts L 1 µ s B m s v 2 ts , (28)
which retains only the contribution coming from the moments involved in the model. By inserting the truncated series (28) into the drift-kinetic Hamiltonian (6) one obtains, with the help of the orthogonality relation (12),

H(g 01 , • • • , g 5N ) = 1 2 s d 3 x n 0 T s 5 i=0 g i 2 s + q s n 0 g 0s φ -q s n 0 v ts c g 1s A . (29)
The functional ( 29) is a conserved quantity for the system ( 22)-( 27) and can therefore be considered as a candidate

Hamiltonian. This functional corresponds, in the appropriate limit, to the total conserved energy functional of the gyrofluid model of Ref. [2].

Concerning the Poisson bracket, the form of the model equations suggests that it be of the form {F, G} = s {F, G} s , with

{F, G} s = - c q s Bn 0 5 i,j,k=0 W ij k s d 3 xg ks [F is , G j s ] - v ts n 0 T s 5 i,j=0 A ij s d 3 xF is ∂G j s ∂z . (30)
In Eq. ( 30), W ij k s and A ij s are constants and the subscripts on functionals indicate functional derivatives, so that, for instance, F is = δF/δg is . Poisson brackets of the form (30) occur very often in Hamiltonian reduced fluid models for plasmas, assuming the presence of a strong magnetic guide field. Examples of such models can be found in Refs.

[ [START_REF] Waelbroeck | [END_REF][29][30][31][32][33][34][35][36]. In order to determine the coefficients W ij k s and A ij s identifying the Poisson bracket, we proceed in the following way. First, we consider the case A ij s = 0, for all i, j = 0, • • • , 5 and look for coefficients W ij k s such that the operation

{F, G} ⊥ s = - c q s Bn 0 5 i,j,k=0 W ij k s d 3 xg ks [F is , G j s ] ( 31
)
is a Poisson bracket of its own and that, when combined with the candidate Hamiltonian (29), it yields the model equations ( 22)-( 27) in the "perpendicular" limit, i.e. when all the terms involving z-derivatives are neglected. Subsequently, inspired by the procedure described in Ref. [29], we identify the coefficients A ij s such that the bracket (30), again combined with the functional (29), yields the full model equations, and verify that the resulting bracket is a Poisson bracket.

To carry out the first part of the procedure, we recall, from Ref. [37], that the bilinear operation ( 31) is a Poisson bracket if and only if the matrices W (j) s , defined as

W (j) s ki = W ij k s , i, j, k = 0, • • • , 5, (32)
with rows labelled by k and and columns by i, are symmetric in the upper indices, i.e. W (j)

s ki = W (j)
s ik , and they all pairwise commute. Such two properties indeed are equivalent to the antisymmetry and the Jacobi identity, respectively, for the bracket (31). We remark that the Leibniz identity is automatically satisfied by a bracket of the form (31). In addition to the constraints related to the properties of a Poisson bracket, we require the coefficients W ij k s (or, equivalently, the matrices W (j) s ) to be such that the Hamiltonian evolution equations

∂g ls ∂t = {g ls , H} ⊥ s = - c q s Bn 0 5 j,k=0 W lj k s [H j s , g ks ], l = 0, • • • , 5 (33)
correspond, when using H given by Eq. ( 29), to the model Eqs. ( 22)-( 27) in the limit when terms involving zderivatives are set to zero. To this scope, it is useful to note that the functional derivatives of the Hamiltonian (29) are given by

H 0s = n 0 T s g 0s + q s n 0 φ, H 1s = n 0 T s g 1s -q s n 0 v ts c A, H is = n 0 T s g is , i = 2, • • • , 5. (34)
With this information, one finds that, in order for Eqs. (33) to match Eqs. ( 22)-( 27), and for the bracket (31) to be antisymmetric, the matrices W (j) s must be of the form

W (0) s =                   
1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0

0 0 0 0 0 1                    , W (1) s
=                    0 1 0 0 0 0 1 0 √ 2 0 0 0 0 √ 2 0 √ 3 0 0 0 0 √ 3 0 0 2α s 0 0 0 0 0 1 0 0 0 2α s 1 0                    , W (2) s
=                    0 0 1 0 0 0 0 √ 2 0 √ 3 0 0 1 0 a 1s b 1s b 2s b 3s 0 √ 3 b 1s b 4s c 1s c 2s 0 0 b 2s c 1s b 7s c 3s 0 0 b 3s c 2s c 3s b 10s                    , W (3) s
=                    0 0 0 1 0 0 0 0 √ 3 0 0 2α s 0 √ 3 b 1s b 4s c 1s c 2s 1 0 b 4s a 2s b 5s b 6s 0 0 c 1s b 5s b 8s c 4s 0 2α s c 2s b 6s c 4s b 11s                    , (35)
W (4) s =                    0 0 0 0 1 0 0 0 0 0 0 1 0 0 b 2s c 1s b 7s c 3s 0 0 c 1s b 5s b 8s c 4s 1 0 b 7s b 8s a 3s b 9s 0 1 c 3s c 4s b 9s b 12s                    , W (5) s
=                    0 0 0 0 0 1 0 0 0 2α s 1 0 0 0 b 3s c 2s c 3s b 10s 0 2α s c 2s b 6s c 4s b 11s 0 1 c 3s c 4s b 9s b 12s 1 0 b 10s b 11s b 12s a 4s                    , where a 1s , • • • , a 4s , b 1s , • • • , b 12s , c 1s , • • • , c
4s are constants to be determined by the constraint of the Jacobi identity.

Indeed we can observe that there exist three types of contributions that can be added to the bracket (31) in such a way that antisymmetry is respected and that, adopting (29) as Hamiltonian, the equations of motion do not get modified. These three types of terms correspond to

d 3 xg is [F is , G is ], i ∈ {2, 3, 4, 5}, (36)
d 3 x[g is [F j s , G j s ] + g j s ([F is , G j s ] + [F j s , G is ])], (37)
i, j ∈ {2, 3, 4, 5}, i = j

d 3 x[g is ([F j s , G ks ] + [F ks , G j s ]) + g j s ([F ks , G is ] + [F is , G ks ]) +g ks ([F is , G j s ] + [F j s , G is ])], i, j, k ∈ {2, 3, 4, 5}, i = j = k. ( 38
)
Because

H is = n 0 T s g is , for i = 2, • • • , 5
, it is not difficult to see that any linear combination of terms of the type ( 36), ( 37) and ( 38), when included in the bracket, is irrelevant for the equations of motion (33) given that [g is , g is ] = 0.

Such terms, however, although transparent for the model dynamics, can be crucial for the Jacobi identity.

The matrices (35) account for the most general linear combinations of terms of the type ( 36)- (38). In particular, As above mentioned, the Jacobi identity for the bracket (31), is satisfied if and only if all the matrices in Eq. (35) pairwise commute. Imposing this condition determines the arbitrary coefficients as follows:

a 1s = 2 √ 2, a 2s = 0, a 3s = 0, a 4s = - 1 + 2α 2 s -8α 4 s √ 6α s , b 1s = 0, b 2s = 0, b 3s = √ 6α s , b 4s = √ 2(1 + 2α 2 s ), b 5s = 0, b 6s = 2 3 α s (1 + 4α 2 s ), b 7s = 0, b 8s = 0, ( 39
) b 9s = - 1 -2α 2 s √ 6α s , b 10s = 2 √ 2α 2 s , b 11s = 0, b 12s = 0, c 1s = √ 2α s , c 2s = 0, c 3s = 0, c 4s = - 2 3 (1 -2α 2 s ).
We remark that the solution is unique.

We have thus shown that the system ( 22)-( 27), in the limit when terms involving the derivatives with respect to z are neglected, possesses a Hamiltonian structure. This structure consists of the Hamiltonian functional ( 29) and of a Poisson bracket of the form (31), with A ij s = 0, for i, j = 0, • • • , 5, and where all the coefficients W ij k s are zero, except for the following ones:

W i0 k s = W 0i k s = δ i k , i = 0, • • • , 5, W 11 0 s = 1, W 21 1 s = W 12 1 s = √ 2, W 11 2 s = √ 2, W 31 2 s = W 13 2 s = √ 3, W 21 3 s = W 12 3 s = √ 3, W 51 3 s = W 15 3 s = 2α s , W 51 4 s = W 15 4 s = 1, W 31 5 s = W 13 5 s = 2α s , W 41 5 s = W 14 5 s = 1, W 22 0 s = 1, W 23 1 s = W 32 1 s = √ 3, W 22 2 s = 2 √ 2, W 25 2 s = W 52 2 s = √ 6α s , W 23 3 s = W 32 3 s = √ 2(1 + 2α 2 s ), W 24 3 s = W 42 3 s = √ 2α s , W 23 4 s = W 32 4 s = √ 2α s , W 22 5 s = √ 6α s , W 25 5 s = W 52 5 s = 2 √ 2α 2 s , W 33 0 s = 1, W 35 1 s = W 53 1 s = 2α s , W 33 2 s = √ 2(1 + 2α 2 s ), W 34 2 s = W 43 2 s = √ 2α s , W 35 3 s = W 53 3 s = 2 3 α s (1 + 4α 2 s ), W 35 4 s = W 53 4 s = - 2 3 (1 -2α 2 s ), W 33 5 s = 2 3 α s (1 + 4α 2 s ), W 34 5 s = W 43 5 s = - 2 3 (1 -2α 2 s ), W 44 0 s = 1, W 45 1 s = W 54 1 s = 1, W 45 3 s = W 54 3 s = - 2 3 (1 -2α 2 s ), W 45 4 s = W 54 4 s = - 1 -2α 2 s √ 6α s , W 44 5 s = - 1 -2α 2 s √ 6α s , W 55 0 s = 1, W 55 2 s = 2 √ 2α 2 s , W 55 5 s = - 1 + 2α 2 s -8α 4 s √ 6α s . ( 40
)
Now we intend to determine the Hamiltonian structure for the full model ( 22)-( 27), including the terms involving z-derivatives. To this scope we follow the procedure of Ref. [29]. According to such procedure, the Hamiltonian functional for the full model corresponds to the one for its perpendicular reduction, that is to the functional (29), in this case. The Poisson bracket for the full model, on the other hand, has the form (30), with the coefficients W ij k s equal to those found for the perpendicular reduction, and with coefficients A ij s such that the the full model equations are retrieved, as equations of motion, and that the properties of a Poisson bracket are respected. Bilinearity and Leibinz identity are clearly satisfied for any choice of A ij s . On the other hand, antisymmetry amounts to the condition A ij s = A ji s , whereas the Jacobi identity implies that the coefficients W ij k s and A ij s must satisfy [29]

W jk r s A ri s = W ki r s A rj s = W ij r s A rk s , i, j, k = 0, • • • , 5, (41)
where the sum over the repeated index r is understood.

With the help of Eq. ( 34), one obtains that the full model equations ( 22)-( 27) are obtained from the Hamiltonian (29) and from a bracket of the form (30), if

A ij s = (W (1) s ) j i , (42)
where the expression for W

(1) s is given in Eq. (35). The matrix W (1) s is symmetric, which guarantees antisymmetry of the Poisson bracket.

Concerning the Jacobi identity, we note that, for any fixed i, one has

W ki r s A rj s = (W (i) s ) k r (W (1) s ) r j = (W (i) s ) j r (W (1) s ) r k = W ij r s A rk s , (43)
which gives the second equality in Eq. ( 41). In the second step of Eq. ( 43) we made use of the fact that the matrix

(W (i) s ) k r (W (1) s )
r j is symmetric, given that all the matrices W (j) s are symmetric and pairwise commute. Similarly, one has

W ij r s A rk s = (W (j) s ) i r (W (1) s ) r k = (W (j) s ) k r (W (1) s ) r i = W jk r s A ri s , (44)
which proves the remaining equality in Eq. ( 41). The Jacobi identity is thus satisfied.

To summarize, we have shown that the model ( 22)-( 27) admits a Hamiltonian structure. The Hamiltonian functional is given by the functional (29). The Poisson bracket is given by Eq. ( 30), with coefficients W ij k s equal to zero except for those indicated in Eq. ( 40), and coefficients A ij s given by Eq. ( 42).

We find it useful, at this point, in order also to help the contact with the existing literature, to reformulate the model equations ( 22)-( 27), in terms of the dimensional physical variables introduced in Eq. ( 14). In terms of such variables the model can be written as

∂n s ∂t = - c B [φ, n s ] + n 0 ∇ u s , (45)
∂ ∂t u s + q s m s c A = - c B φ, u s + q s m s c A + T s m s ∇ T s T s + n s n 0 - q s m s ∂φ ∂z , (46)
∂T s ∂t = - c B φ, T s + 2 n 0 ∇ (q s + n 0 T s u s ), (47)
∂q s ∂t = - c B [φ, q s ] -v ts ∇ √ 6α s q ⊥s - 3 2 n 0 v ts T s , (48)
∂T ⊥s ∂t = - c B [φ, T ⊥s ] + 1 n 0 ∇ q ⊥s , (49)
∂q ⊥s ∂t = - c B [φ, q ⊥s ] -v ts ∇ 2 2 3 α s q s -n 0 v ts T ⊥s , (50)
where 49) and (50) determine the evolution of temperature and heat flux fluctuations due to the motion perpendicular to the magnetic guide field.

∇ f = (1/B)[A, f ] -∂ z f ,
Each dynamical field undergoes the advection by the E × B velocity, represented by the first term on the righthand side of each equations, where the electrostatic potential φ appears. In addition to this effect, each field evolves according to the gradients of neighbouring moments, that form along the magnetic field. For instance, temperature evolutions are affected by the presence of velocity and heat flux gradients. On the other hand, parallel/perpendicular heat fluxes are affected by parallel/perpendicular temperature gradients but also by perpendicular/parallel heat fluxes, through the coupling related to the coefficients α s .

From the inspection of the expression for the coefficients W 45 4 s , W 54 4 s , W 44 5 s , W 55 5 s in Eq. ( 40) it emerges that, in the Hamiltonian model, the constants α s cannot be zero. Therefore, the existence of a Hamiltonian structure for the model, with a Poisson bracket of the form (30), requires the existence of the terms in Eqs. ( 25) and ( 27), which introduce a direct coupling between the parallel and perpendicular heat fluxes. Although, to the best of our knowledge, its Hamiltonian structure is not known, the model in the case α s = 0, is nevertheless energy conserving.

V. NORMAL FIELDS AND CASIMIR INVARIANTS

In this section we show that there exists an alternative set of dynamical variables, in terms of which the model ( 22)-( 27) takes a more compact form, unveiling the existence of some conservation laws. In particular, the Poisson bracket in terms of such new variables takes a remarkably simpler form, which greatly simplifies the identification of Casimir invariants.

In order to identify the new appropriate set of dynamical variables, it is important to note that the generic equation of the model can be written as

∂g ms ∂t = - c B [φ, g ms ] + v ts B (W (1) s ) mn [A, g ns ] -v ts (W (1) s ) mn ∂ ∂z g ns -δ m1 v ts q s T s ∂φ ∂z + √ m!(δ m0 + δ m2 ) v 2 ts c q s T s ∂A ∂z , 0 ≤ m ≤ 5. (51)
The matrices W (1) s are real and symmetric, therefore each of such matrices possesses six real eigenvalues λ 0s , λ 1s • • • , λ 5s . Moreover, for each species s, there exists an orthogonal matrix U s , such that U s T W (1) s U s = Λ s , where Λ s = diag(λ 0s , λ 1s , • • • , λ 5s ). If we perform the change of variables (g 0s , g 1s ,

• • • , g 5s ) → (G 0s , G 1s , • • • , G 5s ),
where

G is = U T s im G ms , (52)
the model equations (51) transform into

∂G is ∂t = - c B φ -λ is v ts c A, G is -v ts λ is ∂G is ∂z -v ts √ m!U T s im δ m1 q s T s ∂φ ∂z -v ts (δ m0 + δ m2 ) q s T s ∂A ∂z , i = 0, • • • , 5. (53)
The form (53) for the model equations presents some advantages with respect to the original formulation (51). Indeed, the evolution equation for each G is depends only on G is (apart from the implicit dependence on the other dynamical variables coming from the potentials φ and A), whereas the evolution equations for g is in general involved explicitly also the dynamical variables g i-1 s and g i+1 s . Also, the form of Eq. ( 53), in the two-dimensional limit assuming z as ignorable coordinate, reveals that the fields G is are Lagrangian invariants. Indeed, in the two-dimensional limit, the generic equations ( 53) can be written as

∂G is ∂t + v is • ∇G is = 0, ( 54
)
where

v is = c B ẑ × ∇ φ -λ is v ts c A (55)
is an incompressible velocity field advecting the dynamical field G is . The two contributions in such velocity field are related to the E × B velocity field and to the free streaming along the "poloidal" magnetic field, respectively, in the original drift-kinetic description. The existence of such advecting field becomes transparent when the model is formulated in terms of the G is variables.

The explicit expressions for the eigenvalues λ 0s , • • • , λ 5s can be found in A. The columns v s of the matrices U s , on the other hand, can be found from W (1) s v s = Λ s v s .

For the sake of our analysis, we point out that the first three rows of the matrices U s have the following structure:

U s =                v (0) s v (1) s • • • • • • v (5) s λ 0s v (0) s λ 1s v (1) s • • • • • • λ 5s v (5) s λ 2 0s -1 √ 2 v (0) s λ 2 1s -1 √ 2 v (1) s • • • • • • λ 2 5s -1 √ 2 v (5) s . . • • • • • • . . . • • • • • • .                , (56)
where v (0) s , • • • , v (5) s are constants determined by the normalization condition. By virtue of the expression (56), the model equations (53) can also be written as

∂G is ∂t = - c B φ -λ is v ts c A, G is -v ts λ is ∂G is ∂z -λ is v ts v (i) s q s T s ∂ ∂z φ -λ is v ts c A , i = 0, • • • , 5. (57)
The Hamiltonian functional (29), in terms of the new variables, reads

H(G 01 , • • • , G 5N ) = 1 2 s d 3 x n 0 T s 5 i=0 G 2 i s +q s n 0 U 0l s G ls φ -q s n 0 v ts c U 1l s G ls A , (58)
and consequently

His = n 0 T s G is + q s n 0 v (i) s φ -λ is v ts c A . (59)
The advection form of Eq. ( 57), and the expression for the functional derivatives (59) suggests that, in terms of the new variables, the Poisson bracket for the model is given by

{F, G} = s 5 i=0 - c q s Bn 0 v (i) s d 3 xG is δF δG is , δG δG is -λ is v ts n 0 T s d 3 x δF δG is ∂ ∂z δG δG is . ( 60
)
It is straightforward to verify, for instance by using the procedure adopted in Sec. IV, that the bracket (60) is a

Poisson bracket and that in particular it satisfies the Jacobi identity. The bracket (60) has a direct sum structure, consisting of a linear combination of Poisson brackets labelled by the indices i and s. Poisson brackets with such a structure are common in reduced fluid models for plasmas and further examples can be found for instance in Refs. [START_REF] Waelbroeck | [END_REF]31,38]. By direct substitution into the expression ∂ t G is = {G is , H}, using the bracket (60) and the Hamiltonian (58), it is easy to see that the model equations (57) are retrieved.

As above noticed, the alternative set of variables G 0s , • • • , G 5s provides a considerable simplification the structure of the Poisson bracket, with respect to the Poisson bracket introduced in Sec. IV in terms of the variables g 0s , • • • , g 5s .

We refer to such variables, reducing a Poisson bracket to its simplest form, as to normal fields. The main advantage of such simplification is that, from Eq. (60), it is staightforward to identify Casimir invariants for the model. We recall that, Casimir invariants of a Poisson bracket { , } (see, e.g., Refs. [25,[START_REF] Marsden | Introduction to Mechanics and Symmetry[END_REF]) are functionals C, such that

{C, F } = 0, (61)
for any functional F . Because Casimirs commute in particular with the Hamiltonian functional, they are also conserved by the dynamics.

From the definition (61), it is easy to see that, for the bracket (60), Casimirs are given by

C is = d 3 x G is , i = 0, • • • , 5. (62)
The conservation of such integrals therefore constraints the dynamics of the model. In particular, in the twodimensional limit, with z as ignorable coordinate, Casimirs correspond to

C is = d 3 x C is (G is ), i = 0, • • • , 5, (63)
where C is are arbitrary functions. Therefore, in the presence of such spatial invariance of the system (which is a sensible approximation for laboratory and solar plasmas where a strong component of the magnetic field makes the dynamics independent, at the leading order on one coordinate) one obtains that the evolution of the model is constrained by an infinite number of conservation laws. These are associated with the role, discussed in Sec. IV, of the fields G is as Lagrangian invariants in the two-dimensional limit.

We remark that if, for a given species s, one has α s = 0, then the W (1) s becomes a block diagonal matrix with eigenvalues given by

λ 0s = 1, λ 1s = -1, λ 2s = -3 + √ 6, λ 3s = 3 + √ 6, λ 4s = -3 - √ 6, λ 5s = 3 - √ 6.
The corresponding matrix U s , on the other hand, reads

U s =                    0 0 -3 √ 3+ √ 6+ √ 6(3+ √ 6) 3 √ 2 3 √ 3+ √ 6- √ 6(3+ √ 6) 3 √ 2 3 √ 3- √ 6+ √ 6(3- √ 6) 3 √ 2 -3 √ 3- √ 6- √ 6(3- √ 6) 3 √ 2 0 0 1 √ 2 1 √ 2 -1 √ 2 -1 √ 2 0 0 -3+ √ 6 6 3+ √ 6 6 -3- √ 6 6 3- √ 6 6 0 0 1 √ 2 1 √ 2 1 √ 2 1 √ 2 1 √ 2 -1 √ 2 0 0 0 0 1 √ 2 1 √ 2 0 0 0 0                    . (64)
From the expression (64) one sees that, in this case, v (0) s = v (1) s = 0. Consequently, the Poisson bracket (60), in this case does not exist, because v (0) s and v (1) s would vanish at the denominator. Thus, in spite of the fact that, also for α s = 0, the system can be cast in the form (57), a Hamiltonian formulation with Hamiltonian functional given by (59) and Poisson bracket of the form (60) is not possible, thus confirming the conclusions drawn in Sec. IV with regard to the case α s = 0. By imposing the additional constraints of antisymmetry and the Jacobi identity the expression for the Lie-Poisson bracket of the model in the 2D limit is thus obtained. The extension to 3D is then carried out following a procedure adopted in Ref. [29]. From the expression of the Poisson bracket it emerges that the proportionality coefficients α s in the closure relation cannot be zero, thus implying the presence of the terms coupling parallel and perpendicular heat fluxes in the model equations, in order to assure a Hamiltonian structure with a Lie-Poisson bracket.

An alternative set of dynamical variables (normal fields) is then introduced, in terms of which the Poisson bracket of the model simplifies considerably. This made the identification of Casimir invariants particularly simple. In the 2D limit, 6N infinite families of Casimir invariants are found. These express the property of the normal fields of being Lagrangian invariants of the model. This unveils a remarkable feature of the model in the 2D limit: the model can indeed be reformulated as a system of 6N advection equations for 6N Lagrangian invariants transported by incompressible generalized velocity fields. Such velocity fields, corresponding to those indicated in Eqs. (55), are associated with stream functions corresponding to linear combinations of electrostatic and magnetic potential. These account for the E × B velocity and free stream along the fluctuating part of the magnetic field. A similar feature occurs in other Hamiltonian reduced plasma models, e.g. in Refs. [START_REF] Grasso | [END_REF][START_REF] Waelbroeck | [END_REF]31,38] but none of such models accounted also for perpendicular temperature and heat fluxes.

Finally, it is worth emphasizing a remark about the adopted closure in the case α s = 0. This corresponds to eliminating the terms directly coupling parallel and perpendicular heat fluxes for a given species and is the closure that appears to be adopted in the models of Refs. [2,7], when taken in the limit of negligible finite Larmor radus effects, background inhomogeneities and dissipation. The Hamiltonian formulation found for the present model fails in the case α s = 0. This special case turns out to be energy conserving but no Poisson bracket of the form ( 30) is able to yield the model equations when combined with the assumed Hamiltonian functional. This suggests that, for the particular closure α s = 0, the Hamiltonian structure requires a different form of Poisson bracket or that ingredients neglected in the present analysis (e.g. finite Larmor radius effects) are needed to restore a Hamiltonian structure with the adopted form of Poisson bracket.

The author acknowledges useful discussions with the Nonlinear Dynamics Team of the Centre de Physique Théorique. Financial support was received from the CNRS through the PEPS project GEOPLASMA2.

Appendix A: Determining the eigenvalues of the matrices W (1) s For a fixed species, the eigenvalues of the matrix W (1) s correspond to the roots of the characteristic polynomial Following a standard procedure to determine the roots of an algebraic equation of third degree, we operate the transformation u = y + (7 + 4α

a

  1s , • • • , a 4s are coefficients of the terms of the type (36), whereas b 1s , • • • , b 12s and c 1s , • • • , c 4s are coefficients of terms of the type (37) and (38), respectively.

  for a generic function f . Such operator represents a normalized gradient along the magnetic field B = ∇A × ẑ + B ẑ. Eqs. (45) and (46) express particles number and canonical momentum conservation. Eqs. (47)-(48) govern the evolution of the temperature and heat flux fluctuations in the direction parallel to the constant component of the magnetic field, and similarly (

  We have presented a Hamiltonian anisotropic electromagnetic fluid model evolving six field variables for each particle species, accounting for the evolution of density, canonical fluid momentum, as well as parallel and perpendicular temperature and heat flux fluctuations. The model is derived by taking moments with respect to Hermite-Laguerre polynomials of a generalized distribution function and imposing as closure a proportionality relation between the moments g 40s and g 21s and the perpendicular and parallel heat flux fluctuations, respectively. Such closure amounts to introduce terms coupling directly the evolutions of the parallel and perpendicular heat fluxes, for each species.As Hamiltonian for the fluid model a conserved functional obtained self-consistently form the Hamiltonian of the parent drift-kinetic model is assumed. The Poisson bracket is then constructed by assuming it to possess the form of a Lie-Poisson bracket, which is very common for reduced fluid plasma models. A general form of bilinear functional satisfying the Leibniz condition and yielding, from the assumed Hamiltonian, the model equations in 2D, is considered.

λ 6 -

 6 (7 + 4α 2 s )λ 4 + (9 + 12α 2 s )λ 2 -3,(A1)which we can rewrite asu 3 -(7 + 4α 2 s )u 2 + (9 + 12α 2 s )u -3, (A2)where u = λ 2 .

  2 s )/3, which casts the polynomial (A2) in the form

	The roots of (A3) are then given by
	y 1 =	2 3	√	2 8α 4 s + 10α 2 s + 11 cos		1 3	Arg	4 27	(16α 6 s + 30α 4 s + 12α 2 s + 25)
					+	2 3	i √	2 24α 8 s + 40α 6 s + 36α 4 s + 30α 2 s + 1 ,
	y 2 =	2 3	√	2 8α 4 s + 10α 2 s + 11 cos		1 3	Arg	4 27	(16α 6 s + 30α 4 s + 12α 2 s + 25)
				y 3 -	16α 4 s + 20α 2 s + 22 3	y -	32α 2 s + 80α 4 s 9	-	128α 6 s + 200 27	.	(A3)