Ha Dong Hwang
email: hadong.hwang@insa-lyon.fr

Kerem Ege

Laurent Maxit

Nicolas Totaro

Jean-Louis Guyader

A METHODOLOGY FOR INCLUDING THE EFFECT OF A DAMPING TREATMENT IN THE MID-FREQUENCY DOMAIN USING SmEdA METHOD

An additive damping treatment is an effective tool to control the dynamic response of builtup structures, and it is widely utilized through industrial applications. By applying a viscoelastic layer on a given structure, the vibratory energy is dissipated through shear and in-plane motions at the layer interface. Modeling the effect of such a treatment in a complex mechanical system for the mid frequency domain is of interest. Statistical modal Energy distribution Analysis (SmEdA) has been developed as an alternative approach to Statistical Energy Analysis (SEA) for describing subsystems with low modal overlap. This technique is developed from the knowledge of the uncoupled subsystem modes. In this paper, one proposes to extend SmEdA by including the effect of a damping treatment. A damped subsystem consisting of a composite layer is modeled with the equivalent modulus of a single layer, which gives the same transverse displacement as a multi-layered system. The modal loss factor of a partially damped structure is estimated by the Modal Strain Energy method (MSE), and the results are well agreed with the Complex Eigenvalue Method (CEM). Finally, energy transmission between the damped structure and a coupled cavity can be deduced from SmEdA modeling, knowing the modeshapes and modal loss factors of the equivalent single layer and of the cavity. This method is applied for modeling a rectangular plate partially damped with an unconstrained viscoleastic layer coupled to a small acoustic cavity.

Introduction

Deterministic FEM analysis and SEA method are typical means to predict dynamic response of a complicated mechanical system, yet neither method is appropriate for a mid-frequency domain. Over the years, a number of alternative/hybrid techniques have been developed to cover this gap between low and high frequency domains. SmEdA (Statistical Modal Energy Distribution Analysis) [START_REF] Maxit | Estimation of SEA coupling loss factors using a dual formulation amd FEM modal information, part I: Theory[END_REF][START_REF] Guyader | Estimation of SEA coupling loss factors using a dual formulation amd FEM modal information, part II: Numerical applications[END_REF][START_REF] Maxit | Extension of SEA model to subsystems with non-uniform energy distribution[END_REF][START_REF] Totaro | SEA coupling loss factors of complex vibroacoustic systems[END_REF][START_REF] Totaro | Extension of statistical modal energy distribution analysis for estimating energy densityin coupled subsystems[END_REF] is among them, where the basic SEA theory is modified on the energy flow exchanged by two oscillators coupled through a gyroscopic element and use of the dual modal formulation for describing coupled subsystems. With modal information of uncoupled subsystems deduced from FEM, individual modal energies between coupled resonant modes can be apprehended rather than ensemble energy of modes that dictates classical SEA theory. Since eigenvectors and natural frequencies of uncoupled subsystems are deduced from FEM, the SmEdA technique can evaluate a system that has a complicated geometry and consists of composite materials.

This paper aims at analyzing energy transmission of a structure-fluid coupled system for a mid-frequency domain. SmEdA can reveal energy distribution of coupled modes in a given frequency band for a rectangular steel plate coupled to an acoustic cavity. This is extended to a case where a plate subsystem is treated with an unconstrained viscoelastic damping layer. The influence of additive damping on the energy transmission is of interest. A rectangular steel plate partially damped with a viscoelastic patch is characterized by equivalent modulus of a single layer plate [START_REF] Guyader | Viscoelastic properties of singly layer plate material equivalent to multi-layer composites plate[END_REF][START_REF] Guyader | Acoustic transmission through orthotropic multilayered plates, part I: Plate vibration modes[END_REF] that gives the same transverse displacement as a multi-layered plate. In this approach, multi-layers are replaced by a single homogeneous layer with equivalent mechanical properties calculated by the mixture of classical plate theories. The frequency-dependent properties of a viscoelastic material are considered in the equivalent plate modeling which can compute Young's modulus and damping loss factors at resonant modes for an infinite plate. In order to calculate damping levels of a finite plate, MSE (Modal Strain Energy) method [START_REF] Koruk | Assesment of the complex eigenvalue and the modal strain energy methods for damping predictions[END_REF][START_REF] Assaf | Vibration and damping analysis of plates with partially covered dampinf layers[END_REF][START_REF] Fronk | Finite element modeling of damping in constrained layer composite structures induced by inplane loads using ADINA[END_REF] is employed. The fundamental assumption of MSE is that the modeshapes of a damped structure is the same as those of undamped. Then, a ratio of a damped plate loss factor to the viscoelastic material loss factor can be estimated from a ratio of the elastic strain energies stored in the entire structure to the viscoelastic layer. Finally, with the modal loss factors of a damped plate calculated by series of mentioned techniques, the energy transmission between subsystems can be deduced.

This paper presents numerical analysis of SmEdA on undamped/damped rectangular steel plate coupled to an acoustic cavity. When a plate subsystem is partially damped with a viscoelastic patch, the composite panel is modeled as an equivalent single layer, and modal damping levels of total subsystem are evaluated and compared to those obtained from CEM (Complex Eigenvalue Method) [START_REF] Koruk | Assesment of the complex eigenvalue and the modal strain energy methods for damping predictions[END_REF][START_REF] Werner | Intergrated vibro-acoustic calculation procedure of composite structures using coupled finite and boundary element method[END_REF] . The modal coupling energy loss factors between subsystems and their energy relation is evaluated.

SmEdA Method For The Mid-Frequency Domain

SmEdA is dictated by the basic linear equation of power exchanged between coupled subsystems as the SEA principle. However, the method differentiates by dealing with energy levels of individual modes of subsystems rather than ensemble energy within a considered frequency band Δf. This compensates rather disputable SEA assumption of the energy equipartition among subsystem modes. Two subsystems with natural mode p and q, and powers injected into mode p of subsystem 1 and mode q of subsystem 2 can be written as .

(

) 1
where Π inj p and Π inj q are the time-averaged power injected into mode p and q of subsystem 1 and 2, E p and E q are the energy of mode p and q of subsystem 1 and 2, η p and η q are the damping loss factor of mode p and q of subsystem 1and 2, η pq is the coupling loss factor between mode p and q of subsystem 1 and 2. Consider a rectangular plate coupled to an acoustic cavity illustrated in Fig. 1, the dual modal formulation uses the modal bases of uncoupled subsystems described by displacement and pressure with blocked and free boundary conditions accordingly.

Modeshapes of each subsystem are extracted from FEM, and interactions between coupled modes are evaluated by considering the discretized modes at a coupling surface. The modal coupling loss factor is given by 20th International Congress on Sound and Vibration (ICSV20), Bangkok, Thailand, 7-11 July 2013 ICSV20, Bangkok, Thailand, 7-11 July 2013 3

(2)

where β ijkmn , β W and β ω are the modal coupling loss factor between mode (m,n) of the plate and mode (i, j,k) of the cavity, the modal spatial coupling factor and the modal spectral coupling factor respectively. W ijkmn M mn and M ijk are the inter-modal work between mode (m,n) of the plate and mode (i, j,k) of the cavity and the modal mass of mode (m,n) of the plate and mode (i, j,k) of the cavity respectively. η ijk and η mn are the modal damping loss factors of each subsystem. If a plate is treated with a viscoelastic damping material, the complex modulus of elasticity is considered for the modal damping level. Deducing this term, η mn is rendered by MSE method and is presented in Section 3.3. The inter-modal work (W ijkmn ) is given by .

where p ijk (Q) and W mn (Q) are modeshapes of the cavity and the plate. This term represents the modal interaction between displacement modeshape (m, n) of the plate and pressure modeshape (i, j,k) of the cavity.

Analysis of A Structure-Fluid Coupled System

The damping loss factor of the plate subsystem

The properties of a steel plate are experimentally determined. Young's modulus and the loss factor are measured from modal analysis of a plate with a free-free boundary condition. The plate is driven by an impact hammer at one corner of the plate, and the response is measured with an accelerometer at another corner. The loss factors are obtained from well-separated modes by the high-resolution modal analysis method based on the ESPRIT algorithm [START_REF] Ege | High-resolution modal analysis[END_REF] on the recorded time signal. Table 1 shows the loss factors averaged over frequency bands of 400 Hz. The total (global) average loss factor is considered for the equivalent single plate modelling as presented in Table 2.

Equivalent single layer modelling of a partially damped plate

Developed by Guyader et al. [START_REF] Guyader | Viscoelastic properties of singly layer plate material equivalent to multi-layer composites plate[END_REF][START_REF] Guyader | Acoustic transmission through orthotropic multilayered plates, part I: Plate vibration modes[END_REF] and implemented in the software MOVISAND, the equivalent single layer modelling simplifies a multi-layered panel with a single homogenous layer of equivalent viscoelastic moduli. The method is based on transverse waves motion computation (Love-Kirchhoff thin plate theory), the equivalent viscoelastic moduli of a single layer are determined in order to give the same transverse displacements of multi-layered panel. The method considers bending, membrane and shear motions in each layer and uses the continuity condition on displacements and shear stresses at layer interfaces to obtain an equation of motion that is independent of a number of layers. In order to characterize an equivalent single layer material, equivalent moduli e.g. Young's modulus, mass density and Poisson's ratio must be determined. There parameters are given by where D and h are bending stiffness of the equivalent layer and thickness of each damping layer respectively. Derivation of the analytical expression of the equivalent bending stiffness D can be found in the literature [START_REF] Guyader | Acoustic transmission through orthotropic multilayered plates, part I: Plate vibration modes[END_REF] . A rectangular steel plate subsystem damped with a rectangular PVC (polyvinyl chloride) patch is modelled as an exact multi-layered plate and an equivalent plate as seen in Fig. 2. The PVC patch bears 33% of the base plate, and its location is arbitrarily chosen over the plate surface. Young's modulus of the steel plate presented in Table 2 is obtained from the averaged relative natural frequency differences between numerical FE calculation and experimental results. MOVISAND software calculates equivalent parameters for an infinite plate model according to Eq. 4; Young's modulus, density, Poisson's ratio and damping loss factor. Figure 3 shows the variations of Young's modulus and damping factors with respect to vibration frequency for the equivalent infinite plate model. In fact, neither quantity varies significantly since complex moduli (frequency-dependent) of PVC are neglected in the calculation. Once equivalent parameters are determined for an infinite plate model, they can be taken into a computation of the dynamic behaviour of a finite plate partially damped with a PVC patch. To evaluate the equivalent modelling method, the frequency response function of two plate models shown in Fig. 2 is computed. Seen from the equivalent plate model in Fig. 2, the composite core is defined by the equivalent parameters presented Table 2 and Fig. 3, and the rest is the steel plate whose parameters are also given in Table 2. Since complex moduli of the equivalent plate shown in Fig. 3 do not vary significantly in the frequency range of interest (up to 2 kHz), constant averaged values E eqv = 2.95 GPa and η eqv = 0.0096 are considered for the computation. Both plate models with clamped boundaries are driven with a harmonic transverse point force at the same position, and the point response is calculated with NASTRAN FE software. The results are in good agreement between two models, demonstrating a validity of the equivalent single layer modelling.

Numerical estimation of modal loss factor of the finite equivalent plate

