Alexey Aleksandrov
email: alexey.aleksandrov@polytechnique.edu

A molecular mechanics model for flavins

Keywords: flavin, riboflavin, flavin adenine nucleotide, flavin mononucleotide, flavoprotein, CHAMM, force field, molecular dynamics

Flavin containing molecules form a group of important cofactors that assist a wide range of enzymatic reactions. Flavins use the redox-active isoalloxazine system, which is capable of one-and two-electron transfer reactions and can exist in several protonation states. In this work, molecular mechanics force field parameters compatible with the CHARMM36 all-atom additive force field were derived for biologically important flavins, including riboflavin, flavin mononucleotide and flavin adenine dinucleotide. The model was developed for important protonation and redox states of the isoalloxazine group. The partial charges were derived using the CHARMM force field parametrization strategy, where quantum mechanics water-solute interactions are used to target optimization. In addition to monohydrate energies and geometries, electrostatic potential around the compound was used to provide additional restraints during the charge optimization. Taking into account the importance of flavin-containing molecules special attention was given to the quality of bonded terms. All bonded terms, including stiff terms and torsion angle parameters, were parametrized using exhaustive potential energy surface (PES) scans. In particular, the model reproduces well the butterfly motion of isoalloxazine in the oxidized and reduced forms as predicted by quantum mechanics in gas-phase. The model quality is illustrated by simulations of four flavoproteins. Overall, the presented molecular mechanics model will be of utility to model flavin cofactors in different redox states.

INTRODUCTION

Flavoproteins participate in a wide range of biological functions and perform redox reactions and nonredox reactions, such as transferase, lyase, isomerase, and ligase reactions. [START_REF] Macheroux | Flavogenomics-a genomic and structural view of flavin-dependent proteins[END_REF] The flavoproteins have been extensively studied by experiments and represent one of the best-characterized families of proteins. [START_REF] De Colibus | New frontiers in structural flavoenzymology[END_REF] Flavoproteins have either flavin mononucleotide (FMN) or flavin adenine dinucleotide (FAD) as a cofactor, and more rarely both cofactors. [START_REF] Hefti | Deflavination and reconstitution of flavoproteins[END_REF] Apart from FAD and FMN, another biologically important flavin is riboflavin (vitamin B2) that serves as the precursor for FMN and FAD in almost all organisms. [START_REF] Fischer | Biosynthesis of vitamin B2: Structure and mechanism of riboflavin synthase[END_REF] FMN and FAD utilize in enzymatic reactions the redox-active isoalloxazine system, capable of one-electron and two electron transfer reactions, which can be coupled to the proton transfer. Flavins can also act in synergy with other cofactors such as iron-sulfur clusters, heme, molybdopterin, or thiamine diphosphate, increasing the diversity of reactions. [START_REF] Macheroux | Flavogenomics-a genomic and structural view of flavin-dependent proteins[END_REF] Time-resolved spectroscopy allowed gaining a wealth of information on the protein dynamics of flavoproteins, on time-scales directly accessible by Molecular Dynamics simulations. [START_REF] Sundstrom | Femtobiology[END_REF][START_REF] Laptenok | Ultrafast real-time visualization of active site flexibility of flavoenzyme thymidylate synthase ThyX[END_REF][START_REF] Nag | Identification of the TyrOH(*+) Radical Cation in the Flavoenzyme TrmFO[END_REF] Molecular modeling techniques based on force fields have been applied to study different aspects of flavins in solvent [START_REF] Van Den Berg | Dynamic Conformations of Flavin Adenine Dinucleotide: Simulated Molecular Dynamics of the Flavin Cofactor Related to the Time-Resolved Fluorescence Characteristics[END_REF] and in protein complexes. [START_REF] Sattelle | Calculating Chemically Accurate Redox Potentials for Engineered Flavoproteins from Classical Molecular Dynamics Free Energy Simulations[END_REF][START_REF] Schneider | A molecular dynamics simulation of the flavin mononucleotide-RNA aptamer complex[END_REF][START_REF] Dittrich | When Light Falls in LOV: A Quantum Mechanical/Molecular Mechanical Study of Photoexcitation in Phot-LOV1 of Chlamydomonas reinhardtii[END_REF][START_REF] Freddolino | Signaling mechanisms of LOV domains: new insights from molecular dynamics studies[END_REF][START_REF] Lans | The Dimer-of-Trimers Assembly Prevents Catalysis at the Transferase Site of Prokaryotic FAD Synthase[END_REF] The major requirement for such computer simulations is the existence of an adequate molecular mechanics force field model. Moreover, the underlying force field models determine the quality of the results of simulation studies relying on such models. A recent study compared the performance of different force field models including CHARMM36, Amber and OPLS in simulations of protein complexes with flavin mononucleotide. [START_REF] Robertson | Performance of Protein-Ligand Force Fields for the Flavodoxin-Flavin Mononucleotide System[END_REF] It was found that the performance of the CHARMM36 force field was similar to OPLS in reproducing experimental 3J couplings, but worse in reproducing the effect of amino residue substitutions by alanine on FMN binding. [START_REF] Robertson | Performance of Protein-Ligand Force Fields for the Flavodoxin-Flavin Mononucleotide System[END_REF] However, it should be noted that for the CHARMM simulations, atom types and force field parameters were assigned by the MATCH Web server. [START_REF] Yesselman | MATCH: an atom-typing toolset for molecular mechanics force fields[END_REF] Thus, parameters for similar compounds existing in the CGenFF force field were used, which were not necessarily optimal for flavin. Moreover, such "similar" compounds may not be available in CGenFF for the unique tricyclic heterocycle isoalloxazine structure of flavin.

While the CHARMM force field model exists for a wide range of molecules including proteins [START_REF] Mackerell | An all-atom empirical potential for molecular modelling and dynamics study of proteins[END_REF][START_REF] Mackerell | 3rd, Improved treatment of the protein backbone in empirical force fields[END_REF][START_REF] Best | Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ(1) and χ(2) dihedral angles[END_REF] and nucleotides, [START_REF] Foloppe | All-atom empirical force field for nucleic acids: I. Parameter optimization based on small molecule and condensed phase macromolecular target data[END_REF] flavins, despite their importance, have been omitted in the force field development. The present study represents a systematic development of a force field model for flavins in important protonation and redox states. Using the standard CHARMM parametrization protocol, the force field is developed to be compatible with the CHARMM36 force field for proteins and nucleic acids, [START_REF] Mackerell | An all-atom empirical potential for molecular modelling and dynamics study of proteins[END_REF][START_REF] Mackerell | 3rd, Improved treatment of the protein backbone in empirical force fields[END_REF][START_REF] Best | Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ(1) and χ(2) dihedral angles[END_REF] and CGenFF force field for small molecules. [START_REF] Vanommeslaeghe | CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields[END_REF] To summarize, the model developed in this work is suitable to investigate interactions of flavins with a wide range of flavoprotein and can be used as to parametrize other flavins, such as covalently linked to proteins. [START_REF] Mewies | Covalent attachment of flavin adenine dinucleotide (FAD) and flavin mononucleotide (FMN) to enzymes: the current state of affairs[END_REF][START_REF] Heuts | What's in a covalent bond? On the role and formation of covalently bound flavin cofactors[END_REF]

RESULTS

Set of parametrized molecules

We parametrized three redox forms of flavin shown in Figure 1. For each redox form of flavin the most important protonation states were considered at the physiological pH of 6.5. In particular, the flavin group in the reduced form can be protonated and deprotonated on the N1 atom. The associated pKa is 6.3 pH units suggesting that both forms coexist in solution at the physiological pH. [START_REF] Kyte | Mechanism in Protein Chemistry[END_REF] The pKa of the N5 nitrogen of flavin is 8.3 pH units in the radical form. [START_REF] Kyte | Mechanism in Protein Chemistry[END_REF] When bound to a protein and in the absence of a proton donor, the anionic radical form can exist for a significant period of time after the oxidized form of flavin is reduced by the transfer of one electron. [START_REF] Nag | Identification of the TyrOH(*+) Radical Cation in the Flavoenzyme TrmFO[END_REF] Thus, this protonation state of the flavin radical was also parametrized in this work. We considered the flavin group in the following important cofactors: lumiflavin, riboflavin, flavin mononucleotide (FMN), and flavin adenine dinucleotide (FAD). The flavin group common to these molecules was first parametrized in lumiflavin, and then parameters were transferred to other molecules in accord with the development of the standard CHARMM force field.

Charge optimization

The CHARMM atomic charges were derived targeting water-compound minimum interaction energies and geometries along with the dipole moment magnitude and orientation, and electrostatic potential. In what follows the charge optimization is discussed in detail only for lumiflavin in the oxidized and neutral reduced forms.

The statistics for water-compound interactions for all compounds is given in Table 1. Empirical and ab initio interaction energies and distances are given in Table S1-S5 in the Supplementary Information. The Root Mean Square (RMSD) deviation for interaction energies with the optimal set of charges is 0.34 kcal•mol -1 and 0.42 kcal•mol -1 for lumiflavin in the oxidized and the neutral reduced forms respectively. With the initial charges provided by ParamChem the RMSD is 2.70 kcal•mol -1 and 0.61 kcal•mol -1 for the oxidized and neutral reduced forms respectively. The RMS deviation for interaction distances is 0.20 and 0.17 Å with the optimal set of charges and 0.95 and 0.18 Å with the initial set for the oxidized and reduced forms of lumiflavin respectively. Overall, the optimization of charges improved water-compound interactions for all compounds. Empirical and ab initio dipole moments and electrostatic potentials are given in Table 2. To be consistent with the protocol used to derive the CHARMM force field for proteins and small molecules the dipole moment was included only for neutral compounds. The resulting CHARMM dipole moment for lumiflavin in the oxidized form of 11.48 Debye is in good agreement with the overestimated QM dipole moment of 11.77 Debye demonstrating the expected overestimation. The dipole moment is also significantly improved relative to the initial estimate of 3.85 Debye with the ParamChem charges. The orientation of the dipole moment is also improved from 105.6º to 0.7º with the initial and optimal set of charges respectively. For the reduced form of lumiflavin the CHARMM dipole moment of 4.88 Debye with the optimal set of charges is smaller than the corrected QM dipole of 6.39 Debye and is not improved in comparison with 5.00 Debye computed with the ParamChem charges.

However, the direction of the dipole moment is improved relative to the initial ParamChem dipole moment from 20.9º to 0.1º with the optimal set of charges. The dipole moment for charged compounds was not targeted in the charge fitting as in the previous force field determination works, since it depends on the system of coordinates. [START_REF] Xu | Additive CHARMM force field for naturally occurring modified ribonucleotides[END_REF] In principle, if the molecule is considered in the same system of coordinate and in the same orientation, the ab initio and empirical dipole moments can be nevertheless compared. The dipole moment for charged compounds computed with the empirical and QM models is given in Table 2. In both ab initio and empirical calculations, the molecule was positioned in the same "standard" orientation used by the Gaussian software. Though it was not explicitly included in the charge fitting, the CHARMM model reproduces well the unscaled ab initio dipole moment for the charged molecules parametrized in this work. For example, in the case of the reduced deprotonated form of lumiflavin the QM and MM dipole moment with the optimal set of charges is 12.25 Debye and 12.12 Debye respectively. However, the dipole moment of 21.90 Debye is significantly overestimated with the initial ParamChem set of charges. The direction of the dipole moment is also more consistent with the optimal charge set with an angle between the QM and MM dipoles equal to 1.1º, in contrast to 9.8º obtained with the initial charges.

Electrostatic potential was included as an additional restraint during charge fitting to provide a better charge distribution in the model compound in accord with the previous study. [START_REF] Huang | Automated Force Field Parameterization for Nonpolarizable and Polarizable Atomic Models Based on Ab Initio Target Data[END_REF] In all cases, including anionic model compounds, electrostatic potential is significantly improved relative to the initial values. For the oxidized form of lumiflavin the RMS deviation between MM and QM electrostatic potential is 3.27 and 15.87 kcal•mol - 1 •Å -1 with the optimal and initial set of charges respectively. For the neutral reduced form, RMSD of electrostatic potential is 1.00 and 3.58 kcal•mol -1 •Å -1 with the optimal and initial set of charges respectively. Targeting the QM electrostatic potential was found particularly important for charge compounds, since the number of probe water interactions was fewer than for neutral compounds, due to the dominant contribution of the net charge to water interactions. Finally, to test the dependence of the results on the QM theory employed to optimize geometries and compute the dipole moment, geometry optimization and dipole moment calculations were also performed using the MP2 functional. The results are given in Table S10 in the Supplementary Information. Compound geometries optimized using the B3LYP and MP2 functionals are very similar, with the RMS deviation lower than 0.1 Å.

Consequently, the dipole moments computed using the B3LYP and MP2 optimized structures also very similar, demonstrating that the B3LYP functional is appropriate to optimize non-bonded parameters of flavins.

Optimization of bonded terms

Stiff degrees of freedom

In this work, all bonded terms including stiff terms were parametrized based on the reproduction of energies and geometries of potential energy surface (PES) scans. We first present results for the parametrization of stiff terms. Selected results are discussed in detail only for lumiflavin in the oxidized form and in the fully reduced neutral form. The results for empirical and ab initio structures and conformation energies are summarized in Table 3. To parametrize bonded terms of lumiflavin in the oxidized form a total of 324 optimizations of PES scans were performed. The agreement between QM and MM energies of PES scans for the oxidized and reduced forms of lumiflavin is illustrated in Figure S1 in the Supplementary Information. Overall, a good agreement between QM and MM energies was achieved with the RMS deviation for energies is 0.11 kcal•mol -1 and 0.50 kcal•mol -1 for the oxidized and reduced forms respectively. The equilibrium values for bonds and angles for the structure optimized with the MM model are also in good agreement with the QM structure. Generally, the deviation for bonds is less than 0.02 Å; for covalent angles the RMSD is lower than 2.4º and for torsions the RMS angle deviation is in the range of 0.2 to 6.2º. The RMS deviation between ab initio and empirical geometries is 0.04 Å and 0.10 Å for lumiflavin in the oxidized and reduced forms respectively. Overall, with the optimized bonded parameters, the empirical model reproduces the QM equilibrium geometries very well. b RMS deviation between QM and MM optimized equilibrium structures for all atoms and heavy atoms of the compound; c RMS deviation for bond length, angles and dihedrals.

The agreement for energies is demonstrated in Figure 2 for selected degrees of freedom in the oxidized form of lumiflavin with the optimal set of bonded parameters. All energies are within 2.0 kcal•mol -1 of the minimum energy, which we achieved by using an additional calculation to determine the range of values (Equation 1and 2 in the Method section) for each bonded term. Overall, the model reproduces well the QM equilibrium conformations of the model compounds as well as QM energies of various deformations along parametrized degrees of freedom.

Soft dihedral angles

These terms include flexible dihedrals along which energy variations are within a range that is sampled in MD simulations, hence accurate treatment of these dihedral PES is important to describe adequately the conformational space of molecules. The parameters of these flexible dihedrals were determined based on points of PES scans. Here we illustrate the quality of the parameters that define the most important conformational changes. In particular, the torsion angle that defines rotations of 10-ethyl group in 10-ethyl-flavin was parametrized to reproduce the complete rotation of 360º. This torsion defines orientation of the flavin group relative to other groups, such adenine in FAD. The results of the PES scan are shown in Figure 3. In the oxidized form, the 10ethyl group has two symmetric minimum positions -90º and 90º relative to the plane of the isoalloxazine group separated by barriers at 0º and 180º. In the protein crystal structures with FMN and FAD in the oxidized form (PDB access codes: 2XOD 26 and 5GV7 [START_REF] Takaba | Distribution of valence electrons of the flavin cofactor in NADHcytochrome b5 reductase[END_REF] ) this torsion angle is -83.1º and 92.3º respectively, in agreement with the QM and MM PES scans in vacuum. In the reduced form, the torsion of the 10-ethyl group in vacuum has two asymmetric minima at -67.8º and 90.0º. The PES near the minimum at 90º is shallow, suggesting that the reduced form of flavin is more flexible around this torsion in protein complexes than in the oxidized form. Indeed, in the protein crystal structures with FMN and FAD in the reduced form (PDB access codes: 4UTK [START_REF] Werther | Redox-dependent substratecofactor interactions in the Michaelis-complex of a flavin-dependent oxidoreductase[END_REF] and 4U2S [START_REF] Golden | High-resolution structures of cholesterol oxidase in the reduced state provide insights into redox stabilization[END_REF] ) this torsion is -66.2º and 102.9º. The first value is close to -67.8º, but the second torsion angle is different from 90.0º, in agreement with the shallow region of PES. Finally, we examine the bending (also called the butterfly) motion of the isoalloxazine group of flavin.

The bending motion of flavin was measured by the torsion defined by atoms C4, N5, N10 and C9 of the isoalloxazine group. It is worth to note that the deformation along this degree of freedom is contributed by several dihedral energy terms and there is no particular dihedral angle term associated with this torsion. Thus, it illustrates the overall quality of the force field parameters. Figure 4 shows ab initio and empirical energies of PES scans along the torsion for the oxidized and fully reduced form of lumiflavin. The oxidized form is characterized by a single minimum at 180.0º and high energies for relatively small deformations. In contrast, the reduced form has two shallow minima connected by a small barrier, which is explained by the electron-rich π molecular orbital system of the fully reduced form. [START_REF] Kyte | Mechanism in Protein Chemistry[END_REF] The force field model reproduces well positions of the two minima in the reduced flavin. The empirical energy barrier of 1.9 kcal•mol -1 is in good agreement with the QM energy barrier of 1.8 kcal•mol -1 , which is important to reproduce structural transitions in simulations with the reduced form of flavin.

Molecular Dynamics simulations of protein complexes

To illustrate the quality of the model, Molecular Dynamics (MD) simulations of protein complexes with bound flavins were performed. The protein structures with a high-to-medium resolution were chosen for MD simulations. MD simulations were performed for 50 nanoseconds. The resulting RMS deviation between simulation and experimental structures is given in Table 4. Figure S2 in the Supplementary Information demonstrates the evolution of RMSD in the MD simulations. The RMS deviation for heavy atoms of flavin cofactors in all simulations is in the range of 0.29 Å to 0.53 Å and in all cases lower than the RMS deviation for unrestrained protein backbone atoms, demonstrating that the model performs as good as the standard CHARMM force field for proteins in protein simulations. As an additional test, the most important soft dihedral angles parametrized in this work that define the geometry of the flavin group were further investigated. The torsion angles observed in MD simulations are given in Table 6. The angle between atoms C4, N5, N10, and C9 defining the butterfly motion of the isoalloxazine group is well reproduced in MD simulations with flavin in the oxidized and reduced forms. For example in MD simulations with FMNred this angle is 153.8º in the experimental structure 4UTK and 156.6±6.6º in MD simulations, which is within the standard deviation from the experimental value. The torsions between atoms C2', C1', N10, and C9A, and C1', N10, C10, and N1 defining the position of the flavin group relative to the rest of molecule are also reproduced within the standard deviation. Overall, the force field model reproduces well the torsions between atoms, which define the overall orientation of the isoalloxazine group relative to the rest of molecule.

Discussions

The present study represents a systematic development of a force field model for flavin cofactors, including riboflavin, flavin mononucleotide and flavin adenine dinucleotide in the most important protonation and redox states. The parametrization was performed to be consistent with the standard method used to develop the CHARMM additive force field. The model is, thus compatible with the other components in the CHARMM force field. The parameters were optimized targeting QM data and validated against experimental structural data. The initial guess for atomic charges was obtained from the CGenFF force field. The charges were then adjusted to reproduce interaction energies and distances of a large number of solute-water interactions, which was important to parametrize the large heterocycle isoalloxazine group. Targeting solute-water interactions was also important to maintain the balance between flavin interactions with solvent and protein residues. In addition, the model reproduces the magnitude and direction of the ab initio dipole moment as well electrostatic potential. Though the dipole moment for charged compounds was not targeted, the model reproduces the QM dipole moment for the anionic molecules. Including the QM electrostatic potential in the charge optimization, in accord with the previous study, [START_REF] Huang | Automated Force Field Parameterization for Nonpolarizable and Polarizable Atomic Models Based on Ab Initio Target Data[END_REF] was found helpful to obtain a better charge distribution in the large isoalloxazine group.

A special emphasis was given to the quality of bonded parameters, with all bonded parameters, including soft torsions and stiff harmonic terms adjusted using computationally intensive PES scans.

METHODS AND MATERIALS Parametrization Scheme

