On the reproducibility of ultra-intense terahertz experiments: Comment on "Sub-cycle insulator-to-metal transition in vanadium dioxide by terahertz-field-driven tunneling" by Giorgianni, Vicario,.. and Hauri (arXiv:1706.00616) Mostafa Shalaby* SwissFEL, Paul Scherrer Institut, Villigen 5232, Switzerland

The manuscript by Giorgianni et al [1], and co-authored by myself, claims and concludes on a purely electronic insulator to metal transition (IMT) in VO 2 using ultra-intense THz fields without any lattice interaction. The underlying mechanism behind IMT in VO 2 is a highly debatable subject with the "electron-only" IMT concept being increasingly resisted in the field especially under intense excitations. If Giorgianni/Hauri's claim were sustained, it would be an important step towards fundamental understanding of complex systems such as the model VO 2 [2][3][4] and may have implications on low power ultrafast switching applications. The experiment was enabled by the uniquely-intense THz source I previously developed at the SwissFEL [5,6] and is thus not possible to reproduce elsewhere under similar excitation conditions. However, some of the measurements under the reported experimental settings showed clear phonon oscillations on the short time scale, and sample damage using the maximum reported field intensities, which may be contradicting the purely electronic/non-phonon paper claims/conclusions. These measurements were known to all PSI authors (Feb./March 2017) before the paper was drafted. However, the principal authors (Giorgianni, Vicario, and Hauri) 1 decided to exclude these measurements from the paper, disapproved the call to further repeat and verify the paper measurements, and submitted it. They then publically published the text on the arXiv without the consent of all the authors. The external co-authors (from U.C. Berkeley, USA, and Tsinghua U., China) requested the addition of the note "We (Kai Liu, Junqiao Wu and Kevin Wang) were unaware of the situation, of the hidden data, and of the PSI disagreements, and we have never been contacted or informed of any of these issues. We merely provided the samples." † The experiment was performed at a SwissFEL laboratory, but not the experimental user laboratories. This comment should not lay responsibility on the SwissFEL facility. PSI attempted to withdraw the submission, but the arXiv is non-retractable. ]. Generally speaking, phonons can also be measured through ∆" but the effect is negligible when compared to that of IMT [7]. The arXived text [1] shows only ∆" In the original text [1], these oscillations are not obvious because the paper shows only the transmission and focuses on long wavelength probe to magnify the change in the transmission [8]. Both ∆" and ∆" were taken under the same conditions (experimental setup [6], sample, 2x 20 THz + 1x 6 THz QMC LPF excitation filters). Lower cut-off filters (4.2 THz) have also shown trace of phonon excitations. These oscillations match the previously reported phonons in VO 2 (Fig. b) [9,10].

Sample damage

The measurements in the paper [1] were done under very high field excitation of up to 18 MV/cm in the sub-6 THz range. When the experiment was repeated (Fig. a), gradual sample damage was observed at such fields and exactly same experimental conditions in ref. [1] (such a damage at such a low repetition rate typically leads to irreproducible transition traces and reduction in the unmodulated transmission of the optical probe). Therefore, to take the scans in Fig. a, the maximum field was reduced (< 9 MV/cm) keeping the same spectral contents. This observation is consistent with previous report on damage around 4 MV/cm in [2]. However, the comparison between the two reports is complex. Ref. [2] uses higher (1 kHz) repetition rate (that has lower damage threshold) but the excitation frequency was low (~ 1 THz), this is stronger tunneling effect. In the present experiment, the repetition was 100 Hz but the excitation frequency was high and VO 2 absorption is higher. Nevertheless, sample damage changes the surface properties, transmission, and complicates the overall structural features being measured.

Conductivity calculation (minor)

The used formula "Tinkham" is an approximation that is valid only under specific conditions (when the film is so thin that the probe does' not experience phase change across the film). This is clearly not satisfied as the transmission change is > 30%. The fit on a logarithmic (tunneling) effect does not indicate much. Unless you have a very large dynamic range (not the case here), any high order function will probably give good matching.

Conclusion:

The published text on the arXiv by Giorgianni et al [1] does not reflect the opinion of all the listed authors. It does not fully represent the undertaken measurements. The excluded phonon measurements represent an integral part of the overall structural dynamics and may contradict the paper's conclusion. This case may raise concerns on the reproducibility of results in ultra-intense experiments in the THz community.

Science is based on peer-review and reproducibility. Outside the SwissFEL laboratories, the most powerful existing THz sources (in the sub-6 THz range) are 100 times less intense (in W/m 2 ) than the one used for this experiment. Such intensities have shown to me nonlinear response from every sample I used in the recent years. However, extensive repetition and in-depth experimental analysis have proven that many of them were misrepresentations of results in an unprecedented regime of ultrastrong-low frequency fields. Most of these issues are not known in the field to the wide scientific (reviewers) community. Nor is there a way of getting such field intensities to reproduce such experiments elsewhere. In this regard, the underlying theory/simulations papers typically focus only on what they want to prove and do not deal with the associated unwanted other correlated dynamics. The author predicts this theory-experiment deviation to be a major issue in the coming years in the field of ultra-strong THz-induced structural dynamics, beyond the perturbative regime.

Sub

Transmission/ reflection

Very weak (The reason why it is not clear in ref. [1])

Tunneling is extremely strong (reported in ref. [1])

Birefringence/ polarization rotation.

Very strong (excluded from [1]) strong varies MOKE configuration + ! direction.

2 1 .

 1 Phonon excitationsIn transmission-type THz pump/optical probe spectroscopy, the basic experimental techniques are transmission modulation ∆" and birefringence/polarization rotation ∆" . Understanding of the structural dynamics in VO 2 strongly depends on the ultrafast electron and lattice dynamics which are coupled and difficult to disentangle using an optical probe. While ∆" reflects the conductivity change and IMT, ∆" is the tool to trace phonon dynamics (excitations) [Table

  Fig. (a) THz-induced ∆" on a collinear 800 nm probe in 90 nm-VO 2 . Different curves test on different THz excitation levels and 800 nm probe polarization. They all reveal traces of the 5.5-6 phonon. Exactly the same sample and experimental THz pump conditions as used for ∆" in [1].Fig. (b) Phonon spectrum in VO 2 taken from ref. [9] and measured in a degenerate pump-probe experiment 1.55 eV. The spectrum corresponds to the THz oscillations in (a).

  Sensitivity could vary significantly depending on the probe photon energy, incidence angle, etc.

			-system		Phonon	Electron		Spin
		THz excitation	resonant	non-resonant	instantaneous	delayed	resonant (magnon)	non-resonant
	POWER LAW	Material damage occurs here	! (IR active phonons) ! ".$%& (ISRS); coupling to ! " (ISRS), (Raman & IR active phonons) other phonons, change phonon bandwidth/frequency	! " (nonlinear optics Kerr), any material Exponential response, easily seen ! " (thermalization), mainly conductors and semiconductors in semiconductors	Same as resonant	phonon excitation	∓" low field ! " high field till ! " Coherent Magnetic switching (debatable)
		Optical probe
	Signal / noise	Rep. rate dep.

In connection with [1] submission: Affiliations: C.V.: Staff scientist, SwissFEL, PSI; ext. funding: IZLRZ2_164051(Swiss SNF); C.P.H.: Head of the SwissFEL laser group, PSI/ EPFL (Lausanne); ext. funding:

200021_146769 (Swiss SNF).2 PSI: was informed on June 26 th , 2017; requested from Nature Publishing Group to put the manuscript under review on hold on July 6 th , 2017; attempted to remove the arXived version[1] on Sept. 1 st , 2017. *most.shalaby@gmail.com ; †wuj@berkeley.edu

Table: Crib sheet for THz-pump structural dynamics experiments:

A rough and very generalized picture of how experimental nonlinear THz experiments could be designed and judged. (mainly for optical probe)


Abdelghani Chahmi

Abdelhak Djoudi

Diagnosis of the Induction Machine by the Kalman Filter

Keywords: Induction Motor, Modeling, Rotor Defect, Stator Defect, Extended Kalman Filter

come

I. INTRODUCTION

arge industrial systems are widely confronted to the challenges of reliability and availability of the production devices. Electric drives containing an induction machine are largely used in the industrial applications thanks to their low costs, high performances and robustness [START_REF] Chahmi | Fault Detection In Electrical Drives -Approach Signal[END_REF]. As known, a system can realize the assigned task, only under conditions ensuring the security. An early detection of anomalies in electrical motors may help avoiding downtimes [START_REF] Ferradj | Synthesis of Electrical Quantities Applied to Squirrel Cage Induction Motor Bearing Faults Detection[END_REF]. Thus, in order to detect an incipient fault, we should give a special attention to the spectral analysis of stator currents [START_REF] Boudinar | Three Phase Induction Motor Incipient Rotor's Faults Detection Based on Improved Root-MUSIC Approach[END_REF], [START_REF] Lebaroud | Classification of induction machine faults by optimal time-frequency representations[END_REF]. Most of recent researches on defects detection of induction machines use the stator currents analysis. Several condition-monitoring methods [5] have been developed and were used to help the detection of electric defects in induction machine. Among these methods, the authors in [17] consider the vibrations monitoring, torque control, the temperature monitoring, etc. In this paper, the rotor resistance is estimated using a method of Extended Kalman Filter (EKF). This state observer is proved to be a useful tool to determine the parameter variations [6], [7]. In this paper, it is shown that EKF can be employed to identify the disturbance in rotor resistance. Recently, the Kalman filter was applied successfully for defect forecasting in [8]. For example it is used to interpolate the trends of a signal-learned g from the recorded data and then estimate the evolution of defects in the proposed diagnosis methods in [9]. In this study, the proposed approach considers the value of rotor resistance as fixed for condition monitoring. This value in the diagnostic tools which one uses is not fixed contrary to the classical approaches of control of machine. Hence, we propose a method that considers a variation of rotor resistance in order to give a good estimate as well as a good robustness for the control of the induction machine. This paper is organized as follows: part II describes the induction machine model under simplifying assumptions. Induction machine models have been developed under healthy and faulty conditions in the (abc) reference frame. Part III, presents a state of the art review of Extended Kalman Filter. Finally, part IV gives digital simulations of the proposed faulty model. The obtained results prove the effectiveness of the proposed monitoring method based on an Extended Kalman Filter.

II. INDUCTION MOTOR DRIVES

A. Overview of the Drive

The drive under study is composed of a machine with induction, inverter of width of impulse, ordering of rotor flow -directed, a loop of measurement of the current and a loop of speed control. Figure 1 presents the total diagram of the unit.

Asynchronous machi ne Inverter

B. MODEL OF THE INDUCTION MOTOR

The objective of this section is to present the different models used for our study. This is first of all to model the induction machine finely enough to introduce a relevant modeling of defects. Then we present how the default templates are not used in their simplified and balanced version diagnostic purposes.

1) Models Simulation Proposed in the Reference Frame (abc):

The induction machine is a known system-well [10]; models exist and allowing, in the majority of the cases and under certain simplifying assumptions, a representation of the machine without fault. These assumptions are the following ones [11-12]:

-The air-gap is of constant width without effect of notch (not of eccentricity, of unbalance, reduction amongst harmonics in the development of inductances expression).

-Assimilation of the cage to a short-circuit having the same number of phases as static winding (simplification of the

Diagnosis of the Induction Machine by the Kalman Filter

Abdelghani Chahmi, Abdelhak Djoudi L equations of the machine by taking into account the average effect of the cage).

-Sinusoidal distribution, along the air-gap, of the magnetic fields of each winding (simplification of the equations of the machine by taking into account of the expressions in complex of the currents, allows to use a vectorial model).

-Absence of saturation in the magnetic circuit (identical behavior of each phase, allows use of a vectorial model).

-The influences of the effect of skin and the heating of the drivers are neglected (not variation of the parameters intervening in the machine under normal functioning).

-Star coupling of windings (the currents form a system balanced under normal functioning).

Under the conditions pointed out previously, the equations of the electrical circuits utilize clean and mutual inductances making it possible to define flows according to the currents [11].

On the level of the stator, the tensions check:

                                  sc sb sa dt d isc isb isa R Vsc Vsb Vsa s    . ( 1
)
Where the resistances matrix is written according to the supposition of symmetry

             s s s s r r r R 0 0 0 0 0 0 (2)
On the level of the rotor, under -assumption-that -rotor-is comparable to a rotor with three-phase rollings up, the voltage are written

                                  rc rb ra dt d irc irb ira R Vrc Vrb Vra R    . ( 3
)
Matrix RR has the same form that RS, flux vector are expressed according to the inductances matrix and stator and rotor currents

                                    irc irb ira M isc isb isa L sc sb sa SR s . .    (4)                                     irc irb ira L isc isb isa M rc rb ra R RS . .    (5)
In order to be able to give an account of the interactions between the phases, we expressed each resistance and inductance matrix according to the number of whorls of rollings up of the machine.

                                               
                                                   
                    ) cos( ) 3 / 2 cos( 3 / 2 cos( ) 3 / 2 cos( ) cos( ) 3 / 2 cos( ) 3 / 2 cos( ) 3 / 2 cos( ) cos( ) ( )                sr rc
m n n m n n m n n m n n m n n m n n m n n m n n m n n M (8)     T RS SR M M ) ( ) (   

2) Cage Electric Failures in the (abc) Model:

