Chunhua Zhang

Xiaoyu Zhang

Yuheng Fu

Lin Zhang

Alexander Kuhn1
email: kuhn@enscbp.fr

Metal-Organic Framework functionalized bipolar electrodes for bulk electroenzymatic synthesis

Keywords: NADH regeneration, bipolar electrochemistry, metal organic frameworks, electrocatalysis, electroenzymatic synthesis

Electrochemical regeneration of the cofactor nicotinamide adenine dinucleotide (NADH) is crucial for the enzymatic synthesis of a variety of valuable chemicals. [Rh(Cp*)(bpy)Cl]+ has been reported as one of the most efficient non-enzymatic catalysts for electrochemical regeneration of NADH, but the efficiency of the existing systems is still far from satisfying. Here we propose to suspend metalorganic framework functionalized glassy carbon beads in a bipolar electrochemical cell under mechanical stirring for bulk electroenzymatic synthesis. The beads are modified with a metal-organic framework NU-1000 shell, in which [Rh(Cp*)(bpy)Cl]+ catalyst is incorporated via the solvent-assisted ligand incorporation approach. Enzymatically active 1,4-NADH is electrochemically regenerated, and finally bipolar electrochemical bulk conversion of pyruvate into L-lactate is performed in the presence of NADH dependent L-lactate dehydrogenase as a model system. The very high turnover frequency (TOF) of 3100±106 h-1 opens up promising perspectives for employing this concept also for various other electroenzymatic bulk processes.

Introduction

Enzymes are biocatalysts with high specificity and catalytic efficiency, among which NADHdependent oxidoreductases constitute one of the largest class of enzymes for the synthesis of many valuable chiral products [START_REF] Macedo | Three-dimensional catalysis and the efficient bioelectrocatalysis beyond surface chemistry[END_REF][START_REF] Sheldon | The Hitchhiker's guide to biocatalysis: Recent advances in the use of enzymes in organic synthesis[END_REF]. In this context, the regeneration of the cofactor 1,4-NADH is essential for large scale enzymatic synthesis, since 1,4-NADH is always required as an electron donor for the oxidoreductase to catalyze the enzymatic reaction. However, the high price of NADH imposes a recycling of the oxidized form NAD+ back to NADH [START_REF] Mordhorst | Round, round we go -strategies for enzymatic cofactor regeneration[END_REF][START_REF] Wang | Cofactor NAD(P)H Regeneration Inspired by Heterogeneous Pathways[END_REF]. Among different regeneration methods, the electrochemical regeneration of 1,4-NADH in combination with the enzymatic reactions has been considered as a green and sustainable approach, and different electroenzymatic systems have been developed for the synthesis of a variety of chiral products in recent years [START_REF] Steckhan | Electroenzymatic Synthesis[END_REF][START_REF] Chen | The progress and outlook of bioelectrocatalysis for the production of chemicals, fuels and materials[END_REF][START_REF] Lee | Advances in electrochemical cofactor regeneration: enzymatic and non-enzymatic approaches[END_REF][START_REF] Wu | Exploring Ferredoxin-Dependent Glutamate Synthase as an Enzymatic Bioelectrocatalyst[END_REF][START_REF] Wu | Enzymatic electrosynthesis as an emerging electrochemical synthesis platform[END_REF][START_REF] Wu | Methods for the regeneration of nicotinamide coenzymes[END_REF]. Since the direct electrochemical reduction of NAD+ requires high overpotentials, enzymatically inactive side products such as NAD2 dimers and other isomers of NADH are formed. Thus, it is necessary to use a catalyst for the regeneration to obtain selectively enzymatically active 1,4-NADH [START_REF] Steckhan | Electroenzymatic Synthesis[END_REF][START_REF] Gorton | Electrochemistry of NAD(P)+/NAD(P)H[END_REF]. Up to now, in addition to the possibility to use an enzymatic catalyst (eg. diaphorase) [START_REF] Lee | Advances in electrochemical cofactor regeneration: enzymatic and non-enzymatic approaches[END_REF][START_REF] Yuan | Efficient NADH Regeneration by a Redox Polymer-Immobilized Enzymatic System[END_REF], [Rh(Cp*)(bpy)Cl]+ and its derivatives are reported to be among the most efficient and robust non-enzymatic catalysts for NADH regeneration [START_REF] Steckhan | Analytical Study of a Series of Substituted(2,2'-Bipyridyl) (pentamethylcyclopentadienyl)rhodium and -iridium Complexes with Regard to Their Effectiveness as Redox Catalysts for the Indirect Electrochemical and Chemical Reduction of NAD (P)+[END_REF][START_REF] Ganesan | Efficient Nicotinamide Adenine Dinucleotide Regeneration with a Rhodium-Carbene Catalyst and Isolation of a Hydride Intermediate[END_REF]. The classical strategy to achieve efficient electrochemical regeneration of NADH is to immobilize the [Rh(Cp*)(bpy)Cl]+ catalyst on porous electrodes, in order to increase the loading of catalyst. Different porous electrodes (carbon fibers [START_REF] Zhang | Covalent immobilization of (2,2'-Bipyridyl) (Pentamethylcyclopentadienyl)-Rhodium Complex on a Porous Carbon Electrode for Efficient Electrocatalytic NADH Regeneration[END_REF][START_REF] Zhang | Molecular and Biological Catalysts Coimmobilization on Electrode by Combining Diazonium Electrografting and Sequential Click Chemistry[END_REF], carbon nanotubes [START_REF] Zhang | Electrocatalytic Biosynthesis using a Bucky Paper Functionalized by [Cp*Rh(bpy)Cl]+and a Renewable Enzymatic Layer[END_REF][START_REF] Tan | Regeneration of the NADH Cofactor by a Rhodium Complex Immobilized on Multi-Walled Carbon Nanotubes[END_REF], metal organic frameworks [START_REF] Li | Fine-Tuning the Electrocatalytic Regeneration of NADH Cofactor Using [Rh(Cp*)(bpy)Cl]+-Functionalized Metal-Organic Framework Films[END_REF]) were functionalized with [Rh(Cp*)(bpy)Cl]+ catalyst and applied for NADH regeneration. In a recently reported work, a NU-1000 metal-organic framework film was used as supporting matrices on a working electrode, and [Rh(Cp*)(bpy)Cl]+ entities were regularly positioned on the Zr nodes of NU-1000. The so-obtained electrode has been shown to be very efficient for NADH regeneration with a high TOF of 1400 h-1 for 1,4-NADH regeneration, indicating that the metal-organic framework film is an interesting support for designing efficient electrocatalytic systems [START_REF] Li | Fine-Tuning the Electrocatalytic Regeneration of NADH Cofactor Using [Rh(Cp*)(bpy)Cl]+-Functionalized Metal-Organic Framework Films[END_REF].

In addition to normal electrochemistry, bipolar electrochemistry provides the possibility to further improve the efficiency of the system. According to the principles of bipolar electrochemistry, multiple electrodes can be placed between two feeder electrodes and polarized under the influence of an electric field. In this case, electrochemical reactions will be triggered at the two extremities (polarized anode δ+ and cathode δ-) of each individual electrode, if the electric field exceeds a certain threshold value [START_REF] Bouffier | Bipolar electrochemistry[END_REF][START_REF] Shida | Bipolar Electrochemistry: A Powerful Tool for Electrifying Functional Material Synthesis[END_REF][START_REF] Koefoed | Bipolar electrochemistry-A wireless approach for electrode reactions[END_REF][START_REF] Rahn | Recent Advancements in Bipolar Electrochemical Methods of Analysis[END_REF][START_REF] Fosdick | Bipolar electrochemistry[END_REF]. In our previous work, several [Rh(Cp*)(bpy)Cl]+ functionalized carbon beads were suspended in the reaction compartment of a bipolar electrochemical cell with the help of viscous cellulose as electrolyte. The 3D assembly of these electrodes in the bulk solution allows a fine-tuning of the efficiency of the electrochemical regeneration of 1,4-NADH [START_REF] Zhang | Bulk electrocatalytic NADH cofactor regeneration with bipolar electrochemistry[END_REF]. However, if this bipolar electrochemical regeneration system is coupled with an enzymatic follow-up reaction, there are still two important problems which need to be solved. First, as the electron transfer between the enzymes and the Rh catalyst is achieved by shuttling 1,4-NADH, the diffusion of this cofactor can become the rate limiting step in a viscous electrolyte. Second, it has been reported that some functional groups (eg.-SH, -NH2), present at the surface of enzymes, can react with the metal catalyst, leading to a degradation of the Rh complex activity [START_REF] Walcarius | Factors affecting the electrochemical regeneration of NADH by (2,2′bipyridyl) (pentamethylcyclopentadienyl)-rhodium complexes: Impact on their immobilization onto electrode surfaces[END_REF]. Thus a physical separation between the Rh complex and the enzymes is necessary during long-term electroenzymatic synthesis [START_REF] Hildebrand | Stable Electroenzymatic Processes by Catalyst Separation[END_REF].

Taking into account the above mentioned two limitations, we propose in this work a bipolar electroenzymatic system based on multiple glassy carbon beads, functionalized with a NU-1000 shell containing [Rh(Cp*)(bpy)Cl]+ catalyst, incorporated at its Zr nodes. On the one hand, the NU-1000 layer is expected to increase the local concentration of Rh catalyst on the surface of the beads, leading to an improvement of the bipolar electrocatalytic efficiency for NADH regeneration. On the other hand, the NU-1000 shell can be used as a physical separation matrix between the Rh catalyst and the enzymes. The pore size of NU-1000 being 3.4 nm, the large NADH-dependent dehydrogenases will be kept outside the pores of NU-1000, while the NADH molecules (1.3 1.1 1.9 nm) are small enough to diffuse into the pores and shuttle electrons between the Rh catalyst and the enzymes. Thus, the degradation of the Rh catalyst and enzyme during long-term electrochemical synthesis can be avoided.

The general concept of this work is shown in Scheme 1. In order to enhance mass transport, the beads are kept in suspension in the reaction compartment of a bipolar cell by simple mechanical stirring. The bipolar electrocatalytic NADH regeneration, catalyzed by immobilized Rh catalyst, is achieved at the cathodically polarized extremity (δ-), and the sacrificial electrooxidation of ascorbic acid occurs at the anodically polarized side (δ+). The efficiency of the bipolar electrocatalytic NADH regeneration scales with the number of suspended beads. Finally, L-lactate dehydrogenase (with dimensions of 13.6 8.6 6.0 nm [START_REF] Ei-Zahab | Enabling multienzyme biocatalysis using nanoporous materials[END_REF]) is added to the solution, and the bipolar electrochemical bulk conversion of pyruvate into L-lactate is carried out under stirring with this model system. Scheme 1. Illustration of (A) the experimental set-up for the bulk bipolar electroenzymatic synthesis and (B) the catalytic cascade with NU-1000-Rh for the bipolar electroenzymatic synthesis of lactate in the presence of NADH-dependent l-lactate dehydrogenase.

Results and discussion

The glassy carbon beads with a NU-1000 shell are prepared following the process illustrated in Scheme 2. The glassy carbon beads (CB) are first polished with sand paper in order to obtain a rough surface, then the beads are functionalized with carboxyl groups via chemical diazonium grafting. The NU-1000 shell is grown on the surface of CB via a solvothermal synthesis step, followed by a solventassisted ligand incorporation process for the [Rh(Cp*)(bpy)Cl]+ immobilization at the nodes of the NU-1000 shell [START_REF] Li | Fine-Tuning the Electrocatalytic Regeneration of NADH Cofactor Using [Rh(Cp*)(bpy)Cl]+-Functionalized Metal-Organic Framework Films[END_REF]. The crystal structure and the chemical analysis of the obtained CB-NU1000-Rh are supported by XRD (Figure S1) and XPS (Figure S2) measurements. Furthermore, EDS measurements (Figure S2D) indicate an atom% of Zr and Rh of 85.6% and 14.4%, respectively. This corresponds to an atomic ratio of 5.94 for Zr to Rh, which means that, within error bars, there is statistically one Rh catalyst per Zr6 node. The CB-NU1000-Rh catalyst particles are then placed in the reaction compartment for bipolar electrochemical NADH regeneration, and the amount of generated enzymatically active 1,4-NADH is measured. In this work we employed a combination of UV-vis spectroscopy and biological assays in order to correctly measure the produced 1,4-NADH. UV/Vis spectroscopy for characterizing 1,4-NADH has been already reported in the literature. 1,4-NADH has a typical absorption peak at 340 nm, whereas possible side products such as NAD2, 1,2-NADH, and 1,6-NADH have different characteristic absorption peaks. To further verify that the absorption peak at 340 nm is originating from the enzymatically active 1,4-NADH, rather than inactive 1,6-NADH or NAD2, we used a direct enzymatic assay, which consists of the conversion of acetaldehyde by an alcohol dehydrogenase (ADH). This reaction can only occur in the presence of enzymatically active 1,4-NADH. Thus, the complete disappearance of the absorption peak at 340 nm indicates that the formed product is 1,4-NADH (details are described in the experimental part). First, forty CB-NU1000-Rh beads are placed at the bottom of the bipolar reaction compartment (assembled in a 2D array), and electric fields with different strengths are applied to trigger the bipolar electrochemical reduction of NAD+ at the δextremity of each bead. The product generated after 1.5 h is characterized by UV/vis spectroscopy (Figure S3). As shown in Figure 2A, the concentration of produced 1,4-NADH is almost zero for an electric field of 10 V/cm, however the concentration of 1,4-NADH starts to be measurable for an electric field of 12 V/cm and further increases for 12.5 V/cm, whereas a decrease is observed for even higher electric fields (13.5 V/cm). This trend can be understood by calculating for the different electric fields the fractional surface area of the carbon beads where the bipolar electrochemical reaction can take place due to a sufficient polarization (see details in supporting information). The conversion barely starts at an electric field of 10 V/cm with a very small cathodically polarized surface area, and gradually increases when the electric field is raised to 12.5 V/cm. For even higher fields (13.5 V/cm) the direct non-catalytic electrochemical reduction of NAD+ starts, leading to the partial formation of enzymatically inactive side product. Since the polarization potential difference is proportional to the effective length of bipolar electrodes, 12.5 V/cm has been chosen as the optimum electric field for the electrochemical regeneration of 1,4-NADH for beads with a characteristic size of 960±20 µm. A long-term bipolar electrochemical transformation experiment has been carried out as shown in Figure 2B. The typical absorption peak of 1,4-NADH at 340 nm is increasing as a function of time, corresponding to a linear increase in concentration up to 6 h (the inset of Figure 2B). The same peak completely disappears during the enzymatic test (Figure S4A) indicating the high selectivity of this system during long-term bipolar electrochemical synthesis. A blank experiment with forty glassy carbon beads, covalently modified by diazonium grafting with Rh catalyst was carried out under the same experimental conditions (Figure S4B). UV/Vis spectra were used to quantify the bipolar electrochemical product in both cases, and a roughly two times higher amount of 1,4-NADH was obtained in the presence of the NU-1000 shell on the beads, indicating that the NU-1000 matrix is playing a significant role for improving the electrocatalytic efficiency for NADH regeneration. Before combining the bipolar NADH regeneration system with enzymatic reactions, we first tested the possibility to suspend the beads in the bulk solution by mechanical stirring. This can be achieved by stirring with a magnetic bar at 1400 rpm (Video S1). The bipolar electrochemical regeneration of NADH is carried out with forty CB-NU1000-Rh microelectrodes in the 3D stirring mode, and the obtained 1,4-NADH is characterized by UV-Vis spectra. As shown in Figure 3A, the absorption peak at 340 nm increases from 0.2 to 0.6, almost three times more compared to the previous experiment with stationary beads in the bipolar cell. In the latter mode, the beads sink to the bottom, where only a limited number (maximum sixty) can be placed without mechanical or electrical interference. A physical interaction between beads should be avoided, because a contact between them will increase their effective length, leading to the formation of side products caused by a higher polarization potential difference. In contrast, under stirring, over one hundred CB-NU1000-Rh beads can be accommodated in the 3D reaction space before physical interactions become statistically significant.

The bipolar electrochemical regeneration of NADH with different numbers of CB-NU1000-Rh electrodes was performed both in the 2D and 3D mode, and the efficiency of both systems was evaluated by characterizing the concentration of 1,4-NADH after 1.5 h. Figure 3B reports the concentration of produced 1,4-NADH as a function of the number of CB-NU1000-Rh electrodes for the 2D and 3D system. The enzymatic activity of 1,4-NADH was again confirmed by addition of alcohol dehydrogenase and acetaldehyde as shown in Figure S5 and Figure S6. For the 2D system, a linear increase of 1,4-NADH concentration is observed up to 40 beads, which is expected due to the increase of total active surface area of the electrodes in the cell. Then a slight decrease is observed with 60 microbeads, which can be explained by diffusion restrictions due to smaller distance between the individual electrodes. Similarly, a linear increase is also observed with the stirred system, and the saturation occurs for around 80 beads, followed by a gradual decrease with an increasing number of electrodes. The slope of the linear range for the 3D system is higher than for the 2D system, indicating a faster reaction kinetics, which can be attributed to an improved mass transport. In addition, the reaction is also promoted by the movement and rotation of the CB-NU1000-Rh beads. When the beads rotate during the bipolar electrochemical process, the effective length of the beads with respect to the electric field remains constant due to their isotropic shape, but the region on their surface where the reduction takes place changes continuously. Thus, the educt and the product have enough time to diffuse in/out of the catalytic layer during the rotation to a new position where fresh catalyst surface becomes activated by the electric field. Compared to our previously reported NADH regeneration system with a gellified medium [START_REF] Zhang | Bulk electrocatalytic NADH cofactor regeneration with bipolar electrochemistry[END_REF], that we used to maintain bipolar beads suspended in the 3D reaction space, the here employed liquid medium allows continuous stirring and therefore a considerably improved mass transport. This gives the possibility to couple this coenzyme regeneration system with an electroenzymatic follow-up reaction in order to perform the conversion of a substrate with high efficiency. Finally, 100U L-lactate dehydrogenase (dimensions of 13.6 8.6 6.0 nm), 2 mM pyruvate, and 0.5 mM NAD+ have been added to the solution for the bipolar electroenzymatic bulk conversion from pyruvate to L-lactate as a model reaction. 5 mM ascorbic acid has also been added as a sacrificial reactant for the electrooxidation at the anodically polarized face of the beads. The bipolar electrochemical bulk conversion is carried out under stirring with eighty CB-NU1000-Rh beads with an electric field of 12.5 V/cm for 4 h, and the product is analyzed by HPLC (Figure S7). According to the calibration curve of L-lactate, (2.23 0.08) 10-6 mol L-lactate is produced, the quantity of [Rh(Cp*)(bpy)Cl]+ catalysts on each CB-NU1000-Rh microbead is estimated to be 2.25 10-12 mol, and thus the corresponding turnover frequency (TOF) for [Rh(Cp*)(bpy)Cl]+ can be calculated to be 3100 106 h-1 (see calculation details in supporting information). To the best of our knowledge, this is the highest TOF value reported so far for a non-enzymatic catalyst in electroenzymatic synthesis, very close to what has been obtained when using an enzyme as a catalyst for NADH regeneration (3680 h-1 with diaphorase [START_REF] Yuan | Efficient NADH Regeneration by a Redox Polymer-Immobilized Enzymatic System[END_REF], and 3600 h-1 with MeFDH1-β [START_REF] Lee | Communication-Highly Efficient Electroenzymatic NADH Regeneration by an Electron-Relay Flavoenzyme[END_REF]). Last-but-not-least, the catalytic layer on the beads seems to be rather stable, as indicated by XRD and EDX measurements performed after the catalysis (Figure S1).

Conclusions

In conclusion, we propose a new strategy for electroenzymatic synthesis, based on bulk bipolar electrochemical NADH regeneration, carried out with suspended functionalized carbon beads. A MOF shell (NU-1000), grown on the surface of glassy carbon microbeads, is post-functionalized with [Rh(Cp*)(bpy)Cl]+ at the zirconium nodes. The resulting beads are then exposed in solution to an electric field, and their cathodically polarized face δ-can be used for triggering the electrocatalytic regeneration of 1,4-NADH catalyzed by the immobilized Rh catalyst. Electrochemical bulk regeneration of NADH is achieved when suspending the CB-NU1000-Rh beads by stirring, and the efficiency of the system as a function of the number of beads is evaluated. The proposed approach shows a much higher efficiency compared to other systems (2D ensembles of beads or their suspension in a viscous medium). Coupling this NADH regeneration strategy with L-lactate dehydrogenase as an enzymatic model system to transform in situ pyruvate into L-lactate, reveals a very high TOF of 3100 106 h-1, thus making this approach very attractive for the bulk electroenzymatic synthesis of other high added-value products.

Scheme 2 .

 2 Scheme 2. Elaboration process of a [Rh(Cp*)(bpy)Cl]+ functionalized NU-1000 layer on a glassy carbon bead. The process includes a chemical diazonium grafting step for carboxyl group functionalization, a solvothermal step for the growth of NU-1000, and a solvent-assisted ligand incorporation step for [Rh(Cp*)(bpy)Cl]+ immobilization.

Figure 1 .

 1 Figure 1. SEM images of a CB-NU1000-Rh electrode recorded with different magnifications and the corresponding elemental mapping.

Figure 2 .

 2 Figure 2. (A) Bipolar electrochemically produced enzymatically active 1,4-NADH concentration as a function of the applied electric field after 1.5 h; (B) UV-Vis absorption spectra of the bipolar electrochemical product obtained with an electric field of 12.5 V/cm after (a) 1 h, (b) 2 h, (c) 3 h, (d) 4 h, (e) 5 h, and (f) 6 h. The inset indicates the corresponding concentration of 1,4-NADH as a function of time. Experiments were carried out with 40 CB-NU1000-Rh electrodes (960±20 µm) placed at the bottom of the bipolar cell containing 10 mL Tris buffer at pH 7.2 in the presence of 5 mM ascorbic acid and 1 mM NAD+ under nitrogen atmosphere. Error bars refer to the standard deviations of three measurements.

Figure 3 .

 3 Figure 3. (A) UV-Vis absorption spectra of the bipolar electrochemical product obtained with fourty CB-NU1000-Rh electrodes (a) placed at the bottom of the bipolar cell (b) suspended in the bulk phase of the bipolar cell under stirring (1400 rpm). (B) Bipolar electrochemically produced 1,4-NADH concentration as a function of the number of CB-NU1000-Rh beads in the 2D mode (black) and in the 3D mode (red). Experiments were carried out with an electric field of 12.5 V/cm for 1.5 h in 5 mL Tris pH 7.2 in the presence of 5 mM ascorbic acid and 1 mM NAD+ under nitrogen atmosphere. Error bars refer to the standard deviations of three measurements.

Acknowledgement

The work has been funded by the National Natural Science Foundation of China (No. 22272044, No.21902045) and the Henan center for Outstanding Overseas Scientist (GZS2020005, HNGD2022051). The project has also been supported by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no 741251, ERC Advanced grant ELECTRA). L. Z. acknowledges Prof. Wensheng Yang for very fruitful scientific discussions and support via the Postgraduate Cultivating Innovation and Quality Improvement Action Plan of Henan University (SYLYC2022157).


Jérémie Labasse

Uwe Ehrenstein
email: uwe.ehrenstein@univ-amu.fr

Guillaume Fasse

Frédéric Hauville

Thrust scaling for a large-amplitude heaving and pitching foil with application to cycloidal propulsion

Keywords: Pitching and heaving foil, Thrust, Chimera method, Cycloidal propulsion

A numerical solution procedure using the mesh-superposition approach, known as the Chimera method, together with the OpenFOAM toolbox environment is used to compute the forces generated by large amplitude heaving and pitching foil. The possibility of fitting thrust prediction laws, based on classical potential flow theories, with the numerically computed forces is explored, for a Reynolds number of 5 10 4 . It is shown, first for a pure heaving motion and subsequently by adding a harmonic pitching motion, that theoretical scaling may be fitted to numerical time-averaged thrust data, even in the case of large amplitude motions. The thrust-prediction law is shown to still apply to pitching-rotating motions, such as those of blades in cycloidal propulsion devices, the mean pressure correction due to the additional surging motion being small. The synchronized rotation-pitching of three foils typical of a cross-flow propeller configuration is addressed as well. The numerical global thrust results are shown to be in general agreement with the theoretical prediction, but also with blade-embedded load cell measurements for an exper

Introduction

Oscillating bodies in moving fluids are intimately related to swimming and flying and the question of thrust performance in this context has received considerable attention in the literature [START_REF] Triantafyllou | Review of hydrodynamic scaling laws in aquatic locomotion and fishlike swimming[END_REF]. Pitching and heaving motions are often considered as archetypes of aquatic locomotion and it is the reason why they have widely been considered for thrust scaling investigations. It has been recognized that flapping wings' performance may be inferred from the wake structure of the wings, which exhibits in the case of thrust production jet-like average velocity profiles, associated with a reversed Kármán vortex street [START_REF] Von Kármán | Aerodynamic Theory: General aerodynamic theory : Perfect fluids[END_REF]. Applying for instance the integral momentum theory to a control volume in which the oscillating body is embedded, the generated thrust can be predicted by velocity measurements in the wake [START_REF] Koochesfahani | Vortical patterns in the wake of an oscillating airfoil[END_REF][START_REF] Bohl | Mtv measurements of the vortical field in the wake of an airfoil oscillating at high reduced frequency[END_REF], the reliability of this approach having been recently addressed in [START_REF] Moubogha Moubogha | Forces on a pitching plate: an experimental and numerical study[END_REF][START_REF] Brousseau | Fluid-structure interaction effects on the deformable and pitching plate dynamics in a fluid flow[END_REF], comparing measurements with numerical simulation results. Theoretical models for thrust prediction based on potential flow theory have been proposed since the seminal work by Garrick [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF], for swimming plates in [START_REF] Wu | Swimming of a waving plate[END_REF] and more recently for example in [START_REF] Ramananarivo | Passive elastic mechanism to mimic fish-muscle action in anguilliform swimming[END_REF][START_REF] Piñeirua | Resistive thrust production can be as crucial as added mass mechanisms for inertial undulatory swimmers[END_REF]. The formulas derived in [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF] have been compared with numerical thrust force results, see for example [START_REF] Jones | Numerical computation of flapping-wing propulsion and power extraction[END_REF][START_REF] Moored | Inviscid scaling laws of a self-propelled pitching airfoil[END_REF][START_REF] Labasse | Numerical exploration of the pitching plate parameter space with application to thrust scaling[END_REF], as well as with direct measurements (e.g. [START_REF] Mackowski | Direct measurement of thrust and efficiency of an airfoil undergoing pure pitching[END_REF]). Scaling laws for pitching and heaving foils have also been compared with experimental data in [START_REF] Anderson | Oscillating foils of high propulsive efficiency[END_REF][START_REF] Floryan | Scaling the propulsive performance of heaving and pitching foils[END_REF], among others. But despite the numerous numerical and experimental investigations on the hydrodynamics of flapping foils available (see [START_REF] Wu | A review on fluid dynamics of flapping foils[END_REF] for a recent review), there is no general consensus about thrust models, in particular when moderate or large flapping amplitudes are required, as for example in the field of renewable power generation [START_REF] Abbaszadeh | Experimental investigation of an optimised pitch control for a verti cal-axis turbine[END_REF] and naval propulsion [START_REF] Bartels | The voith schneider propeller: current applications and new developments[END_REF]. Indeed, large-amplitude oscillations of foils appear to be necessary to achieve efficient thrust magnitudes (e.g. [START_REF] Scherer | Experimental and theoretical investigation of large amplitude oscillating foil propulsion systems, HYDRONAUTICS, Incorporated[END_REF]). The blades of a cycloidal propeller for instance undergo large-amplitude pitching-rotating motions and the associate performance of such devices has been addressed numerically in [START_REF] Hu | Two-dimensional simulation of the hydrodynamic performance of a cycloidal propeller[END_REF], while a low-order hydrodynamic model has been proposed and experimentally validated in [START_REF] Halder | Hydrodynamic modeling and experimental validation of a cycloidal propeller[END_REF]. Cycloidal fluid flow engines are also designed for tidal energy conversion and a Darrieu-type cross flow turbine has for instance been addressed in [START_REF] Paillard | Uranse simulation of an active variable-pitch cross-flow darrieus tidal turbine: Sinusoidal pitch function investigation[END_REF]. Concerning scaling laws, there have been attempts (see e.g. [START_REF] Chopra | Large amplitude lunate-tail theory of fish locomotion[END_REF]) to extent the linear potential theory to foils undergoing large amplitude oscillations, experimental results having been compared with linear and nonlinear theory [START_REF] Anderson | Oscillating foils of high propulsive efficiency[END_REF]. In [START_REF] Floryan | Large-amplitude oscillations of foils for efficient propulsion[END_REF], scaling relations have been proposed in particular for heaving amplitudes of the order of magnitude of the foil's chord. When addressing numerically large amplitude oscillations of a foil, a key issue is to efficiently update the computational mesh. An alternative which has become widely used is to immerse the moving geometry on a simple Cartesian grid using more or less ad hoc momentum sources to model the kinematic condition at the interface between the velocity field and the body's moving boundary. This approach is known as the immersed boundary method and goes back to Peskin's seminal work [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF]. Even though this approach may be used for different Navier-Stokes discretization procedures as well as all types of structures (see e.g. [START_REF] Pinelli | Immersedboundary methods for general finite-difference and finite-volume navier-stokes solvers[END_REF][START_REF] Goza | A strongly-coupled immersed-boundary formulation for thin elastic strucures[END_REF]), a difficulty associated with this method is however to accurately resolve the very vicinity of the body surface, in particular in the case of thin structures [START_REF] Kallemov | An immersed boundary method for rigid bodies[END_REF]. When large amplitude motions are addressed, the so-called overset grid approach, also known as the Chimera method, appears to be particularly appropriate as a body-fitted discretization procedure. The idea is to superimpose moving subdomains to a fixed meshed background global computational domain. The independently moving domains contain the structures in motion and the communication with the fixed mesh is achieved through interpolation within transmission zones. This method has first been documented in [START_REF] Steger | On the use of composite grid schemes in computational aerodynamics[END_REF] and more recent analyses are given in [START_REF] Tang | An overset-grid method for 3d unsteady incompressible flows[END_REF][START_REF] Houzeaux | A chimera method for the incompressible navier-stokes equations[END_REF][START_REF] Roget | Robust and efficient overset grid assembly for partitioned unstructured meshes[END_REF][START_REF] Ye | Benchmark computations for flows around a stationary cylinder with high reynolds numbers by ransoverset grid approach[END_REF], for instance. The open-source OpenFOAM toolbox has become very popular in academic research [START_REF] Tabor | Inlet conditions for large eddy simulation: A review[END_REF][START_REF] Meldi | Quantification of errors in largeeddy simulations of a spatially evolving mixing layer using polynomial chaos[END_REF][START_REF] Lysenko | Modeling of turbulent separated flows using OpenFOAM[END_REF] as well as industrial flow simulations [START_REF] Gao | Numerical simulation of turbulent flow past airfoils on OpenFOAM[END_REF][START_REF] Selma | Optimization of an industrial heat exchanger using an open-source cfd code[END_REF], the Chimera-type approach having been implemented in several version of this numerical simulation environment and in particular in [40] (for a description of the overset-option, see for instance [START_REF] Tisovska | Description of the overset mesh approach in esi version of openfoam[END_REF]). The OpenFOAM toolbox and the overset-option is used in the present investigation, addressing the question of thrust-scaling for a large amplitude heaving foil as well as for a heaving-pitching foil. Heaving, pitching and surging are the three motions into which the displacement of pitching foils in cross-flow propellers can be decomposed. Considering such a propeller, which is part of an academic experimental project at the French Naval Academy Research Institute (IRENav) [START_REF] Fasse | An experimental blade-controlled platform for the design of smart cross-flow propeller[END_REF], simulations are performed for rotatingpitching foils for parameters close to those in the experiment. The paper is organized as follows. A heaving foil is addressed in Section 2 using the OpenFOAM simulation environment, the numerical procedure using the overset-approach being presented and the numerical convergence is assessed as well. The forces generated by the foil's motion are computed for a large set of heaving parameters, allowing a successful fit with a thrust scaling law. Subsequently in Section 3 a pitching motion is added to heaving, the aim being again to propose a thrust prediction law, as function of the parameters. In Section 4 the foil's motion is made more complex adding a surging back-and-forth movement, such as to mimic motions in cross-flow propeller devices. In particular, the reliability of a mere heaving-pitching thrust law in this context is examined. Finally, dynamic force computations are compared with experimental measurements for the cycloidal propulsion device described in [START_REF] Fasse | An experimental blade-controlled platform for the design of smart cross-flow propeller[END_REF]. The results are discussed in the concluding Section 5.

Heaving foil

We consider a foil with chord L in a two-dimensional domain (x, y) immersed in an incoming incompressible flow of velocity U ∞ , heaving in the direction normal to the incoming flow direction, the heaving motion being

h = h 0 L sin(2πf t) (1)
with h 0 the (dimensionless) ratio between the heaving amplitude and the chord, f being the heaving frequency. The heaving velocity is written ḣ and the forces exerted on the foil are equivalent to those that would experience a fixed foil immersed in a flow with an angle of attack α and relative flow velocity U such that (see also [START_REF] Floryan | Large-amplitude oscillations of foils for efficient propulsion[END_REF])

α = arctan - ḣ U ∞ , U = U 2 ∞ + ḣ2 . (2)
The foil therefore is subject to a lift force f l as well as a form drag f d respectively in the direction e l and e d , forming en angle α with the streamwise x-direction, as sketched in figure 1. In the forthcoming analysis the flow quantities are made dimensionless with the foil's chord L and the incoming uniform flow velocity U ∞ . Also, the convention will be used that thrust is positive and drag negative.

In the following, we consider what is known as the reduced frequency

k = πf L U ∞ (3)
as dimensionless quantity and accordingly the dimensionless heaving law is

h * = h L = h 0 sin(2kt * ) and ḣ U ∞ = 2kh 0 cos(2kt * ). (4)
with t * = tU ∞ L the time made dimensionless. Accordingly, the angle of attack due to the heaving velocity can be written

α = arctan (-2kh 0 cos(2kt * )) . (5)
Here dimensionless frequencies and heaving amplitudes up to the order of one are to be considered which means that the time-varying angle of attack as experienced by the foil may be as large as ±60 ○ . Reliable lift formula are however only known for relatively small angles of attack in the limit of inviscid flow, the lift force for a slender symmetrical body being known to be approximately π sin(α)U 2 Lρ ( [START_REF] Batchelor | An Introduction to Fluid Dynamics[END_REF]) with U the incoming flow. Of course, there is no reason that such a formula quantitatively applies to the present heaving foil at large amplitudes. Considering nevertheless tentatively this formula (which would correspond to instantaneous lift f l for varying α in the e l direction as sketched in figure 1), its projection on the x-axis yields (with the convention that thrust is positive) π sin 2 (α)U 2 Lρ, U being given by (2). According to (5)

sin 2 (α) = 4k 2 h 2 0 cos 2 (2kt * ) 1 + 4k 2 h 2 0 cos 2 (2kt * ) while U 2 U 2 ∞ = 1 + 4k 2 h 2 0 cos 2 (2kt * ) (6)
and introducing a hypothetical proportionality factor C, a possible lift-induced dynamic thrust formula could be

f * thrust = C 4πk 2 h 2 0 cos 2 (2kt * ). (7)
While uncountable data are of course available for drag coefficients as function of angles of attack for NACA foils, it seems however that general analytic drag laws are not available, besides the evidence that pressure drag increases with the angle of attack. Making the ad-hoc assumption in the present case that the time-dependent drag is qualitatively proportional to sin(α) U 2 Lρ, then the (dimensionless) drag in the e d direction (see figure 1) would vary as

D 2k h 0 cos(2kt * ) 1 + 4k 2 h 2 0 cos 2 (2kt * ),
given [START_REF] Brousseau | Fluid-structure interaction effects on the deformable and pitching plate dynamics in a fluid flow[END_REF], D being a hypothetical factor of proportionality. Projection on the x-axis (multiplying by cos(α) with the convention that drag is negative) leads to the formula

f * drag = -D 2k h 0 cos(2kt * ) (8)
and adding this negative force to [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF] would yield a theoretical total force f * thrust + f * drag in the streamwise x direction. The flow motion is governed by the Navier-Stokes equations

∇ ⋅ u = 0 , ρ ∂u ∂t + ρ(u ⋅ ∇)u = -∇p + µ∇ 2 u , (9)
where u = (u, v) is the velocity field and p is the pressure, ρ and µ being the constant density and dynamic viscosity. The dimensionless (written with an asterisk) time-dependent pressure force (per unit span) in the fixed coordinate system (x, y) along the heaving foil is

f * = (f * x , f * y ) = - 1 ρU 2 ∞ L ∮ (-p)n ds (10)
where n = (n x , n y ) is the outward unit vector normal to the foil, the integration being performed along the foil's upper and lower face (the minus sign being applied such that a propulsive force is positive, which is the convention generally used). Also, the streamwise component f * x will be the key quantity in the forthcoming analysis. At the flow conditions considered, the viscous drag

f * µ,x = - 1 ρU 2 ∞ ∮ µ 2 ∂u ∂x n x + ∂u ∂y + ∂v ∂y n y ds (11)
(with the minus sign for the drag to be negative) will not contribute significantly to the overall force balance.

In the following, the question will be addressed to which extend the thrust formula [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF] and the drag prediction (8) are retrieved in the numerical pressure force data for a heaving foil. Prior to this we summarize in the nest paragraph the numerical environment used for the simulations and we briefly address the numerical convergence issue.

Numerical procedure and convergence

Direct numerical simulations of equations ( 9) are performed using the OpenFOAM toolbox (the v1906-version [40]); the solver is inherently threedimensional and only one void cell is used in the z direction for the present two-dimensional flow setting. The methodology used in OpenFOAM is based on a finite volume discretization of the differential operators together with a mesh-handling using structured or unstructured grids. The main features of the discretization, based on Gauss' theorem applied to the volume integrals, are briefly summarized. For the approximation of the resulting integrals over the control volume boundaries, the values of the flow quantities at the volume's faces have to be interpolated from the values at the computational nodes. Interpolation practices for finite volume methods are described for instance in [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF]. In the forthcoming simulations, for the convective term an upwind interpolation known as upwind differencing scheme is used, where the direction of the flow determines which node provides the cell face value. For the diffusive term a linear interpolation between the faces' two nearest nodes is performed, which is of second order. Linear interpolation is also used for the computation of the gradients normal to the surfaces. Time discretization is based on the first order implicit Euler scheme. In order to deal with the quadratic non-linearity in the Navier-Stokes system, the equations to be solved at each time step are linearized about the result at the previous time step, in the way that the neglected terms are of second order in time, that is smaller than the error of the time discretization. The time step is adapted during the time integration via the Courant number, well known in computational fluid dynamics. All linear systems are solved using the built-in pimpleFoam predictor-corrector solver, that merges the PISO and SIMPLE algorithms widely used in computational fluid dynamics [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF].

In the present analysis we address large amplitude motions of foils and the so-called overset grid approach (also called Chimera approach) appears to be particularly suitable in this context. This approach consists of meshed subdomains into which the objects are embedded, these meshes being free to move independently one to each other. These subdomains are superim- posed to a fixed meshed background domain. For each independent mesh the discretized Navier-Stokes system is advanced in time and at each time step the coupling with the overall solution domain and the background mesh is achieved by transmitting, using interpolation techniques, flow conditions, that is fluxes in the finite volume context used in OpenFOAM. General considerations concerning the coupling of the subdomains can be found in [START_REF] Houzeaux | A chimera method for the incompressible navier-stokes equations[END_REF] and the use of the Chimera approach in the OpenFOAM-v1906 environment is explained in [START_REF] Tisovska | Description of the overset mesh approach in esi version of openfoam[END_REF] and [START_REF] Chandar | On overset interpolation strategies and conservation on unstructured grids in openfoam[END_REF], for instance. The trajectory of the moving mesh (or independently moving meshes when several moving objects are considered) is prescribed via a text file which is generated prior to the numerical integration of the system. The coordinates in time, with respect to its initial position, of a specific point of the object is to be provided together with the possibility of rotating the object around the axis normal to the computational domain, allowing the combination of translational and pitching motions. The discrete time subdivision of the prescribed motion is to be chosen as small as (or smaller than) the time-step, the value of which being automatically adapted via the Courant number of the Navier-Stokes time integration. In the forthcoming simulations this latter time step is generally of the order of 10 -3 and the the prescribed motion has been discretized with time-intervals as small as 10 -4 . The moving object is here a NACA0018 foil (which is considered in all the simulations) and figure 2 is a zoom showing the overset mesh superimposed to the background mesh (here the position is associated with combined heavingpitching motion which will be addressed later in this paper). The transmission from the background mesh to the moving mesh is achieved by interpolation (the cellVolumeWeight option [START_REF] Tisovska | Description of the overset mesh approach in esi version of openfoam[END_REF] is used) on the outer border cells of the overset mesh, whereas in return the transmission from the moving mesh back to the background mesh uses the background mesh cells immediately adjacent to the foil's boundary.

For the convergence analysis the Reynolds number Re = ρU ∞ L µ has been set to 5 10 4 . As heaving parameters, the reduced frequency k = 0.5 and the heaving amplitude h 0 = 0.8 have been chosen, these values being representative of the forthcoming heaving parameter exploration. Uniform free-stream u = U ∞ e x is prescribed at inflow, whereas at the upper and lower boundary zero normal-gradient conditions are imposed for the two components of the flow velocity. A zero pressure is imposed on the boundaries, except at inflow where zero normal-gradient pressure applies. When starting a new computation, uniform free-stream velocity is imposed in the domain as initial condition. The table 1 briefly summarizes the OpenFOAM options and boundary conditions used.

In figure 3 a typical snapshot of the streamwise velocity in the vicinity of the foil is shown, together with the pressure distribution along the upper and lower faces. The foil is in upward motion at a position where the heaving velocity is approximately the maximum, that is the associated angle of attack (see figure 1) is at its maximum, which explains the large separation region

M 1 M 2 M 3 M + 1 M ++ 1 L x ∆ x 0.01L 0.0083L 0.0067L 0.01L 0.01L ∆ y 0.01L 0.0083L 0.0067L 0.01L 0.01L ∆ x,2 0.016L 0.0125L 0.01L 0.016L 0.018L ∆ y,2 0.016L 0.0125L 0.01L 0.016L 0.018L max(f * x ) 0.436 0.408 (-6.4%) 0.427 (4.7%) 0.458 (5%) 0.383 (-12%)
Table 2: Mesh parameters for the computational domain sketched in figure 2. Note that ∆ x,2 , ∆ y,2 are the values in the coarse mesh region of the moving mesh. The last line gives the maximum values of the periodic pressure force f * x (taken at t * = 3π) for each mesh, the percentage values in parentheses being those of the differences between the value of the finer to the coarser mesh (M 2 to M 1 and M 3 to M 2 ) and the difference of the values for the larger domains M + 1 and M ++ 1 with respect to M 1 .

at the lower side. Table 2 lists the parameters for different meshes considered. Here L x and L y provide the dimension of the global computational domain with the background mesh. A specific point on the foil at one quarter of the chord from the leading is chosen, which in the forthcoming pitching-heaving motion will be the pitch-pivot point, set at (x 0 , y 0 ) = (0, 0) in the computational domain for the foil at its initial position. For instance L x = [-3L, 7L] means that the overall length is ten times the foil's chord L, with 3L the distance from the inflow boundary to x 0 and 7L the distance from x 0 to the outflow boundary. The vertical dimension L y = [-3L, 3L] means a domain's height of 6L with a distance 3L from the foil's chord line to the bottom as well as to the top. The dimensions of the moving overset-mesh are given by L x,2 and L y,2 respectively (see figure 2). Furthermore ∆ x , ∆ y (respectively ∆ x,2 , ∆ y,2 ) are the grid sizes in both directions, for the background and moving mesh, respectively. While for each domain the background mesh has a constant grid size, on the moving mesh, when approaching the foil, two successive refinements are performed in a transient zone, dividing the grid size by 4. At last, in the very vicinity of the foil and inside the boundary layer three additional refinements are performed, dividing again the grid size by 8.

The meshes labeled M 1 , M 2 , M 3 correspond to the same domain size, but to successive grid-size refinements by an almost constant factor (as advocated for instance in [START_REF] Wilson | Comprehensive approach to verification and validation of cfd simulations-part2: application for rans simulation of a cargo/container ship[END_REF]), the grid size being divided here two times by ≈ 1.25.

The meshes labeled M 1 , M + 1 and M ++ 1 in table 2 have the same density of dis-cretization (besides a somewhat coarser grid in the moving domain for mesh M ++ 1 ), the domain sizes being successively increased by a factor of ≈ √ 2. Considering for instance the mesh labeled M 2 , the grid size of the background mesh is (∆x, ∆y) = (0.0083L, 0.0083L). In the coarse mesh region of the moving mesh the values are (∆ x,2 , ∆ y,2 ) = (0.0125L, 0.0125L) which is refined to (∆ x,2 32, ∆ y,2 32), yielding a grid size as small as 0.0004L in the boundary layer region. Concerning the numerical resolution in the vicinity of the foil, the Blasius formula 5L r L 1 Re for the boundary-layer thickness (r L being the dimensionless distance from the leading edge) could provide an idea about the number of grid points inside the boundary-layer region. This formula yields 0.022L r L and given the grid size 0.0004L close to the foil, one may estimate that roughly 55 r L points are inside the boundary layer (that is still 10 points at a dimensionless distance r L as small as 0.033 from the leading edge). It should also be mentioned that the background mesh has to be chosen sufficiently fine to ensure a reliable transmission of the flow conditions from the moving mesh. Indeed, the discrete contour associated with the background mesh cells the closest to the foil must fit in the best possible way with the foil's contour. The pressure forces generated by the foil's motion is at the heart of the investigation and we focus on this quantity for numerical validation. The dynamic force [START_REF] Piñeirua | Resistive thrust production can be as crucial as added mass mechanisms for inertial undulatory swimmers[END_REF] computed for different meshes are superimposed in figure 4(a) and the curves are seen to be very close one to each other, besides slight differences at the peak values. The viscous drag [START_REF] Jones | Numerical computation of flapping-wing propulsion and power extraction[END_REF] has been computed as well and is depicted in figure 4(b) (for two meshes only, the curves being close) and it is interesting to note that (after a short transient) at t * = π, 2π... the viscous force peak is positive, given that the heaving velocity ḣ * = 2kh 0 cos(t * ) (for k = 0.5) undergoes a minimum or maximum, leading to large separation of the flow field alternately on the upper face and lower face. Note that the drag force is almost two orders of magnitude smaller than the pressure force, what justifies to only consider the pressure force for thrust calculations, at least for this range of Reynolds numbers. After a short transient regime the dynamic pressure becomes periodic with peaks at t * = π, 2π... (see figure 4(a)) and for the convergence analysis, table 2 provides the maximum value of f *

x (taken at t * = 3π) for each mesh. The differences, expressed in percentage, between these values when performing successive grid refinement, that is for M 2 compared to M 1 and M 3 compared to M 2 , are provided as well. The differences slightly decrease and change sign, which indicates oscillatory convergence (see for instance [START_REF] Wilson | Comprehensive approach to verification and validation of cfd simulations-part2: application for rans simulation of a cargo/container ship[END_REF] for general considerations on validation and verification of numerical simulation results). Concerning the influence of the computational domain size, the val- ues for M + 1 and to M ++ 1 are compared to the value for M 1 , the pressure peak being higher for M + 1 and lower for M ++ 1 . It is not surprising that finite-size effects are observed, the differences remaining however relatively small. It has to be emphasized that a definite mesh convergence evidence is hardly to be achieved, for this unsteady flow field. Indeed, as mentioned earlier for instance the time step is automatically adapted via the Courant number and it hence changes when refining the grid. Also, when enlarging the domains' dimensions, the relative size between the moving domain and the background domain has been kept constant, which modifies however the distance from the foil to the domains' interpolation region. Note that the background mesh for the M 2 discretization contains a bit more than 8.7 10 5 points and the moving mesh about 2 10 4 points (due to the refining close to the foil) which yields a total number of almost 9 10 5 points. This M 2 mesh has been chosen in all the forthcoming simulations, as a compromise between precision and computational cost.

Thrust scaling

The formula [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF] for lift-induced thrust has been derived according to classical fluid mechanics, the underlying hypotheses being slender body and potential flow theory as well as small angles of attack. When considering for instance cycloidal propulsion, the foils' rotation can be decomposed into heaving and surging, the heaving and surging amplitude being likely of the same order than the foils' chord length (in such devices the foils undergo in general also additional pitching motions). For instance, the Voith Schneider propulsion addressed in [START_REF] Esmailian | Numerical investigation of the performance of voith schneider propulsion[END_REF] would correspond to h 0 = 0.75 whereas the reduced frequency k is of order one and the angle of attack (2) experienced by the heaving foil takes values during a heaving period as large as 50 ○ . The question is whether a formula such as ( 7) is meaningful for the heaving NACA0018 foil considered here. Numerical data have been generated for frequencies k = 0.25, 0.5, 0.75, 1, the heaving amplitudes being h 0 = 0.4, 0.6, 0.8, 1, that is 16 (k, h 0 ) parameter combinations have been considered, the Reynolds number being Re = 5 10 4 . Of course, the Reynolds number associated with real propulsion in water would at least be one or two orders of magnitude higher. The value considered is however likely to be large enough for the pressure stress computations to provide significant thrust results, viscous effects being negligible for the force outcome at this Reynolds number. The oscillatory pressure force f *

x can be decomposed as

f * x = f * l,x + f * d,x (12)
with f * l,x the lift force and f * d,x the drag force, both projected on the x-axis. For the comparison with the formula ( 7) and ( 8), the factors of proportionality C and D have to be estimated. The choice has been made that the time-averaged values ⟨f * l,x ⟩ and ⟨f * d,x ⟩ (⟨⋅⟩ being the average over one heaving period) coincide with those of the formula and given that ⟨cos 2 (2kt * )⟩ = 1 2 and ⟨ cos(2kt * ) ⟩ = 2 π, C and D have been chosen such that

⟨f * l,x ⟩ = C 2πk 2 h 2 0 , ⟨f * d,x ⟩ = -D 4 π kh 0 , (13)
C and D having been computed for all the heaving parameter combinations. The figures 5 and 6 show the computed thrust and drag contributions to the pressure force over several heaving periods, for the heaving parameter k = 0.5, h 0 = 0.8 as well as k = 0.75, h 0 = 0.6, respectively. The theoretical formula [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF] and ( 8 computed dynamic thrust and drag, but it has to be emphasized that different C, D values are retrieved. Considering the 16 parameter combinations, it appears that C varies between 0.3 and 0.6, whereas the D coefficients are in an even larger range from 0.3 to 1.2. Assuming a linear dependence of these coefficients with respect to k and h 0 , writing range. The theoretical model for the mean pressure force is accordingly

C(k, h 0 ) = c 1 + c 2 k + c 3 h 0 , D(k, h 0 ) = d 1 + d 2 k + d 3 h 0 , (14)
F * theory = 2πC(k, h 0 )k 2 h 2 0 - 4 π D(k, h 0 )kh 0 (16)
(with C(k, h 0 ) and D(k, h 0 ) given by 14). The time-average F * x of the oscillatory pressure force [START_REF] Moored | Inviscid scaling laws of a self-propelled pitching airfoil[END_REF] has been computed and the values are shown in figure [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF] against its model counterpart; it is seen that the theoretical thrust prediction is close to the computed values. Note that classical thrust coefficient formula have been derived in [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF] for small heaving amplitudes, based on potential flow formulas [START_REF] Theodorson | Genneral theory of aerodynamic instability and the mechanism of flutter[END_REF], for an idealized airfoil represented as a straight line, that is without pressure form drag. This theory yields a formula 2πk 2 h 2 0 (F 2 + G 2 ), known as Garrick's theory where F and G are function of the reduced frequency tabulated in [START_REF] Garrick | Propulsion of a flapping and oscillating airfoil[END_REF]. The corresponding thrust values have for instance been compared in [START_REF] Jones | Numerical computation of flapping-wing propulsion and power extraction[END_REF] with numerical values (obtained with a two-dimensional panel code) for plunging NACA foils, however for h 0 values smaller than 0.4 which is the lowest value considered here. It is interesting to note that the thrust part in the model ( 16) is at least roughly of the order of magnitude with what one gets with Garrick's formula: it has been checked that for the (k, h 0 ) combinations considered, the latter formula would however provide values about 20% to 40% higher. Large amplitude oscillations have also been experimentally addressed in [START_REF] Floryan | Large-amplitude oscillations of foils for efficient propulsion[END_REF], considering the Strouhal number formed with the maximum foil excursion 2Lh 0 , that is St = 2f Lh 0 U ∞ , as dimensionless frequency. It is easily seen that St = 2kh 0 π and in [START_REF] Floryan | Large-amplitude oscillations of foils for efficient propulsion[END_REF] some evidence is provided for a thrust dependence on St 2 (that is ∼ k 2 h 2 0 ). The present investigation however suggests that such a simplified law only provides a rough estimate of the thrust performance.

Heaving-pitching foil

In order to put a step forward towards the modeling of propulsion devices, a pitching motion is superimposed to the heaving motion considered so far. The pitching motion considered has the same period as the heaving motion, introducing a phase-shift φ such that the pitching angle is

θ(t * ) = α 0 sin(2kt * + φ), (17)
(the heaving motion being h * (t) = h 0 sin(2kt * )). The pitch-pivot point is set at one quarter of the foil's chord from the leading edge, that is the chordline from the pivot point to the trailing edge forms an angle θ(t * ) with the horizontal inflow direction. The phase-shift φ = -π 2 has been chosen and accordingly θ(t * ) = -α 0 cos(2kt * )), which corresponds to the foil's motion in an experimental device to be discussed in the next section.

Simulations have been performed for three values of the maximum pitching angle, that is α 0 = 10 ○ , 20 ○ , 30 ○ , again at Re = 5 10 4 . The left picture in figure 8 shows a sketch of the foils position during the upward and downward heaving motion. An instantaneous flow vorticity field for the pitching angle α 0 = 20 ○ , the heaving amplitude h 0 = 0.8 and the reduced frequency k = 0.5 is shown as well, the foil being in its upward heaving motion at the instant such that θ(t) ≈ -0.5 α 0 and h * (t) ≈ 0.866 h 0 . For the same heaving amplitude h 0 = 0.8 and frequency k = 0.5, the oscillatory pressure force f * x is shown for the three pitching angles in figure 9(a). As can be seen, the maximum peak is for α 0 = 10 ○ , the result for α 0 = 20 ○ is close to the values for pure heaving, whereas the thrust performance is the lowest for the highest pitching angle α 0 = 30 ○ . For these parameter values the viscous force f * µ,x in the x-direction associated with the viscous stress has been computed as well and it is shown in figure 9(b). It is interesting to note that while for α = 0, 10 ○ the peak at t * = π, 2π... is positive, for higher pitching angles however it is negative. This can be explained by the fact, that when the foil is pitching with a phase shift -π 2 with respect to the heaving motion, the apparent incoming flow direction sketched in figure 1 forms, when pitching is added, an angle α(t)-θ(t) with the chord line. This means that the pitching lowers the angle of attack experienced by the foil and hence the intensity of separation. In all cases the mean viscous force is of course negative, contributing as drag. Note gain that the magnitude of viscous drag is negligible with regard to the thrust. In the previous section it has been attempted to scale the dynamic force in terms of the reduced frequency k and the heaving amplitude h 0 and the model [START_REF] Floryan | Scaling the propulsive performance of heaving and pitching foils[END_REF] proved to fit to some degree the time-dependent thrust force as well as its time-averaged values. Adding pitching to the heaving motion, it may be conjectured that the induced force is a combination of a contribution due to heaving as well as one associated with pitching, together with a coupling component. Pure pitching motions for a plate with vanishing thickness has been considered in [START_REF] Labasse | Numerical exploration of the pitching plate parameter space with application to thrust scaling[END_REF] for relatively small pitching angles and it could be shown that in terms of reduced frequency and maximum pitching angle α 0 the timeaveraged thrust scales scales as

a α 2 0 k 2 + b α 2 0 k, (18)
a and b being coefficients depending on the pitch-pivot point location. Note that a similar scaling has been proposed for instance in [START_REF] Floryan | Scaling the propulsive performance of heaving and pitching foils[END_REF]. When adding the pitching to the heaving motion, the outward unit vector normal to the foil rotates with respect to the vertical axis e y by the pitching angle θ. The resulting acceleration of the fluid in the x-direction may hence be estimated to be ḧ sin(θ). Note that L U 2 ∞ ḧ = -4k 2 h 0 sin(2kt * ), given the expression (1) of the heaving motion, and the pressure force in dimensionless form associated with this acceleration is proportional to

k 2 h 0 sin(2kt * ) sin(α 0 sin(2kt * + φ)), (19)
given the expression [START_REF] Wu | A review on fluid dynamics of flapping foils[END_REF]. Even though the pitching's maximum angle α 0 is not particularly small, one is nevertheless tempted to consider the expression above at the leading order in α 0 which would yield by taking the time-average an expression proportional to k 2 h 0 α 0 . While the proportionality factor will probably depend on the the phase shift φ of the pitching motion with respect to heaving, there is however no reason that it may be inferred from the time average of the product sin(2kt * ) sin(2kt * + φ) (which would be zero for instance for φ = -π 2 considered here). Indeed, the fluid-foil dynamics by itself is likely to introduce phase shifts between the periodic foil motion and the induced flow acceleration.

In the previous section, the time-dependent pressure induced by the pure heaving motion could be reasonably fitted for the set of parameter values considered and the law ( 16) could be derived for the time-averaged force as the sum of lift-induced thrust and drag, with however non constant coefficients given by ( 14). At leading order in the parameters, the lift-induced thrust is proportional to k 2 h 2 0 , whereas the drag is proportional to kh 0 . Retaining these leading order terms, adding terms proportional to α 2 0 k 2 and α 2 0 k (associated with the pitching law ( 18)) as well as a term proportional to k 2 h 0 α 0 which stands for the force component due to the coupling between pitching and the heaving acceleration, the following time averaged pressure force expression

F * theory = B 1 α 2 0 k 2 + B 2 α 2 0 k + B 3 k 2 h 2 0 + B 4 kh 0 + B 5 α 0 k 2 h 0 ( 20
)
is proposed. The coefficients B j , j = 1, ⋯, 5 have been computed by a least-square fit for the 36 parameter values considered, that is all the (α 0 , k, h 0 ) combinations (with α 0 = 10 ○ , 20 ○ , 30 ○ , k = 0.25, 0.5, 0.75, 1 and h 0 = 0.6, 0.8, 1). The values found are

B 1 = -0.9165, B 2 = -1.5409, B 3 = 0.6431, B 4 = 0.0496, B 5 = 2.4108 (21)
and the computed force F * x against its model counterpart ( 20) is shown in figure 10. It is seen that the model fits well for the range of parameter values considered, where small heaving amplitudes or large frequencies are discarded. As already mentioned, the parameter values are however in the range of an experimental cycloidal propulsion device discussed in the next section. 20), with coefficients [START_REF] Hu | Two-dimensional simulation of the hydrodynamic performance of a cycloidal propeller[END_REF], the symbols representing the values for the 36 parameter values (α 0 , k, h 0 ), with α 0 = 10 , 20 ○ , 30 ○ , k = 0.25, 0.5, 0.75, 1 and h 0 = 0.6, 0.8, 1; the solid line is the bisector of the first quadrant angle.

B 1 α 2 0 k 2 + B 2 α 2 0 k + B 3 k 2 h 2 0 + B 4 kh 0 + B 5 α 0 k 2 h 0 -0.

Application to cycloidal propulsion

The most widely known cycloidal propulsion device is the Voith-Schneider propeller the performance of which has for instance been numerically addressed in [START_REF] Esmailian | Numerical investigation of the performance of voith schneider propulsion[END_REF]. These type of propellers consist of several blades rotating around a vertical axis, the individual plates themselves undergoing pitching motions around their own pivot points. These propulsion devices are capa- ble to generate thrust forces in all directions. They are very efficient when maneuverability is required and are hence frequently used for the propulsion of tugs and ferries, for instance. A cross-flow propeller has recently been subject of experimental investigations at the Naval Academy Research Institute (IRENav) in France, where such a device has been designed, the pitching motions being controlled electrically by servo-motors enslaved by the main motor ( [START_REF] Fasse | An experimental blade-controlled platform for the design of smart cross-flow propeller[END_REF]), numerical simulations having been performed for this configuration too [START_REF] Fasse | Numerical study of a sinusoidal transverse propeller[END_REF]. The device with three foils is sketched in figure 11, where the solid foils' contours are at their positions at t 0 when foil labeled 1 has the pitching amplitude α 0 . The circle drawn corresponds to the trajectory of the foil and the hatched foils' contours stand for their respective positions at t 0 + T 6. The angle between two foils is 120 ○ . The pitch-pivot point is at L 4 from the leading edge, the incoming flow arriving from the left. The configuration, for which measurements have been performed using Ifremer facilities at Boulogne-sur-Mer, is that of a NACA0018 foil with chord length L = 0.35 m, whereas the radius of the circle of motion is R = 0.4 m. The period of circular rotation and pitching is the same for each foil and the motion can be decomposed into heaving, pitching and a surge motion parallel to the x-direction of incoming flow. The motions are synchronizes such that the dimensionless heaving motion h * and surge motion s * are (h * being the distance from the axis's center in the vertical direction y and s * the distance in the horizontal direction x). The pitching angle around the individual pivot points is

h * = h 0 sin(2kt * ), s * = h 0 cos(2kt * ) (22)
θ(t) = -α 0 cos(2kt * ). (23)
that is the phase shift with the heaving motion is -π 2 as in the previous section. The position of foil labeled 1 in figure 11 would correspond to t * = nπ k (with n an integer).

In the experimental investigation the incoming flow velocity is of the order of 1 m s in the Ifremer water tank which yields a Reynolds number (based on the foil's chord) of the order of 3. these parameters are within the range of those considered in the previous section for the heaving-pitching foil. In section 2, for the pure heaving motion, the influence of the mesh as well as the size of the computational domain has been assessed for this latter Reynolds number and the mesh labeled M 2 (for a computational domain 10L in the x-direction and 6L in the y-direction) proved appropriate to get reliable force results.

Simulations have been performed for these mesh parameters, at Re = 5 10 4 and considering the complete configuration with the three rotating foils, with h 0 = 1.143, α 0 = 22 ○ and k = 0.6. A snapshot of the flow field is shown in figure 12. The question naturally arises to what extent the force generated by one foil rotating with the two others differs from what one would get for a single foil in an equivalent rotating motion. This comparison has been performed for the same parameter values h 0 = 1.143, α 0 = 22 ○ k = 0.6 as those corresponding to the snapshot in figure 12. The results are plotted as function of angle β, β = 0, 360 ○ ... being the maximal vertical position (see x are shown as the symbols for a pitchingrotating foil against its model counterpart (20) for pitching-heaving foil with coefficients B j , j = 1, ⋯, 5 given by ( 21). The solid line being the bisector of the first quadrant angle. figure 11) of the foil around which the pressure force is computed. It is seen that when the foil's trajectory is on the left semicircle (0 < β < 180 ○ ) the pressure force components in the x direction almost superimpose, whether a single rotating foil is considered or a foil in the device. As expected, when the foil's trajectory crosses the wake structure of the device (at the right semicircle of the trajectory) the peaks of the pressure force are somewhat shifted between the two configuration. On the whole however, and at least for a rotation radius h 0 close to one, the results are rather close and in particular when time-averages are considered. This supports the idea that to some degree the thrust for one single rotating flow may be extrapolated to recover the global thrust for a configuration like that in figure 11 (by superimposing three times the oscillatory pressure force for one foil with shifts of 120 ○ in angle β). The rotating-pitching foil differs from the heaving-pitching foil by the additional horizontal back and forth motion and the question is to what degree this latter surge motion affects the thrust results. The oscillatory pressure force for a single heaving-pitching foil and a rotating-pitching foil has been computed, for different parameter combinations. The result for k = 0.5, h 0 = 1, α 0 = 10 ○ is shown in figure 14, the two force components being depicted as function of β (it is recalled that β = 0 ○ , 360 ○ ... are the highest vertical position for each period). The curves are close, the force in the x direction exhibiting however a somewhat higher peak when the rotating foil moves toward the incoming flow, and a lower peak when moving in the opposite direction. One may wonder whether the heaving-pitching force model [START_REF] Scherer | Experimental and theoretical investigation of large amplitude oscillating foil propulsion systems, HYDRONAUTICS, Incorporated[END_REF] may still be relevant when rotation instead of heaving is considered. The time-averaged thrust force for the rotating and pitching foil has been computed for 6 parameter combinations and the values are compared with this thrust law in figure 15. The computed values are shown as the symbols and it is seen that the model prediction is still rather reliable, which means that in average the additional surging motion does not contribute significantly to the global force. For the experimental device at IRENav, measurements have not be performed at a Reynolds number of 5 10 4 , but for instance at Re = 2.8 10 5 , that is for an incoming flow velocity of 0.8 m s in a water tank, the blades' chord being 0.35 m. Of course, to perform the numerical simulation at such a high Reynolds number is questionable, even in the two-dimensional case, at least if one aims at capturing complex flow structures. Concerning the pressure force at the heart of the present investigation, it is however conjectured that it is not strongly affected by the Reynolds number (the value 5 10 4 considered so far being already in a relatively high Reynolds number range). To give some evidence for this assumption, for the same computational mesh M 2 as used so far, the oscillatory pressure force of one rotating-pitching foil has been computed for Re = 2.8 10 5 and the comparison with the result at Re = 5 10 4 is shown in figure 16, for k = 0.6 and h 0 = 1.143 (this latter amplitude is that of the experimental device). It is seen that the thrust values f *

B 1 α 2 0 k 2 + B 2 α 2 0 k + B 3 k 2 h 2 0 + B 4 kh 0 + B 5 α 0 k 2 h 0 -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 F * x -0.2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1 k = 0.18, h 0 = 1.143, α 0 = 30 • 2 k = 0.3, h 0 = 1, α 0 = 10 • 3 k = 0.7, h 0 = 0.6, α 0 = 20 • 4 k = 0.5, h 0 = 1, α 0 = 10 • 5 k = 1, h 0 = 0.6, α 0 = 30 • 6 k = 0.6, h 0 = 1.143, α 0 = 30 • 7 k = 0.6, h 0 = 1.143, α 0 = 15 • 8 k = 0.6, h 0 = 1.143, α 0 = 22 • 9 k = 0.9, h 0 = 0.85, α 0 = 30 • 10 k = 1, h 0 = 0.8, α 0 = 30 • 11 k = 0.95, h 0 = 1, α 0 = 10 • 12 k = 1, h 0 = 1, α 0 = 20 •
x are rather close, besides a slightly stiffer peak for the high Reynolds number case at the most downstream position (with respect to the incoming flow) of the foil on its circular trajectory (at β = 270 ○ , 630 ○ ...).

In the experimental thrust performance investigation at the Boulogne-sur-Mer Ifremer facilities, the time-dependent hydrodynamic loads are measured with a blade embedded load cell. Note that the blades have a rather large span-to-chord aspect ratio of 2.8 and three-dimensional effects are possibly limited regarding the overall thrust generated by each foil. A detailed description of the experimental propeller design is given in [START_REF] Fasse | An experimental blade-controlled platform for the design of smart cross-flow propeller[END_REF]. The hydrodynamic force is measured on one blade and the total force is reconstructed by shifting the loads by respectively 120 ○ and 240 ○ , which means that the reconstructed total force is T 3 periodic with T the period of rotation. The flow structure generated by the rotating plates might however not exhibit exactly this T 3-periodicity. In the numerical simulation the instantaneous force values for each blade have been computed and added to form the total force. Measurements for a pitching angle α 0 = 30 ○ have been performed and figure 17 shows the comparison with the simulation results, for both the force components, at k = 0.29. It is seen that the results coincide fairly well, but the T 3 periodicity, or equivalently the 120 ○ periodicity, is not exactly observed in the numerical simulation results. Indeed, the individual force values depicted in figure 18 for each blade do not exactly coincide at identical positions (β = 0 ○ , 360 ○ ⋯ being the highest vertical position, four rotation periods being shown), in particular in the downstream half of the circular trajectory. This is to be attributed to slightly different flow structures encountered by the individual plates at identical positions. The comparison between the measurements and the simulations at the higher frequency k = 0.49 is shown in figure 19. Now, the simulation results are not far from exhibiting a T 3 period: for this frequency the curves differ a little more, but the general trends are however retrieved in the simulations.

Regarding the time-averaged values, for the frequency k = 0.29 one finds F * x = -0.0094 in the experiment whereas the simulation result is F * x = -0.0118. At the higher frequency k = 0.49, one gets F * x = 0.6008 and F * x = 0.6359 for the measurements and the simulations, respectively, thus the measured dynamic forces as well as the time averaged quantities are approximately retrieved in the numerical simulations.

Concluding discussion

The present investigation has shown that the overset-option and the general modelization capabilities of OpenFOAM are suitable for the simulation of a rigid structure undergoing large scale motions in an incoming flow field. The advantage of such a body-fitted mesh approach is to provide accurately the dynamic forces generated in this system. A theoretical model, which proved reliable for the time-dependent pressure stress as well as for the global time averaged thrust, could be derived for a heaving foil in the form

F * theory = 2πC(k, h 0 )k 2 h 2 0 - 4 π D(k, h 0 )kh 0 .
The first term of the right-hand side is associated with the lift-induced thrust, the second term standing for form-drag due to a virtual angle of attack during the heaving motion. The expressions C(k, h 0 ) and D(k, h 0 ) could not reliably be taken as constants, but rather as linear functions of the heaving amplitude h 0 and the reduced frequency k. A scaling for the thrust merely proportional to k 2 h 2 0 , as proposed for instance in [START_REF] Floryan | Large-amplitude oscillations of foils for efficient propulsion[END_REF], should hence be used with caution. The present investigation has been guided by the objective to numerically simulate fluid-structure interactions for motions associated with cycloidal propulsion. For pitching-heaving motions and focusing on time-averaged thrust prediction, a theoretical scaling (20) is proposed. Besides the k 2 h 2 0 and kh 0 terms of heaving, the expressions α 2 0 k 2 and α 2 0 k associated to pitching (see [START_REF] Labasse | Numerical exploration of the pitching plate parameter space with application to thrust scaling[END_REF][START_REF] Floryan | Scaling the propulsive performance of heaving and pitching foils[END_REF]) enter, as well as the coupling term α 0 k 2 h 0 . This law with 5 scaling parameters has been proven to fit with the numerical simulation data, as shown in figure 10. It appears that this law is still reliable when adding a surging motion, or equivalently, when addressing pitching-rotating foils typical for cross-flow propellers, at least for the simulation data at a Reynolds number of 5 10 4 . The experimental device designed at the IRENav operates however at higher Reynolds numbers, but still below those of industrial devices. Also, this cycloidal propulsion device has a plate's span-to-chord aspect ratio of 2.8 and the resulting quasi-two-dimensional measurements have been shown to be in good agreement with the two-dimensional simulations, for the same Reynolds number Re = 2.8 10 5 as in the experiment. The values [START_REF] Hu | Two-dimensional simulation of the hydrodynamic performance of a cycloidal propeller[END_REF] for the proportional parameters B 1 , ⋯, B 5 of the thrust law [START_REF] Scherer | Experimental and theoretical investigation of large amplitude oscillating foil propulsion systems, HYDRONAUTICS, Incorporated[END_REF] have been computed by least-square fit with numerical simulation data at Re = 5 10 4 and for a pitching-heaving foil. Addressing in this concluding discussion briefly the reliability of the law for pitching-rotating motions when changing the Reynolds number, figure 20 for h 0 = 1.143 corresponding to the experimental device and for α 0 = 30 ○ (the three frequencies chosen being in the parameter range of the experimen-tal measurements). The computed force values are relatively close to the law and they only slightly differ with the Reynolds number, at least for this motion-parameter range. One may conjecture that a law such as (20) might be helpful, in the sense that typical pitching-rotation parameters may be inferred for particular thrust regimes, such as zero thrust for a given flow velocity or maximum thrust, for instance. It has to be emphasized that of course the proportional parameters B 1 , ⋯, B 5 in (20) will depend on the specific foils considered (here a NACA0018 profile, which was used in the experimental device) and more importantly on the phase-shift between the pitching and rotating motions. Of course, the propulsive efficiency also plays an important role for the performance of moving foils. This efficiency is in general measured by the ratio of thrust power output to the power input to the foil. The time-averaged power input for a heaving-pitching foil is P i = ⟨f y ḣ⟩ + ⟨m z θ⟩ (e.g. [START_REF] Mackowski | Direct measurement of thrust and efficiency of an airfoil undergoing pure pitching[END_REF], [START_REF] Floryan | Scaling the propulsive performance of heaving and pitching foils[END_REF]), with m z the moment about the pitch-pivot point and f y the pressure force in the direction of heaving motion. The efficiency can be written as the ratio of dimensionless quantities and one gets the expression for the present foil's motion with a phase shift -π 2 between the pitching and heaving

Figure 21 shows the efficiency η as function of thrust F * x , the different solid curves being line-fits at constant k through the values for different pitching angles (the heaving amplitude is set to h 0 = 1). It is seen that for a given frequency k the efficiency is the lowest for pure heaving (α 0 = 0) and it increases with the pitching angle. Efficiency of relatively large amplitude heaving and pitching foils has been addressed for instance in [START_REF] Floryan | Large-amplitude oscillations of foils for efficient propulsion[END_REF], these authors pointing out discrepancies in the literature concerning efficiency measurements. The present computations provide efficiency values in the range between 0.1 and 0.7. The issue of the present investigation was thrust-prediction and a detailed efficiency analysis is not attempted. One may however remark that the efficiency globally decreases with the frequency (for a fixed heaving amplitude) which is in general line with the measurements reported in [START_REF] Floryan | Large-amplitude oscillations of foils for efficient propulsion[END_REF]. In this latter investigation the phase shift between pitching and heaving is φ = 3π 2 which means that in this case the pitching increases the angle of attack associated with the heaving motion. In our case however (that is for a phase shift φ = -π 2) the angle of attack decreases with pitching (see section 3). Therefore, the fact that the efficiency in the present case increases with the angle of pitching corresponds to the known fact that large angles of attack are detrimental to efficiency (see [START_REF] Floryan | Large-amplitude oscillations of foils for efficient propulsion[END_REF]). Finally, it can be inferred from figure 21 that while thrust globally increases with the oscillation frequency, there is however for each k-value an optimal (with regard to thrust) pitching angle.

Figure 1 :

 1 Figure 1: Sketch of the heaving foil immersed in an incoming flow: lower figure shows U as seen in the foil's coordinate system as well as the sketch of the lift f l and the form drag f d .

Figure 2 :

 2 Figure 2: Zoom on the mesh structure, showing the moving overset mesh with dimensions (L x,2 , L y,2 ) superimposed to the background mesh with dimensions (L x , L y ).

Figure 3 :

 3 Figure 3: Snapshot of the streamwise velocity component (a) and the corresponding pressure distribution along the upper face and lower face of the foil as function r L from the leading edge (b). The foil is in the upward motion at approximately the highest heaving velocity. The heaving parameters are k = 0.5 and h 0 = 0.8, the Reynolds number being Re = 5 10 4 .

Figure 4 :

 4 Figure 4: Streamwise component f * x of the pressure force (a) and viscous drag (b) as function of time t * and for different meshes (for the mesh-labels M 1 , M 2 , M 3 , M ++ 1 see table2), for a heaving amplitude h 0 = 0.8 and frequency k = 0.5, at Re = 5 10 4 .

Figure 5 :

 5 Figure 5: Dimensionless lift-induced thrust f * l,x (a) and pressure drag f * d,x (b) as function of time t * , for a heaving amplitude h 0 = 0.8 and frequency k = 0.5, at Re = 5 10 4 . The theoretical formulas f * thrust = C 4πk 2 h 2 0 cos 2 (2kt * ) and f * drag = -D 2kh 0 cos(2kt * ) are superimposed as the broken lines, for C = 0.477, D = 0.728.

the coefficients (c 1

 1 , c 2 , c 3 ) = (0.745, -0.277, -0.150), (d 1 , d 2 , d 3 ) = (-0.387, 0.640, 0.923)[START_REF] Anderson | Oscillating foils of high propulsive efficiency[END_REF] have been obtained performing a least-square fit with the 16 values of C and D as given by the computations for the 16 parameter configurations considered. Note that C(k, h 0 ) and D(k, h0) are always positive in the parameter

Figure 6 :

 6 Figure 6: Dimensionless lift-induced thrust f * l,x (a) and pressure drag f * d,x (b) as function of time t * , for a heaving amplitude h 0 = 0.6 and frequency k = 0.75, at Re = 5 10 4 . The theoretical formulas f * thrust = C 4πk 2 h 2 0 cos 2 (2kt * ) and f * drag = -D 2kh 0 cos(2kt * ) are superimposed as the broken line, for C = 0.419, D = 0.637.

Figure 7 :

 7 Figure 7: Time averaged pressure force F *x against its model counterpart (16), the symbols representing the values for the 16 parameter values (k, h 0 ), with k = 0.25, 0.5, 0.75, 1 and h 0 = 0.4, 0.6, 0.8, 1, the solid line being the bisector of the first quadrant angle.

Figure 8 :

 8 Figure 8: Sketch (left) of the heaving-pitching foil, in its upward motion (solid contour) and downward motion (broken contour). The snapshot (right) is the vorticity field for the parameters (k, h 0 , α 0 ) = (0.5, 0.8, 20 ○ ), the foil being in its upward motion at θ(t) ≈ -0.5 α 0 and h * (t) ≈ 0.866 h 0 .

Figure 9 :

 9 Figure 9: Oscillatory pressure force f * x (a) and viscous force f * µ,x (b) for the foil undergoing heaving-pitching motion, for the frequency k = 0.5, heaving amplitude h 0 = 0.8 and four different pitching angles α 0 = 0, 10 ○ , 20 ○ , 30 ○ .

Figure 11 :

 11 Figure11: Sketch of the experimental cycloidal propulsion device designed by the Naval Academy Research Institute in France. Three NACA0018 foils labeled 1, 2, 3 rotate in a circular motion around an axis. The foils' trajectory is marked as the circle, the angle between two foils being 120 ○ . The foils sketched with solid contours are at their positions at t 0 (when foil 1 has the highest pitching amplitude α 0 ) and the hatched contours indicate their positions T 6 later, T being the period of rotation.

Figure 12 :

 12 Figure 12: Snapshot of the vorticity field for three rotating and pitching foils, for the parameter values h 0 = 1.143, k = 0.6 and α 0 = 22 ○ . This configuration corresponds to that sketched in figure 11, the Reynolds number being Re = 5 10 4 .

Figure 13 :

 13 Figure 13: Comparison between the pressure force for one pitching-rotating foil and the pressure force around one foil pitching and rotating together with two others. Component in the x-direction (a) and in the y-direction (b), as function of angle β (by convention, β = 0 ○ , 360 ○ , ... is the highest vertical position for the successive periods of rotation, see figure11). The parameters are k = 0.6, h 0 = 1.143, α 0 = 22 ○ .

5 10 5 .

 5 The dimensionless heaving amplitude is h 0 = R L = 1.143, the reduced frequencies considered varying between k = 0.25 and k = 0.45. The pitching angle α 0 range is between 10 ○ and 40 ○ . Besides the Reynolds number, which in the present simulations is Re = 5 10 4 ,

Figure 14 :

 14 Figure 14: Comparison between the pressure force for one pitching-heaving foil and for a pitching-rotating foil. Component in the x-direction (a) and in the y-direction (b), as function of angle β (by convention, β = 0 ○ , 360 ○ , ... is the highest vertical position for the successive periods of rotation, see figure 11). The parameters are k = 0.5, h 0 = 1, α 0 = 10 ○ .

Figure 15 :

 15 Figure 15: Time averaged pressure force F *x are shown as the symbols for a pitchingrotating foil against its model counterpart (20) for pitching-heaving foil with coefficients B j , j = 1, ⋯, 5 given by[START_REF] Hu | Two-dimensional simulation of the hydrodynamic performance of a cycloidal propeller[END_REF]. The solid line being the bisector of the first quadrant angle.

Figure 16 :

 16 Figure 16: Comparison between the pressure force components in the x-direction (a) and y-direction (b), at two different Reynolds numbers Re = 5 10 4 and Re = 2.8 10 5 , for one pitching-rotating foil for the parameter values k = 0.6, h 0 = 1.143, α = 22 ○ .

Figure 17 :

 17 Figure 17: Comparison between measurements for the cycloidal device and the simulation results: pressure component in the x-direction (a) and in the y-direction (b), at Re = 2.8 10 5 , with pitching angle α 0 = 30 ○ and frequency k = 0.29. The dimensionless radius of the circular trajectory is h 0 = 1.143. The results for one period of rotation is shown, β = 0 corresponding to the instant when one of the three foils is at the top (cf. figure 11).

Figure 18 :

 18 Figure 18: Computed individual force values f * x (a) and f * y (b), for each blade (labeled 1, 2, 3) in the device, β = 0 ○ , 360 ○ , ⋯ corresponding to the highest vertical position of the individual blade's trajectory (four periods of rotation are shown). The parameters are the same as in figure 17.

Figure 19 :

 19 Figure 19: Comparison between measurements for the cycloidal device and the simulation results: pressure component in the x-direction (a) and in the y-direction (b), at Re = 2.8 10 5 , with pitching angle α 0 = 30 ○ and frequency k = 0.49. The dimensionless radius of the circular trajectory is h 0 = 1.143. The results for one period of rotation is shown, β = 0 corresponding to the instant when one of the three foils is at the top (cf. figure 11).

1 2 3Figure 20 :

 220 Figure 20: Time averaged pressure force F *x are shown as the symbols for a pitchingrotating foil against its model counterpart (20) for pitching-heaving foil with coefficients B j , j = 1, ⋯, 5 given by (21), at Re = 10 4 ( * ); Re = 5 10 4 (△) and Re = 2.8 10 5 (○). The solid line being the bisector of the first quadrant angle.

Figure 21 :

 21 Figure 21: Efficiency parameter η computed as the formula (24) versus thrust F * x , for different pitching-heaving parameters with h 0 = 1. The symbols represent the different pitching angles and the solid curves are line-fits through the values at the different frequencies k. The Reynolds number of the computations is Re = 5 10 4 .

  * y 2k cos(2kt * )⟩ + ⟨α 0 m * z 2k sin(2kt * )⟩.

Table 1 :

 1 Main numerical simulation options.

	Options	OpenFOAM v1906, overset, cellVolumeWeight, pimpleFoam
	Reynolds numbers 5 10 4 (2.8 10 5 for comparison with experiments)
	Velocity boundary	inlet: u = U ∞ e x
	and initial	other boundaries: zero normal gradient
	conditions	initial condition: u 0 = U ∞ e x
	Pressure boundary	inlet: zero normal gradient
	conditions	other boundaries: zero pressure

  Figure 10: Time averaged pressure force F *x against its model counterpart (

		1.4
		1.2
		1
		0.8
	F * x	0.6
		0.4
		0.2
		0
		-0.2	2	0	0.2	0.4	0.6	0.8	1	1.2	1.4

Acknowledgment

Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resources.

imental device developed by the French Naval Academy Research Institute