Amongst several techniques to predict dynamic behaviours of a damped structure, MSE method has been proved effective and widely used in practice. It relies on the linear dynamics of undamped systems in which a viscously damped linear system is assumed to posses the same real normal modes as undamped systems. In other words, a ratio of the modal composite loss factor to a viscoelastic loss factor can be approximated as a ratio of the elastic strain energy stored in the viscoelastic layer to the total elastic energy stored in the total structure for a given undamped mode ψ. This is given by where K and V are a real part of the stiffness matrix and strain energy. This method is applied to calculate the modal damping levels of a finite damped plate. Note that MSE calculation is performed on the equivalent plate model only. The equivalent parameters calculated for an infinite plate model are taken into FE normal mode calculation on the equivalent plate model with clamped boundaries in order to give the modal damping levels. As the same manner presented in the frequency response function computation given in Section 3.1, constant averaged values of E eqv = 2.95 GPa and η eqv = 0.0096 are considered for the damped part of the finite plate. The results are evaluated with those obtained from the CEM (Complex Eigen Method) [START_REF] Koruk | Assesment of the complex eigenvalue and the modal strain energy methods for damping predictions[END_REF][START_REF] Werner | Intergrated vibro-acoustic calculation procedure of composite structures using coupled finite and boundary element method[END_REF] and are in good agreement.

The damping loss factor of the cavity subsystem

The loss factor of a cavity is experimentally determined. A dimension of the cavity is 0.5 * 0.6 * 0.7 (m) as illustrated in Fig. 1. The cavity is driven with a white noise at a random position, and the pressure decay inside the cavity is recorded after the noise source is turned off. Two methods are used to derive the loss factors. The method presented in the literature [START_REF] Ege | High-resolution modal analysis[END_REF] based on modal analysis and the standard ISO method [START_REF] Iso | Measurement of room acoustic parameters -Part 2: Reverberation time in ordinary rooms[END_REF] are used. For the latter, the decay curve of measured pressure is fitted to a straight line by the linear regression approach in order to calculate the reverberation time. Then the modal loss factor of an acoustic space can be obtained as where A, S, ᾱ, V and T R are total absorption area, surface area of the cavity, angle averaged absorption coefficient, cavity volume and the reverberation time (T 30 ) respectively. Table 3. Frequency band averaged damping loss factor of an acoustic cavity obtained from laboratory measurement with the high-resolution method [START_REF] Iso | Measurement of room acoustic parameters -Part 2: Reverberation time in ordinary rooms[END_REF] . The bandwidth (Δf) is 400 Hz. Center Frequency (Hz) 200 600 1k 1.4k 1.8 Total average η ijk 0.0047 0.0037 0.0025 0.0026 0.0011 0.0029

Energy relation between subsystems

SmEdA computations are rendered on two different coupled systems: undamped steel plate coupled to a cavity (model 1) and partially damped steel plate with a rectangular PVC patch coupled to a cavity (model 2). Dimensions of the plate, cavity and PVC patch are 0.5*0.6*0.001 (m), 0.5*0.6*0.7 (m) and 0.25*0.4*0.003 (m) respectively. The total (global) average loss factor of the cavity shown in Table 3 is considered in the SmEdA computation for both models. For model 2, an average value (η mn = 0.0054) of the modal damping levels of the damped plate subsystem is considered. The plate subsystem is driven with a normal harmonic point force at the same position randomly chosen for both models. The modeshape at the excitation point is considered for the injected power, which is used to compute subsystem energy levels from Eq. 1. SmEdA reveal the energy transmission between coupled resonant modes in a given frequency band. As seen in Fig. 6 (c1), the coupled modes (11 th plate mode at 1097.13 Hz, 10 th cavity mode at 1096.33 Hz) of model 1 dominate the energy exchange rate within a frequency band centered at 1 kHz, and this strong coupling occurs when their resonant frequencies coincide. This introduces a considerable modal energy disparity as overlapping between subsystem modes is low. In such case, SEA is not able to describe uneven modal energy distribution and overestimates the energy transfer.

The influence of damping in model 2 can be clearly detected from the modal coupling loss factor (β ijkmn ) of model 2 as presented in Fig. 6 (c2). For the same frequency band considered (1 kHz), the resonant coincidence that conducts the highest level of energy transfer is shifted to a different mode couple (19 h plate mode at 993.2 Hz, 15 th cavity mode at 993.11 Hz) compared to model 1. While an additive damping layer surely reduces the vibratory energy of a plate subsystem, it does not necessarily lead to a reduction of energy transmission to a cavity subsystem. In fact, for the frequency band of 1 kHz an energy ratio of the cavity subsystem to the plate subsystem (E2/E1) of model 2 exceeds that of model 1, which indicates damped plate permits more energy to be transmitted to a cavity than the case of undamped. As seen in Fig. 6 (a1) and (a2), the interactions between subsystem modeshapes at the coupling surface do not significantly differentiate between two models as total levels (sum) of their modal spatial coupling factors (β W ) within the frequency band are nearly the same. On the other hand, the damping augments the frequency coupling as a sum of the modal spectral coupling factor (β ω ) of model 2 is higher than that of model 1 as presented in Fig. 6 (b1) and (b2). This results in approximately 40 % higher modal damping loss factor (β ijkmn ) for model 2, which indicates that the influence of damping over frequency coincidence is more responsible for the strong coupling between subsystems than spatial coincidence. However, such influence of the damping is shown to be disproportionate over different frequency bands, which depends on the modal density.

Conclusion

The dynamic behaviour of the structure-fluid coupled system is investigated with SmEdA methodology. Energy transmission between a plate subsystem and an acoustic cavity subsystem is deduced based on energies of individual modes rather than ensemble energy in a given frequency band. When a plate subsystem is partially damped with an additive viscoelastic layer, the composite plate can be modelled as an equivalent single layer which gives the same transverse displacement as a multi-layered plate. Such a method has been demonstrated computationally efficient and proved to yield comparable results. The modal loss factor of a damped plate is obtained from the MSE method. A ratio of the composite system loss factor to the viscoelastic loss factor can be approximated from a ratio of the elastic strain energy store in viscoelastic core and total structure. This provides a rapid estimate of modal damping levels with a reasonable accuracy. With the modal loss factors deduced by MSE, the influence of a damping layer over the vibratory energy flow can be evaluated by the SmEdA technique. The advantage of such technique is that it can reveal the energy transmission between individual coupled modes and evaluate the influence of damping on the modal energies for a frequency band under study. It is shown that it is possible for a damping mechanism to allow high energy transmission to a cavity subsystem largely due to the frequency coincidence of subsystems rather than the spatial coincidence in a given frequency band. However, this is a particular case where a damping layer geometry, thickness and location can yield different outcome.

Acknowledgement

This work was co-funded by the French government (FUI 12 -Fonds Unique Interministériel) and European Union (FEDER -Fonds européen de développement régional ). It was carried out in the framework of the LabEx CeLyA ("Centre Lyonnais d'Acoustique", ANR-10-LABX-60) and the research project CLIC ("City Lightweight Innovative Cab") labelled by LUTB cluster (Lyon Urban Truck and Bus), in partnership with Renault Trucks, Arcelor-Mittal, ACOEM, CITI Technologies, FEMTO-ST (Univ. de Franche-Comté) and LVA (INSA de Lyon).

Figure 1 .

 1 Figure 1. Top: a rectangular plate coupled to a cavity via a coupling surface (red area). Bottom: subsystem descriptions with their boundary conditions of the SmEdA model.

Figure 2 .

 2 Figure 2. Equivalent single layer modelling. Boundaries are clamped.

Figure 3 .

 3 Figure 3. Equivalent moduli of the infinite steel-PVC composite panel: (a) Young's modulus, (b) Loss factor

Figure 4 .

 4 Figure 4. Frequency response function (a single response point) of the equivalent and multi-layered plate models, clamped at their boundaries.

Figure 5 .

 5 Figure 5. A steel plate coupled to a cavity (model 1) and a steel plate damped with a PVC patch coupled to a cavity (model 2).

Figure 6 .

 6 Figure 6. Modal coupling factors of model 1 (left) and model 2 (right) in a frequency band centered at 1 kHz. (a1) (a2) The modal spatial coupling loss factor (β W ). (b1) (b2) The modal spectral coupling factor (β ω ). (c1) (c2) The modal coupling loss factor (β ijkmn ).

Table 1 .

 1 Frequency band averaged damping loss factor of a steel plate with a free-free boundary condition obtained from laboratory measurement. The bandwidth (Δf) is 400 Hz.

	Center
	Frequency	200	600	1k	1.4k	1.8	Total
	(Hz)						average
	η	0.0027	0.0017	0.0014	0.0019	0.0009	0.0017

Table 2 .

 2 MOVISAND calculation results: equivalent moduli of the infinite steel-PVC composite panel. Frequency-dependent E(f) and η(f) of the equivalent plate model are presented in Fig.3.

		Steel	PVC	Equivalent
	h (m)	0.001	0.003	0.004
	E (Pa)	2.03e11	2.4e7	E(f)
	(kg/m 3 )	7523.1	1200	2780.77
		0.33	0.45	0.42
	η	0.0017	0.35	η(f)


S Guisset

D Aregba

S Brull

B Dubroca

The M 1 angular moments model in a moving reference frame for rarefied gas dynamics applications

Keywords: Angular moments models, entropy minimisation closure, Galilean invariance, HLL schemes AMS: 76P, 76N, 65D, 65C

published or not. The documents may come

The M1 angular moments model in a moving reference frame for rarefied gas dynamics applications S Guisset, D Aregba, B Brull, B

Introduction