This model that we have just described in the previous part allows in particular creating asymmetrical faults in both stator and rotor by the variation amongst whorl of the rolling up of the phase at fault. The defect can be localized on one or more phases. It should be noted that this modeling can be criticized: the defects of type short-circuit are represented in the model like removals of whorls i.e. openings partial of phase; that does not correspond completely to reality since it would also be necessary to take account of the shading rings [12]. This model constitutes a progress compared to the balanced models or /and considering only one variation of resistance usually used to validate the algorithms [13].

A. Model "The Good Performance" for Control and of Monitoring

Models are necessary to design and improve the controllers and the observers. They are commonly balanced. One passes from the reference frame (abc), fixes, with the fixed reference frame (α, β) by the transformation of Concordia [14].

      X C Y U B X A pX c c c With X= [isα isβ 𝜑r𝛼 𝜑r𝛽] T                                           1 0 1 0 0 0 m m c L L K A T s s c L L B              0 0 1 0 0 0 0 1          0 0 1 0 0 0 0 1 c C and        1 s s L R r r R L   r s m L L L    r s m L L L 2 1    D.

Induction Machine Discrete Model

We develop in this paragraph a discrete model of an induction machine in order to be able to use these equations in the development of Kalman Filter. We consider indeed that the equations of the observer at ends of diagnosis will be inevitably used in a system containing microprocessor or other and thus entered in the form of a discrete model. The discrete version is obtained by using the formula of Euler:

          ) ( ) ( ) ( ) ( ) ( ) 1 ( k x C k y k U B k x A T k x k x c c c s Of or         ) ( ) ( ) ( ) ( ) ( ) 1 (

III. INDUCTION MACHINE OBSERVER-BASED EKF

This part presents the proposed method to supervise the application. We will develop the method by using the observer model previously described.

A. The Extended Kalman Filter

The basic equations of the extended Kalman filter are deduced from that of the discrete model of the machine in good fonctionnement.la discrete matrix of state. One introduces noises of state and noises of measurement [15] - [16]. The function F is defined so that one to extend filter of Kalman has all the parameters of the machine. The estimate of rotor resistance is considered like means of monitoring.

         ) ( ) ( ) ( ) ( ) ( ) ( ) 1 ( k v k x C k y k w k U B k x A k x d d d Where   ) ( ) ( ) ( ), ( k U B k x A k U k x f d d  
The process noise w (k) is characterized by

    Q k w k w E k w E T   ) ( ) ( 0 ) ( Where Q is definite positive. The measurement noise v (k) respects:     R k v k v E k v E T   ) ( ) ( 0 ) ( R is definite positive too.
If we call P the matrix of covariance of the error in estimation, this matrix is obtained to leave:

     T k x k x k x k x E k k p ) ( ) ( . ) ( ) ( ) / (   
The notation k/k corresponds to a prediction at time K based on data up to and including time K (a priori covariance) and k+1/k corresponds to a prediction at time k+1 starting from the data until the time K (a posteriori covariance). In the following equations, t represents the Jacobian matrix of the nonlinear function f is used

[6]. ) ( ), ( ) ( k U k x k x f J   
Applied to the system under study, J can be expressed as

                                    1 0 0 0 0 0 1 0 0 1 0 1 1 0 1 0 1              
Q J k k JP k k P T    ) / ( ) / 1 ( (17)
3/Computation of the outputs of the filter

  ) / 1 ( ) 1 ( ). 1 ( ) / 1 ( ) 1 / 1 ( ˆk k x C k y k G k k x k k x d d          (18)
4/Computation of the feedback gain matrix of the Kalman Filter

    1 ) ( ) / 1 ( ) ( ) 1 ( ) / 1 ( 1        R k C k k P k C k C k k P k G T d d T d d (19)
5/ Computation of the a priori covariance matrix P(k+1/k+1)

) / 1 ( ) 1 ( ) 1 ( ) / 1 ( ) 1 / 1 ( k k P k C k G k k P k k P d d         (20)
We consider that Q and R matrix are constant but they can be time variant. The behavior and the stability of the Kalman Filter are functions of Q and R.

IV. Simulation Results and Discussion

The test conditions are the following:

-From t=0 to t=0.3s, flux application -From 0.3 to 1.8s, liner increase of the speed reference -From 1.8s constant speed -At t=3s, fault in the machine The observers are enabled from beginning.

A. Good performance

Our approach in this part of simulation is to see the behavior of the induction machine in healthy mode and with electric defect (rotor, stator, inverter) in the mode transitory and permanent, after from there will introduce an observer on the outlet side of the machine (extended Kalman filter) is to see its reaction. The results result from rotor resistance, the stator current FFT.

The following figure takes the form of the current in a phase of the induction machine. One can distinguish there the various phases related to the order from the system: phase of fluxing, phase of application of the instruction of couple related to the slope of speed, then stabilization the speed and instruction of couple. Up to 0.3 second, the purpose of the applied voltages is to position the vector flow in the machine, as from this moment, an instruction of couple is applied to the order and speed believes. The following figure takes the form of stator flux following the axis alpha. One finds well the phase of fluxing as well as the evolution of flow according to the increase number of revolutions of the machine and thus of the frequency of the supply voltages. The following figure presents the evolution of the produced torque by the induction machine. One finds the phase of acceleration (with a going beyond related to the speed control) and the phase of stabilization speed and thus couple. The following figure takes the form the speed of the system. One finds well there the various phases, which one describes previously. The following figure takes the form of rotor resistance at exit of the filter of Kalman (for a healthy system). One notes a drift of rotor resistance according to speed. This derivative was corrected: according to speed, the estimated parameter is compensated for so that it is closer to the expected value. This compensation has the form of a linear relation between the drift of rotor resistance and speed.

B. Defect with the rotor

Rotor faults Figure 8 takes the form of the stator current for an unbalanced defect (40% of whorls removed in the phase at fault), note the appearance of the lines in the spectrum of the current. The defect with the rotor is characterized by the appearance of two frequencies (1±2g) or g represents the slip (92,2Hz et107Hz) around the electric frequency f electric (100Hz). If the FFT of this signal is made, this defect will be very difficult to detect in the spectrum of the current because the amplitude of the line at the fundamental frequency is very important compared to that of the required line and these lines are at a frequency frame the line with fundamental with a difference related to the slip (in general very weak). The following figure takes the form of resistance estimated for a defect unbalanced to the rotor (40% of whorls removed at fault on a phase). Noting however that the estimate of rotor resistance is assigned to the moment t=3 second (time of the application of defect) from where appearance of the low frequency oscillations. This existence of the frequency allows then on the basis of resistance rotor estimated characterizing that it is well a defect with the rotor Fig. 9 Estimation resistance rotor for an asymmetrical fault witch the rotor C. Defect with the Stator Figure 10 takes the form of the stator current for a defect unbalanced to the stator (40% of whorls removed at fault on a phase). Let us note the appearance of the frequency, which characterizes the defect, which is 3felectric (300Hz). The appearance of the lines in the spectrum of the current is due to the application of the defect to the moment t= 3s, it is easy to detect the frequency, which characterizes the defect in low frequency. (See figure below after the stator phase current FFT). The following figure presents the current of the phase at exit of the machine and the current for the same phase at exit of the observer, it should be specified that the two currents (measured and estimated) are similar. The observer thus modified the rotor size resistance so that the two currents are -identical‖. The error, which exists between the two currents, explains by the fact why the observer is a discrete system: he thus does the calculations of the sizes on the basis of measure taken at previous time. This explains the light error of hauling. Figure 12 illustrates the application of this estimator on the exit of the machine leads to the convergence of resistance up to 3 second (application of the defect). It is noted whereas the observer does not manage to make converge the resistance estimated following the application of the defect. This is explained because the model used for the equations of the observer is a balanced model. However the defect results in an unbalanced behavior of the machine. The observer has of another choice to only make evolve dynamically rotor resistance so that the other states are coherent compared to measurements. The following figure has the amplitude of the line associates with the frequency of defect (2fe). The characteristic is well detected by this monitoring system. The defect is applied to 3 seconds.

D. Inverter fault

Figure 16 takes the form of the stator current of the machine at fault of inverter. We chose to give an account of the effects of an open component. For that, we included a block allowing to block the control commands arriving at the switch considered as failing. After the period of fluxing (0s-0.3s), one notices that the size is sinusoidal until the application of the defect to the moment t=3seconde. The occurrence of the defect thus leads to the appearance of the frequencies, which characterize this defect in low and high frequency (see figure 17 which presents the FFT of the stator current). The study presented in this paper is a contribution to the diagnosis of the defects of electric machinery, with special focus on squirrel-cage induction machine. Indeed, induction machines are appropriate for applications where high power and speed are required. This is due, to multiple advantages (simplicity of construction, robustness, performances…).However, for certain applications requiring high degree of reliability, it is necessary to carry out a fast and early detection of possible failures. The presented study focus on the contribution of the analytical in this article exposes the contribution of the analytical method containing models of the system to the diagnosis of the electric defects of the asynchronous machines. The vectorial control of the induction machine requires observers (possibly extended). We first of all showed that the use of an extended observer for diagnosis purpose made it possible to detect the principal defects in the electrical machine. The parameter used for the diagnosis is the estimated rotor resistance. We showed that by coupling a spectral analysis of the estimated parameter that the technique could appear useful. The Kalman filter uses dynamics of the rotor resistance, which defines the evolution its value in time domain to obtain reliable data in the case of healthy machine (convergence of the estimate). We are aware of having studied only certain aspects of a very vast subject of study, which requires the comprehension of the system and the use of suitable tools. The prospects with this work are multiple, from the point of view of the monitoring, the use of other parameters (stator resistance, stator inductance) to evaluate the performance of the wide observer. Then, from the point of view the detection of the defects, it will be interesting to study the monitoring of this kind of machine used as a generator of power (of wind turbines type for example) connected to the electrical. APPENDIX Rs = 0.0142Ω; Ls = 0.001314 H; Rr = 0.0126Ω; Lr = 0.001286H; Lm = 0.00123H; p = 3; σ = 0.1047; Q = 0, 1; R= 0, 01

Fig. 1 .

 1 Fig.1. Schematic diagram of the system under study

  the number of whorls of the phase a has stator.

2 /

 2 Computation of the a posteriori covariance matrix P(k+1/k)

Fig. 2 .

 2 Fig. 2.Stator current

Fig. 3 .

 3 Fig. 3.Stator flux

Fig. 4 .

 4 Fig.4. Machine couples

Fig. 5 .

 5 Fig. 5.Speed of the machine

Fig. 6 .Fig. 7 .

 67 Fig. 6.Estimate of rotor resistance by the filter of Kalman without compensation

Fig. 8 .

 8 Fig. 8. Stator current under rotor fault

Fig. 10 .

 10 Fig. 10. Stator current of the machine at fault to the stator

Fig. 12 .

 12 Fig.12. Phase current (Isa), Observed current (Isakal)

Fig. 13 .

 13 Fig.13.Estimation rotor resistance for an asymmetrical fault witch the stator 1) Synchronous Demodulation Applied to Rotor Resistance: Each defect has a frequency which characterizes it, the idea of this method consists with research the amplitude of the line associated with this frequency to detect the evolution of this one. The principle of this demodulation is exposed to the following figure.

