Prof Eric Markiewicz
email: eric.markiewicz@univ-valenciennes.fr

G Haugou

J Fabis

G Gary

Contribution to the definition of a partial overlapping plastic strain rates domain for moderate loadings -application to tensile testing on metallic materials

Keywords: Sheet specimens, moderate strain rates, no-direct tensile loading, Hopkinson bars device, overlapping domains, metallic materials

eq equivalent plastic strain (.) g plastic shear strain (.) s eq equivalent plastic stress (MPa) s conventional stress (MPa)

INTRODUCTION

The identification of materials properties is based on the measurement of the evolution of the plastic stress/strain relations up to the fracture of the specimen over a large plastic strain rates range [0.0001;1000] s -1 . This range of strain and rates of strain is commonly observed during crash or impact of transportation frameworks. Consequently, thermoviscoplastic parameters of material laws, supplying FE models, have to be identified for the highest plastic strains through large experimental testing programmes. However, no universal facility allows to overlap this large range of strain and strain rates. The use of various testing facilities, bounded to partial strain rates ranges for a given plastic strain, is then necessary. Quasistatic machines, high-speed hydraulic machines and split Hopkinson bars are typically used as recurrent facilities. Technical limitations of high-speed hydraulic machines concern the exploitation capabilities of the collected raw data for a maximal plastic strain rate depending on the natural frequency of the assembly unit [START_REF] Beguelin | Approche expérimentale du comportement mécanique des polymères en sollicitation rapide[END_REF][START_REF] Fabis | Caractérisation expérimentale de la loi de comportement d'un acier XES en dynamique[END_REF][START_REF] Rusinek | Comportement viscoplastique des tôles en traction et cisaillement[END_REF] and the closed loop control system. Moreover, the split Hopkinson bar technique is available for a limited time duration t where the maximum strain is depending on the length and the speed of the striker [START_REF] Rota | Application de méthodes inverses au dépouillement de l'essai aux barres de Hopkinson[END_REF][START_REF] Gary | Étude expérimentale du comportement dynamique des matériaux[END_REF]. Deconvolution techniques have been recently developed [START_REF] Bussac | An optimisation method for separating and rebuilding onedimensional dispersive waves from multi-point measurements -Application to elastic or viscoelastic bars[END_REF][START_REF] Othman | Extension du champ d'application du système des barres de Hopkinson aux essais à moyennes vitesses de déformation[END_REF] for lowest strain rates but

are not yet commercially available. The paper proposes a tensile assembly unit composed of threaded sleeves for sheet specimens and placed between a split Hopkinson bars apparatus. The interest of the configuration lies in the minimisation of the impedance mismatch between the bars and the tensile assembly unit (see Fig. 3). The tensile loading is obtained indirectly by the impact of a long cylindrical projectile on the input bar in contact with the output bar. The dimensions of the bars and the location of the strain gauges are established with the help of a Lagrangian x/t diagram [START_REF] Nicholas | Tensile testing of materials at high rates of strains[END_REF]. The influence of the fixing ends is evaluated to quantify the impedance mismatch. A tensile tests programme on two metallic alloys is carried out at moderate strain rates and compared to others tests results for the same metallic alloys obtained on a high-speed hydraulic machine [START_REF] Haugou | Iterative experimental/numerical procedure for improvement of dynamic experimental facilities[END_REF] and a SHPB apparatus using hat-specimens [START_REF] Mouro | Dynamic tensile testing of sheet metal[END_REF].

TESTING TECHNIQUE

Description of the no-direct tensile testing configuration

The aim of the ductile materials characterization is to have a maximum duration time of the output bar's loading.

For instance, a typical low carbon steel admits an ultimate engineering plastic strain close to 0.4 before fracture.

Consequently, for a constant strain rate of 200 s -1 , the needed duration time t is closed to 2 ms. To ensure such duration time pulses, a 5 m long projectile and a 10 m long input bar to ensure a correct separation of the pulses without considering launch system limitations would be theoretically required. In order to ensure a maximum time duration of the incident pulse e I (t) for our set-up -required for ductile materials -the length of the projectile is closed to 1.1 m. For the corresponding duration time, thus 420 s -1 , an engineering plastic strain of 0.1 is obtained for a plastic strain rate equal to 250 s -1 which is sufficient for the identification of behaviour laws parameters [START_REF] Langrand | Contribution à la caractérisation numérique et expérimentale d'assemblages structuraux rivetés sous sollicitation dynamique[END_REF].

As represented in Figure 1, a Lagrangian x/t diagram illustrates the three elastic pulses propagation along the bars [START_REF] Nicholas | Tensile testing of materials at high rates of strains[END_REF]. Considering the length of the projectile, the output bar (OB) is 2.7 m long which is necessary to avoid any superposition of the first reflected pulse with the incident one. The initial elastic pulse (bolded line) propagates along the joined bars and through the assembly unit. After being reflected at the free face of the output bar, the initial compressive incident wave becomes a tensile incident wave. The specimens are then loaded and, e I being known, usual complementary elastic pulses are then obtained (e R , e T ). The location of the strain gauges was determined to avoid superposition of the reflected pulse e R at the free ends with the incident pulse e I after its path along the two bars.

The no-direct tensile set-up is composed of two different bars lengths (L 1 and L 2 ) and a cylindrical projectile named Lproj (see Fig. 2).

The two bars are made of 25 mm diameter maraging steel rods. The projectile is made of stainless steel with a 22.4 mm diameter. The Hopkinson bars are supported on 3 point-shapes spaced out 500 mm. The initial pulse is obtained from the impact of a cylindrical projectile launched by a gas gun. The mechanical properties of the set-up is summarized in Table 1. The difference in diameter between the projectile and the bars induces a correction of the measured strain [START_REF] Gary | Some aspects of dynamic testing with wave-guides[END_REF] in order to take into account the lack of energy brought back by the projectile. minimization of the impedance mismatch [START_REF] Mouro | Dynamic tensile testing of sheet metal[END_REF][START_REF] Nguyen | On spurious reflection waves in Hopkinson bar tensile tests using a collar[END_REF][START_REF] Yokoyama | Impact tension testing of sheet metals for automobile structural uses[END_REF][START_REF] Leblanc | Dynamic tensile testing of sheet material using the Split Hopkinson bar technique[END_REF].

The speed of the striker bar has been limited to 10 m.s -1 to avoid strain gauges deterioration. Moreover, the lowest speed is closed to 3 m.s -1 considering the limitations of the launching system. The gauge length for sheets specimen is chosen close to 11.3 mm. It allows to reach the moderate strain rates domain [200;700] s -1 for tensile tests. The full section is equal to four gauge sections of a specimen in order to obtain a significant amplitude of the transmitted signal and to avoid its amplification. The cross section of the sleeves is slightly decreased (2%) by the four grooves tooled on the cylindrical sleeves surfaces for the positioning of the specimens (see Fig. 3). The gauge lengths of the specimens lies between 6.8 and 25 mm to have access to a theoretical strain rates range from 120 to 1400 s -1 . For the lowest rate of strain, plastic strains close to 0.05 are measured. In the same way, the current set-up allows to reach plastic strains up to fracture for higher strain rates according to [START_REF] Beguelin | Approche expérimentale du comportement mécanique des polymères en sollicitation rapide[END_REF].

e e t

= ˙◊ [START_REF] Beguelin | Approche expérimentale du comportement mécanique des polymères en sollicitation rapide[END_REF] Figure 3 illustrates the progress of the 4-glued sheet specimens on the sleeves. An epoxy adhesive is laid on the adhesive areas of the specimens. Then, a pre-load is done at the free faces of the unit assembly to avoid any gap between the sleeves after the polymerisation of the epoxy adhesive (12 hours at room temperature).

Figure 4 shows the typical records on each bar using the no-direct tensile configuration on sheet specimens.

XES steel and 2024 T3 aluminium alloy material characterization

Tensile specimens have been machined by wire electroerosion from a cold-rolled sheet of XES steel (Wt% C = 0.05-Cr < 0.01-Cu = 0.013-Mn = 0.018-Al = 0.064-Ni = 0.02) and 2024T3 aluminium alloy (Wt% Al = 93.5-Zn = 0.22-Mg = 1.4-Cu = 4.2-Mn = 0.48-Fe = 0.17-Si = 0.07). The XES steel specimens were machined from a 1.17 mm thick plate and oriented in the rolling direction. The 2024 T3 aluminium alloy specimens were machined from a 1 mm thick plate and oriented at 90∞ regarding the rolling direction. Conventional quasi-static tests (QS) on the same geometry have been performed using an Instron universal testing machine (Model 1195) at a constant cross speed of 5 mm.min -1 . Oscillations are observed for the highest plastic strain rate close to 420 s -1 and perturbs the stress vs. time calculation (see Fig. 5). The oscillations The length of the strain gauge allows a 145 kHz bandwidth whereas a typical spectral signature of an impact using bars is closed to 50 kHz. The length of the cylindrical bars, the location of the strain gauges and the gauge dimensions of the specimen have been determined using classical relations [START_REF] Gary | Some aspects of dynamic testing with wave-guides[END_REF][START_REF] Zhao | Analyse de l'essai aux barres d'Hopkinson -Application à la mesure du comportement dynamique des matériaux[END_REF].