A force field of the additive all-atom CHARMM form was adopted in this work. [START_REF] Vanommeslaeghe | CHARMM additive and polarizable force fields for biophysics and computer-aided drug design[END_REF] The CHARMM force field form is well described in other studies and reviews, [START_REF] Vanommeslaeghe | CHARMM additive and polarizable force fields for biophysics and computer-aided drug design[END_REF] and hence is not discussed in detail here. The atom types were adopted from the CGenFF force field. [START_REF] Vanommeslaeghe | CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields[END_REF] The ParamChem web server (https://cgenff.paramchem.org) was used to assign existing atomic types and to obtain initial guesses for the parameters of model compounds. [START_REF] Vanommeslaeghe | Automation of the CHARMM General Force Field (CGenFF) I: bond perception and atom typing[END_REF][START_REF] Vanommeslaeghe | Automation of the CHARMM General Force Field (CGenFF) II: assignment of bonded parameters and partial atomic charges[END_REF] Parameters that were assigned zero "penalty" by ParamChem, i.e. already optimized in CGenFF, were not considered for optimization in the present study. Parameters of the Lennard-Jones potential were taken directly from the CGenFF force field and were not optimized. The atom names from the Protein Data Bank were used throughout this work. In particular, the atoms C4A and C5A on the isoalloxazine group in riboflavin and FMN were renamed in FAD to C4X and C5X respectively, to be unique and different from the C4A and C5A atoms of the adenine group.

Determination of the intermolecular force field parameters

The intermolecular energy is contributed from Coulomb and Lennard-Jones terms. Consistent with the development of the CHARMM force field, atomic charges were optimized targeting interactions between the model compound and individual water molecules, and the dipole moment of the model compound. Quantum mechanical (QM) electrostatic potential was also included as additional target data in the charge fitting similar to the other work. [START_REF] Huang | Automated Force Field Parameterization for Nonpolarizable and Polarizable Atomic Models Based on Ab Initio Target Data[END_REF] The charge optimization was performed on the compound structures optimized with the B3LYP functional [START_REF] Becke | Density functional thermochemistry. III. The role of exact exchange[END_REF] and 6-31G(d) basis set. [START_REF] Krishnan | Self-consistent molecular orbital methods. XX. A basis set for correlated wave functions[END_REF] Gaussian09 (Gaussian, Inc, Wallingford CT, 2009) was used for all QM calculation.

Atoms of the model compound that can participate in hydrogen bonds were probed by individual water molecules placed in idealized linear orientations. [START_REF] Mackerell | An all-atom empirical potential for molecular modelling and dynamics study of proteins[END_REF] Different rotations of the water molecule were considered around the interaction axis: for polar atoms, the complex was calculated every 45º or 90º of the water probe rotation, and only one or two orientations for non-polar atoms. Each water-compound structure was then optimized by varying the interaction distance, to find the energy minimum for the water position. The B3LYP/6-31G(d) geometry was used for the compound and TIP3P model geometry for the water molecule. The single angle defining the orientation of the water molecule was held fixed during the optimization. Following the same protocol that was used to parametrize the CHARMM force field, calculations were done at the HF/6-31G(d) level. [START_REF] Mackerell | An all-atom empirical potential for molecular modelling and dynamics study of proteins[END_REF][START_REF] Vanommeslaeghe | CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields[END_REF] The interaction energy was calculated for the optimal distance. In accord with the standard CHARMM parametrization protocol, [START_REF] Mackerell | An all-atom empirical potential for molecular modelling and dynamics study of proteins[END_REF] the ab initio interaction energies were scaled by an empirical factor of 1.16 only for neutral polar compounds and the HF/6-31G(d) minimum interaction distance was corrected by subtracting 0.2 Å for all polar interactions involving neutral compounds. To parametrize radical states of flavin, water interaction calculations

were performed with the B3LYP potential and unrestricted open-shell wavefunctions using the same 6-31G* basis set. No correction of QM interaction energies and distances was made for the radical forms of flavin.

The molecular dipole moment, which is defined by the charge distribution, was used to provide additional target data for the optimization of the atomic charges. The dipole moment was used only for the neutral compounds in the charge fitting; for charged molecules, the dipole moment was not targeted as in the previous works. [START_REF] Xu | Additive CHARMM force field for naturally occurring modified ribonucleotides[END_REF] The dipole moment was calculated in vacuum at the HF 6-31G(d) level using the B3LYP optimized conformation. As in the previous studies, to account the molecular polarizability implicitly, the MM optimization targeted dipole moments overestimated by 30% with respect to the QM values for polar compounds. [START_REF] Vanommeslaeghe | CHARMM general force field: A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields[END_REF] During the charge optimization, both the magnitude and direction of the dipole moment were targeted. [START_REF] Xu | Additive CHARMM force field for naturally occurring modified ribonucleotides[END_REF] Since we consider only probe water interactions with a few proton donors and acceptors, QM water interaction data may not be sufficient to define partial charges on all atoms for large compounds, such as isoalloxazine parametrized in this work. This is especially problematic for charged molecules, for which the dipole moment is not targeted, and not all proton acceptors and donors can be probed due to the dominant contribution of the net charge to the water electrostatic interaction. For this reason, we included electrostatic potential to provide additional target data in charge fitting. [START_REF] Huang | Automated Force Field Parameterization for Nonpolarizable and Polarizable Atomic Models Based on Ab Initio Target Data[END_REF] To remain consistent with the standard protocols used for the CHARMM and AMBER force fields, the QM electrostatic potential calculations were performed at the HF/6-31G* level. [START_REF] Huang | Automated Force Field Parameterization for Nonpolarizable and Polarizable Atomic Models Based on Ab Initio Target Data[END_REF] At each iteration during the charge optimization, the Root Mean Square deviation between QM and MM electrostatic potential was evaluated and added with the corresponding weight to the target function.

The charge optimization was performed with the C++ program that was used to parametrize a large number of modified nucleotides in the previous work. [START_REF] Xu | Additive CHARMM force field for naturally occurring modified ribonucleotides[END_REF] The following terms were included with different weights in the target function: the RMS deviation between empirical and ab initio minimum interaction energies, the RMS deviation between ab initio and empirical minimum interaction distances, the absolute difference between the norms of the empirical and ab initio dipole moments, the angle between the empirical and ab initio dipole moments, the RMS deviation between ab initio and empirical electrostatic potential and a term associated with restraints on the charges. The latter term was introduced to prevent large deviations from the starting guess for charges. Charges of symmetrical atoms had identical values during the charge optimization. The initial partial charges were obtained from the ParamChem online server. Charges that were already optimized in the CGenFF force field, for example for adenosine, were not modified in this work. Finally, charges of aliphatic groups were not optimized, in accord with the standard CHARMM method with all methyl groups having a charge of 0.09e on protons. The LJ parameters were not considered to optimization.

Optimization of flexible dihedral parameters

Dihedrals within a molecule can be considered in two classes, stiff and rotatable. Terms associated with stiff dihedrals characterized by a single minimum and high energy for small deformations. Rotatable dihedrals are characterized by a shallow energy surface, and thus may undergo large fluctuations during simulations and may have several minima. Since the molecule can perform large conformational motions along these degrees of freedom, the treatment of this class of dihedral terms is crucial. To parametrize these terms potential energy surface (PES) scans were performed on the torsions, in which torsion angles were scanned in the range from -180° to 180° in 10° increments. Each conformation for the MM calculations was extracted from the QM scan and minimized with a harmonic restraint force constant of 5•10 4 kcal•mol•radian -2 on the target torsion. The MM dihedral parameters were optimized to achieve a minimum deviation between the QM and MM surfaces only in the lower energy regions, PES points with the QM energy of more than 7 kcal•mol above the minimum energy were not considered.

Determination of bonded harmonic energy terms

Parameters for the bonded terms described by harmonic potentials, i.e. bonds, angles, Urey-Bradley, and improper dihedrals, as well as for the stiff dihedrals were optimized using the following protocol. Only parameters not optimized in the CGenFF force field, i.e. with the non-zero "penalty", were optimized. The initial guess was provided by the ParamChem online server. In this work, we used the same method used to parametrize the terms associated with soft dihedrals in CHARMM described in the previous section. In particular, we perform an adiabatic PES scan for each degrees freedom that has adjustable parameters in the force field. It is worth to note that a similar method is also used in the original CHARMM force field to determine force constants by three-point potential energy scans, when the assignment of contributions of the internal coordinates to the vibrations is ambiguous. [START_REF] Xu | Additive CHARMM force field for naturally occurring modified ribonucleotides[END_REF][START_REF] Vanommeslaeghe | Robustness in the fitting of molecular mechanics parameters[END_REF] This method, in contrast to the original method based on the Normal Mode Analysis, [START_REF] Vanommeslaeghe | Robustness in the fitting of molecular mechanics parameters[END_REF] departs from the local minimum, and is more computationally expensive. For example, to parametrize lumiflavin in the oxidized form in this study, 324 optimizations were performed of lumiflavin deformed along the parametrized degrees of freedom.

Since the PES scans were performed by varying one stiff degree of freedom, the potential energy can become very high, even with small deformations. Such high-energy points on PES are not sampled during typical MD simulations and should not be considered in the optimization of boded parameters. To use geometries of the compound within a range of energies that can be reached during MD simulations we propose and use the following protocol. First, using initial values of distortions, force constants of energy terms are estimated. The initial values of the distortion, ∆𝑥 are 0.06 Å, 4.0º, 25.0º and 25.0º for the bond length, valence angle, dihedral and improper angle respectively. The compound geometry is optimized with a chosen internal coordinate constrained to 𝑥 = 𝑥 0 -∆𝑥, where 𝑥 0 is the equilibrium value. The optimization is performed at the B3LYP/6-31G(d) level. The initial force constant of each energy term being parametrized is estimated using the following equation:

𝑘 = 2 (𝐸(𝑥) -𝐸 0 )/∆𝑥 2 , [Eq. 1]
where 𝐸 0 and 𝐸(𝑥) is the minimum energy and energy of the deformed structure. Using this force constant, the new range ∆𝑥′ of x values was estimated with the maximum energy of |𝐸(𝑥) -𝐸 0 | < ∆𝐸 𝑚𝑎𝑥 relative to the minimum energy 𝐸 0 :

∆𝑥′ = √2 ∆𝐸 𝑚𝑎𝑥 𝑘 ⁄ , [Eq. 2]
To optimize each bonded term we use seven points on PES equally spaced in the range of 𝑥 ∈

[𝑥 0 -∆𝑥′, 𝑥 0 + ∆𝑥′], including the minimum energy structure at 𝑥 = 𝑥 0 . In this work 2.0 kcal•mol -1 was used for ∆𝐸 𝑚𝑎𝑥 in Equation 2. All PES scans were performed at the B3LYP/6-31G(d) level.

Once points on PES were computed, the equilibrium values of the MM parameters and force constants were adjusted simultaneously using a C++ program based on the Powell and Amoeba minimization algorithms from Numerical Recipes. [START_REF] Press | Numerical Recipes 3rd Edition: The Art of Scientific Computing[END_REF] Each conformation for the MM calculations was extracted from the QM scan and minimized with a harmonic restraint force constant of 5•10 4 kcal•mol•Å -2 or 5•10 4 kcal•mol•radian -2 on the target bond and angle respectively. At each optimization iteration of bonded parameters, PES scans were performed with CHARMM using a new set of MM parameters. The target function included RMS deviation between ab initio and empirical PES energies and geometries; and restraints to the initial set of parameters provided by the ParamChem server. The MM parameters were adjusted until the RMSD between the MM and the optimized QM geometries and energies could not be reduced significantly.

Molecular Dynamics simulations

To evaluate the quality of the force field model of flavins, Molecular Dynamics (MD) simulations of protein complexes were performed. The protein complexes used in this work are summarized in Table S9. The crystal structures with a high to medium resolution were obtained from the Protein Data Bank (PDB). The simulation setup was similar to previous studies. [START_REF] Aleksandrov | Molecular Dynamics Simulations of the 30S Ribosomal Subunit Reveal a Preferred Tetracycline Binding Site[END_REF][START_REF] Aleksandrov | Tet repressor induction by tetracycline: a molecular dynamics, continuum electrostatics, and crystallographic study[END_REF] In brief, the simulations included protein residues within a 24 Å sphere, centered on the N5 atom of the bound FAD or FMN cofactor in the crystal structures. Protonation states of histidines were assigned by visual inspection and ideal stereochemistry; protonation states of other residues were assigned using PROPKA. [START_REF] Olsson | PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa Predictions[END_REF][START_REF] Dolinsky | PDB2PQR: an automated pipeline for the setup of Poisson-Boltzmann electrostatics calculations[END_REF] Protein atoms between 20 and 24 Å from the sphere's center were harmonically restrained to their experimentally determined positions. In the case of XenA, simulations also included residues and crystal waters from the second monomer of the dimeric structure within 24 Å of the N5 atom of the flavin group. [START_REF] Werther | Redox-dependent substratecofactor interactions in the Michaelis-complex of a flavin-dependent oxidoreductase[END_REF] In addition to crystal waters, a 75 Å cubic box of water was overlaid, and waters overlapping the protein, ligands and crystal water molecules were removed. Periodic boundary conditions were assumed; i.e. the entire 75-Å box was replicated periodically in all directions. All long range electrostatic interactions were computed efficiently by the particle mesh Ewald method, [START_REF] Darden | Treatment of Long-Range Forces and Potential[END_REF] and the appropriate number of potassium counterions were included to render the system electrically neutral. A smooth switching function was used to truncate all van der Waals interactions at the distance of 11 Å. Long range electrostatic forces were evaluated every 4 steps, while short-range non-bonded interactions were computed at each step. MD simulations were performed at constant room temperature and pressure, after 200 ps of thermalization. Constant pressure was maintained using the Berendsen pressure bath coupling [START_REF] Berendsen | Molecular dynamics with coupling to an external bath[END_REF] with the relaxation of 500 fs the compressibility parameter of liquid water.

Constant temperate was maintained by simulating temperature coupling to a heat bath with a room temperature by correcting forces as implemented in the NAMD program. [START_REF] Phillips | Scalable molecular dynamics with NAMD[END_REF] The CHARMM36 force field was used for the protein [START_REF] Best | Optimization of the Additive CHARMM All-Atom Protein Force Field Targeting Improved Sampling of the Backbone ϕ, ψ and Side-Chain χ1 and χ2 Dihedral Angles[END_REF][START_REF] Huang | CHARMM36 all-atom additive protein force field: validation based on comparison to NMR data[END_REF] and the TIP3P model for water. [START_REF] Mackerell | An all-atom empirical potential for molecular modelling and dynamics study of proteins[END_REF][START_REF] Jorgensen | Comparison of simple potential functions for simulating liquid water[END_REF][START_REF] Neria | Simulation of activation free energies in molecular systems[END_REF] The flavin cofactors were modeled using the force field model specifically developed in this wok. Calculations were done with the NAMD program. [START_REF] Phillips | Scalable molecular dynamics with NAMD[END_REF] MD simulations of the protein complexes were continued for 50 nanoseconds.

Supporting Information

Tables with water-compound interactions; selected distances observed in molecular dynamics simulations of flavins in protein complexes; experimental protein structures used for Molecular Dynamics simulations.

Figure 1 .

 1 Figure 1. Three oxidation forms of the flavin group parametrized in this work. The most important conjugate acidbases of the flavin forms are shown.

Figure 2 .

 2 Figure 2. PES scans for selected degrees of freedom in lumiflavin in the oxidized form. The PES scan was performed for a bond (top), for a covalent angle (middle), and for a torsion (bottom). The corresponding term is shown by an arrow.

Figure 3 .

 3 Figure 3. PES scan of the 10-ethyl group rotation in 10-ethyl-flavin in the oxidized form (A) and in the fully reduced form (B). Ab initio energies were computed at the B3LYP/6-31G(d) level. The arrow shows the dihedral angle.

Figure 4 .

 4 Figure 4. PES scan of the butterfly motion in lumiflavin in the oxidized form (A) and in the fully reduced form (B). Ab initio energies were computed at the B3LYP/6-31G(d) level. The torsion angle shown by the arrow is defined between atoms in circles.

  simulations are within the standard deviation of the corresponding distances observed in the experimental structures. Overall, the interactions between flavins and the proteins are well reproduced. The superposition of the experimental structures (PDB access codes 2XOD and 4UTK) with the average position of atoms from 50 ns MD simulations is shown in Figure5. To compute the average position, 2,500 structures saved every 10 ps from 50 ns MD simulations were superimposed on the experimental structure based on the backbone atoms and were averaged. Overall, the model reproduces well the structure of FMN and FAD in protein complexes in the oxidized and reduced forms.

Figure 5 .

 5 Figure 5. The experimental structure (in gray) vs the average structure (in color) observed in MD simulations with (A) FMNox in complex with the NrdI protein (PDB access code 2XOD); and (B) with FMNred in complex with xenobiotic reductase A (PDB access code 4UTK).

  Initial parameters were provided by ParamChem, with parameters previously not existing in the CGenFF force field assigned values by analogy. The bonded parameters were then optimized to reproduce PES energies and geometries. The model reproduces the QM geometry for flavin in the oxidized and fully reduced and important protonation states. It is worth to note the quality of the model to reproduce PES for soft torsions. In particular, the force field model reproduces well the bending motion in the oxidized and reduced forms of flavin. Model validation was based on molecular dynamics simulations of four protein complexes containing the flavin group as a cofactor. The result demonstrate that the model reproduces well interactions between the flavin group and proteins as well as the flavin conformation. Importantly, the force field model reproduces the flavin plasticity and conformational difference between the oxidized and reduced forms as observed in the experimental crystal structures.

  Figure with empirical vs ab initio energies of lumiflavin conformations; Figure with RMS deviation observed in MD simulations of the protein complexes; the developed force field model in the CHARMM format.

Table 1 .

 1 Statistics of water-compound interactions. All values are given for the optimal and initial sets of atomic

	charges.
		energy (kcal•mol -1 )	distance (Å)
	Compound	RMSD max |error| RMSD max |error|
	lumiflavin/HFlox 0.34/2.70 0.96/5.61 0.20/0.95 0.39/2.71
	lumiflavin/H3Flred 0.42/0.61 0.69/1.30 0.17/0.18 0.34/0.38
	lumiflavin/H2Flred -0.24/1.60 0.54/2.73 0.27/0.28 0.89/0.93
	lumiflavin/H2Fl• 0.81/3.82 2.98/9.73 0.31/0.29 0.56/0.57
	lumiflavin/H2Fl• ̶ 0.49/4.21 1.18/5.88 0.18/0.37 0.28/1.45

Table 2 .

 2 Empirical and ab initio dipole moments and electrostatic potentials. The corrected dipole moment was used as an additional restraint only for the neutral compounds during the atomic charge optimization.

	Compound	dipole moment, (Debye)	a angleº	b RMSD esp
		ab initio/optimal/initial optimal/initial optimal/initial
	lumiflavin/oxidized	11.77/11.48/3.85	0.7/105.6	3.27/15.87
	lumiflavin/reduced	6.39/4.88/5.00	0.1/20.9	1.00/3.58
	lumiflavin/ H2Flred -	12.25/12.12/21.90	1.1/9.8	1.48/17.77
	lumiflavin/H2Fl•	11.61/9.04/18.72	0.0/46.3	0.72/31.98
	lumiflavin/HFl• ̶ 24.33/24.57/9.33	0.3/73.1	1.28/60.77
	a angle between the ab initio and empirical dipole moment vectors; b RMS deviation in kcal•mol -1 •Å -1 between ab
	initio and empirical electrostatic potential.

Table 3 .

 3 Comparison between empirical and ab initio optimized geometries for equilibrium structures and energies

	of PES scans
			RMSD
	compound	a energy (kcal•mol -1 ) b geometry (Å) c bond (Å)/angle (º)/dihedral (º)
	lumiflavin/HFlox	325/0.11	0.05/0.04	0.02/2.0/0.2
	lumiflavin/H3Flred	244/0.50	0.19/0.10	0.02/1.9/6.2
	lumiflavin/H2Flred -	39/0.18	0.17/0.13	0.02/2.4/4.6
	10-ethyl-flavin/ HFlox	122/0.40	0.07/0.07	0.02/1.5/1.6
	10-ethyl-flavin/ H3Flred	125/0.33	0.10/0.06	0.02/1.6/2.4

a number of points on PES used to optimize bonded parameters and RMS deviation between QM and MM energies;

Table 4 .

 4 Root Mean Square (RMSD) deviation in molecular dynamics simulations. Standard deviations are given RMSD was computed for unrestrained backbone atoms after superposition on the experimental structure; b RMSD was computed based on heavy atoms of the cofactor.The distances between heavy atoms participating in hydrogen bonds are given in Table5for MD simulations of flavin mononucleotide FMNox bound to the NrdI protein (PDB access code 2XOD). The distances for MD simulations with FADox, FMNred, FADred are given in TableS6-S8. All distances observed in MD

	in parenthesis.
	cofactor PDB ref. code a Protein RMSD (Å) b Cofactor RMSD (Å)
	FMNox	2XOD	0.53 (0.10)	0.45 (0.10)
	FADox	5GV7	0.54 (0.05)	0.53 (0.11)
	FMNred	4UTK	0.51 (0.10)	0.29 (0.05)
	FADred	4U2S	0.37 (0.02)	0.32 (0.05)

a

Table 5 .

 5 Selected distances (Å) in molecular dynamics simulations of flavin mononucleotide FMNox in complex

	with the NrdI protein (PDB entry 2XOD)
	a atom pair X-Ray str. MD simulation
	NG46…O4Fl	2.83	3.16 (0.27)
	OM77...N3Fl	3.05	2.98 (0.18)
	NG79...O2FL	2.84	3.05 (0.23)
	NN71...O2FL	2.92	3.34 (0.27)
	NN71...N1FL	3.23	3.48 (0.23)
	OT42...O2'Fl	2.70	2.82 (0.15)

a protein atoms (left) are labeled by the amino acid to which they belong

Table 6 .

 6 Flavin dihedrals observed in molecular dynamics simulations and experimental structures. Standard deviations are given in parenthesis.

	C4-N5-N10-C9	C2'-C1'-N10-C9A	C1'-N10-C10-N1

Acknowledgements

We thank Prof. Hannu Myllykallio and Dr. Marten Vos for helpful discussions. This work was supported by grant ANR-18-CE44-0002. This work was performed using HPC resources from GENCI-CINES (Grant 2018-A0040710436).


Andrea Palumbo
email: palumbo.and@gmail.com

Onofrio Semeraro

Jean-Christophe Robinet

Luigi De Luca

Receptivity to synthetic jet actuation in boundary layer flows

The flow field generated by the interaction between a zero-net-mass-flux (synthetic) jet and a boundary layer crossflow is studied in the present work. The investigation consists in Direct Numerical Simulations (DNS) of the resulting unsteady field, for various boundary layer Reynolds numbers Re δ * 0 and jet reduced frequencies F + . In particular, the actuation frequency is chosen in order to match the unstable (or less stable) frequency of the incoming crossflow. The analysis is focused on the near-field behaviour of the flow, i.e. next to the jet exit orifice, and in the downstream region. It has been found that the presence of the jet deeply modifies the boundary layer transition. Indeed, the instantaneous flow field is characterized by the presence of unsteady, hairpin-like vortical structures, similar to the ones generated by localized roughness elements. Moreover, streaky structures can be detected in the time-averaged flow field. The overall behaviour of the flow field (investigated via time-averaged boundary layer integral parameters) allows to state that early transition can be achieved by using SJ devices. This work is part of a wider project which aims at realizing a reduced-order modeling of a piezo-driven synthetic jet in crossflow, based on simulations and lumped element modeling of the actuator.

I. Nomenclature

II. Introduction

A great effort has been paid in the last years to reduce the pollutant emissions produced by civil aircrafts. As a general issue, the environmental impact of the airplane may be reduced by working on the external aerodynamics of the airplane, its electronic systems, structures and engines. In particular, the research on the external aerodynamics focused mainly on drag reduction techniques, based on laminar wings/fuselages, riblets and active flow control (AFC) systems. AFC techniques are commonly employed to suppress/enhance turbulence and reduce the extent of the separated zones on the lifting surfaces, which are key factors of the vehicle performances. In this framework, synthetic jet (SJ) actuators have been proven to be useful devices, due to their unique ability to transfer momentum to the incoming crossflow with a null-net-mass-flow-rate (per actuation cycle) across the jet boundary. For this reason they do not require any pumping system (namely, they are synthetized within the fluid flow to be controlled), and in many cases they are preferred to other active flow control systems, satisfying also the request of low weight and reduced size. Several techniques to produce a synthetic jet can be found in literature; however two methods have been extensively studied in the last years, namely the piezo-driven [START_REF] Glezer | Synthetic jets[END_REF] and plasma synthetic jet actuators (PSJ) [START_REF] Zong | Plasma synthetic jet actuators for active flow control[END_REF]. In the former one, the jet is created by periodic volume variations of a resonant cavity, produced by the deformation of a piezo-electric element located on the bottom wall of the cavity. In the latter one, a high-speed and high-temperature jet is obtained from an electrical discharge between two or three electrodes, enclosed in a small cavity. The high-pressure plasma formed by the discharge exhausts through an orifice, converting the internal energy of the air into the jet kinetic one. In the end, fresh air is ingested within the cavity, refilling the device for the next cycle. The choice of the actuator to be used for a certain application depends on the desired exit velocity and applied frequency, the available energy and space for the actuator, and safety/reliability considerations. Many works have been recently devoted to the application of SJ actuators to control aerodynamic flows, both in industrial applications and research projects. The former studies are typically devoted to evaluate the effectiveness of such devices to improve the performances of an aircraft. In this framework, the works are also intended to solve practical problems and give design guidelines for the implementation of the actuators on aircraft components, such as wings [START_REF] Arena | Feasibility studies for the installation of Plasma Synthetic Jet Actuators on the skin of a morphing wing flap[END_REF] and vertical tails [START_REF] Jansen | Interaction of a Synthetic Jet with Separated Flow over a Vertical Tail[END_REF]. The potential impact of this technology on the performances of the aircraft (for example, in terms of reduced fuel consumption) has been also discussed [START_REF] Crowther | An evaluation of the mass and power scaling of synthetic jet actuator flow control technology for civil transport aircraft applications[END_REF]. The latter studies, instead, generally investigate the interaction between the jet and canonical flows such as flat-plate laminar boundary layers, circular cylinder in crossflow [START_REF] Ma | Vortex formation and evolution for flow over a circular cylinder excited by symmetric synthetic jets[END_REF], liquid or gas free jets [START_REF] Smith | Jet vectoring using synthetic jets[END_REF][START_REF] Tamburello | Active control of a free jet using a synthetic jet[END_REF]. The main aim of these studies is to identify the vortical structures of the unforced flow field and to understand how they are affected by the actuation of the SJ. The findings are useful to shed light on the flow physics when more complex phenomena are involved. For instance, the application of a SJ in a boundary layer can give information on the control of a lifting surface at low angles of attack, whereas the interaction between a SJ and the wake behind a bluff body can mimic the separation control in stall and post-stall configurations. The analyses have been generally carried out in the frequency domain, often with low-order reconstructions of the controlled flow field. It has been found that synthetic jets are able to interact with the crossflow by annihilating or enhancing some modes, or introducing new vortical structures into the flow field. Among the various crossflows analysed in these investigations, a special attention has been devoted to the case of boundary layer; early works were concerned with continuous jets in crossflow, carried out using experimental measurements, flow visualization techniques, direct or large eddy simulations and stability analyses. Various fluid-dynamics instabilities have been detected using global stability analyses. Indeed, the complex flow field generated by a jet issuing in a boundary layer flow is dominated by the interaction of three common types of instabilities, a shear layer Kelvin-Helmholtz instability in the jet region, a (possible) elliptic instability of the counter-rotating vortex pair, and a vortex shedding mechanism, similar to the von Kármán vortex street behind a cylinder, in the near-wall region [START_REF] Regan | Global linear stability analysis of jets in cross-flow[END_REF][START_REF] Iyer | A numerical study of shear layer characteristics of low-speed transverse jets[END_REF]. These recent analyses strengthen the findings of earlier experimental works [START_REF] Kelso | An experimental study of round jets in cross-flow[END_REF][START_REF] Megerian | Transverse-jet shear-layer instabilities. Part 1. Experimental studies[END_REF], exploring a wider range of velocity ratios (between the jet and the unperturbed boundary layer) and crossflow Reynolds numbers. On the other hand, a smaller number of works concerned the interaction of synthetic jets with a boundary layer, both experimentally and via simulations. The vortical structures generated by a piezo-driven SJ issuing in a boundary layer have been investigated in detail in [START_REF] Jabbal | The near wall effect of synthetic jets in a boundary layer[END_REF][START_REF] Zhou | Coherent structures produced by the interaction between synthetic jets and a laminar boundary layer and their surface shear stress patterns[END_REF][START_REF] Xia | Flow characterization and modeling of strong round synthetic jets in crossflow[END_REF] for laminar boundary layers and in [START_REF] Milanovic | Synthetic jets in cross-flow[END_REF][START_REF] Schaeffler | The interaction of a synthetic jet and a turbulent boundary layer[END_REF][START_REF] Wu | Large-eddy simulations of circular synthetic jets in quiescent surroundings and in turbulent cross-flow[END_REF][START_REF] Berk | Trajectory of a synthetic jet issuing into high-Reynolds-number turbulent boundary layers[END_REF][START_REF] Berk | Effects of vortex-induced velocity on the development of a synthetic jet issuing into a turbulent boundary layer[END_REF] for turbulent boundary layers. Similar studies have been carried out also for plasma SJs [START_REF] Zong | Effect of velocity ratio on the interaction between plasma synthetic jets and turbulent cross-flow[END_REF]. Besides, a few studies focused on the desired characteristics to be obtained in the controlled flow field. Synthetic jets were typically employed to obtain the so-called virtual aeroshaping, [START_REF] Mittal | On the virtual aeroshaping effect of synthetic jets[END_REF]: SJs can generate time-averaged recirculation regions that are larger than the actuator diameter, which can be used to introduce a "virtual morphing" of the aerodynamic surfaces. This behaviour can be achieved when the actuation frequency is much larger than the characteristic frequency of the uncontrolled flow field (which, nonetheless, is generally an unknown quantity): such a technique has been successfully used to control the flow over airfoils [START_REF] Chatlynne | Virtual aero-shaping of a Clark-Y airfoil using synthetic jet actuators[END_REF], bubbles separation [START_REF] Dandois | Numerical simulation of active separation control by a synthetic jet[END_REF], and bluff body wakes. Notwithstanding all the cited papers, the optimized design of a new actuator for a fixed crossflow condition is still a difficult task, since its optimization involves many geometric and structural characteristics of the actuator. Moreover, the device should be able to control the incoming crossflow in a wide range of crossflow conditions. For these reasons, a huge number of experimental tests or numerical simulations should be planned to properly design an actuator, which should lead to intractable computational and experimental costs. However, reduced-order modeling (ROM) has gained an increasing interest in the last decades, in order to reduce the required resources for the design of the actuator and its implementation [START_REF] Rowley | Model reduction for flow analysis and control[END_REF]. Physical modelling of piezo-driven and plasma synthetic jet actuators, based on lumped approximations of the thermo-fluid-dynamic quantities in the cavity and at the exit orifice, has been originally developed for the case of a single SJ actuator, communicating with the quiescent environment through a single orifice or slot [START_REF] Gallas | Lumped Element Modeling of Piezoelectric-Driven Synthetic Jet Actuators[END_REF][START_REF] Chiatto | LEM characterization of synthetic jet actuators driven by piezoelectric element: A review[END_REF][START_REF] Zong | Analytic model and frequency characteristics of plasma synthetic jet actuator[END_REF][START_REF] Chiatto | A Calibrated Lumped Element Model for the Prediction of PSJ Actuator Efficiency Performance[END_REF]. Further developments led to dedicated lumped-element models (LEM) of complex devices, with more cavities and/or orifices [START_REF] Chiatto | Numerical and experimental characterization of a double-orifice synthetic jet actuator[END_REF][START_REF] Chiatto | Design approach to predict synthetic jet formation and resonance amplifications[END_REF]. The flow field generated by unconventional SJ actuators, with noncircular orifices or multi-slot and multi-cavity characteristics are still object of current research [START_REF] Lindstrom | Effect of Orifice Geometry on Synthetic Jet Evolution[END_REF][START_REF] Shi | Evolution of elliptic synthetic jets at low Reynolds number[END_REF][START_REF] Palumbo | Measurements versus Numerical Simulations for Slotted Synthetic Jet Actuator[END_REF], and these studies could be the basis for new LEMs. On the contrary, there are no LEMs for the synthetic jets in crossflow, which led the researchers to rely to empirical models or trial-and-error techniques for the design of the actuators. The present study approaches the problem of realizing a reduced-order modeling of a piezo-driven synthetic jet in crossflow, based on simulations and lumped element modeling of the actuator. Direct numerical simulations of the flow field generated by the interaction between the synthetic jet and the crossflow have been carried out, in order to recognize the vortical structures arising in the flow field. In the present investigation the actuator is not fully included in the CFD simulations, but it is simply represented by an oscillating inlet condition, enforced at the orifice lower end. However, this approach is suitable for integration in a LEM-based design of the actuator: the jet velocity at the orifice exit is obtained via existing LEM models for the jet in quiescent conditions, and changes in geometry, actuation frequency and supplied voltage can be easily considered with the LEM formulation. This procedure allows also to include the compressibility effects of the flow in the cavity and the coupling mechanisms between the diaphragm motion and the fluid flow (Helmholtz resonance). Differently from the virtual aeroshaping studies, in this case the actuation frequency is next to the characteristic frequency of the crossflow, in order to deepen the understanding of how SJs modify the modes of the uncontrolled flow. The characteristic frequency of the Blasius boundary layer is evaluated via standard linear stability theory. The paper is organized as follows: the numerical setup (CFD code, numerical domain, boundary conditions, flow parameters and dimensionless quantities) is described in Sec. III. The main results of the investigations, including time-averaged and instantaneous flow features, are reported in Sec. IV. Finally, conclusions are drawn in Sec. V, with a special reference to flow control applications.

III. Numerical approach A. Numerical setup

The direct numerical simulations of the synthetic jet in crossflow are carried out using the open-source CFD code Nek5000 [START_REF] Fischer | nek5000 Web page[END_REF]. The present code was deemed to be adequate for the present investigation, since it has been previously employed in several studies concerning continuous jets in crossflow [START_REF] Peplinski | Global stability and optimal perturbation for a jet in cross-flow[END_REF], boundary layer suction devices [START_REF] Brynjell-Rahkola | On the stability of a Blasius boundary layer subject to localised suction[END_REF], flow around roughness elements [START_REF] Bucci | Roughness-induced transition by quasi-resonance of a varicose global mode[END_REF] and SJs issuing in quiescent environment [START_REF] Capuano | Comparative study of spectral-element and finite-volume solvers for direct numerical simulation of synthetic jets[END_REF]. This code employs the spectral element method (SEM) to solve the unsteady, incompressible Navier-Stokes equations. The computational domain is partitioned in hexahedral spectral elements: pressure and velocity fields are defined in each element as the combination of the tensor product of 1-dimensional Lagrange polynomials of degree N for the velocity and N -2 for the pressure field. For all the simulations shown in the present work it is chosen N = 9. The time integration, instead, is based on a third-order splitting method. The linear systems arising from the pressure and the velocity discretizations are solved using, respectively, GMRES and conjugate gradient methods, with tolerances of 10 -7 and 10 -9 . The numerical domain is shown in Fig. 1: it includes both the external flow field and the jet pipe. All the lengths and the velocities have been made dimensionless against the orifice diameter D and the crossflow velocity U ∞ . The dimension of the external domain are L x = 55D, L y = 10D, L z = 10D, where x denotes the direction of the crossflow, y the direction of the jet and z the spanwise one. The decision of including the inlet pipe in the simulations derives from a recent literature output: in fact, previous studies on continuous jets in crossflow have shown that it is necessary to include the inflow pipe to correctly obtain the flow field at the orifice exit plane [START_REF] Iyer | A numerical study of shear layer characteristics of low-speed transverse jets[END_REF][START_REF] Peplinski | Global stability and optimal perturbation for a jet in cross-flow[END_REF][START_REF] Brynjell-Rahkola | On the stability of a Blasius boundary layer subject to localised suction[END_REF]. The boundary condition at the crossflow inlet, located 10D upstream of the jet, is a standard Blasius profile, whose displacement thickness is chosen to have a prescribed value of Re δ * 0 at the nozzle centerline location (according to the Blasius' similarity solution). The inlet condition for the jet is an oscillating Dirichlet condition: the velocity profile is obtained from Womersley's flow theory, eq. ( 1):

V(r, t) = Real 1 A 1 - J 0 (i 3/2 W or) J 0 (i 3/2 W o) e iωt , A = V a 1 - 2 i 3/2 W o J 1 (i 3/2 W o) J 0 (i 3/2 W o) (1)
where J n is n-th order Bessel function, r the axial distance from the pipe symmetry axis (made dimensionless against the pipe radius), W o = ωD 2 4ν the Womersley number and V a the velocity amplitude on the symmetry axis. The pipe length is chosen large enough to have a fully-developed pulsating jet flow at its lower end; after a preliminary work H = 5D has been chosen. The present velocity profile (1) has been validated by means of simple axisymmetric computations, which are briefly summarized in Appendix.

No-slip conditions have been applied at the plate and on the nozzle boundary, and periodic boundary conditions are enforced at the spanwise ends z = ±L z /2. Such a boundary condition is often employed to simulate a synthetic jet array with a reduced computational cost. In this framework L z is the spanwise spacing between the axis of two consecutive orifices. This distance comes out as an important project parameter of several active and passive flow control systems, as localized roughness elements, micro-vortex generators and fluidic actuators. Indeed, the control efficiency of an array of actuators is strongly dependent on this distance, being used to modulate the control output in the spanwise direction.

The interaction between synthetic jets issuing from orifices in close proximity is widely investigated in literature, both in quiescent conditions and in presence of a crossflow, employing several actuator configurations [START_REF] Palumbo | Numerical study on the flow field generated by a double-orifice synthetic jet device[END_REF][START_REF] Watson | Contribution to the understanding of flow interactions between multiple synthetic jets[END_REF][START_REF] Mcglynn | Large Eddy Simulation of Flow Interactions of Segmented Synthetic Jets on an Airfoil[END_REF]. However, for L z = 10D it is possible to neglect the interaction between adjacent jets, at least in the near-field region.

A natural outflow boundary condition -pn + ν∇U • n = 0 has been applied at the outlet, while a mixed boundary condition u = U ∞ , v y = w y = 0 is used on the top edge of the boundary layer. The boundary conditions for the pressure field are automatically enforced by the code to ensure the incompressibility constraint of the velocity field.

B. Flow parameters

The flow field depends on several crossflow and jet non-dimensional parameters: the boundary layer Reynolds number at the nozzle centerline station Re δ * 0 = U ∞ δ * 0 /ν, the velocity ratio between the jet and the crossflow R = V j /U ∞ , the reduced frequency F + = ω/ω r e f , the jet Reynolds number Re j = V j D/ν and the pipe dimensionless length H/D. In the previous definitions V j is the jet mean velocity amplitude, ω = 2π f its circular frequency, ω re f the value of the circular frequency which maximizes the growth rate of the boundary layer instability at the nozzle position, U ∞ and δ * 0 the boundary layer velocity and thickness at the streamwise position of the pipe axis, ν the kinematic viscosity. The influence of H/D on the flow behaviour can be omitted if the flow is fully developed at the pipe inlet, so it will be not investigated in the present study. The non-dimensional parameters chosen in the simulations are summarized in Table 1. Two crossflow Reynolds numbers have been selected, Re δ * 0 = 500 and 550, the smaller being slightly below the critical Reynolds number for the Blasius linear instability, whereas the larger one is above this threshold. These flow conditions can be also expressed in terms of local Re x 0 = U ∞ x 0 /ν, where x 0 is the virtual position of the nozzle center, thus obtaining Re x 0 ≈ 84500, and 102000. Variations in the crossflow Reynolds number have been obtained by changing the (virtual) distance of the jet from the leading edge of the plate. From an experimental point of view, this procedure is equivalent to move the actuator upstream, for a fixed value of the freestream velocity. The velocity ratio is chosen equal to 0.1, in order to give a small perturbation on the boundary layer and avoid triggering instability mechanisms just downstream the orifice exit. The jet Reynolds number is obtained as a combination of the velocity ratio and the ratio between the boundary layer displacement thickness and the pipe diameter. The latter ratio is chosen as 0.25 for the S2-S4 cases, and 0.227 for the S1 case. As regards the reduced frequency, the instability frequencies, obtained performing spatial, linear stability analyses at the streamwise location of the nozzle centerline (when the jet is off), for each Reynolds number, are chosen as reference frequencies. The actuation frequency has been chosen equal to the less stable frequency of the boundary layer for the S1 case, and equal to the most unstable frequency for the S3 case, whereas it lies outside the instability range for S2 and S4 simulations. Furthermore, a preliminary simulation (not shown herein), in which the SJ is off (a no-slip condition is applied at the lower end of the pipe) was carried out. The flow field, after a short transient phase, reached a steady state, as expected for the present computational setup. The effect of the inlet pipe was confined in the near field only, consisting in a steady, recirculation region of small extension. It is worth noting that three different parameters are commonly used to identify the SJ in a quiescent environment, the Reynolds number based on a proper jet average velocity, the Strouhal number and the Stokes number, defined in eq.( 2)

Re 0 = V 0 D ν , St 0 = ωD V 0 , S = ωD 2 ν = 2W o (2)
In the Reynolds number definition V 0 is the averaged jet velocity, given by eq. ( 3)

V 0 = 1 T ∫ T /2 0 ∫ A V(x, y = 0, z)dxdz dt (3)
and T = 1/ f is the actuation cycle period. The stroke length L 0 is then defined as L 0 = TV 0 . The Stokes number is also the main parameter characterizing the behaviour of the oscillating pipe flow within the pipe and its possible instability and transition to turbulence. The latter parameters are not independent from the aforementioned ones and are summarized in Table 1. Finally, in several flow control applications the momentum coefficient C µ is considered in place of the velocity ratio:

C µ = ρ jet U 2 jet ρ ∞ V 2 ∞ ( 4
)
where ρ jet and ρ ∞ are, respectively, the density of the jet and the crossflow. In the present case Eq. ( 4) simply reduces to

C µ = R 2 .

IV. Results and discussion

The results of the present investigation will be presented in terms of instantaneous and time-averaged velocity fields, for all the configurations introduced in Table 1. In particular, simulation S3 is considered as the reference case, and discussed thoroughly in Sec. IV.A, whereas the effects of changing the actuation frequency (cases S2 and S4) and boundary layer Reynolds number (case S1) are investigated, respectively, in Sec. IV.B and IV.C.

A. Baseline case (Re δ * 0 = 550, ω = ω r e f ) A view of the vortical structures generated by the interaction between the jet and the boundary layer for case S3 is depicted in Fig. 2, through isosurfaces of the instantaneous λ 2 field [START_REF] Jeong | On the identification of a vortex[END_REF]. Hairpin-like structures, characterized by a varicose symmetry, are generated at each actuation cycle and advected downstream by the incoming boundary layer flow. The distance between two consecutive vortices in the near field is ≈ 6D, and it is close to the wavelength predicted by spatial linear stability for the uncontrolled Blasius proflle at the jet position. The vortices slow down during their motion, and eventually meet each other and merge; the ultimate structure of the flow field in very complex, as vortex interactions, subsequent instabilities and dissipation rise through the domain, however the flow is not fully turbulent at the end of the plate. This feature is confirmed by the fact that the varicose symmetry of the flow is preserved until the end of the domain, which in turn can be related to the short streamwise length of the domain, the low jet-to-BL velocity ratio, and the low Reynolds number of the uncontrolled boundary layer at the jet position. The effect of the synthetic jet on the boundary layer can be investigated by analyzing the time-averaged controlled flow field. The streamwise velocity field distribution in the symmetry plane is depicted in Fig. 3(a), whereas the time-averaged vortical structures are represented if Fig. 3(b). Since the instantaneous flow field is varicose-symmetric with respect to this plane, its maximum deviation from the uncontrolled flow field lies in this plane. The effect of the jet is restricted to a very thin region in the wall normal direction (0 ≤ y/D ≤ 2) , basically consisting in a boundary layer growth enhancement. The time-averaged structures include a horseshoe vortex, located upstream the jet, and two streamwise streaks on the leeward side, similarly to low-Reynolds-number continuous jets in crossflow. The control action of the synthetic jet on the crossflow is better appreciated in Fig. 4, where the difference U p = U -U B between the controlled time-averaged flow field and the corresponding Blasius solution U B (x, y) is depicted. Finally, the time-averaged boundary layer integral parameters, namely the displacement thickness δ * , the momentum thickness θ and the shape factor H = δ * /θ, are shown in Fig. 5. It is possible to distinguish three main regions, while moving downstream from the orifice exit: a transition region, where all the integral parameters grow much faster than a laminar boundary layer. In the following region their growth rate reduces, whereas in the far field the growth rates of the boundary layer parameters increase again, similarly to high-Reynolds-number turbulent boundary layers. It is worth noting that the time-averaged shape factor in the far field is much below the laminar value of 2.6, being comparable with the typical value assumed for turbulent boundary layers of 1.4 [START_REF] Schlichting | Boundary-layer theory[END_REF]. Moreover, also the streamwise distribution of the skin friction coefficient, shown in Fig. 6, collapse onto the turbulent correlation curve C f (x) = 0.455/log 2 (0.06Re x ) [START_REF] White | Viscous fluid flow[END_REF], after a transition route from their laminar solutions. These facts allow to conclude that a quasi-turbulent state is reached at the end of a domain.

B. Effect of the actuation frequency

Fig. 7 depicts the evolution of the vortex structures in the streamwise direction for Cases S2 and S4. The topology of the flow field is very similar to Case S3: the vortical structures are still characterized by a varicose symmetry, the vortices spread during their motion and the flow remains laminar till the end of the numerical domain. The main differences in the flow development are related to the streamwise and spanwise dimensions of the controlled flow. Firstly, the characteristic streamwise distance of the vortices is reduced for Case S4, owing to the higher actuation frequency; on the contrary, this distance increases for Case S2. Moreover, the spread of the hairpin vortices is still influenced by the actuation frequency: indeed, Case S4 reveals a higher spreading rate the case S3. Time-averaged results, again, can be used to assess the capability of the SJ actuation to modify the boundary layer development. Fig. 8 shows the streamwise distribution of the boundary layer parameters: the strongest deviation from the laminar solution towards a quasi-turbulent flow is observed for the S4 case, which rapidly departs from the laminar solution, as shown by the shape factor diagram. Similar conclusions can be drawn by analyzing the time-averaged distribution of the skin friction coefficient (Fig. 9).

C. Effect of the Reynolds number

Fig. 10 depicts the evolution of the vortical structures at the end of the tenth actuation cycle for Cases S1 and S3. It can been seen that reducing the boundary layer Reynolds number below its critical value does not modify the behaviour of the flow. Indeed, the vortical structures are still advected dowstream by the crossflow, leading again to a complicated, transitional flow (as shown by the C f streamwise distribution in Fig. 11). However, the spacing between the vortices, the initial vortex celerity and their spreading rate are slightly different from the baseline case. In fact, these differences can be attributed to two conflicting causes: it must be kept in mind that the S1 case is obtained for F + = 1, but the frequency which maximizes the TS growth rate is decreasing with the crossflow Reynolds number. For this reason, S1 case is characterized by a smaller δ * 0 and a higher applied frequency F. These parameters act in opposite ways on the transition process, yielding to similar C f distributions for the considered cases.

D. Discussion

The flow field topology generated by the jet-boundary layer interaction perfectly resembles the unsteady flow field generated by an isolated cylindrical roughness immersed in a boundary layer flow. The reader can compare, for instance, Fig. 2 of the present work with the numerical campaign carried out in [START_REF] Bucci | Roughness-induced transition by quasi-resonance of a varicose global mode[END_REF] for similar crossflow conditions. It is well known nowadays that roughness elements are able to produce streamwise streaks, which, if stable, attenuate the growth of the Tollmien-Schlichting waves and, as a consequence, delay the transition to turbulence. Nonetheless, in many practical situations the streaks generated by the cylindrical element can be unstable themselves, leading to self-sustained hairpin structures and, as a consequence, early transition to turbulence. This result is also well known in literature [START_REF] Klebanoff | Mechanism by Which a Two-Dimensional Roughness Element Induces Boundary-Layer Transition[END_REF], indeed roughness elements are often used as passive flow control devices to promote transition. From a flow control perspective, however, the performances of such devices drop in off-design conditions, as they work properly only within a small range of boundary layer Reynolds numbers. Synthetic jet actuators, on the contrary, assure good performances in a wider range of freestream conditions, owing to their flexibility and wide applied frequency interval. In particular, the actuation frequency and applied voltages can be changed during their operation in order to handle different crossflow conditions.

V. Conclusions and future work

This paper deals with direct numerical simulations of the flow field generated by the interaction between a synthetic jet (produced by a zero-net-mass-flux device) and a boundary layer crossflow, for the purpose of recognizing the vortical structures arising in the flow field. Differently from previous virtual aeroshaping studies, the original aspect of this work consists in the fact that the actuation frequency is next to the characteristic frequency of the crossflow, namely it is equal to the unstable (or less stable) frequency, so as to deepen the understanding of how the synthetic jet modifies the modes of the uncontrolled flow. A baseline flow configuration has been initially chosen, and additional simulations have been carried out in order to ascertain the effects of the boundary layer Reynolds number and applied frequency on the controlled flow fields. The instantaneous controlled flow field is dominated by hairpin structures for all the considered numerical setups, characterized by a varicose symmetry, similarly to the unsteady flow behind a cylindrical roughness. The time-averaged flow field, instead, is characterized by streaky structures, generated at the spanwise ends of the orifice, and a horseshoe vortex, located upstream of the jet. The analysis of the boundary layer integral parameters revealed that the overall effect of the jet is to promote transition, since the time-averaged shape factor and the skin friction coefficient values are in close agreement with turbulent correlation ones. Moreover, the transition region is smaller if higher applied frequencies are chosen. As a final comment, it is worth noting that the present study approaches the problem of realizing a reduced-order modeling of a piezo-driven synthetic jet in crossflow to suppress separated flows, based on simulations and lumped element modeling of the actuator. In particular, the analysis of the boundary layer shape factor H shows that SJ actuators is capable of making flows more resistant to separation. Moreover, their application could be preferred to passive flow control devices, as roughness elements or micro-vortex generators (MVG), since changes in the BL characteristics can be tackled by changing the actuation frequency and the applied voltage on-the-fly.

NF

  = order of the velocity field polynomials L x , L y , L z = dimensions of the (external) computational domain D = pipe diameter δ * 0 = boundary layer displacement thickness at the nozzle center (without the jet) x 0 = virtual distance of the nozzle center from the leading edge of the plate Re δ * 0 = boundary layer Reynolds number at the nozzle center position (based on δ * 0 ) Re x 0 = boundary layer Reynolds number at the nozzle center position (based on x 0 ) V j = space-averaged jet exit velocity at peak blowing R = jet-to-crossflow velocity ratio F + = reduced actuation frequency Re j = jet Reynolds number (based on the maximum velocity amplitude) C µ = jet momentum coefficient V 0 = averaged jet velocity L 0 = stroke length Re 0 = jet Reynolds number (based on V 0 ) St 0 = jet Strouhal number (based on V 0 ) ω * = non-dimensional frequency of the boundary layer based on U ∞ /δ * 0 = non-dimensional frequency of the boundary layer based on U 2 ∞ /ν -α i = spatial growth rate (from linear stability theory) δ * = time-averaged boundary layer displacement thickness θ = time-averaged boundary layer momentum thickness H = time-averaged boundary layer shape factor C f = time-averaged boundary layer skin friction coefficient U p = time-averaged deviation of the controlled flow from the Blasius profile

Fig. 1

 1 Fig. 1 Computational domain and main geometrical parameters of the simulations.

Fig. 2 Fig. 3 ( 4 U p /U ∞ 4

 2344 Fig. 2 Instantaneous view of the vortical structures, highlighted by λ 2 = -0.1 isosurfaces, and coloured by the instantaneous streamwise velocity. Case S3, tenth cycle, t/T = 1 (end of suction phase).

Fig. 5

 5 Fig. 5 Streamwise, time-averaged distributions of the boundary layer integral parameters δ * (a), θ (b) and H (c) in the symmetry plane for Case S3 (continuous line). The laminar solution is depicted with dashed lines for comparison; in (c) also the value H = 1.4 is shown (dot-dashed line) for reference.

Fig. 6

 6 Fig. 6 Time-averaged distribution of the skin friction coefficient C f along the streamwise direction in the symmetry plane for Case S3. The laminar (Blasius) solution (dashed line) and the turbulent correlation C f (x) = 0.455/log 2 (0.06Re x ) are shown for comparison.

Fig. 7

 7 Fig. 7 Instantaneous view of the vortical structures, highlighted by λ 2 = -0.1 isosurfaces, and coloured by the instantaneous streamwise velocity. (a) Case S2, (b) case S4, tenth cycle, t/T = 1 (end of suction phase).

Fig. 8

 8 Fig. 8 Streamwise, time-averaged distribution of the boundary layer integral parameters δ * , θ and H in the symmetry plane for Case S2 (red line), S3 (black), S4 (blue). The laminar solution is depicted with dashed lines for comparison; in (c) also the value H = 1.4 is shown (dot-dashed line) for reference.

Fig. 9 Fig. 10

 910 Fig. 9 Time-averaged distribution of the skin friction coefficient C f along the streamwise direction in the symmetry plane for Case S2 (red line), S3 (black), S4 (blue). The laminar (Blasius) solution (dashed lines) and the turbulent correlation C f (x) = 0.455/log 2 (0.06Re x ) are shown for comparison.

Fig. 11

 11 Fig. 11 Time-averaged distribution of the skin friction coefficient C f along the streamwise direction in the symmetry plane for Case S1 (green line) and S3 (black). The laminar (Blasius) solutions (dashed lines) and the turbulent correlation C f (x) = 0.455/log 2 (0.06Re x ) are shown for comparison.

Table 1 Non-dimensional parameters used in the simulations.

 1

	Parameter		Simulation acronym
		S1	S2	S3	S4
	Crossflow Reynolds number, Re δ * 0	500	550	550	550
	Crossflow Reynolds number, Re x	≈ 84500	≈ 100000	≈ 100000	≈ 100000
	Velocity ratio, R	0.1	0.1	0.1	0.1
	BL thickness/jet diameter, δ * 0 /D	0.227	0.25	0.25	0.25
	Jet Reynolds Number, Re j	220	220	220	220
	Reduced Frequency, F +	1	0.5	1	2
	Jet Reynolds number, Re 0	70	70	70	70
	Strouhal number, St	16.8	7.41	14.8	29.2
	Stokes number, S	34.4	22.8	32.2	45.6
	Jet momentum ratio, C µ	0.01	0.01	0.01	0.01

Downloaded by 87.18.24.238 on January 5, 2020 | http://arc.aiaa.org | DOI: 10.2514/6.2020-0099

Downloaded by 87.18.24.238 on January

5, 2020 | http://arc.aiaa.org | DOI: 10.2514/6.2020-0099

Acknowledgments

Andrea Palumbo was supported by Università degli Studi di Napoli Federico II and Compagnia di San Paolo under the STAR Program 2018. Part of this work was supported by a grant of HPC time on CINES supercomputer Occigen under the project A0052A062 and on CINECA supercomputer MARCONI under the ISCRA-C project INJECT.

Appendix

A. Validation of the inlet profile

Preliminary simulations of the oscillating pipe flow have been carried out to validate the analytical velocity profile used in the present work. Two-dimensional (axysimmetric) simulations of the pipe in cylindrical coordinates have been performed, enforcing a top-hat, oscillating velocity profile at the pipe inlet. The velocity profiles, after a short development region, converge to the analytical solution of Eq. ( 1), as shown in Fig. 12, where D is the pipe diameter, and r the radial coordinate; open circles represent the numerical data, and the continuous lines the analytical solution. Both the Stokes layer (of order D/W o) and the central region are well modelled by the analytic solution.

B. Choice of the actuation frequency via Linear Stability Theory

The actuation frequency used in the simulation have been obtained performing a spatial, linear stability analysis of the unactuated boundary layer at the streamwise position of the centerline of the nozzle. In this framework the velocity field u and the pressure p perturbation fields are decomposed using the ansatz in Eq. ( 5

e iαx e iβz e -iω * t + c.c.

which is introduced in the linearized Navier-Stokes equations. In the present work the spatial stability analysis for the Blasius boundary layer has been carried out by means of a Chebyschev collocation method; the results, in terms of marginal stability curve (α i = 0) and spatial growth rate -α i at four different Reynolds numbers are shown in Fig. 13. The nondimensional frequency F = ω * /Re δ * 0 • 10 6 has been introduced to scale the frequency along the streamwise direction and eliminating the dependence of the displacement thickness on the Reynolds number. As is known, this allows to have a constant frequency for various downstream positions [START_REF] Schmid | Stability and transition in shear flows[END_REF]. The ω r e f values to be used in the simulations are obtained as follows. For each Reynolds number Re δ * 0 the value of ω * which maximizes the spatial growth (namely, that of the most amplified wave) is first evaluated; this frequency is then made dimensional by multiplying it by the reference quantity U ∞ /δ * 0 .


Michele Alessandro Bucci

Onofrio Semeraro

Alexandre Allauzen

Guillaume Wisniewski

Laurent Cordier

Lionel Mathelin

Control of chaotic systems by deep reinforcement learning

come

• Laminar to turbulent transition

• Separation zone • Transonic buffet • Fluid structure interaction • …
To this end, we test a Deep Reinforcement Learning (DRL) algorithm by applying it for the control of a nonlinear, chaotic system governed by the Kuramoto-Sivashinsky (KS) equation [1]. DRL uses reinforcement learning principles for the determination of optimal control solutions and deep Neural Networks for approximating the value function and the control policy [1,2].

-INTRODUCTION

(*) Michele Alessandro BUCCI (Email: bucci.malessandro@gmail.com) FLOWCON PROJECT https://flowcon.cnrs.fr/ DECIPHER @LIMSI https://www.limsi.fr/en/research/aero/menuitem-aero-topics-en/aero-topics-decipher-en We tested Deep reinforcement learning for the optimal non-linear control of a chaotic system [1], using an actor-critic algorithm, the Deep Deterministic Policy Gradient.

-METHODS II: REINFORCEMENT LEARNING

-REFERENCES

• Full knowledge of the system is not required.

• In principle, the policy is a global optimum if the cost function is solution of the Bellman equation.

• Current work: 2D numerical cavity flow and 3D experimental counterpart (fig. 5).

The critical parameter for the KS equation is the domain extent; we choose , corresponding to a chaotic regime. The equation reads and it is numerically solved in a periodic domain. This regime is characterised by 4 equilibria (Fig. 2).

The controller is designed by solely relying on the measurements taken from 8 localised sensors (Fig. 3); the controller acts introducing a forcing ( ) on the system by means of 4 localised actuators.

The entire control design is meant to be as realistic as possible for mimicking realistic conditions (i.e. localised actuations and measurements).

L = 22 ∂x ∂t = -∇ 4 x -∇ 2 x - 1 2 | ∇x | 2 +g g
Three policies are designed, each of which minimising the distance between the current state and the target state, represented by one of the nontrivial equilibria (Fig. 2). The reward is defined using the 2-norm and reads with the discount factor equals .

The training of the control policies takes on average one hour of computation on a standard CPU.

In Fig. 4 we show how in each of the three cases the policies are capable at driving and stabilising the dynamics around the unstable points. The polices are robust with respect of the initial conditions as shown in [1].

r(t) = -∥x(t) -E i ∥ γ = 0.99
Reinforcement learning (RL) is based on the approximation of the Hamilton-Jacobi-Bellman (HJB) equation is the value function, the reward, the system, the control action, and is the discount factor. The discrete version, the Bellman equation leads to different approaches.

Actor-only: one policy is tested on a long-time trajectory .

The solution satisfies the Pontryagin maximum principle, a necessary solution of optimality.

is related to the discount factor.

Critic-only: the problem is decomposed in local problems, each associated with a policy/action

.

ρ𝒥(x(t), t) = max u {r(x, u) + ∇𝒥(x(t), t)f(x, u)} 𝒥 r f u ρ > 0 π 𝒥 π (x n ) = r(x, u) + γ𝒥 π (x n+1 ) γ = e -ρt u Q(x n , u n ) = r(x, u) + γQ(x n+1 , u n+1 )
The function is the value function when is function of the state and the action . In this case the Bellman principle is satisfied, which is a necessary and sufficient condition for the optimality.

We use the Deep Determinist Policy Gradient (DDPG) combining the actor and critic strategies (Fig. 3). The control policy and the function are approximated by using Neural Networks (NN).

The optimisation of the NN is done by stochastic gradient descent with adaptive moments (ADAM) and the partial-observed Markov Decision Process (PO-MDP) is iteratively stacked in memory and used to reduce the temporal difference .

Q-

x u π Q- TD = Q n -(r + γQ n+1 )

[ 1 ]

 1 Bucci, M.A. et al., (2019) arxiv.org/abs/1906.07672 [2] Silver, D. et al., (2014, June), In ICML. [3] Sipp and Schmid, (2016), Appl. Mech. Rev. 68 -20801

Fig. 3 -Fig. 1 -Fig. 5 -

 315 Fig. 3 -In the RL framework, an agent interacts with an environment by making observations and performing actions . In return, the agent receives a reward that depends directly on the changes of the environment induced by the action. The control policy and the temporal difference are approximated by means of neural networks.u

Fig. 2 -Fig. 4 -

 24 Fig. 2 -Dynamics of the Kuramoto-Sivashinsky (KS) equation, for the domain extent . The null equilibrium and the three nontrivial solutions are shown at (leftmost in the insets). L = 22 t = 0


G De Cillis

S Cherubini

O Semeraro

S Leonardi

P De Palma

Harnessing wind energy and converting it in electric energy require the design and installation of large wind farms, constituted by hundreds of turbines. When grouped together, a great part of the wind turbines operates in the wake of upwind turbines, so that the velocity deficit and the turbulence level of the incoming flow induce additional power losses and fatigue blade loading. Thus, understanding the dynamics of wind turbine wakes is crucial for the design and efficiency improvement of wind turbines and farms. The wake of a wind turbine is characterized by the coexistence of several coherent structures, such as tip and root vortices, tower's vortex shedding etc. High turbulence levels promote wake recovery,

••

  , however how individual coherent structures contribute to wake recovery, is still an open problem.INTRODUCTION OBJECTIVES Large eddy simulation (LES)-Smagorinsky moded with 𝐶 𝑆 = 0.09.-Actuator line model (ALM)[START_REF] Sørensen | Proc. of European Wind Energy Conference EWEC '99[END_REF] for rotor blades.-Immersed boundaries[START_REF] Orlandi | [END_REF] for tower and nacelle.-Domain dimensions: 12.5𝐷 × 5𝐷 × 3𝐷 in the streamwise (x), vertical (y) and transverse (z) directions.-Computational grid: 2048 × 512 × 512 grid-points in x,y and z directions.-Krogstad model wind turbine[START_REF] Krogstad | Renewable Energy[END_REF] at tip-speed ratio 𝜆 = 3.-Diameter-based 𝑅𝑒 = 6.5 × 10 5 .-Uniform, laminar inlet velocity profile . • Proper orthogonal decomposition (POD) -M = 2361 snapshots of the velocity field. -Sampling frequency: 1 snapshot every 10°rotation. -SVD method. The most energetic POD modes are associated with the tip and root vortices and their harmonics. Then, two pair of modes connected with the von Karman vortices shed by the tower follow. Other energetic modes show low-frequency oscillations that could be related tothe wake-meandering phenomenon. • Positive values of M.K.E. flux were found for the POD modes linked to the tower's vortex shedding, to the Kelvin-Helmholtz instability and to the wake meandering, which have, therefore, a beneficial effect on the wake recovery. • POD modes related to the tip vortices show negative M.K.E. fluxes, therefore they sustain the wake and slow down its recovery, despite producing a high turbulence intensity. • The present analysis suggests that neglecting the presence of tower and nacelle in the simulation of clustered turbines can lead to an underprediction of the energy production. The M.K.E. flux entering a cylindrical surface 𝑆 of radius 𝑅 = 0.5 provided by each mode is computed as: where ത 𝑢 𝑖 is i-th component of the time-averaged velocity, 𝜆 𝑘 is the eigenvalue related to the k-th mode, 𝜙 𝑖 𝑘 is the i-th component of the k-th POD mode.

CONTACT

Giovanni De Cillis, DMMM, Politecnico di Bari Via Re David 200, 70125 Bari, Italy E-mail: giovanni.decillis@poliba.it

• Identifying and isolating the main flow features characterizing the wake dynamics of a threebladed wind turbine using a three-dimensional Proper Orthogonal Decomposition (POD) of numerical data obtained by large eddy simulation (LES). • Evaluating how each POD mode contributes to wake recovery through mean kinetic energy entrainment.