Kinetic descriptions are known to be very accurate to describe the transport of particles in rarefied gas dynamics [START_REF] Grad | On the kinetic theory of rarefied gases[END_REF], neutron transport [START_REF] Sanchez | Review of neutron transport approximations[END_REF], plasma physics [START_REF] Landau | On the vibration of the electronic plasma[END_REF][START_REF] Lifchitz | Kinetic theory[END_REF] or radiative [START_REF] Berthon | Mathematical models and numerical methods for radiative transfer[END_REF]. However, they are also known to be computationally expensive to describe most realistic physical applications. An alternative way consists in considering fluid descriptions based on averaged physical quantities. However, such macroscopic descriptions are often not sufficiently accurate. The studied particles may have an energy distribution far from the thermodynamic equilibrium so that the fluid description is not adapted. Moreover kinetic effects can be important over time scales shorter than the collisional time so that fluid simulations are insufficient and kinetic codes have to be considered to capture the physical processes. Kinetic approaches are usually limited to times and lengths much shorter than those studied with fluid simulations. It is therefore an important challenge to describe kinetic effects using reduced kinetic codes operating on fluid time scales [START_REF] Filbet | A hierarchy of hybrid numerical methods for multi-scale kinetic equation[END_REF][START_REF] Helluy | Reduced Vlasov-Maxwell modeling[END_REF]. The angular moments models represent an alternative method situated in between the kinetic and the fluid models. They require computational times shorter than kinetic models and provide results with a higher accuracy than fluid models. They originate from an angular moments average [START_REF] Minerbo | Maximum entropy Eddigton Factors[END_REF][START_REF] Pomraning | Maximum entropy Eddington factors and flux limited diffusion theory[END_REF] of the kinetic equations. The idea is to keep the velocity modulus (denoted ζ in this work) as a variable. That allows to consider the particle distributions in energy far from equilibrium, while using a simplified description of particle angular distribution. Such models are obtained by integration of the kinetic equation in angle (integration on the unit sphere). Thus a hierarchy of moments equations can be obtained. There exist several moment models whose differences come from the choice of the closure relation. In this document we consider the angular moments models [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF] based on an entropy minimisation principle. The entropy minimisation problems have been widely studied in [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Minerbo | Maximum entropy Eddigton Factors[END_REF][START_REF] Muller | Rational Extended Thermodynamics[END_REF][START_REF] Struchtrup | Macroscopic Transport Equations for Rarefied Gas Flows[END_REF][START_REF] Alldredge | High-order entropy-based closures for linear transport in slab geometry II: A computational study of the optimization problem[END_REF]. The underlying distribution function is given by an exponential of a polynomial function depending on the particle energy and it is therefore non negative. Moreover, these closures verify fundamental mathematical properties [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Schneider | Entropic approximation in kinetic theory[END_REF][START_REF] Junk | Maximum entropy moment systems and Galilean invariance[END_REF] such as hyperbolicity and entropy dissipation. However, their solutions could be rather different from the solution of the full kinetic equation. Moreover, from the numerical point of view, even if the closure is well defined, computational challenges remain. In particular, the resolution of the entropy minimisation problem can be very computationally costly and we refer to [START_REF] Alldredge | High-order entropy-based closures for linear transport in slab geometry II: A computational study of the optimization problem[END_REF] for a specific treatment. The angular M 1 model is largely used in the context of radiative transfer [START_REF] Turpault | Multigroup half space moment appproximations to the radiative heat transfer equations[END_REF][START_REF] Berthon | An HLLC Scheme to Solve The M1 Model of Radiative Transfer in Two Space Dimensions[END_REF][START_REF] Dubroca | Entropic moment closure hierarchy for the radiative transfert equation[END_REF][START_REF] Turpault | A consistent multigroup model for radiative transfer and its underlying mean opacity[END_REF][START_REF] Charrier | Multigroup model for radiating flows during atmospheric hypersonic re-entry[END_REF][START_REF] Ripoll | An averaged formulation of the M1 radiation model with presumed probability density function for turbulent flows[END_REF][START_REF] Ripoll | A factored operator method for solving coupled radiation-hydrodynamics models[END_REF] and starts to be considered for others application fields such as plasma physics [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF][START_REF] Guisset | Classical transport theory for the collisional electronic M1 model[END_REF][START_REF] Guisset | Limits of the M1 and M2 angular moments models for kinetic plasma physics studies[END_REF][START_REF] Guisset | Asymptotic-preserving scheme for the Fokker-Planck-Landau-Maxwell system in the quasi-neutral regime[END_REF] or radiotherapy [START_REF] Pichard | Relaxation schemes for the M1 model with space-dependent flux: application to radiotherapy dose calculation[END_REF]. Here, the dynamic of neutral particles is investigated. To achieve such an issue the velocity framework is centred on the particles mean velocity. The choice of the velocity framework origin on the particles mean velocity enables a reduction of the velocity grids and plays an important role on the Galilean invariance property of the M 1 angular model. This point is presented in details in this study. The approach presented here is significantly different with the previous studies, since the M 1 angular model [START_REF] Dubroca | Étude théorique et numérique d'une hiéarchie de modèles aux moments pour le transfert radiatif[END_REF][START_REF] Guisset | Limits of the M1 and M2 angular moments models for kinetic plasma physics studies[END_REF][START_REF] Guisset | Classical transport theory for the collisional electronic M1 model[END_REF] has always been considered in the laboratory framework. In addition to our knowledge it is the first time that the M 1 angular moments model is considered for rarefied gas dynamic applications.

In order to derive the M 1 angular moments model in the mean velocity frame, a velocity change is considered to derived the kinetic equation in a moving frame. Note, that rescaled velocity approaches are largely used in different context see [START_REF] Filbet | A Rescaling Velocity Method for Kinetic Equations: the Homogeneous Case[END_REF][START_REF] Bobylev | On some properties of kinetic and hydrodynamic equations for inelastic interactions[END_REF][START_REF] Mischler | Cooling process for inelastic Boltzmann equations for hard spheres, Part II: Self-similar solutions and tail behavior[END_REF] for example. However, the numerical treatment of the additional terms which appear when considering such a procedure on the kinetic equation can be challenging. In [START_REF] Filbet | A Rescaling Velocity Method for Dissipative Kinetic Equations -Applications to Granular Media[END_REF], in the context of granular flows, a numerical algorithm based on a relative energy scaling is proposed. Then, a clever de-coupling with the hydrodynamics equation is used to avoid the problems related to the change of scales in velocity variables. Significant results on the derivation of Galilean invariant minimum entropy systems have been obtained in [START_REF] Junk | Maximum entropy moment systems and Galilean invariance[END_REF]. It has been shown that polynomials weight functions growing super-quadratically at infinity lead to unusable hyperbolic moment systems. Indeed, in this case equilibrium states are boundary points of the admissible set with possible singular fluxes [START_REF] Junk | Maximum entropy moment systems and Galilean invariance[END_REF]. In addition it has also been shown that non-polynomial weight function can not be used without losing the Galilean invariance property of the moments systems. Therefore, one understands here the difficulty in deriving adapted Galilean invariant reduced models. We also mention [START_REF] Groth | Towards physically-realizable and hyperbolic moment closures for kinetic theory[END_REF][START_REF] Mcdonald | Towards realizable hyperbolic moment closures for viscous heat-conducting gas flows based on a maximum-entropy distribution[END_REF] where this issue is addressed.

The present study is original for two main reasons. First of all, contrarily to the previous works, angular moments are considered. The integration of the kinetic equation is only performed in angle (on the unit sphere) while the velocity modulus is kept as a variable. This partial integration leads to angular moments models which do not suffer of the restrictions presented in [START_REF] Junk | Maximum entropy moment systems and Galilean invariance[END_REF]. However, it can be shown that these resulting angular moments models are not Galilean invariant since they are not invariant by translational transformations. Therefore, the second original idea presented here to surpass this major drawback, is to work in the framework of the particles mean velocity. In this work, it is shown that this choice of velocity origin framework enables to recover the invariance property when considering translational transformations. Therefore, in this paper, a Galilean invariant reduced model is presented in the context of rarefied gas dynamics.

The plan of this study is the following. First of all, the derivation of the angular M 1 moments model in the mean velocity frame is introduced. The choice of the mean velocity framework in order to enforce the Galilean invariance property of the model is highlighted. In addition, it is shown that the model rewritten in terms of the entropic variables is Friedrichssymmetric. Also, the derivation of the associated conservation laws and the zero mean velocity condition are detailed. Secondly, a suitable numerical scheme, preserving the realisability requirement of the numerical solution for the angular M 1 moments model in the mean velocity frame is proposed. Thirdly, some numerical results obtained considering several test cases in different collisional regimes are displayed. Finally, some conclusion and perspectives are given.

Derivation of the model

Velocity change of variables procedures are used in various contexts (see [START_REF] Filbet | A Rescaling Velocity Method for Kinetic Equations: the Homogeneous Case[END_REF][START_REF] Filbet | A Rescaling Velocity Method for Dissipative Kinetic Equations -Applications to Granular Media[END_REF][START_REF] Bobylev | On some properties of kinetic and hydrodynamic equations for inelastic interactions[END_REF][START_REF] Mischler | Cooling process for inelastic Boltzmann equations for hard spheres, Part II: Self-similar solutions and tail behavior[END_REF] for example) and can enable the simplification of a collisional operator form or the reduction of the velocity grid size used for numerical applications. In the context of angular moments models [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF], it will be seen in the next section that moving frame formulations play an important role in enforcing the Galilean invariance property. In order to explain in details this point, in this section we introduce the kinetic formulation in a moving frame from which the M 1 angular moments model studied is derived.

The kinetic equation in a moving velocity frame

We start considering the following kinetic equation written in the laboratory framework ∂f (α)

∂t + div x (vf (α)) = C(f (α)), (1)
where f represents the particle distribution function and

α = (t, x, v) ∈ R + t × R 3 x × R 3 v .
The form of the collisional operator C is not detailed here but only the properties used in this study will be detailed.

Proposition 1. The kinetic equation ( 1) written in a moving velocity framework writes

∂ t g(t, x, c) + div x ((c + u)g(t, x, c)) (2)
-div c (∂ t u + ∂ x u(c + u))g(t, x, c) = C(g(t, x, c)),
where the new velocity c is defined by

c = v -u(t, x),
with u a relative velocity which needs be specified.

Proof. We consider the vector

F (α) ∈ R + t × R 3 x × R 3 v defined as F (α) =   f (α) vf (α) 0   ,
then the kinetic equation (1) rewrites under the following form

div α (F (α)) = C(f (α)). Introducing φ ∈ C ∞ c (R + t × R 3 x × R 3 v ; R)
, we consider the associated weak formulation

α∈R + t ×R 3 x ×R 3 v div α (F (α)) -C(f (α)) φ(α)dα = 0, which rewrites α∈R + t ×R 3 x ×R 3 v ∂ α φ(α)F (α) + C(f (α))φ(α)dα = 0. (3)
In order to derive the kinetic equation in a moving velocity frame the following set of coordinates is considered

t β = (τ, y, c) = (t, x, v -u(t, x)),
where u is a relative velocity which depends of time and space. In particular, we define the C 1 -diffeomorphism Φ as

β = Φ(α).
We remark here that only the velocity coordinates are transformed while time and space coordinates are kept unchanged. Eq. ( 3) considered in the new set of variables rewrites

β∈R + t ×R 3 y ×R 3 c ∂ α φ(Φ -1 (β))F (Φ -1 (β)) + C(f (Φ -1 (β)))φ(Φ -1 (β)) (4)
| det J Φ -1 |dβ = 0,
where det J Φ -1 is the determinant of the Jacobian matrix of the transformation. For this transformation det J Φ -1 = 1, indeed the Jacobian matrix J Φ reads

J Φ =           1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 -∂ t u 1 -∂ x 1 u 1 -∂ x 2 u 1 -∂ x 3 u 1 1 0 0 -∂ t u 2 -∂ x 1 u 2 -∂ x 2 u 2 -∂ x 3 u 2 0 1 0 -∂ t u 3 -∂ x 1 u 3 -∂ x 2 u 3 -∂ x 3 u 3 0 0 1           .
The following quantities expressed in the mobile framework are introduced

G(β) = F (Φ -1 (β)), g(β) = f (Φ -1 (β)) and Ψ(β) = φ(Φ -1 (β)).
In order to express ∂ α φ(Φ -1 (β)), one remarks that

∂ β Ψ(β) = ∂ α φ(Φ -1 (β))∂ β (Φ -1 )(β).
Then by derivation of a reciprocal function

∂ β Φ -1 (β) = (∂ α Φ(Φ -1 (β))) -1 , it follows that ∂ α φ(Φ -1 (β)) = ∂ β Ψ(β)(∂ α Φ(Φ -1 (β))).
Using the two previous equations, equation [START_REF] Berthon | An HLLC Scheme to Solve The M1 Model of Radiative Transfer in Two Space Dimensions[END_REF] rewrites

β∈R + t ×R 3 y ×R 3 c ∂ β Ψ(β)∂ α Φ(Φ -1 (β))G(β)) + C(g(β))Ψ(β)dβ.
Finally for all Ψ(β)

∈ C ∞ c (R + τ × R 3 y × R 3 c ; R) β∈R + t ×R 3 y ×R 3 c div β (∂ α Φ(Φ -1 (β))G(β)) -C(g(β)) Ψ(β)dβ = 0, it follows that div β (∂ β Φ(Φ -1 (β))G(β)) = C(g(β)).
In the case of the present change of variables

∂ β Φ(Φ -1 (β))G(β) =    g(β) vg(β) -( ∂u ∂t + ∂u ∂x v)g(β)    .
Finally, using the fact that v = c + u(x, t), one obtains the kinetic equation in a moving velocity frame

∂ τ g(β) + div y ((c + u)g(β)) -div c (∂ τ u + ∂ y u(c + u))g(β) = C(g)(β).
Since time and space are unchanged by the change of variable

∂ t g(t, x, c) + div x ((c + u)g(t, x, c)) -div c (∂ t u + ∂ x u(c + u))g(t, x, c) = C(g(t, x, c)).
Eq.( 2) is used in the next sections to derive the angular M 1 moments model in a moving reference frame. Of course, an additional evolution equation is required to compute the velocity u. In this work, the velocity u is chosen as the particles mean velocity in the fixed frame (laboratory frame). This choice enables the reduction of the size of the velocity grids and plays an important role in enforcing the Galilean invariance property of the angular M 1 model. This point is presented in details in the next sections.

In order to derive the evolution equation for u, the kinetic equation ( 1) is integrated in velocity. This leads to the following conservation laws ∂n ∂t

+ div x (nu) = 0, (5)
∂(nu) ∂t + div x ( v f v ⊗ vdv) = 0.
Injecting the following expansion into (5)

v ⊗ v = (v -u) ⊗ (v -u) + (v -u) ⊗ u + u ⊗ (v -u) + u ⊗ u,
and by using the following identities

v u ⊗ uf dv = nu ⊗ u, v u ⊗ (v -u)f dv = u ⊗ v (v -u)f dv = 0, v (v -u) ⊗ uf dv = v (v -u)f dv ⊗ u = 0,
one obtains the evolution equation for u expressed in the new frame quantities

∂(nu) ∂t + div x (nu ⊗ u) + div x ( v g(c)c ⊗ cdc) = 0, (6)
where n = c g(c)dc.

The equations ( 5) and (6) will be used to compute the relative velocity u at each time step.

M 1 angular moments model in a moving frame

The M 1 angular moments model in a moving frame is derived by performing an angular moments extraction of the kinetic equation [START_REF] Benzonie-Gavage | Multi-dimensonal Hyperbolic Partial Differential Equations[END_REF]. One defines the following three first angular moments of the distribution function g

g 0 (ζ) = ζ 2 S 2 g(Ω, ζ)dΩ, g 1 (ζ) = ζ 2 S 2 g(Ω, ζ)ΩdΩ, g 2 (ζ) = ζ 2 S 2 g(Ω, ζ)Ω⊗ΩdΩ,
where S 2 is the unit sphere. The complete derivation of the M 1 angular moments model in a moving frame is presented in appendix. Removing the collisional operators contribution, this model reads

               ∂ t g 0 + div x (ζg 1 + ug 0 ) -∂ ζ du dt .g 1 + ζ∂ x u : g 2 = 0, ∂ t g 1 + div x (ζg 2 + u ⊗ g 1 ) -∂ ζ g 2 du dt + ζg 3 ∂u ∂x + g 0 Id -g 2 ζ du dt + ∂u ∂x g 1 -g 3 ∂u ∂x = 0, (7)
where du dt is defined as

du dt = ∂u ∂t + ∂u ∂x u,
and the third order moments g 3 as

g 3 (ζ) = ζ 2 S 2 g(Ω, ζ)Ω ⊗ Ω ⊗ ΩdΩ. ( 8
)
The evolution law [START_REF] Charrier | Multigroup model for radiating flows during atmospheric hypersonic re-entry[END_REF] expressed in terms of the angular moments rewrites

∂(nu) ∂t + div x (nu ⊗ u) + div x ( +∞ 0 g 2 (ζ)ζ 2 dζ) = 0. (9)
In order to close the system (7), the higher order moments g 2 and g 3 must be expressed in terms of g 0 and g 1 . Following [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF] and considering the M 1 closure, the distribution function from which the angular moments are derived writes

g(t, x, ζ, Ω) = exp(a 0 (t, x, ζ) + a 1 (t, x, ζ).Ω), (10)
where a 0 is a scalar function and a 1 a vector valued function. The coefficients a 0 and a 1 are the Lagrange multipliers associated to entropy minimisation problem [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF][START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF]. Then extending the ideas of [START_REF] Dubroca | Étude théorique et numérique d'une hiéarchie de modèles aux moments pour le transfert radiatif[END_REF][START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF][START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF] one can show that the closure relation for g 2 is given by

g 2 = g 0 3χ(α) -1 2 g 1 |g 1 | ⊗ g 1 |g 1 | + 1 -χ(α) 2 Id , (11)
where

χ(α) = 1 + |α| 2 + |α| 4 3 , α = g 1 /g 0 . (12)
Similarly the higher order moment g 3 reads

g 3 = 3|g 1 | -χ 2 g 0 2 g 1 |g 1 | ⊗ g 1 |g 1 | ⊗ g 1 |g 1 | + χ 2 g 0 -|g 1 | 2 g 1 |g 1 | ∨ Id , (13)
with

χ 2 (α) = 3|α| -|α| 3 + 3|α| 5 5 ,
and

g 1 |g 1 | ∨ Id = g 1 |g 1 | ⊗ e 1 ⊗ e 1 + e 1 ⊗ g 1 |g 1 | ⊗ e 1 + e 1 ⊗ e 1 ⊗ g 1 |g 1 | + g 1 |g 1 | ⊗ e 2 ⊗ e 2 + e 2 ⊗ g 1 |g 1 | ⊗ e 2 + e 2 ⊗ e 2 ⊗ g 1 |g 1 | + g 1 |g 1 | ⊗ e 3 ⊗ e 3 + e 3 ⊗ g 1 |g 1 | ⊗ e 3 + e 3 ⊗ e 3 ⊗ g 1 |g 1 | .
Before studying the models properties, the realisability conditions associated to the model ( 7) are introduced

A = (g 0 , g 1 ) ∈ R 2 , g 0 ≥ 0, |g 1 | ≤ g 0 . (14)
Since the distribution function g is a nonnegative quantity the realisability conditions ( 14) naturally needs to be satisfied. In addition, these conditions are related to the existence of a nonnegative distribution function from which the angular moments can be derived [START_REF] Pomraning | Maximum entropy Eddington factors and flux limited diffusion theory[END_REF].

Model properties

In this section the main properties of the angular M 1 model in a moving frame ( 7)-( 9) are presented. It is first proved that the choice of working in the mean velocity frame enables to ensure the Galilean invariance property of the model. Secondly it is shown that this model, rewritten in terms of entropic variables, is Friedrichs-symmetric. Finally, the derivation of the conservation laws is detailed.

Galilean invariance property

Galilean invariance is a fundamental feature of the Boltzmann equation. Following [START_REF] Junk | Maximum entropy moment systems and Galilean invariance[END_REF], we start defining translational and rotational transformations. For any vector s ∈ R d and any rotation matrix R ∈ SO(d)

(T s f )(v) = f (v -s), (T R f )(v) = f (Rv), v ∈ R d .
The following translational and rotational invariance properties of the collisional operator C are considered

T s C(f ) = C(T s (f )), T R C(f ) = C(T R (f )). ( 15
)
Note that the Boltzmann collision operator or the BGK collision operator [START_REF] Gross | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF] satisfy such properties. Using equation [START_REF] Grad | On the kinetic theory of rarefied gases[END_REF] the Galilean invariance property of the kinetic equation ( 1) can be shown. Indeed, the reference coordinates system (t, x, v) and a new set of coordinates (t, x, ṽ) can be linked by the following relations