TENSILE TESTS RESULTS

Material and specimens preparation

Many modifications for the fixation of the specimens onto the bars are generally provided to ensure the holding of both undeformed parts on the specimen. For this, an adaptive assembly must be designed to produce a solicitation as pure as possible. The holding of the specimen is usually made with pin attachments (for plate specimens) or threads (for axi-symetric specimens). Those modifications may induce an impedance mismatch because of cross section variations in the bars [START_REF] Nicholas | Tensile testing of materials at high rates of strains[END_REF][START_REF] Chocron-Benlolou | Dynamic tensile testing of aramid and polyethylene fiber composites[END_REF][START_REF] Nguyen | On spurious reflection waves in Hopkinson bar tensile tests using a collar[END_REF][START_REF] Yuan | Tensile response of ductile atitanium at moderately high strain rates[END_REF] or in the assembly [START_REF] Kammerer | Modélisation du comportement plan d'un n verre E/polyester sous faibles et fortes vitesses de déformation -Application au cas de l'impact[END_REF][START_REF] Yokoyama | Impact tension testing of sheet metals for automobile structural uses[END_REF][START_REF] Radford | The mechanical and constitutive behaviour of as-received and irradiated Zr-2.5Nb pressure tube material under high rates of tensile strain[END_REF]. Any section variation must be taken into account for the waves analysis [START_REF] Chocron-Benlolou | Dynamic tensile testing of aramid and polyethylene fiber composites[END_REF]. For the proposed assembly unit, the two sleeves have been machined from the same steel alloy as the cylindrical bars by wire electroerosion. The choice of the principle is based on the period corresponds to the run inside the sleeves of a partial reflected pulse of the incident wave. These interfaces may be responsible of the various reflections in the system [START_REF] Haugou | Nodirect tensile testing on Hopkinson bars apparatus : assessment of an original configuration by means of 4-glued sheet specimens[END_REF]. So, the comparison is made considering filtered data (low frequency pass filtering) for impact tests using Hopkinson bars [START_REF] Gary | Essai de cisaillement plan appliqué à des tôles minces[END_REF]. The input/output forces balance has been previously cheeked [START_REF] Haugou | Nodirect tensile testing on Hopkinson bars apparatus : assessment of an original configuration by means of 4-glued sheet specimens[END_REF]. For further comparisons between the various dynamic facilities, the tests made with the proposed configuration on the split Hopkinson bars apparatus are called "No-direct SHB".

Comparison for different testing devices

Two others testing facilities were used to complete the considered range in terms of strain and rates of strain. This approach is generally led to encompass a large strain rate domain [START_REF] Lindolhm | Effect of strain rate on yield strength, tensile strength and elongation of three aluminum alloys[END_REF][START_REF] Nicholas | Dynamic tensile testing of structural materials using a split Hopkinson bar apparatus[END_REF][START_REF] Townsend | Material Tests -Tensile properties of aluminium alloys 7010T651 and AU4G over a range of strain rates[END_REF]. A hydraulic Schenck machine uses a 50 kN capacity linear actuator (see Fig. 6).

Input bar Output bar

No-direct tensile raw data on 4-glued XES steel sheet specimens using a For the present tests, the servo system was programmed with a ramp function to provide a constant relative jack speed available up to 4 m.s -1 . Previous works have permitted to through out the natural frequency of the assembly. This kind of device is naturally subjected to mechanical vibrations depending on the weight of the added mass and the stiffness of the piezo-electric load cell [START_REF] Beguelin | Approche expérimentale du comportement mécanique des polymères en sollicitation rapide[END_REF][START_REF] Fabis | Caractérisation expérimentale de la loi de comportement d'un acier XES en dynamique[END_REF][START_REF] Rusinek | Comportement viscoplastique des tôles en traction et cisaillement[END_REF][START_REF] Haugou | Iterative experimental/numerical procedure for improvement of dynamic experimental facilities[END_REF]. The natural frequency of the testing set-up is closed to 8 kHz which allows the access to theoretical plastic strain rates close to 200 s -1 according to material properties of the considered alloy [START_REF] Haugou | Iterative experimental/numerical procedure for improvement of dynamic experimental facilities[END_REF][START_REF] Haugou | Moyens d'essais et de caractérisation de lois de comportement matérielles en dynamique moyennes vitesses[END_REF]. The mechanical oscillations correspond to the natural frequency of the set-up due to the dynamic loading of the specimen and remains when fracture occurs [START_REF] Lindolhm | Effect of strain rate on yield strength, tensile strength and elongation of three aluminum alloys[END_REF]. For the tests on the hydraulic testing machine, an optical extensometer (Zimmer XH 200) was employed to measure higher plastic strains [START_REF] Haugou | Moyens d'essais et de caractérisation de lois de comportement matérielles en dynamique moyennes vitesses[END_REF]. The gauge length L1 of the specimen was taken equal to 6.8 mm for a higher strain rates range [START_REF] Nicholas | On the determination of the mechanical properties of materials at high shear-strain rates[END_REF][START_REF] Hamouda | Testing of composites materials at high rates of strain : advances and challenges[END_REF].

The second facility is composed of a SHPB apparatus described in details by Mouro [START_REF] Mouro | Dynamic tensile testing of sheet metal[END_REF] (see Fig. 7). In that configuration, a hat-specimen -machined by wire electroerosion and folded by a quasi static machine -is placed into a specific assembly unit between two cylindrical maraging bars (see Fig. 8).

The behaviour laws obtained using the high-speed hydraulic machine are named "H" and presented in details by Haugou [START_REF] Haugou | Moyens d'essais et de caractérisation de lois de comportement matérielles en dynamique moyennes vitesses[END_REF]. In the same way, the tests results obtained with the second split Hopkinson pressure bars are named "SHPB-LMS". The results are plotted for the different testing devices and measurement systems (see Figs. 9 and 10). The comparison is done for the XES low carbon steel at average strain rates of 230 and 520 s -1 . In the same way, the comparison of the true relations is done for the 2024 T3 aluminium alloy and the same testing devices. It is verified that the considered alloy is not strain rate sensitive from quasi-static rates of strain to 500 s -1 (see Fig. 10).

OVERLAPPING DOMAIN -DISCUSSION

The comparison is made with results of previous tests programmes carried out for the 2024 T3 aluminium alloy [START_REF] Langrand | Contribution à la caractérisation numérique et expérimentale d'assemblages structuraux rivetés sous sollicitation dynamique[END_REF][START_REF] Nicholas | On the determination of the mechanical properties of materials at high shear-strain rates[END_REF]. Nicholas and Lawson [START_REF] Nicholas | On the determination of the mechanical properties of materials at high shear-strain rates[END_REF] have performed shear tests from 10 -4 to 4500 s -1 using three different testing devices (see Fig. 11). A quasi-static machine has covered the range of strain rates [10 -4 ;10 -2 ] s -1 . A pneumatic machine has been used for the next domain from 10 -2 to 100 s -1 . Stress/strain relations from 500 to 4500 s -1 have been obtained with a torsional Hopkinson bars device. Plastic shear stress have been converted according to classical relations (2) and (3) based on the von Mises flow rules [START_REF] François | Comportement mécanique des matériaux[END_REF]. Langrand [START_REF] Langrand | Contribution à la caractérisation numérique et expérimentale d'assemblages structuraux rivetés sous sollicitation dynamique[END_REF] has made quasi-static and dynamic tensile tests from 0.003 to 55 s -1 in order to identify the parameters of Johnson-Cook model for the 2024 T3 aluminium in order to supply FE riveted aeronautic frameworks models. The present comparison shows a good agreement with results obtained with the no-direct tensile configuration. The present testing process is then validated.

s eq = s = t • 3 [2]
e eq = e = g

For the XES steel, previous tests programmes have been prospected under double shear [START_REF] Gary | Essai de cisaillement plan appliqué à des tôles minces[END_REF], compression [START_REF] Gary | Étude expérimentale du comportement dynamique des matériaux[END_REF] and tensile [START_REF] Mouro | Dynamic tensile testing of sheet metal[END_REF] tests on sheet specimens within split Hopkinson bars devices. The domain of plastic strain rates was about 100 to 900 s -1 . Materials are revealed to be strain rate sensitive for these various tests programs and the considered steel alloy, particularly for strain rates up to 100 s -1 (see Fig. 12). However, the use of different testing facilities raises the problem of the interpolation of the mechanical properties identified close to the technological limits of each device, such as oscillations on high-speed hydraulic machines and limited time duration on Hopkinson bars.