Fig. 14 .

 14 Fig.14. Principe of the demodulation

Fig. 15 .

 15 Fig.15. Compensated resistance rotor demodulated

Figure. 16 .Fig. 17 .

 1617 Figure.16.Stator current of the machine at fault of inverter


Francesco Capuano

Donato Vallefuoco

Effects of discrete energy and helicity conservation in numerical simulations of helical turbulence

Helicity is the scalar product between velocity and vorticity and, just like energy, its integral is an inviscid invariant of the three-dimensional incompressible Navier-Stokes equations. However, space-and time-discretization methods typically corrupt this property, leading to violation of the inviscid conservation principles. This work investigates the discrete helicity conservation properties of spectral and finitedifferencing methods, in relation to the form employed for the convective term. Effects due to Runge-Kutta time-advancement schemes are also taken into consideration in the analysis. The theoretical results are proved against inviscid numerical simulations, while a scale-dependent analysis of energy, helicity and their non-linear transfers is performed to further characterize the discretization errors of the different forms in forced helical turbulence simulations.

Introduction

This work is concerned with the incompressible Navier-Stokes (NS) equations, which can be written as follows

∇ • #» u = 0 , (1)
∂ #» u ∂t + # » N ( #» u ) = -∇p + ν∇ 2 #» u , (2)
where # » N ( #» u ) is the non-linear convective term and ν is the kinematic viscosity. In the inviscid case the NS equations are replaced by the Euler equations, which possess linear as well as quadratic invariants, i.e., linear or quadratic functions of the velocity vector whose integral over the domain remains constant in time throughout the flow evolution. The identification of the invariants becomes particularly straightforward upon recognizing that the inviscid NS equations are a Hamiltonian system and can be thus written in the canonical skew-gradient formulation

∂ #» y ∂t = D ( #» y ) ∇H , (3)
where D is a skew-symmetric operator and H is the Hamiltonian functional [START_REF] Olver | A nonlinear hamiltonian structure for the euler equations[END_REF]. Recasting the inviscid version of Eqs. ( 1)-( 2) in the skew-gradient formulation of Eq. ( 3) involves either eliminating the pressure (upon projection of the velocity field onto a divergence-free space) [START_REF] Arnold | Hamiltonian nature of the euler equations in the dynamics of a rigid body and of an ideal fluid[END_REF], or taking the curl of Eq. (2) [START_REF] Olver | A nonlinear hamiltonian structure for the euler equations[END_REF]. Once a system is known to be Hamiltonian, the Noether theorem can be applied to yield an association between the symmetries of the Euler system and conserved quantities; for instance, conservation of linear momentum corresponds to the invariance of the Euler equations under space-time translations [START_REF] Olver | On the hamiltonian structure of evolution equations[END_REF]. Here we are particularly interested in quadratic invariants. One follows directly from Eq. ( 3), since the Hamiltonian functional is seen to be the mean kinetic energy of the flow,

K e = 1 2 #» u • #» u , (4)
where the angular brackets denote the average operator over the entire domain. Another inviscid integral arises from the nontrivial kernel of the operator D, and corresponds to the quantity named helicity, whose average value over the domain reads

K h = #» u • #» ω , (5)
where #» ω = ∇ × #» u is the vorticity vector. Helicity has thus been termed a Casimir invariant of the Euler equations [START_REF] Holm | Geometric mechanics: Part I: Dynamics and symmetry[END_REF].

The discovery of helicity conservation is relatively recent and dates back to 1961 [START_REF] Moreau | Constantes dun ilot tourbillonnaire en fluide parfait barotrope[END_REF]. Since then, helicity has been found to play an important role in laminar and turbulent flows, and has raised significant interest both in terms of fundamental understading and repercussions in engineering applications and flows of practical interest [START_REF] Moffatt | Helicity in laminar and turbulent flow[END_REF]. Helicity is known to inhibit the transfer of energy towards smaller scales, since the statistical alignment of velocity and vorticity leads to partial suppression of the nonlinear term [START_REF] André | Influence of helicity on the evolution of isotropic turbulence at high reynolds number[END_REF][START_REF] Kraichnan | Depression of nonlinearity in decaying isotropic turbulence[END_REF]. On the other hand, the self-similar energy decay rate has been shown to be the same as in the nonhelical case, unless rotation is present [START_REF] Teitelbaum | Effect of helicity and rotation on the free decay of turbulent flows[END_REF]. The dynamics of the joint cascade of energy and helicity constitutes an interesting field of research which is very active nowadays [START_REF] Chen | The joint cascade of energy and helicity in three-dimensional turbulence[END_REF][START_REF] Biferale | Split energy-helicity cascades in three-dimensional homogeneous and isotropic turbulence[END_REF][START_REF] Kessar | Non-kolmogorov cascade of helicity-driven turbulence[END_REF]. Helicity is also important in atmospheric and geophysical flows [START_REF] Lautenschlager | Subgrid parametrization in helical flows[END_REF], while an helicity-based index is used in biomedical research to quantify swirling motions in cardiovascular flows [START_REF] Morbiducci | Helical flow as fluid dynamic signature for atherogenesis risk in aortocoronary bypass. a numeric study[END_REF].

Upon discretization of the NS equations in space and time, the invariant character of energy and helicity is generally lost. A significant research effort has been carried out over the last years to develop numerical algorithms that preserve invariants also in a discrete sense, with the aim of obtaining stable computations and physically relevant solutions. Particular attention has been paid to the development of energy-preserving numerical methods, that have ultimately allowed stable long-time integrations and realistic representations of the energy cascade [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flows[END_REF][START_REF] Verstappen | Symmetry-preserving discretization of turbulent flow[END_REF][START_REF] Capuano | Energy preserving turbulent simulations at a reduced computational cost[END_REF]. On the other side, invariance of helicity has been very seldom considered in the derivation of numerical methods. Notable exceptions include the works by Liu and Wang [START_REF] Liu | Energy and helicity preserving schemes for hydro-and magnetohydrodynamics flows with symmetry[END_REF], for axisymmetric flows, and by Rebholz and coworkers (see, e.g., [START_REF] Rebholz | An energy-and helicity-conserving finite element scheme for the navier-stokes equations[END_REF]), in the framework of finite-element methods.

The aim of the present work is to investigate the discrete helicity-conservation properties of commonly used semi-discrete algorithms based on spectral, finite-difference or finite-volume methods and Runge-Kutta time-advancement schemes. Particularly, the discrete behaviour of the different formulations of the nonlinear term # » N is discussed. The effects of discretization in time are also taken into consideration in the analysis.

The paper is organized as follows. Section 2 deals with the energy and helicity balance equations in a continuous setting, along with the rules of calculus required to derive them. The spatial and temporal discretization of the Navier-Stokes equations are presented in Section 3, in which a fully discrete equation for the helicity evolution is also derived. In Section 4 the numerical results are reported and discussed. Concluding remarks are given in Section 5.

Conservation of energy and helicity

In this section, we derive the conservation laws for the two quadratic invariants of the Euler equations. Clearly, in the viscous case, these quantities are actually not conserved but are subject to a source term which is either strictly disspative (for the energy) or has no definite sign (for the helicity).

For the sake of clarity, we will assume that periodic boundary conditions apply, but this does not come at a loss of generality. In deriving the conservation laws, it will be useful to list a number of identities, which can be readily proven using the standard rules of calculus

∇ • ( #» u #» u ) = #» u • ∇ #» u if ∇ • #» u = 0 , (6)
∇u = #» 0 , (7)
(∇u) v = -u (∇v) , (8)
#» u • (∇ × #» v ) = #» v • (∇ × #» u ) . (9)
All the functions appearing above are supposed to be periodic. The energy evolution equation can be derived by scalarly multiplying Eq. ( 2) by #» u and taking the average over the domain, yielding

dK e dt + #» u • # » N = - #» u • ∇p + ν #» u • ∇ 2 #» u . ( 10
)
The contribution from the pressure term vanishes due to use of the integration-by-parts (8) together with the incompressibility condition [START_REF] André | Influence of helicity on the evolution of isotropic turbulence at high reynolds number[END_REF]. Similarly, for helicity one also has to take the curl of the momentum equation ( 2) and then scalarly multiply the resulting vorticity equation by #» u , leading to

1 2 dK h dt + #» u • ∇ × # » N = ν #» u • ∇ 2 #» ω . ( 11
)
The contributions from the nonlinear term in both Eq. ( 10) and Eq. ( 11) also vanish, but the properties to be invoked to prove their cancellation depend on the formulation adopted for the nonlinear term. This can be indeed expressed in several, analytically equivalent forms:

(Adv.) # » N = #» u • ∇ #» u , (12)
(Div.) # » N = ∇ • ( #» u #» u ) , (13)
(Skew.) # » N = 1 2 ( #» u • ∇ #» u + ∇ • ( #» u #» u )) , (14)
(Rot.) # » N = #» ω × #» u , (15)
which are named the advective, divergence, skew-symmetric and rotational forms respectively [START_REF] Capuano | Energy preserving turbulent simulations at a reduced computational cost[END_REF]. Switching from one form to another requires use of the product rule ( 6) and of Eq. ( 1). Note that when the rotational form is used, the kinematic pressure is subsituted by the dynamic pressure

P = p+| #» u | 2 /2.
The properties required to prove cancellation of the nonlinear term contribution in the energy and helicity equations are reported in Table 1. Note that while the advective and divergence forms always require the use of the product rule, only the integration-by-parts property is needed to prove energy conservation for the skew-symmetric form. However, the skew-symmetric form still necessitates the product rule to verify helicity invariance. On the other hand, the rotational form is automatically energy conserving (due to orthogonality between the velocity and its cross product with vorticity), and requires only the symmetry of the curl operator to demonstrate conservation of helicity.

The energy and helicity equations can be then integrated in time between two generic instants t and t + ∆t, to yield the finite time-variation of the two invariants

∆K e ∆t = ν #» u • ∇ 2 #» u ∆t , (16)
∆K h ∆t = 2ν #» u • ∇ 2 #» ω ∆t , (17)
where • ∆t is the time-average operator over the interval ∆t. In deriving Eqs. ( 16)-( 17), which will be useful in the following, we have implicitly employed the integration-by-parts rule for the time-derivative operator. Indeed, given two differential equations ∂ t f (t) = F and ∂ t g(t) = G, where f (t) and g(t) are two continuous and differentiable functions of time, and ∂ t is a continuous time-derivative operator, the integration-by-parts rule gives

∆( f • g) ∆t = f • G ∆t + g • F ∆t , (18)
which is needed to prove Eqs. ( 16)-( 17) upon substituting f and g with #» u , and with #» u and #» ω respectively. The properties ( 6)-( 9) and ( 18) do not necessarily hold on a discrete level when the continuous space and time operators are substituted by their finite-differencing counterparts, as will be outlined in the following section.

Fully discrete energy and helicity evolution

The ultimate aim of this section is to derive a fully discrete analogue of Eqs. ( 16)- [START_REF] Kerr | Higher-order derivative correlations and the alignment of small-scale structures in isotropic numerical turbulence[END_REF]. In doing this, we also wish to understand which of the properties listed in Table 1, as well as Eq. ( 18) are retained by the discretization procedure. To this end, we employ a semi-discrete approach (also called method of lines) which is routinely used to treat the incompressible Navier-Stokes equations numerically.

Spatial discretization

The first step of the algorithm is to derive a spatially discretized version of Eqs. ( 1)-( 2), which can be expressed as follows

Mu = 0 , ( 19
) du dt + C(u)u = -Gp + νLu , ( 20
)
where u is a discrete vector containing the components of velocity on the three-dimensional mesh, u = [u x u y u z ] T , the matrices G ∈ R N u ×N p and M ∈ R N p ×N u are the discrete gradient and divergence operators, respectively, while L ∈ R N u ×N u is the block-diagonal Laplacian, diag(L, L, L). N u and N p are the number of unknowns on the mesh for velocity and pressure respectively. The spatial discretization reported in Eqs. ( 19)-( 20) possibly encompasses finite-difference, finite-volume, pseudo-spectral as well as finite-element approaches, although the details of the various operators depend on the specific method employed. In what follows, it is assumed that the discretization is built upon a regular mesh with uniform step size, and that gradient and divergence operators are discretized consistently, in such a way that the relation G T = -M holds. The uniform mesh assumption is made here only for the sake of clarity; extension to more general cases is discussed in Section 3.3. We are particularly interested in the convective term, which can be expressed as the product of a linear block-diagonal convective operator C(u) and u:

C(u)u =   C x (u) C y (u) C z (u)     u x u y u z   . (21)
The operator C(u) is obtained by discretizing the nonlinear term starting from one of the possible expressions in which this can be written. As mentioned in Section 2, while all these expressions are equivalent in a continuous framework, this is not the case for the discrete setting. Also, the operator C assumes a specific form which depends on the layout of the variables on the mesh. Typically, the Navier-Stokes equations are either discretized in a colocated arrangement, i.e., with all the variables located at the same points of the mesh, or on a staggered grid, in which the velocity components and the pressures are stored in different locations.

In the subsequent sections we will first analyze the colocated layout of the variables, and then discuss the staggered one.

Colocated layout

In a colocated layout, the operators C x , C y and C z are identical and can assume one of the following forms:

(Adv.) C (u) = (U x D x + U y D y + U z D z ) ≡ A , (22)
(Div.)

C (u) = (D x U x + D y U y + D z U z ) ≡ D , (23)
(Skew.) C (u) = 1 2 (A + D) , (24)
while the rotational form can be expressed as

(Rot.) C(u) = V (Ru) = -V (u) R . (25)
In the expressions above, the matrices D (•) represent the discrete derivative operators along each direction and acting on the whole set of variables on the mesh, while U (•) are the diagonal matrices of the discretized velocity components along the three directions (e.g., U x = diag(u x )). Also, V is a skew-symmetric matrix performing pointwise vector product, while R is the curl operator, which can be expressed as

R =   -D z D y D z -D x -D y D x   . (26)
Generally speaking, the introduced discretization yields truncation errors, due to the finite approximations of derivatives, as well as aliasing-type errors, due to the evaluation of the nonlinear terms in a finite-dimensional space. These errors are likely to produce deviations from the rules of calculus and therefore violation of the conservation of quadratic invariants. As a first step, it is useful to check how the discrete operators behave with respect to the properties listed in Table 1. To this end, we introduce a discrete inner norm, which under the hypothesis of regular mesh simply reads

uv → u T v. (27)
The discrete integration-by-parts rule is easily verified for central-differencing formulas, including compact schemes and Fourier differentiation, which all lead to skew-symmetric derivative matrices. Indeed, in such cases, assuming two discrete functions u and v on a periodic domain, and a generic derivative operator D, one has

u T Dv = v T D T u = -v T Du , (28)
which is the discrete analogue of Eq. ( 8). It is worth to note that Eq. ( 28) holds regardless of the scheme order (i.e., the truncation error), as well as of the presence of aliasing errors. This was already shown by Kravchenko and Moin [START_REF] Kravhcenko | On the effect of numerical errors in large eddy simulations of turbulent flows[END_REF] using Fourier expansion. The divergence property is also easily proven assuming that the discrete derivative operators behave consistently with respect to a constant function.

In this case, one has simply

1 T Du = -u T D1 = 0 , (29)
Table 2. Requirements for a spatial discretization on a colocated grid to satisfy the rules of calculus.

Eq. Property Requirement (6) Product rule Fully de-aliased spectral differentiation (7) Divergence Consistent scheme (D1 = 0) (8) Integration by parts Central-differencing scheme (incl. explicit, compact, spectral) (9) Symmetric curl where 1 is the unit vector. As a consequence of using central-differencing schemes, the curl operator reported in Eq. ( 26) also turns out to be symmetric as in the continuous case, Eq. ( 9), being a skewsymmetric block matrix constituted by skew-symmetric blocks. On the other hand, it is straightforward to verify that the product rule is in general violated by the discrete approximation. Although this property is restored in case of dealiased spectral schemes, for finite-differencing methods of any order the product rule does not hold on a discrete level, even if aliasing errors are removed [START_REF] Kravhcenko | On the effect of numerical errors in large eddy simulations of turbulent flows[END_REF]. A summary of the properties is reported in Table 2.

In light of the above considerations, and upon comparison with Table 1, one could already infer the consequences of using the various formulations of the nonlinear term. Further insight, however, is gained by deriving the semi-discrete, time-continuous energy and helicity equations. The discrete energy is defined as e = u T u/2, whereas the discrete helicity reads h = u T w, where w = Ru is the discrete vorticity vector. The energy equation has been widely reported elsewhere [START_REF] Capuano | Energy preserving turbulent simulations at a reduced computational cost[END_REF][START_REF] Verstappen | Symmetry-preserving discretization of turbulent flow[END_REF] and reads

de dt = -u T C(u)u -u T Gp + νu T Lu . (30)
While the pressure term vanishes as a consequence of the continuity equation (assuming Mu = 0 and the consistency relation G T = -M), the convective term contribution depends basically on the formulation employed for it. It is easy to show that Eq. ( 24) yields a skew-symmetric convective operator provided that the schemes satisfy the discrete integration-by-parts rule, while the expression ( 25) is skew-symmetric a priori. On the other hand, this does not happen for none of the operators given by Eq. ( 22) and Eq. ( 23). As a consequence, only the skew-symmetric and the rotational forms of the convective term preserve energy spatially in the inviscid case. This circumstance was soon recognized [START_REF] Zang | On the rotation and skew-symmetric forms for incompressible flow simulations[END_REF], and since then there has been a long debate regarding the use of these two forms, which allow stable long-time integration [START_REF] Blaisdell | The effect of the formulation of nonlinear terms on aliasing errors in spectral methods[END_REF][START_REF] Horiuti | Comparison of conservative and rotational forms in large eddy simulation of turbulent channel flow[END_REF]. While the rotational form is much more cost-effective than the skew-symmetric one, it is generally agreed that its aliasing errors are higher than those provided by the skew-symmetric form [START_REF] Zang | On the rotation and skew-symmetric forms for incompressible flow simulations[END_REF]. Indeed, aliased simulations carried out with the skew-symmetric form proved to be very similar to the de-aliased ones [START_REF] Kravhcenko | On the effect of numerical errors in large eddy simulations of turbulent flows[END_REF].

Here, we are especially interested in how the helicity balance is affected by the two energyconserving forms. A relation analogous to Eq. ( 30) can be derived by scalarly multiplying Eq. ( 20) by w, then by multiplying the curl of Eq. ( 20) by u, summing the two equations and using the symmetry of R. The result is

1 2 dh dt = -u T RC(u)u - u T RGu + νu T RLu . ( 31
)
The pressure term cancels as long as RG = 0, while the convective contribution depends again on the formulation employed for the convective term. In particular, the rotational form preserves helicity, provided that the curl operator is symmetric. Indeed, upon substitution of the convective operator, C(u) = -V(u)R, the nonlinear contribution reads

u T RV(u)Ru = 0 , (32)
which cancels due to the skew-symmetry of V and the symmetry of R. On the other hand, the skewsymmetric, as well as the advective and divergence forms do not yield cancellation of the nonlinear term contribution in the helicity equation. The above conclusions are in line with those that can be inferred by comparing Table 1 and Table 2. Indeed, the skew-symmetric form requires the product rule to prove conservation of helicity, whereas for the rotational form only the symmetry of the curl operator is needed; the former property is not satisfied discretely (unless aliasing and truncation errors are removed), while the latter holds for any central-differencing method. Although we do not attempt a formal proof, the fact that the discrete skew-symmetric form does not yield conservation of helicity, despite providing a skew-symmetric convective operator, can be intutively linked to the nature of the helicity invariance. As mentioned in the Introduction, this is in fact not related to the skew-gradient nature of the Euler equations, but rather to the degeneracies of the Hamiltonian operator. In other words, both the skew-symmetric and the rotational forms yield a semi-discrete system that can be recast as Eq. ( 3), but only one of them also possesses helicity as an inviscid invariant.

Staggered layout

The colocated arrangement is simple and straightforward to implement, although its use is often complicated by the well-known odd-even decoupling phenomenon, which leads to checkerboard patterns for the pressure field [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF]. Nevertheless, the colocated layout is still employed in many modern codes and special techniques have been developed to avoid the pressure-velocity decoupling [START_REF] Rhie | Numerical study of the turbulent flow past an airfoil with trailing edge separation[END_REF][START_REF] Shashank | A co-located incompressible navier-stokes solver with exact mass, momentum and kinetic energy conservation in the inviscid limit[END_REF][START_REF] Trias | Symmetry-preserving discretization of navier-stokes equations on collocated unstructured grids[END_REF]. On the other hand, a common remedy is represented by the staggered layout, in which the pressure and the velocity components are stored in different locations of the computational cell. The most popular staggered arrangement was proposed by Harlow and Welch [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF], and is outlined in Fig. 1.

A formal matrix analysis, similar to the one performed for the colocated case, becomes particularly involved for the staggered layout and is beyond the scope of this paper. Here, we focus on the scheme proposed by Harlow and Welch (HW) [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF], and test its properties numerically in Section 4. The HW scheme is a second-order method which relies on the divergence form of the convective term, as well as on proper interpolation operators; using the same notation as Morinishi et al. [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flows[END_REF], the nonlinear term is computed as

(Div.) i ≡ δ 1 u 1x i j u 1x j i δ 1 x j , (33)
where δ 1 and • 1x are the half-grid derivative and interpolation operators respectively.

The HW scheme has several interesting properties: it preserves momentum a priori (i.e., regardless of mass conservation) as well as kinetic energy in the inviscid case. Indeed, it can be shown that the discretization in Eq. ( 33) leads to a skew-symmetric convective operator C(u), provided that Mu = 0 [START_REF] Coppola | Energy-preserving discretizations of the Navier-Stokes equations. Classical and modern approaches[END_REF]. On the other hand, to the authors' knowledge the helicity-conservation properties of the HW scheme have not been documented thus far in literature. Note that the definition of helicity on a staggered grid is ambiguous, due to the relative positioning of velocity and vorticity components, see Fig. 1.

Here, we consider two possibilities, namely to compute helicity in the vertices and in the center of the computational cell

h| vertices = u 1x 2 1 1x 3 ω 1 1x 1 + u 1 1x 1 1x 3 ω 2 1x 2 + u 3 1x 1 1x 2 ω 3 1x 3 , (34)
h| centers = u 1 1x 1 ω 1 1x 2 1x 3 + u 2 1x 2 ω 2 1x 1 1x 3 + u 3 1x 3 ω 3 1x 1 1x 2 . (35)
We explicitly note that three interpolations have to be performed for each term to gather all the required variables in the same location. On the other hand, among the various possible definitions, a form of kinetic energy that is preserved by the HW scheme in absence of viscosity is the following [START_REF] Ham | A fully conservative second-order finite difference scheme for incompressible flow on nonuniform grids[END_REF]

e = 1 2 u i u i 1x i , (36)
and is located at the cell center.

Temporal discretization

Once the Navier-Stokes equations have been spatially discretized, one is left with Eqs. ( 19)-( 20), for which a further discretization in time is required. The semi-discrete equations constitute an index-2 Differential Algebraic system [START_REF] Sanderse | Accuracy analysis of explicit Runge-Kutta methods applied to the incompressible Navier-Stokes equations[END_REF], due to the presence of the pressure, which acts as a kinematic constraint to ensure incompressibility. A system of ODE is formally obtained upon introducing a projection operator P,

du dt = f(u)u , (37)
where f = Pf and f (u) = -C (u) + νL, with P = I -GL -1 M and L = MG. Time advancement of Eq. ( 37) is now straightforward and can be carried out by means of any ODE solver. In general, the discretization in time will in turn spuriously contribute to the energy and helicity balance, unless proper care is taken. Here we will focus on the class of Runge-Kutta (RK) time-advancement methods. RK schemes are very popular in the fluid dynamics community due to their favorable properties, such as their self-starting capability and relatively large stability limit.

A general s-stage Runge-Kutta method applied to Eq. ( 37) can be expressed as

u n+1 = u n + ∆t s ∑ i=1 b i f(u i )u i , (38)
u i = u n + ∆t s ∑ j=1 a i j f(u j )u j , (39)
where a i j and b i are the RK coefficients. The RK coefficients are often arranged into the so-called Butcher tableau [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF], and are usually constructed to maximize the temporal order of accuracy of the method p; for s ≤ 4, schemes with p = s can be obtained. Note that the role of the projection operator is to enforce incompressibility at each time instant and at each sub-step through the solution of a Poisson equation for pressure. This is often done in practice through a fractional-step method [START_REF] Capuano | Approximate projection method for the incompressible navier-stokes equations[END_REF].

Similarly to the spatial case, we first proceed to check whether the integration-by-parts rule for the time-derivative operator, Eq. ( 18), is retained by the RK method. The discrete analogous of Eq. ( 18) reads

∆( f • ĝ) ∆t = f • Ĝ ∆t + ĝ • F ∆t , (40)
where • indicates the time-discrete variable. In Runge-Kutta methods, the inner norm is properly defined as

f G ∆t → s ∑ i=1 b i f i G i , (41)
which is reminiscent of the quadrature rule at the base of RK methods. Upon simple manipulations, one has

∆( f • ĝ) ∆t = f • Ĝ ∆t + ĝ • F ∆t -∆t ∑ i, j (b i a i j + b j a ji -b i b j ) Fi Ĝ j . (42)
In practice, upon comparison of Eq. ( 42) with Eq. ( 40), it is evident that RK schemes do not satisfy a discrete integration-by-parts rule unless the identity

g i j ≡ b i a i j + b j a ji -b i b j = 0 ( 43
)
is verified for every i, j. The condition ( 43) is very well known in the community of so-called geometric integration methods, and is the requirement for a RK method to preserve all the quadratic invariants of a system [START_REF] Hairer | Geometric Numerical Integration[END_REF]. Also, for irreducible RK methods, Eq. ( 43) is a necessary and sufficient condition for the method to be symplectic. Schemes satisfying Eq. ( 43) are necessarily implicit.

The derivation of fully-discrete balance equations for energy and helicity is now straightforward. Again, the expression for the kinetic energy variation introduced by Eqs. ( 38)-( 39) has been derived elsewhere [START_REF] Capuano | Energy preserving turbulent simulations at a reduced computational cost[END_REF][START_REF] Sanz-Serna | Runge-Kutta schemes for Hamiltonian systems[END_REF][START_REF] Sanderse | Energy-conserving Runge-Kutta methods for the incompressible Navier-Stokes equations[END_REF] and is reported here for convenience,

∆e ∆t = ν s ∑ i b i u T i Lu i - ∆t 2 s ∑ i, j g i j u T i f T (u i ) f (u j ) u j , (44)
where it has been assumed that the spatial discretization is energy conserving. We focus hereinafter on the fully discrete evolution equation for helicity, which reads

∆h ∆t = 2ν s ∑ i b i u T i RLu i -2 s ∑ i b i u T i RC (u i ) u i -∆t s ∑ i, j g i j f T i R f j . (45)
Note that the equations of both energy and helicity resemble the generic structure given in Eq. ( 42).

Particularly, the three terms appearing in the right-hand side of Eq. ( 45) are respectively the contribution of the discretized (physical) viscous dissipation, the spatial error term due to convection, and the temporal error. The spatial error behaviour follows from the results already obtained in Section 3.1, e.g., it vanishes for the rotational form on a colocated layout. The temporal error stems from the possible lack of the integration-by-parts rule. Therefore, only symplectic RK methods are able to preserve energy and helicity in time in case of inviscid flow. This is not surprising, since these two quantities are inviscid invariants of the Navier-Stokes equations, and, as mentioned above, RK methods satisfying g i j = 0 preserve all the quadratic invariants of a system. When energy-and helicity-preserving methods are employed both in space and time, Eqs. ( 44)-(45) become the discrete counterparts of Eqs. ( 16)-( 17), i.e., the induced balance equations of energy and helicity are correctly enforced, regardless of scheme order, grid size and time step [START_REF] Sanderse | Energy-conserving Runge-Kutta methods for the incompressible Navier-Stokes equations[END_REF]. On the other hand, the discrete dissipation rate is obviously still affected by spatial and temporal discretization errors, when compared to its exact counterpart.

As far as time-integration is concerned, it is worth to note that recently Capuano et al. developed special explicit Runge-Kutta methods for the fluid flow equations, named pseudo-symplectic, with the aim of preserving kinetic energy in time up to an order q, with q > p. Nonetheless, these schemes have actually order q for the conservation of any quadratic invariant, and are therefore applicable also for the purpose of preserving helicity (or enforcing the correct helicity balance) in time to a higher order of accuracy. Expl.

• + • × + × 3 Skew Symp. • + • + × + 4 Skew Expl. • + • × × × 5 Div. Expl. + + × × × × 6 HW Expl. + + • × × ×
By collecting the outcome of the preceding sections with known literature results [START_REF] Morinishi | Fully conservative higher order finite difference schemes for incompressible flows[END_REF][START_REF] Capuano | Energy preserving turbulent simulations at a reduced computational cost[END_REF], the conservation properties (also including mean momentum and energy) of the Navier-Stokes algorithms analyzed in the paper are reported in Table 3. Any explicit RK method produces an error in the conservation of quadratic invariants, although the mean momentum is preserved, since all Runge-Kutta methods preserve linear invariants [START_REF] Rosenbaum | Conservation properties of numerical integration methods for systems of ordinary differential equations[END_REF]. Only the rotational and the skew-symmetric forms on a colocated layout, and the staggered HW scheme are able to preserve kinetic energy, although the condition Mu = 0 (as well as the usual consistency relation G T = -M) must hold. Remarkably, both the convective and the pressure terms in the helicity equation vanish regardless of the discrete enforcement of continuity, and therefore a discretization based on the rotational form preserves helicity a priori (assuming RG = 0). Note that the rotational form, coupled to a symplectic integrator, is able to preserve mean momentum, energy and helicity in the inviscid case.

Extension to nonuniform meshes

The hypothesis of uniform Cartesian grids made in the preceding sections can be easily removed by taking into account a relevant inner product,

uv → u T Sv , ( 46
)
where S is a diagonal matrix containing the metrics of the mesh. In this case Eq. ( 20) is most properly recast as follows

S du dt + C(u)u = -SGp + νLu , (47)
and the requirements for energy and helicity conservation can be easily shown to be the same as in the uniform case, i.e., a skew-symmetric convective operator C(u) for the energy, and a symmetric curl operator R for the helicity. The proper relation between the gradient and the divergence operators in this case is G T S = -M, while RG = 0 must still hold. With a careful choice of the discrete differential operators, extension to general unstructured meshes is also possible in principle, but will not be detailed here.

Numerical results

All the numerical results presented in this section have been produced by a pseudo-spectral code in a three-dimensional periodic box. This choice has allowed very efficient computations and the possibility to analyze the numerical errors in a scale-dependent way through spectra of energy, helicity and their non-linear transfers. It was also possible to analyze the effects of truncation errors by mimicking standard non-dealiased finite-difference schemes, on both a colocated and a staggered arrangement, via the modified wavenumber approach, and to produce accompanying de-aliased spectrally resolved computations for comparison. It is worth to recall that the modified wavenumber approach consists in substituting the wavenumber vectors used for the computation of spectrally-accurate derivatives with their finite-difference counterparts [START_REF] Kravhcenko | On the effect of numerical errors in large eddy simulations of turbulent flows[END_REF]. For instance, for a second-order colocated scheme, the modified wavenumbers for the first-and second-derivative read respectively

k (k) = sin (kh) h , k 2 (k) = 2 (1 -cos (kh)) h 2 , ( 48
)
where h is the grid size. The staggered scheme has been mimicked by introducing proper shift and averaging operators to take into account the necessary interpolation procedures. For all the simulations the time step ∆t is selected based on a CFL-like condition with CFL=0.5.

Inviscid helical dynamics

As a first test, we have simulated an inviscid system with the aim of isolating the errors coming from the discretization of the nonlinear term and the time-advacing scheme. We have thus solved the Euler equations on a grid composed by 32 3 points, and tested the various algorithms listed in Table 3 coupled with various Runge-Kutta methods and both with second-order and spectral accuracy. The initial condition is the superposition of two Arnold-Beltrami-Childress (ABC) flows [START_REF] Childress | Stretch, twist, fold: the fast dynamo[END_REF] at the wavenumbers k 1 = 4 and k 2 = 6. The time evolutions of energy and helicity are shown in the left and right part of Fig. 2 respectively. The colocated algorithms have been used with full spectral accuracy and with the classical fourth-order Runge-Kutta scheme (RK4), except for two runs, in which the rotational form has been coupled to a second-order spatial scheme, and to both a second-order scheme and a pseudo-symplectic RK method. Particularly, the explicit six-stage scheme with p = 4 and q = 7 (named 4p7q) has been tested; the reader is referred to [START_REF] Capuano | Explicit runge-kutta schemes for incompressible flow with improved energy-conservation properties[END_REF] for details. The symplectic scheme in Algorithm #1 is the Gauss midpoint method.

The results fully confirm the theoretical predictions. The rotational form preserves helicity spatially, although it produces significant deviations in both quadratic invariants when used with an explicit timeadvancement scheme. This behaviour is attributed to the prominent accumulation of energy at smaller scales typical of spectrally truncated conservative inviscid systems, which causes the temporal error of the RK method to be significant. Remarkably, the RK error is dissipative for energy (as it is well known), but productive for helicity. A possible explanation for this circumstance might be the fact that the temporal error in the discrete helicity equation does not have a definite sign and therefore helicity can increase. When the rotational form is used in conjunction with a second-order scheme, the errors are much less pronounced, likely due to the attenuation of the high-wavenumber energy content by the finite-difference derivative method. In this regard, note that the group of second-order runs has a different initial helicity content due to truncation errors affecting the computation of initial vorticity. The beneficial impact of the higher-order conservation properties of the pseudo-symplectic RK method is evident in keeping the two invariants practically constant, as done to machine precision by the fully conservative Algorithm 1.

The skew-symmetric form is shown to dissipate the initial helicity content in few characteristic times. On the other hand, on equal spatial accuracy, the skew-symmetric form has a better behaviour on energy than the rotational form, due to its well-known favorable aliasing cancellation properties [START_REF] Kravhcenko | On the effect of numerical errors in large eddy simulations of turbulent flows[END_REF]. When coupled to a lower-order scheme (not shown here), the behaviour of the skew-symmetric form is qualitatively very similar to the one already described for the rotational formulation.

The staggered Harlow-Welch scheme is found to dissipate helicity, with a rate similar to the colocated skew-symmetric form. Actually, both the expressions proposed to compute helicity in a staggered layout, Eqs. ( 34)- [START_REF] Sanderse | Energy-conserving Runge-Kutta methods for the incompressible Navier-Stokes equations[END_REF] have been found to be practically coincident from a numerical point of view. It can thus be concluded that the classical energy-conserving second-order staggered method does not preserve helicity, at least with the definitions of vorticity and helicity given in Section 3.1. The colocated divergence form (Algorithm 5) diverged immediately due to violation of energy conservation.

The analysis reported above has allowed to confirm the theoretical results and gain preliminary insights into the behaviour of the two energy-conserving forms. With the aim of obtaining scale-dependent 3. In all cases, spectral accuracy and a RK4 scheme are used, except for two variants of the Algorithm 2, namely with II order accuracy (square) and with II order accuracy and the 4p7q RK scheme [START_REF] Capuano | Explicit runge-kutta schemes for incompressible flow with improved energy-conservation properties[END_REF] (diamond). t 0 = e -1/2 k -1 1 is the system characteristic time scale.

information, we have computed energy and helicity spectra for the solution of spherically truncated Euler equations, i.e. the system obtained by truncating the Euler equations in Fourier space at a wavenumber |k| = k max . In the absence of numerical errors, this truncated system preserves energy and helicity. In correspondence of the statistically stationary state (absolute or statistical equilibrium), exact expressions for the spectra of both energy and helicity are available [START_REF] Kraichnan | Helical turbulence and absolute equilibrium[END_REF],

E(k) = 4π α k 2 1 -β α 2 k 2 , H(k) = 8πβ α 2 k 4 1 -β α 2 k 2 , (49)
where α and β are constants that depend on the energy and helicity contents. The statistical equilibrium is governed by a single non-dimensional parameter, e.g. the relative helicity H rel = K h /(2k max K e ). In the non-helical case, i.e. H rel = β = 0, the helicity spectrum vanishes and the energy spectrum becomes simply

E(k) = 4π α k 2 . (50)
In the performed simulations k max = 42 and the initial condition is again the sum of two Arnold-Beltrami-Childress (ABC) flows with k 1 = 28 and k 2 = 30, similarly to the test proposed in [START_REF] Krstulovic | Cascades, thermalization, and eddy viscosity in helical galerkin truncated euler flows[END_REF], yielding H rel = 0.687. The same algorithms described in the previous case are tested here with exception of the case with the pseudo-symplectic method.

Energy and helicity spectra are reported in Fig. 3. The fully conservative computation matches very closely the exact solution, while Algorithm 2 (rotational form with spectral accuracy) is slightly dissipative due to the temporal dissipation of the explicit scheme. Importantly, the computation employing Algorithm 4 (skew-symmetric form with spectral accuracy, explicit) drives the system towards the nonhelical equilibrium solution, i.e. the relative helicity vanishes and E(k) is proportional to k 2 , see Eq. (50). The behaviour of Algorithm 2 when the accuracy of the spatial derivatives is reduced to second order deserves attention. Although this method preserves energy and helicity spatially and with overall good accuracy (even better than its spectral version, see Fig. 1), the energy and helicity spectra are distorted with respect to the analytical solution, and present significant errors over almost all the wavenumber range, which can be attributed to the interaction of the truncation error with the mechanism of nonlinear helicity transfer. We conjecture that this distorsion is also responsible for the aforementioned reduced energy and helicity dissipation arising from the temporal error. However, a cross-comparison with a low-order scheme based on the skew-symmetric form cannot be performed in this inviscid framework 3 and the label of Figure 1 for further details. The square symbol refers to Algorithm 2 with second-order spatial accuracy. The grey curves display the exact solution which is reported in Eq. ( 49).

since the helicity is quickly dissipated.

Forced helical turbulence

In this section, the effect of the spatial discretization errors related to the convective term is investigated through viscous helical simulations. In all cases, time-advancement is performed through a fourth-order RK scheme. We choose to force our simulations at large scales in order to achieve a statistically stationary state and then compute spectra by time averaging. The selected forcing scheme is the helical Euler forcing [START_REF] Kerr | Higher-order derivative correlations and the alignment of small-scale structures in isotropic numerical turbulence[END_REF][START_REF] Vallefuoco | Small-scale anisotropy induced by spectral forcing and by rotation in non-helical and helical turbulence[END_REF], which is inspired by the truncated Euler dynamics: the lowest-wavenumber modes, corresponding to wavevectors k such that |k| ≤ k F (k F is the largest forcing wavenumber), obey the threedimensional incompressible Euler equations and are independent of the other modes. Therefore, we set k F = 2.5 and inject both energy and positive helicity in the range k ≤ k F = 2.5. The relative helicity of the Euler system, defined as K E h /(2k F K E e ), is equal to H rel = 0.967. Energy and helicity are then transfered to higher wavenumbers down to the dissipative scales, where they are dissipated. It is worth to recall the Lin equations for energy and helicity spectra, which read

∂ E(k) ∂t = T e (k) -2νk 2 E(k) + Φ e (k) , (51)
∂ H(k) ∂t = T h (k) -2νk 2 H(k) + Φ h (k) , (52)
where Φ e (k) and Φ h (k) are linked to the external forcing and vanish in the unforced range k > k F = 2.5. The non-linear transfers T e (k) and T h (k) represent the ratio of net energy and helicity transfered to scale k through non-linear interactions among modes. For an energy-or helicity-conserving discretization, the sum of T e (k) or T h (k) over all scales vanishes, while a non conservative discretization leads to a net global transfer. Under the assumption of statistically stationary flow, and restricting the balance to the unforced scales, Eqs. (51)-(52) reduce to

T e (k) = 2νk 2 E(k) , (53) T h (k) = 2νk 2 H(k) , (54)
that is, energy and helicity transfers are balanced by viscous dissipation. Three simulations are performed in order to study the behaviour of the rotational and the skew- symmetric forms: a reference simulation, spectrally resolved and fully de-aliased (2/3 rule) with a 192 3 resolution (that is, 128 3 effective resolution); a second-order 128 3 simulation employing the rotational form; and a second-order 128 3 simulation employing the skew-symmetric form. Table 4 reports some relevant parameters of these three runs. Spectra and statistical quantities are obtained by time-averaging over 310 turnover times after the statistically stationary state has been reached. The large-scale resolution is L/(2π) = 0.11 (L is the integral scale and 2π is the periodic domain size), while the small scale resolution is k max η = 1.12 (η is the Kolmogorov scale and k max = 64 is the maximum resolved wavenumber). The Reynolds number is Re = L 2/3e/ν = 173.

The energy and helicity spectra are plotted in Fig. 4. The energy spectra of both the rotational and skew-symmetric algorithms almost superpose with the reference solution spectrum. Only at the smallest resolved scales, in the dissipative range, the discretization error in the rotational simulation causes an energy pile-up, in accordance with previous observations [START_REF] Zang | On the rotation and skew-symmetric forms for incompressible flow simulations[END_REF]. On the other hand, the departure of the skew-symmetric algorithm from the exact energy spectrum is significantly smaller, and of the same order as the slight pile-up of the reference solution. The small departure from the exact solution in terms of energy is confirmed by the values of the global contents of energy reported in Table 4, the relative error for the rotational form being twice the one yielded by the skew-symmetric.

As for the helicity spectra in Fig. 4, the differences among the three simulations are larger. Surprisingly, the rotational form seems to yield a larger departure from the exact helicity spectrum with respect to the skew-symmetric form, despite its helicity-preserving property. In particular, for the skewsymmetric form the helicity spectrum is negligibly larger than the reference one at low wavenumbers, while being significantly lower at smaller scales. On the other hand, for the rotational form, helicity is clearly too large over almost all the resolved range, becoming lower than the reference values only at the smallest resolved scales, due to the slight helicity pile-up in the de-aliased spectrally resolved simulation. These observations are confirmed by the global helicity contents reported in Table 4.

At this point, it is worth to investigate why the discretization error in the (helicity-conserving) algorithm based on the rotational form affects the helicity spectrum more than for the skew-symmetric case, which does not conserve helicity. In this regard, it is worth to preliminarily remind that although energy-and helicity-preserving methods satisfy the global balance equations to machine precision, the discretization error still affects the velocity unknowns, and thus might alter both the scale-by-scale nonlinear transfer as well as the viscous dissipation rate. In light of these clarifications, it appears from Fig. 4 that, while discrete energy conservation leads to accurate energy spectra, the helicity-preserving property is not sufficient to reproduce the correct helicty spectrum. Therefore, we proceeded to carry out a scale-dependent conservation analysis to further characterize the effect of discretization errors on helicity.

We thus compute the energy and helicity non-linear transfers T e (k) and T h (k). Table 4 reports the sums of the transfers over scales 0 ≤ k ≤ k max . For the energy and helicity transfers in the spectrally resolved dealiased algorithm and rotational algorithms, these sums are slightly negative, due to the fact that wavevectors outside the k max -radius sphere have not been included. As expected, the sum of the helicity transfers in the skew-symmetric simulation is significantly different than zero, and in particular it is negative, meaning that helicity is being numerically dissipated. The nonlinear transfers T e (k) and T h (k) are plotted in Fig. 5. For clarity, only the positive values in the unforced range are shown. As expected, for the energy transfers the error is low, the rotational form showing larger values at the largest scales (close to the forcing scale) as well as at the dissipative scales. Importantly, the skewsymmetric helicity transfer almost superposes with the reference one at larger scales, while it becomes significantly lower than the reference computation at smaller scales. The rotational simulation shows helicity transfer values larger than the exact spectrum over almost all the resolved range, and only at the dissipative scales does T h (k) become smaller than in the de-aliased spectrally truncated case.

As a result, while for the skew-symmetric form the discretization error affects mainly the small scales (causing a numerical helicity dissipation) and leaves the helicity budget basically unaltered at large scales, for the rotational form the nonlinear transfer T h (k) results to be significantly distorted by the discretization error, leading to an increased helicity cascade, although its integral over all scales vanishes. These observations explain the good behaviour of the skew-symmetric algorithm and the poor behaviour of the rotational algorithm in reproducing the exact helicity spectrum.

Conclusions

The Navier-Stokes equations possess two inviscid quadratic invariants: the mean kinetic energy #» u • #» u /2 and the mean helicity #» u • #» ω . While the former is intimately related to the conservative structure of the Euler equations, the latter comes instead from the degeneracies of the Euler operator and, like energy, plays a fundamental role in turbulent flows.

Upon discretization, the invariant character of the aforementioned quantities is generally lost due to numerical errors, and particularly due to the lack of the product rule on a discrete level and the contamination of truncation and aliasing errors. Numerical methods capable of preserving energy and helicity might be highly desirable for reliable computer simulations of turbulent flows, although the impact and interference of numerical errors have also to be taken into consideration.

We have characterized the helicity conservation properties of commonly used discretization algorithms for the Navier-Stokes equations, with particular reference to the formulation employed to express the convective term. It is found that, for a colocated layout, spatial conservation of helicity (or, in the viscous case, stric enforcement of the discrete helicity balance) is possible through the use of the rotational form. On the other hand, the energy-preserving and widely used skew-symmetric form fails to preserve helicity. As far as the time integration is concerned, implicit symplectic methods are required to satisfy conservation of both quadratic invariants to machine accuracy, while explicit pseudo-symplectic schemes can ensure higher-order conservation.

Numerical simulations have been carried out in both inviscid and viscous cases. The inviscid cases have fully confirmed the theoretical results and demonstrated that the skew-symmetric form dissipates the initial helicity content in few characteristic times. On the other hand, the rotational form was able to conserve both energy and helicity, although its behaviour was found to be strongly dependent on the spatial accuracy of the schemes (i.e., the truncation error). When used in conjunction with a secondorder method, the rotational form displayed a significant error with respect to the exact energy and helicity spectra predicted by Kraichnan at absolute equilibrium. The inviscid results have also shown that the popular staggered Harlow-Welch scheme does not preserve helicity, although the definition of this quantity is somewhat ambiguous on a staggered grid. As an additional result, it was found that explicit Runge-Kutta schemes, which are typically dissipative for energy, can instead be productive for helicity.

In the viscous case, we performed finite-difference forced simulations of helical turbulence at moderate Reynolds number. The simulations have highlighted that, despite the global conservation properties of the rotational form, the results obtained with the skew-symmetric form are overall in better agreement with the reference dealiased spectral simulation. The explanation for this counter-intuitive behaviour was obtained by performing a scale-by-scale analysis of the nonlinear energy and helicity transfers. It was found that the helicity-conservation error of the skew-symmetric form is primarily biased towards the smaller scales. In particular, a depletion of the helicity transfer is observed, which causes a helicity dissipation affecting only the small-scale part of the helicity spectrum. On the other hand, the nonlinear transfer produced by the rotational form was found to be significantly distorted, leading to helicity spectra in significant disagreement with the reference solution.

Figure 1 .

 1 Figure 1. Locations of the velocity components, vorticity components and pressure in a three-dimensional staggered layout.

Figure 2 .

 2 Figure 2. Time evolution of global energy and helicity for Euler dynamics using the various algorithms listed in Table3. In all cases, spectral accuracy and a RK4 scheme are used, except for two variants of the Algorithm 2, namely with II order accuracy (square) and with II order accuracy and the 4p7q RK scheme[START_REF] Capuano | Explicit runge-kutta schemes for incompressible flow with improved energy-conservation properties[END_REF] (diamond). t 0 = e -1/2 k -1 1 is the system characteristic time scale.

4 Figure 3 .

 43 Figure 3. Energy and helicity spectra in spherically truncated Euler dynamics for various algorithms, see Table3and the label of Figure1for further details. The square symbol refers to Algorithm 2 with second-order spatial accuracy. The grey curves display the exact solution which is reported in Eq. (49).

Figure 4 .

 4 Figure 4. Energy and helicity spectra in helical turbulence for the de-aliased spectrally resolved formulation (red), the rotational formulation (black) and the skew-symmetric formulation (blue).

Figure 5 .

 5 Figure 5. Energy and helicity transfer spectra in helical turbulence for the de-aliased spectrally resolved formulation (red), the rotational formulation (black) and the skew-symmetric formulation (blue).

Table 1 .

 1 Minimum set of properties required to prove cancellation of the nonlinear term contribution in the energy and helicity equations for the various formulations of convection in a continuous setting.

	Eq. Property	Energy Adv. Div. Skew. Rot. Adv. Div. Skew. Rot. Helicity
	(6) Product rule (7) Divergence (8) Integration by parts (9) Symmetric curl	× ×	× ×	×	× ×	× ×	× ×	×

Table 3 .

 3 Conservation properties of mean momentum, energy and helicity for 3D Navier-Stokes discretizations, for various forms of the nonlinear term and symplectic (symp.) or explicit (expl.) Runge-Kutta schemes. Algorithms 1-5 refer to a colocated grid, Algorithm 6 is the Harlow-Welch scheme described in Section 3.1.2. The discrete integration-by-parts rule is assumed to hold. + conservative a priori, • conservative if and only if continuity is discretely satisfied, × non conservative.

	Algorithms # Space Time	Momentum Space Time Space Time Space Time Energy Helicity
	1 Rot. 2 Rot.	Symp. •	+	•	+	+	+

Table 4 .

 4 Parameters of helical turbulence simulations. K E e and K E h are the energy and helicity contents of the Euler system. ε e and ε h are the energy and helicity average dissipation rates.

	Algorithms e/K E e	h/K E h	∑ T e (k)/ε e ∑ T h (k)/ε h
	Deal. Rot. Skew.	2.16 2.22 2.19	1.82 1.95 1.86	-0.012 -0.015 -0.005	-0.021 -0.0032 -0.13