x = Rx -st, ṽ = Rv -s, (16)
for any constant vector s ∈ R d . Distribution function f in the moving frame is defined as f (t, x, ṽ) = f (t, x, v).

Consequently the following relations can be derived

∂ t f (t, x, v) = ∂ t f (t, x, ṽ) -s.∂ x f (t, x, ṽ), ∂ x f (t, x, v) = R∂ x f (t, x, ṽ).
Therefore using (1), it follows that f satisfies

∂ t f (t, x, c) + div x(ṽ f (t, x, c)) = C( f (t, x, c)), (17)
which shows the Galilean invariance of (1).

The same property cannot be directly obtained when considering the M 1 angular moments model ( 7)-( 9). Indeed, when integrating (1) on the unit sphere and applying the change of variables ( 16) lead to unconvenient nonlinear terms and we are not able to show that the form of the M 1 model is invariant by translational transformation. In order to overcome this drawback, in this study, we propose to not derive the M 1 angular moments model from the kinetic equation ( 1) but from the kinetic equation ( 2) which is expressed in a mobile reference frame. In particular, in this work the velocity u used in ( 2) is chosen as the particles mean velocity defined by

u = 1 n R 3 f (v)vdv. ( 18
)
In order to show the advantage in deriving the M 1 angular moments model from the kinetic equation ( 2), the kinetic equation ( 17) is rewritten in its mean velocity frame. This second kinetic equation expressed in a mobile frame reads

∂ t g(t, x, c) + div x((c + ũ)g(t, x, c)) ( 19
)
-div c (∂ t ũ + ∂ x ũ(c + ũ))g(t, x, x) = C(g(t, x, c)), where ũ = 1 n R 3 f (ṽ)ṽdṽ. ( 20
)
The key point which will be useful when considering angular moments models is the relation between the two kinetic equations ( 2) and [START_REF] Guisset | Classical transport theory for the collisional electronic M1 model[END_REF]. Indeed, the two relative velocities u and ũ are linked through the following relation

ũ = Ru -s. (21)
Therefore, we propose to consider the following change of variable

x = Rx -st, c = Rc, ũ = Ru -s. (22)
Indeed, by injecting the change of variable ( 21)-( 22) into (2) and using the following relations

∂ t u = t R(∂ t ũ -(∂ x ũ)s), ∂ x u = t R(∂ x ũ)R, (23)
and

∂ t g = ∂ t g -(∂ x g)s, ∂ x g = (∂ x g)R, ∂ c g = (∂ c g)R, (24)
a direct calculation enables to recover equation [START_REF] Guisset | Classical transport theory for the collisional electronic M1 model[END_REF]. The relationships between the different studied frameworks are summarised on Fig 1 . The starting point is the kinetic equation ( 1) expressed in the fixed frame, denoted A 0 . Since this kinetic equation is Galilean invariant, one obtains (17) denoted B 0 , by using [START_REF] Gross | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF]. Secondly, the kinetic equation in a mobile frame (2) denoted A has been derived. In the present case u is the particles mean velocity defined in [START_REF] Guisset | Asymptotic-preserving scheme for the Fokker-Planck-Landau-Maxwell system in the quasi-neutral regime[END_REF]. The same procedure can be applied on [START_REF] Groth | Towards physically-realizable and hyperbolic moment closures for kinetic theory[END_REF] to obtain [START_REF] Guisset | Classical transport theory for the collisional electronic M1 model[END_REF], denoted B. Finally, one remarks that ( 2) and ( 19), denoted A and B are linked by the change of variable [START_REF] Junk | Maximum entropy moment systems and Galilean invariance[END_REF]. The change of variables [START_REF] Junk | Maximum entropy moment systems and Galilean invariance[END_REF] x = Rx -st ṽ = Rv -s 2) and ( 19) relevant when considering angular moments models. Indeed, this change of variable also enables to link the angular M 1 model derived from the kinetic equation ( 2) to the angular M 1 model derived from the kinetic equation ( 19). This point is detailed in the following result.