On the basis of the present tests programme, an overlapping domain is established for moderate rates of strain using plate specimens made of the same metallic alloys. The classical measurements devices have put into evidence a discontinuous plastic strain domain for moderate strain rates. This is especially clear for strain rate below 500 and 700 s -1 respectively for the 2024 T3 aluminium alloy and XES steel (see Figs. 13 and14). So, for FE modelling of frameworks submitted to large plastic strains and plastic strain rates, an unknown domain remains to correctly identify parameters of viscoplastic laws.

CONCLUSIONS

The paper presents a no-direct tensile testing configuration based on glued sheet specimens on cylindrical sleeves. This assembly unit -containing 4-glued sheet specimens on its external surfaces without any others attachment tools -is threaded on a split Hopkinson bars device. Nevertheless, a problem of oscillations is still to be solved to allow classical SHPB processing of the raw signals.

Moreover, the gluing progress have to be controlled to avoid detachment for high strength materials (i.e. composites, HSS). This involves particular care for the design of the specimens and the related strain rates ranges. The comparison of results obtained using the proposed no-direct tensile solicitation using classical Hopkinson bars apparatus with results from others techniques is shown for a common strain rates range from 200 up to 400 s -1 . A high-speed hydraulic machine and a SHPB apparatus using hat-specimens are used for this purpose in order to have access to the domains of plastic rates of strain [0.002; 500] s -1 and [120; 700] s -1 , respectively. A good agreement is observed for the engineering stress/strain relations provided by the different tests facilities. A restricted overlapping domain has been underlined for the moderate plastic strain rates [120; 500] s -1 and plastic strains up to 0.1 for both metallic alloys (2024 T3 aluminium alloy and XES low carbon steel). This constitutes a limitation for FE analysis of crashworthiness behaviour in particular case where fracture is observed.

Then, direct measurements for the analysis of plastic strains using no-contact measurement devices are considered as optical extensometer techniques, high-speed camera for correlation analysis. In these conditions, a direct tensile configuration on SHB is required using the proposed assembly unit.

1 Figure 2

 12 Figure 2 Schematic diagram of the no-direct tensile loading.

Figure 1

 1 Figure 1 Lagrangian x/t diagram of no-direct tensile configuration on Hopkinson bars.

Figure 3

 3 Figure 3 Progress of the 4-sheet metallic specimens on the cylindrical parts.

Figure 4 1 Figure 5

 415 Figure 4 Raw data under no-direct tensile tests on Hopkinson bars device.

Figure 6

 6 Figure 6 High-speed hydraulic testing machine (Onera-Lille Centre).

Figure 7

 7 Figure 7 Special supports on classical SHPB (LMS-Polytechnique).

  True strain/stress relations for XES steel tensile sheet at an average strain rate of 170 s -160 s -1 No-direct SHB -200 s -1 H -160 s -1True strain/stress relations for XES steel tensile sheet at an average strain rate of 520 s -600 s -1 No-direct SHB -440 s -1 H -580 s -1

Figure 8

 8 Figure 8 Geometry and dimensions of hat specimens (mm).

Figure 9

 9 Figure 9 True relations of XES steel from different testing devices.

Figure 10 Figure 11

 1011 Figure 10 True relations of 2024 T3 aluminium alloy from different devices.

Figure 12

 12 Figure 12 Comparison of the ultimate plastic stress -XES steel.

Figure 13

 13 Figure 13 Overlapping domain for different testing facilities -XES steel.

Figure 14

 14 Figure 14 Overlapping domain for different testing facilities -2024 T3 aluminium alloy.