x = Rx -st c = Rc ũ = Ru -s A 0 : fixed frame ∂ t f + div x (vf ) = C(f ) A : Mobile frame ∂ t g + div x ((c + u)g) -div c (∂ t u + ∂ x u(c + u))g = C(g) B 0 : Uniform translation frame ∂ t f + div x(ṽ f ) = C( f ) B : Uniform translation mobile frame ∂ t g + div x((c + ũ)g) -div c (∂ t ũ + ∂ x ũ(c + ũ))g = C(g) c = v -u(t, x) nu(t, x) = v f vdv c = ṽ -ũ(t, x) nũ(t, x) = ṽ f ṽdṽ

Theorem 1. (Galilean invariance property)

The form of the M 1 angular moments model [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF] expressed in the mean velocity frame is invariant by rotational and translational transformations. Using the fact that

ζ = ζ, Ω = RΩ, (25)
and the equations ( 24), the following relations can be derived

g0 = g 0 , ∂ t g 0 = ∂ t g0 -∂ x g0 s, ∂ x g 0 = ∂ x g0 R, ∂ ζ g 0 = ∂ ζ g0 , (26)
and

g 1 = Rg 1 , ∂ t g 1 = t R(∂ t g1 -∂ x g1 s), ∂ x g 1 = t R∂ x g1 R, ∂ ζ g 1 = t R∂ ζ g1 . (27)
Using the definition of g2 , we remark that

g2 = ζ2 S2 g(t, x, c) Ω ⊗ Ωd Ω = ζ 2 S2 g(t, x, c)R Ω ⊗ Ω t Rd Ω. ( 28
)
Then injecting ( 26)-( 27) into the first equation of (7) and using ( 28) and ( 23) a direct calculation gives

∂ t g0 + div x ( ζ g1 + ũg 0 ) -∂ ζ dũ dt .g 1 + ζ∂ x ũ : g2 = 0.
In order to deal with the second equation of ( 7), we remark that using [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF] the ijk th component of the higher order moments g 3 defined in (8) rewrites

g 3ijk = l,m,n ζ2 S2 t R il t R jm t R kn ( Ω ⊗ Ω ⊗ Ω) lmn d Ω. ( 29
)
Therefore, using ( 28) and (29) a direct calculation gives

div x g 2 = t R div x g2 , (30)
and

g 3 ∂u ∂x = t R g3 ∂ ũ ∂ x . (31)
Consequently injecting ( 26)-( 27) into the second equation of ( 7) and using the relations ( 23)-( 30) and ( 31), one obtains

       ∂ t g1 + div x( ζ g2 + ũ ⊗ g1 ) -∂ ζ g2 dũ dt + ζ g3 ∂ ũ ∂ x + g0 Id -g2 ζ dũ dt + ( ∂ ũ ∂ x g1 -g3 ∂ ũ ∂ x = 0.

Symmetrization property

In this section it is shown that the M 1 model in a moving frame ( 7)-( 9) written in terms of the entropic variables is Friedrichs symmetric. Following [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF], the M 1 model in a moving frame [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF] can be rewritten in terms of the entropic variables a 0 and a 1 . This procedure is sometimes called a Godunov's symmetrisation [START_REF] Goudon | Analysis of the M1 model: well-posedness and diffusion asymptotics[END_REF].

Theorem 2. The M 1 model in a moving frame ( 7)-( 9) written in terms of the variables a 0 and a 1 is Friedrichs symmetric.

Proof. Setting t m = (1, Ω), t α = (α 0 , α 1 ), the distribution function [START_REF] Filbet | A hierarchy of hybrid numerical methods for multi-scale kinetic equation[END_REF] reads

g(t, x, ζ, Ω) = exp(α.m),
and the solution of ( 7)-( 11)-(13) writes

t (g 0 , g 1 ) = ζ 2 exp(α.m)m ,
where the notation . refers to the angular integration on the unit sphere. Consequently, after a direct calculation, the M 1 angular moments model in a moving frame [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF] rewrites

A 0 (α)∂ t α 0 α 1 + j A j (α)∂ x j α 0 α 1 + B(α)∂ ζ α 0 α 1 + S(x, ζ, α) = 0 0 , (32)
where

A 0 (α) =< exp(α.m) 1 t Ω Ω Ω ⊗ Ω >, A j (α) =< (ζΩ j + u j ) exp(α.m) 1 t Ω Ω Ω ⊗ Ω >, B(α) =< -(ζ 2 du dt .Ω + ζ 3 ∂u ∂x : Ω ⊗ Ω) exp(α.m) 1 t Ω Ω Ω ⊗ Ω >, and
S(x, ζ, α) =    (div x u)g 0 - 2 ζ du dt .g 1 -3 ∂u ∂x : g 2 ∂u ∂x g 1 - 2 ζ g 2 du dt -3g 3 ∂u ∂x + g 0 Id -g 2 ζ du dt + ( ∂u ∂x g 1 -g 3 ∂u ∂x )    .
Since A 0 (α) is a positive-definite symmetric matrix and A j (α) and B(α) are symmetric matrices, one obtains that the system (32) is Friedrichssymmetric [START_REF] Friedrichs | Systems of Conservation Equations with a Convex Extension[END_REF][START_REF] Benzonie-Gavage | Multi-dimensonal Hyperbolic Partial Differential Equations[END_REF].

Conservation laws

In this section the derivation of the conservation laws derived from the angular M 1 model in a moving frame ( 7) is detailed. Before deriving the mass and energy conservation equations, we point out that in this work the velocity u is chosen as the particles mean velocity. Therefore, in the considered framework the mean velocity is equal to zero. This point is expressed by the following condition

+∞ 0 g 1 (t, x, ζ)ζdζ = 0. ( 33
)
Multiplying the second equation of ( 7) by ζ and integrating in ζ, one shows using Green's formula that all the terms vanish two by two and that condition [START_REF] Rey | An Exact Rescaling Velocity Method for some Kinetic Flocking Models[END_REF] is preserved over times. The derivation of the mass conservation equation can be directly obtained by direct integration in ζ. Indeed, integrating the first equation of ( 7) in ζ, one obtains

∂ t n + div x (nu) = 0, (34)
where condition [START_REF] Rey | An Exact Rescaling Velocity Method for some Kinetic Flocking Models[END_REF] has been used. In order to derive the energy conservation equation, one starts multiplying the first equation of ( 7) by m 2 ζ 2 and integrate in ζ to obtain the following internal energy equation

∂ t ( 1 2 +∞ 0 g 0 ζ 2 dζ) + div x ( 1 2 +∞ 0 g 1 ζ 3 dζ + u 1 2 +∞ 0 g 0 ζ 2 dζ) (35)
+ (∂ x u : +∞ 0 g 2 ζ 2 dζ) = 0.
One notices that since the mean velocity frame is considered, only an equation on the internal energy is obtained. The kinetic energy equation is derived from the evolution equation ( 9) and writes

∂ t (nu 2 ) + div x ( nu 2 2 u) + u.div x ( +∞ 0 g 2 ζ 2 dζ) = 0. ( 36
)
The energy conservation equation is directly obtained by summing equation [START_REF] Ripoll | A factored operator method for solving coupled radiation-hydrodynamics models[END_REF] with equation (36).

Numerical scheme

In this part an appropriate numerical scheme is proposed for the M 1 model in a moving framework in an one dimensional spatial geometry considering a standard BGK collision operator [START_REF] Gross | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF]. In this case, the collisional operator

C(f ) used in (1) is specified C(f ) = 1 τ (M f -f ), with M f (v) = n (2πT ) 3/2 exp(- (v -u) 2 2T ),
and τ is a collisional parameter which is fixed depending of the collisional regime studied. In this case the M 1 model in a moving framework [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF] writes

               ∂ t g 0 + ∂ x (ζg 1 + ug 0 ) -∂ ζ du dt g 1 + ζ ∂u ∂x g 2 = 1 τ (M g 0 -g 0 ), ∂ t g 1 + ∂ x (ζg 2 + ug 1 ) -∂ ζ du dt g 2 + ζ ∂u ∂x g 3 + du dt g 0 -g 2 ζ + ∂u ∂x (g 1 -g 3 ) = - 1 τ g 1 , (37)
where

M g 0 = 4πζ 2 n (2πT ) 3/2 exp(- ζ 2 2T
).

Derivation of the numerical scheme

In order to derive a suitable numerical scheme for the model (37) which preserves the admissibility of the solution, the different terms of (37) are studied separately. Then the admissibility requirement of the complete scheme is shown under a reduced CFL condition.

Step 1: the first intermediate state is the following

∂ t g 0 + ∂ x (ζg 1 + ug 0 ) = 0, ∂ t g 1 + ∂ x (ζg 2 + ug 1 ) = 0. ( 38
)
In order to derive a numerical scheme preserving the realisability of the numerical solution, we consider an underlying kinetic model from which the system (38) can be derived by direct angular moments extraction

∂ t F (t, x) + ∂ x (a(x)F (t, x)) = 0, (39)
with

F = ζ 2 f , a(x) = ζµ + u(x) and µ ∈ [-1, 1]
. Note that µ is the angular variable in the case of one space dimension.

A natural conservative numerical scheme is proposed for the kinetic equation ( 39)

F n+1 i -F n i ∆t + h n i+1/2 -h n i-1/2 ∆x = 0, ( 40
) with h n i+1/2 = a - i+1/2 F n i+1 + a + i+1/2 F n i ,
and

a ± = 1 2 (a ± |a|).
Rewriting equation [START_REF] Turpault | Multigroup half space moment appproximations to the radiative heat transfer equations[END_REF] as a convex combination

F n+1 i = F n i 1 - ∆t 2∆x (a i+1/2 -a i-1/2 ) -∆t |a i+1/2 | + |a i-1/2 | 2∆x + F n i+1 2∆t ∆x |a i+1/2 | -a i+1/2 (41) + F n i-1 2∆t ∆x |a i-1/2 | -a i-1/2 ,
it follows that the positivity of the numerical distribution function is ensured under the following CFL condition

∆t 1 ≤ ∆x 2||u|| ∞ + ζ . ( 42
)
The numerical scheme (41) rewrites on the following viscous form

F n+1 i -F n i ∆t + a i+1/2 F n i+1 + (a i+1/2 -a i-1/2 )F n i -a i-1/2 F n i-1 2∆x (43) - |a i+1/2 |F n i+1 -(|a i+1/2 | + |a i-1/2 |)F n i + |a i-1/2 |F n i-1 2∆x = 0.
The angular integration can not be directly performed on the scheme (43) because of the angular variable µ which appears in the term |a| in the numerical viscosity. Therefore we modify (43) and consider the following scheme which is suitable for the angular integration.

F n+1 i -F n i ∆t + a i+1/2 F n i+1 + (a i+1/2 -a i-1/2 )F n i -a i-1/2 F n i-1 2∆x (44) -||a|| ∞ F n i+1 -2F n i + F n i-1 2∆x = 0.
Remark: Considering (44), one observes that the numerical viscosity of the scheme is increased in order to enable the angular integration.Therefore the numerical scheme still preserves the nonnegativity of the numerical solution under CFL condition (42).

The angular integration of the scheme (44) leads to a natural discretisation for the intermediate state ( 38) , one can show the scheme (45) preserves the realisability requirement of the numerical solution under the CFL condition (42).

g n+1 0i -g n 0i ∆t + (ζg n 1i+1 + u i+1/2 g n 0i+1 ) + ((ζg n 1i + u i+1/2 g n 0i ) 2∆x - (ζg n 1i + u i-1/2 g n 0i )) -(ζg n 1i + u i-1/2 g n 0i-1 ) 2∆x -(ζ + ||u|| ∞ ) g n 0i+1 -2g n 0i + g n 0i-1 2∆x = 0, g n+1 1i -g n 1i ∆t + (ζg n 2i+1 + u i+1/2 g n 1i+1 ) + ((ζg n 2i + u i+1/2 g n 1i ) 2∆x (45)
- (ζg n 2i + u i-1/2 g n 1i )) -(ζg n 2i + u i-1/2 g n 1i-1 ) 2∆x -(ζ + ||u|| ∞ ) g n 1i+1 -2g n 1i + g n 1i-1 2∆x = 0.
Step 2: the second intermediate step we consider writes

     ∂ t g 0 -∂ ζ ( du dt g 1 + ζ∂ x ug 2 ) = 0, ∂ t g 1 -∂ ζ (g 2 du dt + ζg 3 ∂ x u) = 0. ( 46
)
Following the same procedure than for the first intermediate state, the following underlying kinetic model is proposed

∂ t F (ζ) -∂ ζ (( du dt µ + ζ∂ x uµ 2 )F (ζ)) = 0,
with the following corresponding scheme

F n+1 j -F n j ∆t + b j+1/2 F n j+1 + (b j+1/2 -b j-1/2 )F n j -b j-1/2 F n j-1 2∆ζ (47) -||b|| ∞ F n j+1 -2F n j + F n j-1 2∆ζ = 0, with b = du dt µ + ζ∂ x uµ 2 . The CFL condition associated reads ∆t 2 ≤ ∆ζ 2(|| du dt || ∞ + ζ||∂ x u|| ∞ ) . ( 48
)
The angular integration of (47) leads to the following discretisation for the intermediate state ( 46)

g n+1 0j -g n 0j ∆t + ( du dt g n 1i+1 -ζ j+1/2 ∂ x ug n 2j+1 ) + ( du dt g n 1j -ζ j+1/2 ∂ x ug n 2j ) 2∆ζ - ( du dt g n 1j -ζ j-1/2 ∂ x ug n 2j )) + ( du dt g n 1j-1 -ζ j-1/2 ∂ x ug n 2j-1 ) 2∆ζ -(| du dt | + ||ζ|| ∞ |∂ x u|) g n 0j+1 -2g n 0j + g n 0j-1 2∆ζ = 0, g n+1 1j -g n 1j ∆t + ( du dt g n 2j+1 -ζ j+1/2 ∂ x ug n 3j+1 ) + ( du dt g n 2j -ζ j+1/2 ∂ x ug n 3j ) 2∆ζ (49)
-

( du dt g n 2j -ζ j-1/2 ∂ x ug n 3j ) -( du dt g n 2j-1 -ζ j-1/2 ∂ x ug n 3j-1 ) 2∆ζ -(| du dt | + ||ζ|| ∞ |∂ x u|) g n 1j+1 -2g n 1j + g n 1j-1 2∆ζ = 0.
Remark: The scheme (49) preserves the realisability domain under the CFL condition (48).

Step 3: the third state we consider is the following

   ∂ t g 0 = 0, ∂ t g 1 + g 0 -g 2 ζ du dt = 0.
We choose the following classical scheme for this first model

     g n+1 0ij = g n 0ij , g n+1 1ij = g n 1ij -∆t g 0ij -g 2ij ζ j ( du dt ) i .
Remark: This scheme preserves the realisability conditions under CFL conditions ∆t 3 ≤ ζ

| du dt | 1 + α 1 -χ(α) ,
where α is defined by [START_REF] Filbet | A Rescaling Velocity Method for Kinetic Equations: the Homogeneous Case[END_REF].

Proof. This result is directly obtained by computing g n+1 0i ± g n+1 1i .

Remark:

The term 1 + α 1 -χ(α)
does not tend to zero as α tends to -1. Indeed, using the definition of χ given in [START_REF] Filbet | A Rescaling Velocity Method for Kinetic Equations: the Homogeneous Case[END_REF], one can show that 1 + α 1 -χ(α) tends to 1/2 as α tends to -1.

Step 4: the fourth intermediate step we consider writes

∂ t g 0 = 0, ∂ t g 1 + ∂ x u(g 1 -g 3 ) = 0.
Following the third step we propose

g n+1 0i = g n 0i , g n+1 1i = g n 1i + ∆t(∂ x u) i (g 1i -g 3i ).
Remark: This scheme preserves the realisability conditions under CFL conditions ∆t 4 ≤ 1

|∂ x u| 1 + α α -χ 2 (α)
.

Using the definition of χ 2 , we remark that 1 + α α -χ 2 (α) tends to -1/2 as α tends to -1.

In order to derive a admissible numerical scheme for the complete model [START_REF] Schneider | Entropic approximation in kinetic theory[END_REF], we propose to consider the following time semi-discretisation

U n+1 = U n + ∆t N k=1 F k (U n ), (50)
where

U n+1 = g n+1 0 g n+1 1 .
F k represents the discretisation proposed for the k th intermediate step and N is the number of intermediate step considered. Equation (50) rewrites under the form of a convex combination

U n+1 = N k=1 1 N [U n + (N ∆t)F k (U n )]. ( 51
)
Setting ∆t = N ∆t, one shows that if each intermediate step

Ũ n+1 = U n + ∆tF k (U n ),
preserves the realisability conditions of the numerical solution under CFL condition ∆t ≤ C k .

Therefore the general scheme (50) preserves the realisability conditions of the numerical solution under the following CFL condition

∆t ≤ min k ( C k N ).
The following result is then obtained Theorem 3. The general scheme (50) preserves the realisability conditions under the following CFL condition

∆t ≤ 1 4 min(∆t 1 , ∆t 2 , ∆t 3 , ∆t 4 ). ( 52
)
Proof. Each step preserves the realisability conditions under CFL condition. Therefore, by convexity of the admissible set, considering the convex combination (51) and using the condition (4.1), we directly obtain that the general scheme (50) preserves the realisability conditions under the CFL condition (52).

The discretisation of the collision operator is performed by using a standard implicit scheme. For the numerical test presented in the next section, an usual second order Van Leer's slope limiter [START_REF] Van Leer | Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow[END_REF] is used.

Enforcement of the discrete energy conservation and zero mean velocity condition

In this section, the enforcement of the discrete energy conservation and zero mean velocity condition is discussed. In a recent work [START_REF] Rey | An Exact Rescaling Velocity Method for some Kinetic Flocking Models[END_REF], a numerical scheme has been proposed to enforce the discrete zero mean velocity condition considering a kinetic equation. However, this strategy does not directly apply in the present case since a nonlinear set of equations ( 37) is considered associated to the realisability conditions [START_REF] Goudon | Analysis of the M1 model: well-posedness and diffusion asymptotics[END_REF]. The enforcement of the discrete energy conservation and the zero mean velocity condition while preserving realisability conditions (14) of the numerical solution is particularly challenging and beyond the scope of the present study. However, in order to be able to present numerical results, in this section a correction of the numerical solution is proposed.

In order to enforce the correct energy conservation, we start considering the following conservation laws associated to (37)

     ∂ t ρ + div x (ρu) = 0, ∂ t (ρu) + div x (ρu ⊗ u + p -s) = 0, ∂ t E + div x ((E + p -s)u + q) = 0, ( 53
)
where E is the total energy. The pressure tensor p, the stress tensor s and the heat flux q expressed in terms of the angular moments read

p -s = m 2 +∞ 0 g 2 ζ 2 dζ, q = m 2 +∞ 0 g 1 ζ 3 dζ.
At each time time step, the set of conservation laws (53) is numerically solved. Then the numerical solution is corrected by using

g 0p = α exp(βζ p 2 )ḡ 0p , ∀p ∈ {1; ...; pf },
where g 0 is the corrected solution and ḡ0 the solution which requires a correction computed with the scheme (50). The coefficients α and β are numerically computed such that

m pf p=1 g 0p ∆ζ = ρ, m 2 pf p=1 g 0p ζ 2 p ∆ζ = E - ρu 2 2 ,
where the quantities E, ρu 2 2 and ρ are known at each time step since the set (53) has been numerically solved. This procedure enables the enforcement of the correct energy conservation. As it will be shown in the next section this correction is important for the numerical results, in particular in order to numerically capture shock waves.

In order to enforce the zero mean velocity condition [START_REF] Rey | An Exact Rescaling Velocity Method for some Kinetic Flocking Models[END_REF] at the discrete level, one could think in proposing an adapted discretisation for the source terms which appears in the second equation of [START_REF] Schneider | Entropic approximation in kinetic theory[END_REF]. However, this procedure leads to an unsuitable CFL condition when considering the realisability requirements [START_REF] Goudon | Analysis of the M1 model: well-posedness and diffusion asymptotics[END_REF] for the numerical solution. Therefore the following correction is proposed based on the resolution of the convex optimisation problem min

g 1 ∈R pf 1 2 ||g 1 -ḡ1 || 2 L 2 = 0, under equality constraint pf p=1 g 1p ζ p ∆ζ = 0,
where g 1 is the corrected solution and ḡ1 the solution before correction given by the scheme (50). One observes that this procedure does not enforce the realisable conditions of the numerical solutions. In such unfortunate case, g 1 is simply projected on the realisable set.

Numerical results

In this section, several test cases are presented. Depending on the regime considered, the numerical results obtained with the scheme introduced in the previous part for the angular M 1 moments model in a moving frame, denoted M 1 mobile, are compared either with an exact solution or with a kinetic reference solution. The results are given with and without the correction procedure. In the following, the kinetic solution has been obtained considering a standard kinetic 1D3V BGK model using an usual Lax-Friedrichs scheme with the second order Van Leer's slope limiter [START_REF] Van Leer | Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow[END_REF]. The results obtained with this scheme are denoted BGK 1D3V. In addition, the results obtained considering a second order HLL scheme for the Euler equations using the second order Van Leer's slope limiter are also given. These results obtained using this scheme are denoted Euler.

Test 1: Temperature gradient test case in different collisional regimes.

The first test case we study consists in considering a strong temperature gradient at initial time and studying the temporal evolution of density, velocity and temperature. The initial distribution function is supposed to be a Maxwellian distribution function defined by

f ini (x, v) = n ini (x) (2πT ini (x)) 3/2 exp - (v -u ini (x)) 2 2T ini (x) , with n ini (x) = 1, u ini (x) = 0, T ini (x) = 2 -arctan(x).
The space range chosen is [-40, 40], and the velocity range [-15, 15] 3 . For the present test case, 400 cells in space and 200 3 cells in velocity have been considered for the 1D3V BGK kinetic approach. Also, 400 cells in space and 200 cells in velocity modulus have been considered for the M 1 mobile scheme. Finally, 400 cells in space have been considered for the Euler description.

Neumann boundary conditions are considered, the values in the boundary ghost cells set to the values in the corresponding real boundary cells.

1.a Fluid regime. The first regime we consider is the fluid regime. The collisional parameter τ is set equal to zero. In Figure 2, the density, velocity and temperature profiles are displayed at time t = 10 for the kinetic BGK 1D3V scheme in continuous blue, the M 1 mobile scheme in dashed green, the M 1 mobile scheme with correction in dashed blue and for the Euler scheme in dashedpoint pink. It is observed that all the schemes converge towards the same solution. This behaviour is expected since working in fluid regime the distribution remains a Maxwellian distribution function and the three descriptions give the same solution.

1.b Rarefied regime.

The second regime we consider is a rarefied regime where the collisional parameter τ is set equal to 1. In Figure 3, the density, velocity, temperature and heat flux profiles are displayed at time t = 10 for the kinetic BGK 1D3V scheme in continuous blue, the M 1 mobile scheme in dashed green, the M 1 mobile scheme with correction in dashed blue and for the Euler scheme in dashed-point pink. The Euler scheme gives the same results than in the previous case 1a. This is expected since the description is not able to distinguish the different regimes. In this case the heat flux is equal to zero. One observes that M 1 mobile scheme gives close results to the ones obtained with the kinetic BGK 1D3V scheme. When looking at the heat flux profiles, one observes that the general trends are qualitatively similar with some notable differences in the amplitude reached. Since the heat flux is a high order velocity moment, the differences between the models are particularly visible. The M 1 model is accurate in collisional regimes, however as pointed out in [START_REF] Guisset | Limits of the M1 and M2 angular moments models for kinetic plasma physics studies[END_REF], it can be inaccurate in collisionless regimes. The differences observed here, are due to the inaccuracy to the M 1 model in rarefied regime. 1.c Non-homogeneous collisional parameter. When considering realistic physical applications, the collisional parameter varies according to the gas conditions. Therefore, in the third case we consider that the collisional parameter τ is variable in space and is defined by

τ (x) = 1 2 (arctan(1 + 0.1x) + arctan(1 -0.1x)).
In Figure 4, the density, velocity, temperature and heat flux profiles are displayed at time t = 10 for the kinetic BGK 1D3V scheme in continuous blue, the M 1 mobile scheme in dashed green and for the Euler scheme in dashed-point pink. It is observed that the profiles obtained using the M 1 mobile sheme and the BGK 1D3V scheme are very close. One also observes that even the heat flux profiles are very similar. These results show the interest in using an angular moment model.

f ini (x, v) = n ini (x) (2πT ini (x)) 3/2 exp - (v -u ini (x)) 2 2T ini (x) , with (n ini (x), u ini (x), T ini (x)) = (1.00 • 10 -4 , 0, 4.80 • 10 -3 ) if x < 0, (1.25 • 10 -5 , 0, 3.84 • 10 -3 ) if x > 0.
The space range chosen is [0, 0.6], and the velocity range [-20, 20] 3 . For the present test case, 200 cells in space and 200 3 cells in velocity have been considered for the 1D3V BGK kinetic approach. Also, 200 cells in space and 200 cells in velocity modulus have been considered for the M 1 mobile scheme. Finally, 200 cells in space have been considered for the Euler description. Neumann boundary conditions are considered, the values in the boundary ghost cells are set to the values in the corresponding real boundary cells. For this test case, we consider the fluid regime therefore the collisional parameter τ is set equal to 0. For this test case an exact solution is known, a rarefaction wave, a contact discontinuity and a shock wave appear. In Figure 5, the mass density, velocity and temperature solution profiles are displayed at time t = 7.34 • 10 -2 . It is observed that the rarefaction wave (left side) is correctly captured by the M 1 mobile scheme (solution displayed in dashed green). However, one remarks that the shock amplitude is not correctly captured. It has been observed that this incorrect behaviour is due to the wrong discrete energy conservation. Indeed, by using the correction procedure introduced in the previous part the results in dashed blue are obtained, in this case the correct amplitude is recovered. We notice, the importance of the correct discrete energy conservation for capturing shock waves. This point is highlighted in the next test case.

Test 3: Double shock wave test case

The third test case we study is the double shock wave test case in fluid regime. The initial distribution function is supposed to be a Maxwellian distribution function defined by

f ini (x, v) = n ini (x) (2πT ini (x)) 3/2 exp - (v -u ini (x)) 2 2T ini (x) , with (ρ(x), u(x), T (x)) = (1, 2, 0.4) if x < 0, (1, -2, 0.4) if x > 0.
The space range chosen is [0, 1], and the velocity range [-15, 15] 200 cells in velocity modulus have been considered for the M 1 mobile scheme. Finally, 200 cells in space have been considered for the Euler description. Neumann boundary conditions are considered, the values in the boundary ghost cells are set to the values in the corresponding real boundary cells. For this test case, we consider the fluid regime therefore the collisional parameter τ is set equal to 0. For this test case an exact solution is known, two shock waves are created. In Figure 6, the mass density, velocity and temperature solution profiles are displayed at time t = 0.15. Similarly as remarked in the previous test case, it is observed that the M 1 scheme does not capture the correct amplitude profile nor the correct shock positions (results in dashed green). The results displayed in dashed blue are obtained using the corrected scheme. It is observed that the correction enables to correctly captures the shock profiles. This example confirms the importance of the discrete energy conservation.

Conclusion

In this work, the M 1 angular moments model in the particles mean velocity frame has been derived. Several fundamental properties of the model have been presented. In particular, the importance of working in the mean velocity frame has been highlighted. Indeed, this choice of framework is relevant when considering the Galilean invariance property of angular moments models. The derivation of the associated conservation laws has been detailed in addition to the zero mean velocity condition. A numerical scheme preserving the realisable sets has been proposed and validated with numerical test cases in different collisional regimes. Also, the importance of the correct discrete energy conservation has been emphasized.

As a short term perspective, one needs to derive a numerical scheme enforcing the discrete energy conservation and the zero mean velocity condition. Such an issue is challenging since it should be done preserving the realisable property of the numerical solution. As long term perspective, it would be interesting to study the motion of charged particles. One could consider the electron particle transport working in the ion mean velocity framework. This choice would enable a great simplification of the electronion collisional operator and an important step toward the multispecies particle transport for plasma physics applications.

Appendix: derivation of the angular M 1 model in a moving frame [START_REF] Dubroca | Angular moment model for the Fokker-Planck equation[END_REF] In this section, the derivation of the angular M 1 model in a moving frame is detailed. The kinetic equation ( 2) is considered for the angular integration. Introduce a test function φ, we consider the integral Finally one obtains the second equation of ( 7)

∂

Figure 1 :

 1 Figure 1: Diagram presenting the relations between the different frames

Proof.

  Before showing the Galilean invariance property of the M 1 angular moments model (7), we define the quantities in the new frame. Consider the velocity modulus ζ and Ω the angular direction in the mobile frame ζ = |c|, c = ζ Ω, we define the two first angular moments g0 and g1 in the mobile frame g0 = ζ2 S2 g(t, x, c)d Ω, g1 = ζ2 S2 g(t, x, c) Ωd Ω.

Remark: Computing g n+1 0 +g n+1 1 and g n+1 0 -g n+1 1

 101

Figure 2 :

 2 Figure 2: Test 1a -Solution profiles obtained for the temperature gradient test case with τ = 0 at time t = 10.

Figure 3 :

 3 Figure 3: Test 1b -Solution profiles obtained for the temperature gradient test case with τ = 1 at time t = 10.

Test 2 :

 2 Sod tube test case in fluid regimeThe second test case we study is the Sod tube test case in fluid regime. The initial distribution function is supposed to be a Maxwellian distribution

Figure 4 :

 4 Figure 4: Test 1c -Solution profiles obtained for the temperature gradient test case with variable collisional parameter at time t = 10.

Figure 5 :

 5 Figure 5: Test 2 -Sod tube test case with τ = 0 at time t = 7.34 • 10 -2 .

Figure 6 :

 6 Figure 6: Test 3 -Double shock wave test case with τ = 0 at time t = 0.15.

c∂ 2 (: ζg 2

 22 t g + div x ((c + u)g) -div c ( du dt g + ∂u ∂x cg φ(ζ)dc = 0. By using the Green formulae and c = ζΩ c (∂ t g + div x ((c + u)g))φ(ζ)dc + Ωζgζ)φ (ζ)dc = 0.In spherical coordinates the previous equation leads toζ S ∂ t g+div x ((c+u)g))φ(ζ)ζ 2 dΩdζ+ φ (ζ)ζ 2 dΩdζ = 0.By using the definitions of the angular momentsζ (∂ t g 0 + div x (ζg 1 + ug 0 ))φ(ζ)dζ + )φ (ζ)dζ = 0.Finally by integration by part,ζ ∂ t g 0 + div x (ζg 1 + ug 0 ) -∂ ∂ζ du dt .g 1 + ∂u ∂x : ζg 2 φ(ζ)dζ = 0.This holds true for all test function φ then one obtains the first equation of (7)∂ t g 0 + div x (ζg 1 + ug 0 ) -∂ ∂ζ du dt .g 1 + ∂u ∂x : ζg 2 = 0.Introduce a test function φΩ, we consider the integralc ∂ t g + div x ((c + u)g) -div c ( du dt g + ∂u ∂x cg) φ(ζ)Ωdc = 0.By using the Green formulaec ∂ t g + div x ((c + u)g) φ(ζ)Ωdc + ∂c = φ (ζ)Ω ⊗ Ω + φ(ζ) Id -Ω ⊗ Ω ζ .Then the second term of the left side of the equation gives c ζ)dζ. The third term of the right side leads to c Ω ⊗ Ω du dt φ (ζ)g(c)dc = -

  3 . For the present test case, 200 cells in space and 200 3 cells in velocity have been considered for the 1D3V BGK kinetic approach. Also, 200 cells in space and

								2
					Initial condition		Initial condition
					Exact solution		Exact solution
		4			BGK 1D3V			BGK 1D3V
					Euler M1 mobile			1	Euler M1 mobile
					M1 mobile corrected	M1 mobile corrected
	rho	3						velocity	0
		2						-1
		1						-2
		0	0,5 x			1		0	0,5 x	1
				2
								Initial condition
								Exact solution
								BGK 1D3V
			1,5				Euler M1 mobile
			Temperature	1				M1 mobile corrected
			0,5
				0	0,2	0,4	x	0,6	0,8	1

  t g 1 +div x (ζg 2 +u⊗g 1 )-∂ ζ g 2

	du dt	+ζg 3	∂u ∂x	+	f 0 Id -f 2 ζ	du dt	+(	∂u ∂x	g 1 -g 3	∂u ∂x	= 0.


Jean-Louis Guyader
email: jean-louis.guyader@insa-lyon.fr

Ha Dong Hwang

Kerem Ege

SANDWICH PATCH WITH THERMOVISCOUS FLUID CORE, FOR INCREASING DAMPING OF PANELS

The introduction of a thermo-Viscous Fluid (TVF) between two skins, produces a sandwich with high dissipation. This was studied by Hussain and Guyader 1 leading to equivalent material properties depending on skins and fluid core properties. An important result is that contrary to standard sandwich panels with polymers cores, the damping loss factor of such TVF sandwich is constant with frequency. When applying damping treatment on structures, bonded patches are generally used. In this paper the use of TVF sandwich patches on a panel is described through a model of heterogeneous plate, the predicted vibration response of the structure shows an increase of damping due to the patch. Estimation of modal damping loss factors for the plate with TVF sandwich patches are given depending on the patch size and TVF thickness.

Introduction

A method has been developed by Hussain and Guyader 1 to predict the damping behavior of an equivalent plate formed out of sandwich panels with thermoviscous fluid core, which is based on an asymptotic approach. The technique works on the breakup of the basic physical quantities of the fluid core such pressure, temperature and particle velocities as constants and linear functions of the z coordinate which signifies the perpendicular direction from the plate mid surface. The asymptotic modeling is performed in a very thin and highly thermoviscous fluid layer on the full set of linearized Navier Stokes equations. The use of highly viscous fluid layers yields high damping coefficients in all the frequency range and not strongly variable with temperature. Experimental validation provides a good agreement.

Equivalent thin plate equation

In order to consider the plate and the thermo-viscous fluid layer as a single structure, Hussain and Guyader 1 developed an asymptotic analysis of a double plate system with internal thermoviscous fluid. In the limit of thin fluid film, the transverse vibrations of both plates are identical and the fluid film motion can be expressed from that of the plate in order that a global plate equation can be derived:

( ) { } ( ) { } 0 , , = + y x W L y x W L F P (1)
Where the plate displacement is ( )

y x W ,
(identical for both skins of the sandwich),

( ) { } y x W L P , (resp. ( ) { } y x W L F ,
) is the operator associated to the plates (to the thermo-viscous fluid film).

Considering harmonic motions of angular frequency ω , these operators are expressed as:

( ) { } ( ) ( ) 2 2 1 0 2 2 1 1 2 , ∆ + + + + - = D D h h y x W L P δ ρ ρ ρ ω (2) ( ) { } ∆         + + = δµ χ ω 2 , 2 1 h h j y x W L F ( 3
)
Where ∆ is the Laplacian, the skin plates thicknesses are

3 1 + = B .
The equation (1) corresponds to that of an equivalent thin damped plate which is expressed as:

( ) ( ) { } 0 , 2 2 = ∆ + ∆ + - y x W j D M eq eq eq ωλ ω (4) Where ( ) δ ρ ρ ρ 0 2 2 1 1 + + = h h M eq , ( )
2 1 D D D eq + = and
( ) δµ χ λ + + = 2 2 1 h h eq
The equivalent mass is the sum of that of plates ( )

2 2 1 1 h h ρ ρ +
and of the added fluid mass δ ρ 0 . The equivalent rigidity is the sum of the plates bending stiffnesses ( )

2 1 D D +
. The damping term is associated to a Laplacian operator that does not appear in the standard plate equation. So, the equation for the equivalent sandwich plates with thermo-viscous fluid core is different from the standard Love Kirchhoff bending plate equation. The equivalent plate damping effect is related to the skin plate thicknesses and viscosity of the visco-thermic fluid.

Equivalent structural damping

In order to find out the structural damping factor of the equivalent plate let us consider a one dimensional plate solution of progressive wave form in the x direction whose displacement function given as:

( ) jkx Ae y x W = , (5)
Substituting this expression in the homogeneous equivalent plate equation yields:

0 2 4 2 = + + - k j k D M eq eq eq ωλ ω (6)
20th International Congress on Sound and Vibration (ICSV20), Bangkok,Thailand,[7][8][9][10][11] July 2013 ICSV20,Bangkok,Thailand,7891011Eq. ( 6) can be resolved for k giving the two values. With these values, Eq. ( 6) can be rewritten as:

0 2 4 2 2 2 2 2 2 2 =         + ∆ ∆ + + ∆ - ∆ + - k j D M eq eq eq eq eq eq λ ω ωλ λ ω λ ω ω (8)
Eq. ( 8) supposes that ∆ is real, this is true in the majority of cases however for very light materials with low rigidity, ∆ becomes imaginary and the following calculations must be changed accordingly. Let η be the loss factor due to the viscous layer, then this factor is identified as the ratio of the imaginary and the real stiffness: One important tendency appears here; the damping loss factor is not frequency dependant contrary to sandwich plates with visco-elastic cores, in addition the coefficients of viscosity of thermoviscous fluids are generally not temperature dependant in a wide range of temperature and the loss factor remains identical in a larger range of temperatures.

In the following numerical results are presented in order to show the efficiency of sandwich plates with thermo-viscous fluid core to get highly damped panels. Results are compared to standard sandwich panels with viscoelastic cores and the case of partial and total coverage of panels is discussed.

Modal Loss Factors of Sandwich Plates

Modal damping levels of a simply supported rectangular aluminum plate treated with a constrained layer damping (CLD) patch of a thermoviscous fluid core are computed, and the results are compared to the same plate structure of a viscoelastic polymer core. Dimensions of the base plate and the patch are 0.5*0.6*0.0005 (m 3 ) and 0.16*0.25 (m 2 ) respectively, and the location of the patch is arbitrarily selected. Aluminum and High-Density PolyEthylene (HDPE) are considered for the constraining layer materials, and their thicknesses are varied. The base plate thickness is kept constant (0.5 mm) for all cases. The modal damping loss factor of the composite layer is deduced from the equivalent single layer modeling developed by Guyader et al 2 . This methodology is implemented in MOVISAND software which computes the equivalent moduli of a single layer plate model for an infinite plate, and then obtained parameters are used to compute the modal damping loss factor of finite plate models. Nastran software is used to render the modal analysis of the finite plate models. Note that the complex moduli of viscoelastic polymer, which considers the Williams-Landel-Ferry model i.e. the temperature dependence of viscosity, are taken into the calculations.

Table 1. Mechanical properties of Aluminum and Polymer

Aluminum

High density rigid constraining layer

The CLD patch with an aluminum constraining layer of 0.5 mm applied on a base plate is considered for two plate models with a thermoviscous core and a polymer core. Both cores have the same thickness (δ) of 0.2 mm. For the thermoviscous core, Rhordsil™ 47 500,000 is considered, and its viscosity is 486 kg/s•m. Equivalent moduli (Young's modulus, density, Poisson's ratio and damping loss factor) of the sandwich plate models with both core materials are calculated. The composite layer of a polymer core displays the frequency/temperature dependent characteristics. These parameters are taken for calculating finite plate models, and the results are shown in Fig. 3. At 20 o C, the CLD patch with a polymer core has higher damping levels throughout the frequency range. The polymer core gives almost twice higher damping than the thermoviscous core for the same thickness applied.

Low density elastic constraining layer

The constraining layer of aluminum is replaced with a HDPE (High-Density PolyEthylene) layer of 0.5 mm. Thicknesses of the base plate (h 2 ) and the core (δ) for both plate models are kept 0.5 mm and 0.2 mm respectively. The mechanical properties of HDPE taken from Ege et al. 3 are E = 850 MPa, ρ = 945 kg/m 3 and η = 15%. The damping loss factor of the infinite thermoviscous composite plate is clearly seen higher than the polymer for computed range of frequency, and this is consistent for finite plate models. Since a density of polymer is higher than that of the thermoviscous fluid, a thickness (δ) of the thermoviscous core is increased to 0.31 mm in order to apply the same total mass of core materials to the composite plate models. Rhodorsil™ 47 1,000,000 oil is used for the computation, and its viscosity is 972 kg/s•m. A thickness (δ) of the polymer core is kept constant as 0.2 mm. As seen in figure 5, modal damping levels of the thermoviscous CLD patch are significantly higher for both infinite and finite plate models.

Conclusion

Rectangular plates partially treated with a CLD patch have been investigated for two different core materials: thermoviscous and viscoelastic. Rhodorsil™ 47 oil and a polymer layer are considered for their damping loss factors. It is shown that the damping provided by the thermoviscous fluid is proportional to its viscositiy (Rhodorsil™ 47, 500,000 and 1,000,000). The thermoviscous fluid can yield higher damping than the polymer layer for the same core thickness applied if the constraining layer is a light-weight elastic material. With such a material, the thermoviscous fluid guarantees significantly higher damping loss factors than the polymer layer at different temperatures if the same total mass is applied. The main advantage of the thermoviscous fluid core compared to the viscoelstic one is that the high damping property of the sandwich is not frequency and temperature dependent.

  of the damping loss factor can be derived:

Figure 1 .

 1 Figure 1. Schematic of CLD (Constrained Layer Damping). (a) Infinite sandwich panel with a core material (thermoviscous or viscoelastic). (b) Equivalent single layer modelling of a simply supported rectangular plate damped with a CLD patch.

Figure 2 .

 2 Figure 2. Equivalent moduli of the Aluminum-Polymer-Aluminum sandwich infinite plate model at different temperatures. (a) Young's modulus, (b) Loss factor.

Figure 3 .

 3 Figure 3. The 1/3 octave band averaged modal loss factor of finite plate models damped with a CLD patch at 20 o c: Aluminum-Rhodorshil™-Aluminum and Aluminuim-Polymer-Aluminuim.

Figure 4 .

 4 Figure 4. The 1/3 octave band averaged modal loss factor of finite plate models damped with a CLD patch at 20 o c: HDPE-Rhodorshil™-Aluminum and HDPE-Polymer-Aluminuim.

Figure 5 .

 5 Figure 5. The octave band averaged modal loss factor of finite plate models damped with a CLD patch at different temperatures: HDPE-Rhodorshil™-Aluminum and HDPE-Polymer-Aluminuim. Rhodorshil™ 47 oil is temperature-independent.

Table 2 . MOVISAND calculation: equivalent parameters for infinite plate models. Thicknesses of Aluminum (h 1 and h 2 ), Rhodorsil™ 47 (δ) and Polymer layers (δ) are 0.5 mm, 0.2 mm and 0.2 mm respectively.

 2

		Sandwich plate
		Aluminum	Aluminum
	Equivalent parameters	Rhodorsil™ 47 500,000	Polymer
		Aluminum	Aluminum
	E (Pa) ρ (kg/m 3 )	1.04E+10 2487.2	E(f)

Table 3 . MOVISAND calculation: equivalent parameters for infinite plate models. Thicknesses of HDPE (h 1 ), Rhodorsil™ 47 (δ) and Polymer layers (δ) are 0.5 mm, 0.2 mm and 0.2 mm respectively.

 3

		Sandwich plate
		HDPE	HDPE
	Equivalent parameters	Rhodorsil™ 47 500,000	Polymer
		Aluminum	Aluminum
	E (Pa) ρ (kg/m 3 ) ν	5.34E+10 1718.4 0.31	E(f) 1806.25 0.337
	η	0.1345	η(f)

Table 4 . MOVISAND calculation: equivalent parameters for infinite plate models. Thicknesses of HDPE (h 1 ) and Rhodorsil™ (δ) are 0.5 mm and 0.31 mm respectively.

 4

	Sandwich plate

Acknowledgements

This work was co-funded by the French government (FUI 12 -Fonds Unique Interministériel) and European Union (FEDER -Fonds européen de développement régional ). It was carried out in the framework of the LabEx CeLyA ("Centre Lyonnais d'Acoustique", ANR-10-LABX-60) and the research project CLIC ("City Lightweight Innovative Cab") labelled by LUTB cluster (Lyon Urban Truck and Bus), in partnership with Renault Trucks, Arcelor-Mittal, ACOEM, CITI Technologies, FEMTO-ST (Univ. de Franche-Comté) and LVA (INSA de Lyon).