Table 1 Specific dimensions of Hopkinson bars apparatus

 1 Measurement group) bonded on the Hopkinson bars with a cyanoacrylate adhesive. The output signals are obtained from a full-bridge into a Vishay 2310 amplifier. A 8-channel transient recorder (Nicolet MultiPro -sampling rate 1 MHz) is used to storage tests data (elastic pulses and signals from photo-diodes) afterwards.

		L (m)	∆ (mm)	Section Density r Velocity (mm 2 ) (kg.m -3 ) C (m.s -1 )
	input bar	3.6	25	491	8198	4627
	output bar	2.7	25	491	8171	4599
	projectile	1.1	22.4	394	7702	5288
	The three elastic pulses (e I (t), e T (t) and e R (t)) are
	recorded with a set of 3.2 mm strain gauges CEA-06-
	125-UN-350 (ref.

ACKNOWLEDGMENTS

This research was carried out with the financial support of the CNRS, ONERA-Lille Centre and the Conseil Régional du Nord/Pas de Calais. The present authors wish to thank especially A. Deudon, P. Mahelle and R. Barre for there contributions to the tests performed at Onera-Lille Centre and LMS-Polytechnique, respectively.


Yongliang Feng
email: yongliang.feng@univ-amu.fr

Pierre Boivin

Jérôme Jacob

Pierre Sagaut
email: pierre.sagaut@univ-amu.fr

Hybrid recursive regularized lattice Boltzmann simulation of humid air with application to meteorological flows

Keywords: numbers: 47.55.t, 47.11.j,47.45.-n Keywords: lattice Boltzmann method, atmospheric convection, cloud dynamics, shallow cumulus convection

An extended version of the hybrid recursive regularized Lattice-Boltzmann model which incorporates external force is developed to simulate humid air flows with phase change mechanisms under the Boussinesq approximation. Mass and momentum conservation equations are solved by a regularized lattice Boltzmann approach well suited for high Reynolds number flows, whereas the energy and humidity related equations are solved by a finite volume approach. Two options are investigated to account for cloud formation in atmospheric flow simulations. The first option considers a single conservation equation for total water and an appropriate invariant variable of temperature. In the other approach, liquid and vapor are considered via two separated equations, and phase transition is accounted for via a relaxation procedure. The obtained models are then systematically validated on four well-established benchmark problems including a double diffusive Rayleigh Bénard convection of humid air, 2D and 3D thermal moist rising bubble under convective atmospheric environment as well as a shallow cumulus convection in framework of large-eddy simulation.

I. INTRODUCTION

Due to its central role in weather forecast [1] for being a key element in cloud formation and rain, atmospheric moisture has a significant impact on air quality and consequently on public health [2]. Moist atmospheric dynamics is composed of motions occurring on a wide range of temporal and spatial scales, from 0.1-1km for the convective macro-structures (e.g. clouds), down to the mm-sized water drop governing the atmospheric moisture through liquid-vapor equilibrium. Although all these mechanisms can be captured using the fully compressible Navier-Stokes equations incorporating energy and water conservation laws coupled with a condensation model, the number of degrees of freedom is expected to be reduced by filtering out specific modes of motion. To tackle this multiscale issue, several mathematical approximations based on the compressible Navier-Stokes equations have been developed e.g. Boussinesq equations [3], pseudo-incompressible equations [4] and the anelastic equations [5]. The Boussinesq approximation is widely adopted in numerical weather prediction tools in which the continuity equation is replaced by an incompressibility condition.

Numerical modelling of atmospheric flows under above mathematical approximations is still challenging due to the large-scale motions and circulations and coupling among multi-physical variables. Accurately modelling the interplay between fluid dynamics and the thermodynamics related to reversible and irreversible moist con-densation or evaporation processes is then key in atmospheric flow modelling. There are two popular numerical treatments of moist microphysics in the Navier-Stokes equations framework [6]. In the first approach, often referred to as the conservative variables (liquid potential temperature and total water mass fraction) approach, the equations of motion are defined using appropriate conservative variables such that terms resulting from phase change are eliminated in the governing equations [START_REF] Cuijpers | Large eddy simulation of trade wind cumulus clouds[END_REF]. This condensation model is mainly employed for shallow convection under Boussinesq approximation. In the second approach, which is more commonly adopted, the equations are again defined using the sensible potential temperature, and separate transport equations for mass fractions of liquid water and vapor water. In the second step a saturation adjustment procedure is performed that correctly allocates the water into its two phases based on the Clausius-Clapeyron formula [START_REF] Duarte | [END_REF]. Balancing accuracy, efficiency and the physical realisability is problem specific and rarely straightforward. Consequently, there is a continuing quest for improved hydrodynamic core incorporated with condensation models applied in climate, weather and chemistry-transport models [9].

The lattice Boltzmann method (LBM) is an alternative approach to simulate fluid flows based on the Boltzmann equation. Due to its advantages for massively parallel computing as well as its suitability to handle very complex geometries, the LBM has been quickly extended to large scale and multi-scales applications. These applications cover, among others, aerodynamic studies on fullscale vehicles, turbulent flows in urban areas [10], aerodynamic predictions on airfoils [11]. Furthermore, the LBM also has been extended to complex multi-physics phenomena through the additional distribution function approach proposes additional conservation law to solve multi-physics fields [START_REF] Luo | Encyclopedia of Aerospace Engineering[END_REF][START_REF] Feng | [END_REF]. Extending the application scope of lattice Boltzmann method to large-scale meteorological flows by incorporating cloud dynamics has received some interest [14,15], but it is worth noting that most existing LBM approaches for multiphase flows have been assessed considering laminar flows with resolved liquid/gas interfaces (see [16] for a survey) and are therefore not well suited for atmospheric flow simulations.

When the lattice Boltzmann method is applied in complex large-scale geophysical flows [17], such as moist atmospheric convection and condensation, the special limitation on the range of pressure, physical sound speed, relaxation procedure for phase change and associated time and space resolution should be carefully considered, e.g. the numerical instability induced by the use of large time steps, coarse grid and thermal stratification. Besides, the accuracy and efficiency of the condensation models are also highly related with temporal and spatial resolution as well as the numerical schemes [18]. For instance, the numerical stability of the lattice Boltzmann collision model has long been identified among the key issues [19]. In order to overcome the insufficient stability observed in the Bhatnagar-Gross-Krook (BGK) model [20], several improved collision models with enhanced stability have been proposed, e.g. multi-relaxation-time models [21,22], entropic LB models [23,24], regularized BGK models [25,26], cumulant/cascaded LB models [27,28]. The MRT models were successful in moderately stabilizing the LB method but still remain challenged by high Reynolds numbers [29,30]. In addition, their high number of free parameters have motivated the development of alternative models. For example, entropic LB models were developed for turbulent flows through reestablishing H theorem, which was lost during the discretization [31]. The cumulant and cascaded LB models were developed to avoid the deficits of the MRT models, e.g. violation of Galilean invariance or an insufficient level of hyperviscosity [27,28]. Moreover, the recursive regularized LB models were designed by considering that the relaxation parameters of the ghosts are chosen to annihilate the higher-order non-equilibrium moments in the post-collision state. The RLB models were successful in simulation of high Reynolds flows [32], even exhibited superconvergence properties for the acoustic problems compared with the MRT models in inviscid limit [33]. In particular, regularized BGK model with hybrid recursive procedure presented both fairly good robustness and accuracy in large scale simulations on urban flows [10,34].

In this paper, we focus on developing a hybrid recursive regularized lattice Botlzmann method with condensation scheme for cloud dynamics with application to turbulent meteorological flows. The paper is organized as follows. Section II presents governing equations along with the Boussinesq approximation. Section III reviews the lattice Boltzmann method and the approach used to implement the condensation models and to solve humidity and temperature equations, respectively, and the numerical procedure of the proposed method is also presented. Section IV investigates and discusses simulations on a two dimensional double diffusive Rayleigh Bénard convection, as well as a 2D and 3D thermal moist rising bubble under convective atmospheric environment. In particular, the proposed LB approach for humid air and cloud dynamics is assessed considering a shallow cumulus convection in framework of large-eddy simulation. Finally, Section V draws conclusions and perspectives.

II. MACROSCOPIC PHYSICAL MODEL

A. Governing equations under the Boussinesq approximation

The atmosphere is assumed to be a mixture of dry air, water vapor, and liquid water, with respective mass fractions q d , q v and q l (q v being often referred to as specific humidity). The Boussinesq approximation, widely used for the simulation of atmospheric flows, consists in neglecting density differences everywhere but the gravity terms. Under this assumption, the mass and momentum conservation equations read [6]

∂ u α ∂x α = 0, (1)
∂ u α ∂t + u β ∂ u α ∂x β = - 1 ρ 0 ∂ p ′ ∂x α + ∂ ∂x β ν( ∂ u α ∂x β + ∂ u β ∂x α ) + A g,z (2)
with use made of Einstein's summation notation. u α represents the components of the velocity vector and x α represents the components of the position vector (x, y, z); ν is the kinematic viscosity. The hydrodynamic pressure p ′ = p -p 0 (z) is related to pressure p departing from the height-dependent reference state pressure p 0 (z) and ρ 0 the constant reference state density. A g,z is the buoyancy term due to the gravitational acceleration g, to be detailed in Eq. (9).

Next, let us define the potential temperature θ as a function of temperature T following [35]

Π = T θ = p 0 (z) p 0 R d /cp , (3)
where R d is the air gas constant (per mass unit), c p the average mass heat capacity, assumed to be constant and equal to that of dry air.

Π = 1 - g.z c p .θ 0 (4)
is the Exner function [35], θ 0 = T 0 is the reference temperature and p 0 is the (constant) reference temperature at ground level, not to be confused with p 0 (z). With the above definitions, the energy conservation equation may now be expressed in terms of the potential temperature θ as [6].

∂ θ ∂t + u α ∂ θ ∂x α = ∂ ∂x α D θ ∂ θ ∂x α + L v θ c p T Q, (5)
where D θ is the heat diffusion coefficient and Q is the mass transfer between the liquid and gas phases associated to L v , the mass latent heat of water. Associated conservation equations for the liquid and vapor mass fractions (q l and q v , respectively, the dry air mass fraction being deduced as q d = 1 -q v -q l ) are [6]

∂ q v ∂t + u α ∂ q v ∂x α = ∂ ∂x α D q ∂ q v ∂x α -Q (6a) ∂ q l ∂t + u α ∂ q l ∂x α = ∂ ∂x α D q ∂ q l ∂x α + Q (6b)
with D q being the water fraction diffusion coefficient.

As to derive the buoyancy acceleration in Eq. ( 2), let us now define the virtual temperature θ v as [36]

θ v = θ 1 -1 - 1 ǫ q v -q l , (7)
with ǫ the molecular mass ratio of dry air to that of water:

ǫ = W v /W d = R d /R v . (8)
The buoyancy acceleration required in the momentum conservation equation (2) then reads:

A g = g θ v,H θ v -θ v,H (9)
where θ v,H is virtual temperature of reference base state.

B. Phase transition modelling

In the present work, we assume that the rate of phase transition is infinitely fast, or equivalently, that the liquid and gas phases are in thermo-chemical equilibrium everywhere at all times. Under this assumption, saturation properties provide additional relations between q v and q l . The saturation specific humidity is

q sat v = ǫp sat p 0 (z) -(1 -ǫ)p sat , (10)
in which the saturation pressure is approximated following [37] p sat (T ) = 610.78 exp 17.269

T -273.16 T -35.86 , ( 11
)
where T is in Kelvin and p sat in Pa.

Under the infinitely fast relaxation approximation, the source term Q in Eq. ( 6) can then be computed from Q = -q l /δ t if q v < q sat v and q l < ∆q v ∆q v /δ t otherwise [START_REF] Luo | Encyclopedia of Aerospace Engineering[END_REF] where

∆q v = c p R d (θ 0 Π) 2 (q v -q sat v ) c p R d (θ 0 Π) 2 + ǫq sat v L 2 v ( θ0 θ -θ0 θ R d θ0Π ǫLv ) , (13)
with δ t being time step following [6,38] . Note that in this framework, Eq. ( 13) can readily be adapted to address finite rate phase transition, including for instance the liquid droplets size. In the case where only infinitely fast relaxation towards thermochemical equilibrium is of interest, the saturation relation Eq. ( 10) can be directly included into the conservation equations, thereby reducing the number of partial differential equations, and resulting in a condensation model based on two invariant variables (mass fraction of total water q t and liquid water potential temperature θ l ). Eq. (6a) and Eq. (6b) can be combined to cancel the source term Q in a

q t = q v + q l ( 14
)
conservation equation, becoming

∂ q t ∂t + u α ∂ q t ∂x α = ∂ ∂x α D q ∂ q t ∂x α . ( 15
)
It then becomes more convenient to rewrite the energy conservation equation using an approximate expression of the liquid water potential temperature [36]

θ l ≈ θ - L v c p Π q l , (16)
to read

∂ θ l ∂t + u α ∂ θ l ∂x α = ∂ ∂x α D θ ∂ θ l ∂x α , (17)
as to cancel the source term. Where required, q l can then be computed from q t as q l = max(0, q t -q sat v ).

In that framework, however, q sat v cannot be evaluated directly from Eq. ( 10), as T is not available. It is instead approximated from [6]

q sat v = q * 1 + ηq t 1 + ηq * , η = ǫ L 2 v c p R d T 2 l , (19)
where q * corresponds to Eq. ( 10), evaluated at liquid temperature T l = Π.θ l , rather than Π.θ. Compared to the former approach, only two equations of invariant variables (θ l , q t ) need to be solved in the reduced model, and the relaxation step is replaced by a recovery procedure, which is performed using the approximate expression of saturation humidity following Eq. ( 19). It is worth noting that the set of invariant variables (θ l , q t ) could often be simply derived in the shallow convection.

In the following, the model based on set of governing equations about invariant variables (q t , θ l ) is referred to as 1eq. model, whereas the model based on the variables (q l , q v , θ) is referred to as the 2eq. model, as summarized in Table I. Eq. ( 5) Eq. ( 17) Water conservation Eq. (6a) and Eq. (6b) Eq .( 15)

III. HYBRID RECURSIVE REGULARIZED LATTICE BOLTZMANN MODEL

This section presents the hybrid recursive regularized Lattice Boltzmann model (HRR-LBM) proposed to solve the governing equations. The idea is to solve the mass and momentum equations (1,2) following a LB approach [39], whereas the energy and water conservation(s) equation(s) are solved by a classical finite volume (FV) approach. The present developments have been performed in the ProLB solver.

A. Lattice Boltzmann solver

Reminders about LBM

Lattice Boltzmann methods aim at solving the lattice Boltzmann equation through space, time, and velocity discretization [39,40]. Space and time are classically discretized on a cartesian grid, whereas speeds are discretized on a so-called DnQm lattice (n dimensions and m discrete velocities c i ).

The flow problem is then solved for f i (x, t), the density distribution of particles with velocity c i at (x, t), which can be obtained at time t + δ t through the so-called BGK collision model [40]

f i (x + c i δ t , t + δ t ) = f i (x, t) - 1 τ [f i (x, t) -f eq i (x, t)] + (1 - 1 2τ )δ t F i (x, t), (20)
a succession of a streaming and a collision step [41,[START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF] which is equivalent to a Strang-splitting-based time integration method [START_REF] Dellar | [END_REF]. The characteristic collision time τ is related to the kinematic viscosity ν through ν = c 2 s (τ -1 2 )δ t , and c s is the lattice sound speed, inherent of the lattice DnQm choice [START_REF] Krüger | The Lattice Boltzmann Method: Principles and Practice[END_REF]. Density and momentum can then be recovered as

ρ(x, t) = m i=1 f i (x, t), ρu(x, t) = m i=1 c i f i (x, t) + δ t 2 m i=1 c i F i (x, t), (21)
where F i (x, t) is an external force [44].

The recursive regularized Lattice-Boltzmann solver

Considering that the third-order equilibrium distribution function and recursive regularized collision model exhibit considerably better performance on numerical stability and accuracy compared to the original LBGK model [17,[45][46][47][48][49], a recursive regularized collision model associated with the third-order equilibrium distribution function is employed as the collision model. The lattice Boltzmann equation Eq. ( 20) is rewritten as

f i (x + c i δ t , t + δ t ) = f eq i (x, t) + (1 - 1 τ )f neq i (x, t) + δ t 2 F i (22) where f neq i = f i -f eq i + δ t F i /
f eq i = w i ρ + c iα ρu α c 2 s + a (2) 0,αβ H (2) iαβ 2c 4 s + a (3) 0,αβγ H (3) iαβγ 6c 6 s , a (2)
0,αβ = ρu α u β , H (2)
iαβ = c iα c iβ -c 2 s δ αβ , a (3) 0,αβγ = ρu α u β u γ , H (3) iαβγ = c iα c iβ c iγ -c 2 s [cδ] αβγ . ( 23
)
where H i,αβ and H i,αβγ correspond to the second and third order Hermite polynomials. The corresponding offequilibrium distribution function f neq i is also expended in Hermite polynomial according to Chapman-Enskog technique as

f neq i = a (2) 1,αβ H (2) iαβ 2c 4 s + a (3) 1,αβγ H (3) iαβγ 6c 6 s , (24)
and reconstructed by a Recursive Regularized BGK Lattice-Boltzmann Method (RR-BGK-LBM) with forcing terms as

a (2) 1,αβ = i H (2) i,αβ (f i -f eq i + δ t 2 F i ), (25) a (3) 1
,αβγ = u α a (2) 1,βγ + u β a (2) 1,γα + u γ a (2) 1,αβ . (26)
Lastly, to account for the buoyancy force Eq. ( 9) in the momentum equation ( 2), the external force term F i (x, t) in Eq. ( 22) is prescribed [44,51] as

F i = ρw i A α c iα c 2 s + (u α A β + u β A α )(c iα c iβ -c 2 s δ αβ ) 2c 4 s . ( 27
)
Through the Chapman-Enskog technique [39], it can be shown that the lattice Boltzmann BGK equation (22), appended with the expression for the equilibrium (23) and non-equilibrium (24) functions, and the forcing term (27) recovers the mass (1) and momentum (2) conservation equations in the low-Mach regime.

In the simulations, physical and lattice units for length and time are related through a reference length scale L 0 , a physical reference sound speed c sp for space and time.

Accordingly, the physical variables can be converted from the quantities in lattice units (l subscript) as

δ x = L 0 /N n , δ t = δ x . c s c sp , t = N t .δ t x = N i .δ x , u = u l . δ x δ t , ν = ν l δ 2 x δ t (28)
The physical reference sound speed c sp is kept as a free parameter. This allows to freely accelerate the convergence rate, by varying the time step in the same manner of Courant-Friedrichs-Lewy (CFL) number for conventional Navier-Stokes solvers.

B. Finite Volume solver

Numerically, the multiple distribution function model for scalar transport equation is not optimal from the computational standpoint, even though this numerical inefficiency can be improved somewhat by using some redundant degree of freedom in LB models in diffusion dominated condition [52]. Besides, it is convenient to use an upwind biased scheme to preserve the numerical stability in meteorological flows in which the strong advection effect and stiff terms due to phase transition commonly exist. Energy and water conservation equations (5, 6a, 6b, 15 and 17) of both the 2eq. and 1eq. models (see Tab. I) are solved with a finite volume method. They may all be written under the non-conservative generic form

∂φ ∂t + u α ∂φ ∂x α = ∂ ∂x α (D ∂φ ∂x α ) + Ψ, ( 29
)
where φ is the scalar (e.g. q l , θ l ,...), D is the diffusion coefficient, provided via specification of Prandtl numbers

Pr θ = ν/D θ , Pr q = ν/D q , (30)
for the heat and molecular diffusion, respectively. Ψ is a source term (e.g. related to the phase transfer rate Q), and velocity u α required in the convective term is obtained from the LB solver via Eq. ( 21).

In the present hybrid LB approach, the scalars transport equations e.g. temperature, humidity are solved by a vertex based finite volume (FV) method in the same temporal-spatial discrete coordinates of LB method. The coupling between LB and FV is detailed in Fig. 1. The first-order explicit Euler scheme is adopted as temporal discretization, which is given as

φ n+1 = φ n + δ t RHS C (φ n ) + RHS D (φ n ) + Ψ n (31)
where RHS C , RHS D represents convection term and diffusion term, respectively. The convective flux is constructed using MUSCL scheme [53] with the van Albada limiter [START_REF] Hirsch | Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics[END_REF]. The classical second-order central difference scheme is adopted for the diffusion term.

C. Implementation of boundary conditions

The boundary conditions of distribution function are implemented by a finite difference reconstruction approach along with a regularization procedure [25,[START_REF] Latt | [END_REF].

1. First, the macroscopic velocities u on the boundary nodes are prescribed in the case of a Dirichlet boundary condition or extrapolated in the case of a Neumann boundary condition.

2. Next, the shear stress tensor is computed using the velocity gradients on boundary nodes, which are computed on these nodes using a first order biased scheme.

3. Then, the density distribution function is reconstructed by Eq. ( 23) and off-equilibrium moment as follows

a (2) 1,αβ ≈ -ρτ c 2 s ∂ u α ∂x β + ∂ u β ∂x α , (32) a (3)
1,αβγ = u α a

(2)

1,βγ + u β a (2) 1,γα + u γ a (2) 1,αβ , (33)
It is worth noting that the density values on the bottom surface node are set to constant reference density and densities on other boundary are extrapolated from inside nodes.

4. For humidity and temperature fields which are solved by the finite volume approach, the macroscopic values θ or θ l and q t or (q v , q l ) are simply prescribed or extrapolated on all boundary nodes.

D. Summary of the algorithm

The algorithm consists of a LB and a FV solver running simultaneously. The only data exchanges between the two solvers are :

• the velocity u obtained from Eq. ( 21) in the LB solver and sent to the FV solver to assess convective fluxes,

• θ v , obtained from the energy equation in the FV solver, and sent to the LB solver to update the forcing term.

The complete algorithm is detailed in Fig. 1 in the case of the 2eq. model, the 1eq model merely being a reduction of the former. The 2D (resp. 3D) computational examples presented in the next section were obtained within the D2Q9 (resp. D3Q19) framework in ProLB [START_REF]ProLB: high-fidelity CFD in exceptional turnaround times[END_REF].

Temp erature

Vap or Liquid u, θ, q v and q l at (n-1) th time Compute θ * by Eq. ( 5) Compute q * v by Eq. (6a)

Compute

q * l by Eq. (6b) Calculate T → p sat (T ) → q sat v by Eqs. (10,11) Calculate Q and ∆q v by Eqs. [START_REF] Luo | Encyclopedia of Aerospace Engineering[END_REF][START_REF] Feng | [END_REF] θ v (7) Update q v ,q l and θ θ q v q l FV advection step

IV. RESULTS AND DISCUSSION

A. Double diffusive Rayleigh Bénard convection

Albeit not technically a moist convection condensation problem, a two dimensional Rayleigh-Bénard convection is studied in a square domain with height of H = 1 m to assess the proposed LB method on coupling of temperature, humidity and velocity. The Rayleigh-Bénard convection is regarded as a classical motion of atmospheric flows under convective stratification. In moist atmosphere, the humid air is heated in the hot ground surface and cooled in the cold top boundary. Density variations of humid air caused by the temperature and humidity variations drive the convection flow under gravity acceleration while the viscosity will counteract to equilibrate the motion which is simplified and modeled as a double diffusive Rayleigh-Bénard convection.

For this test case, we use the 2eq. model where the condensation term is turned off and heat specific capacity c p is set to infinity. Thus the Exner function Π is set to unity and the following temperatures are equal: θ = θ l = T . In addition, the virtual temperature ( 7) is modified as

θ v = θ 1 -1 - 1 ǫ q l + q v ( 34
)
in the buoyancy force (9) as to account only for the total humidity q t = q v + q l . The Prandtl number Pr θ and Pr q are both set to 0.71 for air. The enclosure boundary conditions are non-slip. The bottom and top surface is respectively at a constant temperature 299.5 K and 300.5 K. The specific humidity is set to q t = 0.003 (q l = 0.001, q v = 0.002) and q t = 0 (q l = 0, q v = 0) at the bottom and top boundary, respectively. Adiabatic conditions for both temperature and humidity are set on the right and left boundaries. The Rayleigh number is defined as where θ 0 is 300. For Ra = 10 4 , the non-dimensional relaxation time τ is about to 0.589 and 0.678, respectively. For Ra = 10 5 , the value of τ is approximately equal to 0.528 and 0.556. Figure 2 and Fig. 3 respectively display the contours of potential temperature and total water humidity when the steady state has been achieved. It is seen that both the temperature field and water field obtained by the hybrid LB model are in good agreement with results displayed in [START_REF] Ouertatani | [END_REF]. The reference solution [START_REF] Ouertatani | [END_REF] was obtained using a finite-volume method discretized using QUICK scheme in the momentum equation and a second order central differencing one in the energy equation on a non-uniform grid with 256 2 grids. A comparison between the results obtained by the present LB model and the results in [START_REF] Ouertatani | [END_REF] is reported in Fig. 4. Both the horizontal velocity u through vertical middle line of square enclosure and the vertical velocity v through horizontal center line are plotted and compared in Fig. 4. Both components of the velocity are in very good agreement with the data in literature, indicating that the hybrid recursive regularized LB model can be accurately applied in double diffusive convection problems.

Ra = Pr θ gH 3 ∆θ v ν 2 θ 0 (35)
For further validation, the different τ associated with various values of gravity g, e.g., 9.8, 0.98, 0.098 m s -2 in Eq. ( 35) are investigated in the convection at Ra=10 with different values of τ .

N u(x) = L (θ v,H -θ v,0 ) ∂θ v ∂x | bottom , (36)
as obtained with the present and in [START_REF] Ouertatani | [END_REF], showing an excellent quantitative agreement.

B. A two-dimensional rising moist thermal bubble

A rising moist thermal bubble in atmosphere is examined using the present LB model. The simulation is twodimensional, with a domain height of 2.4 km and width of 3.6 km. Periodic boundary conditions are set on either side of the domain and symmetric nonconducting boundary at the top and bottom surfaces where

∂u ∂z = 0, ∂ 2 ∂z 2 (q v , q l , θ) = 0, at z = 0, H (37)
The initial unperturbed environment is calm (zero initial wind everywhere) and defined by a constant stratification effect, which is introduced by the Brunt-Väisälä frequency,

N = g θ Dθ Dz ≈ 0.0113 s -1 . (38)
A constant relative humidity of 20 percent is given as a base state humid environment with a surface temperature θ 0 and pressure p 0 of 283 K and 85,000 Pa. Therefore, a circular perturbation only on relative humidity is centered in the domain at z = 0.8 km. The relative humidity (RH) is set to 100 percent from radius r ≤ 200 m.

For 200 m ≤ r ≤ 300 m, a transition layer is assumed in which

RH = 20% + 80%cos 2 π 2 r -200 100 ( 39
)
The vapor specific humidity q v is converted from

q v = RH ǫ.p sat p 0 -(1 -ǫ)p sat (40)
The kinematic viscosity ν is assigned 1 m 2 s -1 and Prandtl number Pr θ and Pr q are both set to unity in the simulations. The simulations based on the present hybrid recursive regularized lattice Boltzmann method with condensation model (1eq) and (2eq) are carried with δx = 5 m. The time step is calculated by δ t = (c s /c sp )δ x with a sound speed in physical unit c sp = 85 m s -1 . Thus, the non-dimensional relaxation time τ ≈ 0.504 is used in the simulation. Results of the simulations are presented in Fig. 7. Due to the symmetry of the case tested, the streamlines and mass fraction of liquid water obtained by 2eq model are plotted in the left sub-domain (1.1-1.8 km) while the solution calculated from 1eq model is presented in right sub-domain. In excellent agreement with the results in Ref.[see 58, Fig. 6], the warm moist bubble rises and expands over time. The phase transition from vapor to liquid gradually happens, while the bubble rises. The benchmark solution in [58] was obtained with a multidimensional positive definite advection transport algorithm on 2.5 m grid spacing.

In order to quantitatively validate the proposed LB model, the highest vertical position of contour of 20% maximum q l (denoted as H 20 ) and vertical fluid velocity on the top central position of interface (denoted as W f ) are given in Table II. The computed velocities at 3, 5 and 7 min are in good agreement with benchmark solution given in Ref. [58]. The benchmark solution was obtained using an anelastic flow solver using variables (θ, q v , q l ) on 2.5 m grid spacing. The vertical position of contour of 20% maximum q l reported by our LB models is consistent with reference values at 3 min. However, the deviation on the vertical location between our models and reference values gradually increases with physical time. This deviation may be the consequence of Boussinesq approximation which is used in our model rather than the anelastic approximation adopted in Ref. [58].

To further investigate numerical stability of the LB model, random perturbations with a small amplitude are introduce during the simulation as same as calculations in Ref. [58]. An excitation at t = 3 min was provided by adding a perturbation to potential temperature θ for 1eq model θ l for 2eq model which ranged between ±0.05 K on every node.

Figure 8 and Fig. 9 show contours of vorticity and liquid humidity field after excitation of instability on resolution of 5 m at 7 min and 10 min, respectively. The red solid lines represents vorticity from -0. Besides, a slight different of cloud formatting is obtained by 1eq model and 2eq model at 10 min. The 2eq model obtained a little more complex structure of vorticity and cloud when the bubble rose to the high altitude, which is consistent with the solution in [58]. A complex instability was excited in their results which was obtained by two equation condensation model using monotone finite difference method on 2.5 m resolution. The difference between our two results and the difference between our results and their solution could be caused by different numerical characteristics of two methods, which strongly influenced the interactions between numerical noises and excitation instability. Considering the 2eq model is more general in both shallow and deep atmospheric flows and that it can be conveniently extended for involving the precipitation microphysics, only the 2eq model is used in the following study.

C. A three-dimensional rising moist thermal bubble

The simulation of the rising moist thermal bubble is extend to three-dimensional study in this section. The investigated domain is three-dimensional, with a domain height of 2.4 km, length of 3.6 km and width of 3.6 km. Periodic boundary conditions are specified four laterals of the domain and free slip nonconducting boundary at the top and bottom surfaces, which is the same with the setting in previous 2D simulation. The same initial unperturbed environment is adopted and defined by the constant stratification effect. A sphere perturbation only on relative humidity is centered in the domain at z = 0.8 km. The relative humidity (RH) is set to 100 percent from radius r ≤ 200 m. For 200 m ≤ r ≤ 300 m, a similar transition layer is assumed with the one in the previous 2D case. The present lattice Boltzmann method with condensation model (2eq) incorporated with grid refinement technique are assessed in this case. The sub-domains with different grid resolution is used as illustrated in Table III and Fig. 10. The kinematic viscosity ν is assigned 1 m 2 s -1 and Prandtl number Pr θ as well as Pr q are set to unity in the simulations. The time step is calculated with a sound speed in physical unit c sp = 20 m s -1 . The corresponding non-dimensional relaxation time parameters are also given in Tab. III.

Results of the simulations with grid refinement after 2, 4 and 6 min are presented in Fig. 11 In order to quantitatively validate the proposed LB model, the liquid humidity and vertical velocity along mid-line are plotted in Fig. 13. The computed velocities at 2, 4 and 6 min are confirmed closely with benchmark solution in Ref. [59]. The benchmark solution was obtained by an anelastic solver using variables (θ, q v , q l ) on minimum grid spacing of 6.25 m.

D. Shallow Cumulus Convection

At last, a shallow cumulus convection is simulated to evaluate the capability of the present LB method to capture the moist thermodynamics and its interaction with the meteorological flows. The shallow cumulus convection simulations follows the setup of the Barbados Oceanographic and Meteorological Experiment (BOMEX) model inter-comparison case [35]. This is the most prevalent shallow cumulus large-eddy simulation (LES) case.

In this case, an altitude-dependent geostrophic wind U g is given by a linear formula U g = (-10 + 1.8 × 10 -3 z) m s -1 and the Coriolis parameter is set to f = 0.376 × 10 -4 s -1 . The initial conditions for velocity, liquid water potential temperature and total water mixing ratio are linear profiles following the values given on Table IV. The temperature and humidity surface fluxes are 8 × 10 -3 K m s -1 and 5.2 × 10 -5 m s -1 , respectively. The shear stresses are prescribed by u i w = -u 2 * u i /(u 2 1 +u 2 2 ) 1/2 , with u * = 0.28 m s -1 . Moreover, additional terms are added to represent the large-scale forcing which could not be represented directly in the LES. The source terms of momentum conservation equations, temperature equation and water equations are parameterized considering the effects of largescale subsidence, radiative cooling and moisture effects. The details on large-scale forcing were reported in [35].

The computational domain size is 5000m × 5000m × 3000m. Here, the lattice Boltzmann method with condensation scheme 2eq is used to compute that case considering on grid resolution of δ x = 40 m. The time step is equal to 0.27s with c sp being 85 m s -1 in the simulation. The kinematic viscosity of air is 1.5 × 10 -5 m 2 s -1 and Prandtl number P r θ = P r q = 0.71 for potential temperature and water fraction are used in the simulation. Thus the non-dimensional relaxation time τ based on the kinematic viscosity is approximately equal to (0.5 + 7 × 10 -9 ). The classical Smagorinsky model is adopted as the subgrid model in the large-eddy simulation [60] and the surface model of the horizontal momentum components, temperature and humidity is using the Monin -Obukhov similarity theory (MOST) [61].

The incorporating LES within framework of lattice Boltzmann method for atmospheric boundary layer flows are detailed in [62]. The constant of the Smagorinsky subgrid viscosity model is taken equal to C S = 0.15 and turbulent Prandtl numbers of K t = 0.33 and K q = 0.33 are adopted both for potential temperature and water fraction equations. The simulation is carried out over 6 hour of physical time. The statistics are computed averaging the solution over 1 hour. Figure 14 shows the computed vertical mean profiles of velocities, potential temperature, vapor water and liquid water. All of the results are in good agreement with the reference data from [35], proving the accuracy of the proposed LB model with condensation scheme. The mixed region below the altitude of 540 m is well captured, demonstrating the accuracy of the wall model with complex thermophysics and cloud dynamics. Furthermore, the conditionally unstable layer from 540 m to 1500 m, and the inversion layer from 1500 m to 2000 m are also clearly observed in the results.

Figure 15 shows iso-surfaces of liquid water and temperature in the full computational domain. In order to emphasize interaction between temperature field and humidity field as well as instantaneous flow features of LES, cloud (iso-surface of liquid water) is displayed. The results is obtained at time of the 1.5th hour. To highlight the thermal plumes near surface, the potential temperature is chosen equal to 298.

V. CONCLUSION

In this paper, a hybrid LB based condensation models for moist atmospheric convection is proposed. A recursive regularized BGK model with forcing term is employed in LB equation and the finite volume scheme is adopted on temperature and water transport equations. Four benchmark problem have been used to assess the LB method incorporated two condensation models. Both the 1eq model and the 2eq model present accurate solution in a 2D moist rising bubble simulation. In summary, the 2eq model is more general in both shallow and deep atmospheric flows and the 2eq model is preferred as a framework for involving the precipitation microphysics in the future. The proposed LB approach for humid air and cloud dynamics has been assessed considering a 3D moist rising bubble simulation and cloudy convective atmospheric boundary layer with phase change and large scale sources. To further develop a LB model for moist atmospheric deep convection, an appropriate model under the anelastic approximation will be proposed in the future.

FIG. 1 .

 1 FIG.1. The flowchart of LB method with condensation model (2eq). (potential temperature θ, vapor humidity qv and cloud liquid specific humidity q l ).

FIG. 2 .

 2 FIG. 2. Temperature contours obtained by the LB method (2eq. model) for Ra = 10 4 (top), and Ra = 10 5 (bottom).

FIG. 3 .

 3 FIG. 3. Contours of water fraction obtained by the LB method (2eq. model) on Ra = 10 4 (top), Ra = 10 5 (bottom).

  5 K and kinematic viscosity ν is calculated by Ra given in the simulations. Simulation of double diffusive Rayleigh Bénard convection with Ra = 10 4 and 10 5 are presented hereafter, on 50 × 50 and 100 × 100 grids. An equivalent Reynolds number can be computed as Re = Ra/Pr θ , which approximately equals 119 and 375, respectively. The time step is calculated by δ t = (c s /c sp )δ x with a sound speed c sp = 1 m s -1 .

FIG. 4 .

 4 FIG. 4. Profile of velocity along midline of domain (top: velocity u along y; bottom: velocity v along x). The double diffusive Rayleigh Bénard convection with Ra = 10 4 and 10 5 are presented on 50 × 50 and 100 × 100 grid resolutions by the proposed LB method.

4 on 50 ×FIG. 5 .

 505 FIG. 5. Comparison of velocities (u, v) along midlines of domain with different non-dimensional relaxation time τ at Ra=10 4 on 50 × 50 grids.

FIG. 6 .

 6 FIG.6. Local Nusselt number through the hot wall. The double diffusive Rayleigh Bénard convection with Ra = 10 4 and 10 5 are presented on 50×50 and 100×100 grid resolutions by the proposed LB method.

Figure 6

 6 Figure 6 reports a comparison of the local Nusselt number

FIG. 7 .

 7 FIG. 7. Part I: Streamlines and contours of liquid humidity field at 2 min (top) and 3 min (bottom) obtained by the present LB method with 1eq model and 2eq model. In each figure, the left part is obtained by 2eq model while the right part is from 1eq model.

  FIG. 7. Part II: Streamlines and contours of liquid humidity field at 5 min (top) and 7 min (bottom) obtained by the present LB method with 1eq model and 2eq model. In each figure, the left part is obtained by 2eq model while the right part is from 1eq model.

FIG. 8 .

 8 FIG. 8. Contours of vorticity and liquid humidity field obtained by 1eq model (top) and 2eq model (bottom) at 7 min with excitation of instability on resolution of 5 m in domain of [1100 m, 2500 m] × [800 m, 1800 m].

No

  FIG.10. Setup of three-dimensional geometry and subdomains for simulation of 3D rising moist thermal bubble.

  FIG. 11. Contours of liquid humidity field (left) and vapor humidity field (right) at 2, 4 and 6 min (top to bottom)as obtained by the present LB method with 2eq. condensation model.

FIG. 12 .

 12 FIG. 12. Contours of vertical velocity (left) and x component of vorticity (right) at 2, 4 and 6 min (top to bottom) as obtained by the present LB method with 2eq model and grid refinement.

FIG. 13 .

 13 FIG.[START_REF] Feng | [END_REF]. Profile of liquid humidity and vertical velocity along vertical middle line at 2, 4, 6 min. The symbols represent benchmark solution in Ref.[59].

FIG. 14 .

 14 FIG.14. Mean profiles of vapor water (top, left), liquid water (top, right), velocity (bottom, left) and potential temperature (bottom, right) obtained in the convective cumulus convection case. The reference data in symbols are results at 6 hr from[35].

  Figure15shows iso-surfaces of liquid water and temperature in the full computational domain. In order to emphasize interaction between temperature field and humidity field as well as instantaneous flow features of LES, cloud (iso-surface of liquid water) is displayed. The results is obtained at time of the 1.5th hour. To highlight the thermal plumes near surface, the potential temperature is chosen equal to 298.7 K, and the color represents the altitude. The results are also at time of the 1.5th hour with δ x = 40 m in the figure. The thermal plumes are well generated by the coupling between LB core and passive scalars solvers. It can be observed that the instantaneous formation of cloud is well captured by the present lattice Boltzmann model with condensation scheme.

FIG. 15 .

 15 FIG. 15. Iso-surface of liquid water (top) and temperature (bottom) at the 1.5 hour for shallow cumulus convection.

TABLE I .

 I Governing equations for the 2eq. and 1eq. models.

	Governing equations	2eq. model	1eq. model
	Mass conservation	Eq. (1)	Eq. (1)
	Momentum conservation	Eq. (2)	Eq. (2)
	Energy conservation

TABLE II .

 II The highest vertical location of contour of 20% maximum q l

TABLE III .

 III Setup of grid refinement for simulation of 3D rising moist thermal bubble.

TABLE IV .

 IV Initial conditions for shallow cumulus convection.

	0	17.0	298.7	-8.75	0
	520	16.3	298.7
	700			-8.75
	1480	10.7	302.4
	2000	4.2	308.2
	3000	3.0	311.85 -4.61	0

Height (m) qt (g kg -1 ) θ l (K) u (m s -1 ) v (m s -1 )

ACKNOWLEDGEMENTS

The first author would like to gratefully acknowledge to Dr. Johann Miranda for fruitful discussions. This work was supported by the French project CLIMB, with the financial support of BPIFrance (Grant P3543-24000), in the framework of the program "Investissement d'Avenir: Calcul Intensif et Simulation Numérique". This work was performed using HPC resources from GENCI-TGCC/CINES (Grant 2018-A0052A07679). This work was also supported by project of ProLB software (http://www.prolb-cfd.com).


Giulia Lioi

Simon Butet

Mathis Fleury

Anatole Lécuyer

Isabelle Bonan

Christian Barillot

Ef cacy of EEG-fMRI Neurofeedback in stroke in relation to the DTI structural damage: a pilot study

Introduction:

If many interventions for motor rehabilitation after stroke are available, recent studies have shown the potential of neurofeedback (NF) (Wang, Mantini and Gillebert, 2017). The majority of these NF approaches have relied solely on one imaging technique: mostly on EEG recordings. If EEG offers the advantages of practicability and high time resolution, it suffers from a limited spatial resolution and access to deeper areas of the brain. This has prompted the effort of exploring the implementation of NF approaches with other sensing techniques such as fMRI (Sitaram et al., 2012), which allows to more precisely identify cortical targets. Recent study have gone further (Zotev et al., 2014;Perronnet et al., 2017), revealing the potential of integrating complementary techniques such as EEG and fMRI to achieve a more speci c regulation. In this exploratory work, for the rst time multisession bimodal EEG-fMRI NF for upper limb motor recovery was tested in four stroke patients. In particular, the feasibility and functional ef cacy of the NF training were investigated with respect to the integrity of the corticospinal tract (CST), a well-established predictor of the potential for clinical improvement.

Methods:

Four chronic stroke patients (54-76 years, 2 females) with left hemiparesis took part to the study. At inclusion, diffusion tensor imaging (30 directions DTI, b0=1000) was performed to assess the asymmetry between ipsilesional and controlesional CST. To this end, the diffusion tensor model was estimated and the fractional anisotropy (FA) calculated. The CST was then reconstructed (Jong et al., 2005) and an index of FA asymmetry between the affected and unaffected CST was calculated. Disruption of the bers integrity is associated to a FA decrease therefore an index of FA asymmetry gives important indication about the structural de cit. In particular it has been shown that for a FA asymmetry index value greater than 0.25 limited capacity for recovery are expected (Stinear et al., 2007). The experimental protocol for NF training included an alternation of bimodal EEG-fMRI and unimodal EEG NF sessions. We expected that during bimodal NF training the patient, receiving richer and more speci c information, could develop a strategy that could then be "transferred" to unimodal EEG sessions. During each NF session, patients alternated 16 20 seconds blocks of rest and motor imagery of the affected upper limb with NF. Information about their brain activity was given to the patients by means of a visual feedback re ecting the EEG and BOLD activity of regions of interest identi ed in the ipsilesional motor cortex after a preliminary calibration block. Results:

All the patients were able to self-regulate their brain activity in the motor cortex during the NF training. The analysis of the CST integrity revealed that for 3 patients a high degree of symmetry between ipsilesional and controlesional CST was preserved (FA asymmetry equal to 0.03, 0.06 and 0.05). For one patient however, the ipsilesional CST was more affected, giving an index of FA asymmetry of 0.105, closer to the threshold indicating very poor recovery potential. This patient also failed to up-regulate the brain activity in the ipsilesional motor cortex at the end of the training with respect to the rst session (Figure 1). On the other hand, 2 of the 3 patients exhibiting a high degree of integrity of the ipsilesional CST showed signi cant increased activation of the ipsilesional M1 at the end of the training (p<<0.001, Wilcoxon test) and exhibited a larger involvement of the ipsilesional motor and premotor areas. These preliminary ndings con rm the critical role of the CST integrity for stroke motor recovery and indicate that this is importantly related also to functional brain regulation of the ipsilesional motor cortex.

  ww5.aievolution.com/hbm1901/index.cfm?do=abs.viewAbs&subView=1&abs=2634 2/4

Conclusions:

This exploratory study indicates that functional performances are related to CST integrity of and give useful indications for future studies inclusion criteria.

Disorders of the


Simon Butet

Giulia Lioi

Mathis Fleury

Anatole Lécuyer

Christian Barillot

Isabelle Bonan

A multi-target motor imagery training using EEG-fMRI Neurofeedback: an exploratory study on stroke

Introduction:

Upper limb recovery after stroke is a complex process involving a variety of structures and mechanisms, and current rehabilitation techniques have limitations (Langhorne, Bernhardt and Kwakkel, 2011). Recent studies have revealed the potential of neurofeedback training (NF) (Cervera et al., 2018) as an alternative or in aid to traditional therapies. Studies on cerebral plasticity and recovery after stroke indicate that premotor areas (supplementary motor area -SMA, premotor cortex PMC) should be a preferred target for NF in the most severe patients (Plow et al., 2015) while M1 stimulation may be more ef cient for patients with better recovery potential (Favre et al., 2014). Moreover, fMRI-NF studies (also on stroke patients) have shown that SMA is a robust correlate of motor imagery, while the activation of M1 is more dif cult to achieve, especially for short training sessions (Mehler et al., 2019). Based on these results, in this exploratory work we tested a dynamic NF training more strongly rewarding SMA activation in the rst NF training session and then increasing the M1 activation contribution in the nal NF session. We tested this novel approach on four stroke patients in a multisession bimodal EEG-fMRI NF training. To this end, we de ned an adaptive cortical region of interest (ROI) equal to a weighted combination of ipsilesional SMA and M1 activities and then varied the weights in order guide the patient training towards an improved activation of M1.

Methods:

Four chronic stroke patients with left hemiparesis participated to the study. The experimental protocol included an alternation of bimodal EEG-fMRI NF (N=2) and unimodal EEG-only NF (N=3) sessions. For each bimodal NF session, the protocol included a calibration step (motor imagery of hemiplegic hand) and 3 NF training blocks (5 minutes each). Each NF block consisted of epochs of rest (20 s) alternated to period of closed-loop MI training (20 s). MRI imaging was performed using a Siemens 3T Prisma scanner and 64 channels EEG was recorded using an MRI compatible Brain Product solution. The experiment was run using a NF platform (Mano et al., 2017) performing real-time EEG-fMRI processing and NF presentation. The NF metaphor consisted of a ball moving on a gauge proportionally to the average BOLD and EEG activity in regions of interest (ROI) identi ed over the ipsilesional motor cortex during calibration. The fMRI NF score was calculated as a weighted combination of the activity in the SMA and M1 cortical areas. For the rst session, an equal weight of 0.5 was attributed to SMA and M1, while in the second session the activity in the primary motor area was much heavily weighted (0.75) than the SMA (0.25) (Figure 1).

  All the patients were able to upregulate the brain activity in the target cortical areas during the NF training (as compared to rest). The increase in brain activation during the NF task was particularly signi cant for the BOLD signal, but varied depending on the ROI targeted (SMA, M1) and the subject. The EEG activity was harder to modulate than the BOLD activity. Two patients, that had a less severe upper limb de cit, showed signi cant increased activation of the ipsilesional M1 at the end of the training (p<<0.001, Wilcoxon test across training blocks) and exhibited a larger involvement of the ipsilesional motor and premotor areas (with respect to controlesional ones) in the second NF training session (Figure2). These preliminary results, on a short training duration, indicate the potential of a dynamic, multi-target/multimodal NF training approach on chronic stroke patients.

Disorders of the Nervous System: