Si-Min Yuan

A-Li Chen

Liyun Cao

Hua-Wei Zhang

Shi-Wang Fan

Badreddine Assouar
email: badreddine.assouar@univ-lorraine.fr&alchen@bjtu.edu.cn

Yue-Sheng Wang

Tunable multifunctional fish-bone elastic metasurface for the wavefront manipulation of the transmitted in-plane waves

Keywords: Tunable, elastic metasurface, in-plane wave, Generalized Snell's law, simultaneous manipulation, multifunctions

In this work, a tunable fish-bone elastic metasurface is used to manipulate the in-plane waves which are more complex than flexural or out-of-plane waves because of the coupling of the P-and SV-waves. First, the refracted generalized Snell's law (GSL) for mixed in-plane wave is deduced which is also testified by the numerical results. Then according to the GSL, the elastic metasurfaces are designed to realize multifunctions such as directional refraction and wave focusing for different working frequencies. The critical refraction is discussed for the directional refraction and the behavior of SV-wave is predicted when the P-wave is designed to focus at any point by the metasurface.

Introduction

Manipulating elastic waves propagating along a desired way has attracted a lot of substantial attentions. The metamaterials proposed in the recent years provide new ideas and methods for controlling of the elastic waves [START_REF] Ma | Polarization bandgaps and fluid-like elasticity in fully solid elastic metamaterials[END_REF][START_REF] Lai | Hybrid elastic solids[END_REF][START_REF] Wu | Elastic metamaterials with simultaneously negative effective shear modulus and mass density[END_REF]. However, metamaterials cannot meet the increasing requirements of integration and miniaturization in modern devices because of their bulky and heavy characteristics. The recently proposed concept, metasurface, a kind of ultrathin metamaterial, provides broad prospects for elastic wave manipulation. It is compact, easily fabricated, and taking up less physical space than three-dimensional (3D) metamaterials [START_REF] Holloway | An overview of the theory and applications of metasurfaces-the two dimensional equivalents of metamaterials[END_REF]. It can modulate elastic waves by modifying their phase states in a predefined and desired manner based on the generalized Snell's law (GSL) [START_REF] Yu | Light propagation with phase discontinuities: generalized laws of reflection and refraction[END_REF]. For the first time, the conception of the metasurface was proposed in optics [START_REF] Chen | A review of metasurfaces: physics and applications[END_REF][START_REF] Li | Principles and applications of metasurfaces with phase modulation[END_REF] and then extended to acoustics [START_REF] Assouar | Acoustic metasurfaces[END_REF][START_REF] Liang | Wavefront manipulation by wavefront manipulation by acoustic metasurfaces: from physics and applications[END_REF]. Up to now, significant progress has been made in the research of both kinds of metasurfaces. They have exhibited extraordinary wave manipulating properties absent in natural materials, such as anomalous reflection and refraction [START_REF] Pfeiffer | Metamaterial huygens' surfaces: tailoring wave fronts with reflectionless sheets[END_REF][START_REF] Li | Theory of metascreen-based acoustic passive phased array[END_REF][START_REF] Li | Reflected wavefront manipulation based on ultrathin planar acoustic metasurfaces[END_REF], self-accelerating beams [START_REF] Greenfield | Accelerating light beams along arbitrary convex trajectories[END_REF][START_REF] Zhang | Generation of acoustic self-bending and bottle beams by phase engineering[END_REF][START_REF] Zhao | Delivering sound energy along an arbitrary convex trajectory[END_REF][START_REF] Zhu | Implementation of dispersion-free slow acoustic wave propagation and phase engineering with helical-structured metamaterials[END_REF], focusing [START_REF] Chen | Generation of wavelength-independent subwavelength Bessel beams using metasurfaces[END_REF][START_REF] Xie | Coding acoustic metasurfaces[END_REF][START_REF] Zhao | Continuously Tunable Acoustic Metasurface for Transmitted Wavefront Modulation[END_REF][START_REF] Fan | Tunable Broadband Reflective Acoustic Metasurface[END_REF][START_REF] Zhu | Reconfigurable origami-inspired metamaterials for controllable sound manipulation[END_REF][22] and so on. Recently, the extension to elastic metasurface has received wide attention from researchers.

The studies of elastic metasurfaces can be classified according to the types of elastic waves (out-of-plane wave, in-plane wave and surface wave). Zhu and Semperlotti [START_REF] Zhu | Anomalous refraction of acoustic guided waves in solids with geometrically tapered metasurfaces[END_REF] have presented the first numerical and experimental demonstration of elastic metasurfaces with geometric tapers to manipulate the phase profiles of Lamb waves. After then, flexural wave metasurfaces have become a research hotspot due to the wide applications of flexural waves in engineering. Such as the source illusion of flexural waves with a zigzag-shaped metasurface [START_REF] Liu | Source illusion devices for flexural lamb waves using elastic metasurfaces[END_REF], steering flexural waves with the metasurface constituted by the thin composite beams [START_REF] Zhang | Metasurface constituted by thin composite beams to steer flexural waves in thin plates[END_REF], and retroreflecting of a flexural wave with a simplified designed elastic metasurface [START_REF] Ruan | Retroreflection of flexural wave by using elastic metasurface[END_REF]. Cao et al. also did a lot of works on flexural wave regulation, such as realizing flexural wave deflection [START_REF] Cao | Deflecting flexural wave with high transmission by using pillared elastic metasurface[END_REF][START_REF] Cao | Pillared elastic metasurface with constructive interference for flexural wave manipulation[END_REF][START_REF] Cao | Disordered Elastic This is the author's peer reviewed[END_REF], flexural wave asymmetric transmission [START_REF] Cao | Asymmetric flexural wave transmission based on dual-layer elastic gradient metasurfaces[END_REF], and flexural wave absorption [START_REF] Cao | Flexural wave absorption by lossy gradient elastic metasurface[END_REF]. As for the in-plane wave metasurface, Hu's team has designed a composite-plate-based metasurface that can steer the wavefront of SV-waves [START_REF] Zhang | Wavefront steering of elastic shear vertical waves in solids via a composite-plate-based metasurface[END_REF]; and the asymmetric transmission can be achieved by combining this kind metasurface with phononic crystals [START_REF] Zeng | Asymmetric transmission of elastic shear vertical waves in solids[END_REF]. Lee et al. [START_REF] Lee | Mass-stiffness substructuring of an elastic metasurface for full transmission beam steering[END_REF] have designed metasurfaces with the mass-stiffness substructure to realize anomalous steering and focusing of in-plane longitudinal waves. Furthermore, some works have reported on manipulating in-plane P-and SV-waves simultaneously, e.g. splitting of P-wave and SV-wave [START_REF] Su | Elastic metasurface for splitting SV-and P-waves in elastic solids[END_REF][START_REF] Zheng | Non-resonant metasurface for broadband elastic wave mode splitting[END_REF], converting P-wave to SV-wave [START_REF] Kim | Transmodal elastic metasurface for broad angle total mode conversion[END_REF] and manipulating multimodal ultrasonic guided waves [START_REF] Tian | Elastic phased diffraction gratings for manipulation of ultrasonic guided waves in solids[END_REF]. Recently, topology optimization has been used to manipulate longitudinal and shear waves [START_REF] Rong | Multifunctional elastic metasurface design with topology optimization[END_REF]. Besides, the surface wave metasurface with a non-resonant meta-ridge was proposed to deflect Rayleigh surface waves [START_REF] Xu | Deflecting Rayleigh surface acoustic waves by a meta-ridge with a gradient phase shift[END_REF].

However, the aforementioned elastic metasurfaces are untunable, i.e., it can only operate for a fixed or an extremely narrow frequency range, which will limit their application [START_REF] Wang | Tunable and active phononic crystals and metamaterials[END_REF]. Hence, it is of practicale importance to study the continuously tunable elastic metasurface [START_REF] Li | Tunable modulation of refracted lamb wave front facilitated by adaptive elastic metasurfaces[END_REF][START_REF] Chen | A programmable metasurface for real time control of broadband elastic rays[END_REF][START_REF] Xia | In situ steering of shear horizontal waves in a plate by a tunable electromechanical resonant elastic metasurface[END_REF][START_REF] Yuan | Switchable multifunctional fish-bone elastic metasurface for transmitted plate wave modulation[END_REF]. Li et al. [START_REF] Li | Tunable modulation of refracted lamb wave front facilitated by adaptive elastic metasurfaces[END_REF] and Chen et al. [START_REF] Chen | A programmable metasurface for real time control of broadband elastic rays[END_REF] used piezoelectric materials to achieve the tunability of the lamb wave and flexural wave, respectively. And a tunable electromechanical resonant metasurface is proposed to steer SH waves in a plate [START_REF] Xia | In situ steering of shear horizontal waves in a plate by a tunable electromechanical resonant elastic metasurface[END_REF]. In our previous work, a switchable multifunctional fish-bone elastic metasurface was proposed to manipulate the transmitted flexural waves [START_REF] Yuan | Switchable multifunctional fish-bone elastic metasurface for transmitted plate wave modulation[END_REF]. In the present research, the transmitted P-and SV-waves will be manipulated simultaneously at a wide frequency range based on the tunable structure. In other words, the metasurface can be designed at each frequency within a frequency band to realize the manipulation of in-plane waves by changing the screw-in depth for different frequencies. Different from the flexural wave, the coupling between P-and SV-waves needs to be considered in the case of the in-plane waves where the relation between the phase changes for the P-and SV-waves is unknown. We first derive the generalized Snell's law for the in-plane waves and then numerically observe the phase shift and transmittance of the functional units. In addition, the mechanism of phase-shift tunability is revealed analytically. The designed metasurface can be designed to manipulate P-and SV-waves simultaneously.

The refracted generalized Snell's law of the mixed in-plane wave

First, we need to get the generalized Snell's law of the refracted waves in case of the mixed in-plane modes. As shown in Fig. 1, the surface x = 0 between media A and B is the metasurface which provides the phase shifts

) 1
where A is the polarization vector; k is the wavenumber;  is the incident angle; c is the wave velocity; the subscripts P or SV represent for the P-or SV-wave, respectively. After the incident wave passes through the metasurface, the wave mode transformation will occur. Then the wave fields in medium B of the transmitted P-and SV-waves can be expressed as

i k z x c t i k z x c t e e           + - +     + - +  =  = PP SV SV uA uA , (2)
where the superscript "  " means the variables for the transmitted waves. Taking the P-wave incidence as an example, according to the continuous conditions at the interface (i.e., x = 0), we can get ( sin )= ( sin ) = ( sin )

P P P P P P P SV SV SV SV k z c t k z c t k z c t            - - +  - +  . (3) Considering
= P P P P SV SV k c k c k c      ==
with ω being the angular frequency, we can obtain the following relations: sin = sin = sin

P P P P P SV SV SV k z k z k z          +  +  . ( 4
)
Taking the derivation w.r.t z on both sides of Eq. ( 4), we can get the generalized Snell's law for the P-wave incidence

           = + = + . ( 5
)
Similarly, the generalized Snell's law for the SV-wave incidence can be obtained as:

d d sin sin sin dd SV P SV SV P P SV SV k k k zz            = + = + . ( 6
)
Fig. 1 Schematic illustration of the GSL for the transmitted waves of the mixed in-plane modes.

From Eq. ( 5) or [START_REF] Chen | A review of metasurfaces: physics and applications[END_REF]

                - + - +  - + - +  . ( 7
)
Substituting Eq. ( 4) into Eq. ( 7), we can obtain cos cos ( cos cos )

P SV P P SV SV P P SV SV k x k x x k k               - = - = -
which implies that P SV at an arbitrary point (x0) is independent of the variable z, i.e.

( ) 0

P SV z   -=  . ( 9
)
Eq. ( 9) holds in the whole transmision region (x0). Considering the fact that the phase difference between P-and SV-waves (

P SV  -) is equal to P   - SV  
at the side of the metasurface facing the transmission region (i.e. x → 0 + ), we obtain, from Eq. ( 9),

d d d d SV P z z     = , (10)
Substituting Eq. ( 10) into Eq. ( 5) or Eq. ( 6), we can get sin sin

P P SV SV k k       = , (11)
i.e., sin == sin

SV PP SV P SV k c kc        . (12)
The above formula shows that the relations between the refraction directions of the P-and SV-waves are the same as the classical Snell's law and are independent of the incident wave. Eq. ( 12) also shows that the refraction angle of the P-wave is greater than that of the SV-wave.

Therefore, when the refraction angle of the P-wave Eqs. ( 5), ( 6) and ( 12) can be applied to design elastic metasurfaces for wave front modulation. If one steers an incident P (or SV) wave, Eq. ( 5) [or [START_REF] Chen | A review of metasurfaces: physics and applications[END_REF]] should be used; and the propogation direction of the SV (or P) wave can be determined by Eq. [START_REF] Li | Reflected wavefront manipulation based on ultrathin planar acoustic metasurfaces[END_REF].

The metasurface and tunable functional units

According to the above GSL for the mixed in-plane wave modes, the metasurface can be designed to manipulate the transmitted in-plane waves. First, the functional unit providing 2π phase shift should be designed. As shown in Fig. 2, the functional unit with the finite element mesh as designed in Ref. [START_REF] Yuan | Switchable multifunctional fish-bone elastic metasurface for transmitted plate wave modulation[END_REF] is used to manipulate the in-plane wave mode. The values of the geometric parameters are shown in Table I. It can be seen that the shape of the functional unit (Fig. 2(a)) looks like a fish-bone with 6 pairs of screw bolts (bones) of equal length and two symmetrical rows of nuts as vibrators, so the metasurface is named as a "fish-bone metasurface". As we know that the position of the nut on each screw-bolt can be adjusted by screwing in or out according to the screw-and-nut mechanism, then the phase shift and transmittance of the transmitted wave can be tuned continuously. The orange and green regions represent the host (aluminum) and vibrators (steel), respectively. The material parameters of the aluminum are: mass density ρ = 2700kg/m 3 , longitudinal wave velocity cP = 5351m/s and shear wave velocity cSV = 3097m/s, and those of the steel are: ρ = 7670kg/m 3 , cP = 6010m/s and cSV = 3150m/s.

It should be noticed that for the manipulation of the refracted waves, not only the phase shift but also the efficiency of the functional unit needs to be guaranteed at the same time. Then the phase shift and the transmittance of the functional unit varying with the screw-in depths (s) are calculated by using the finite element software COMSOL Multiphysics 5.3. A line source with in-plane excitation is set on the left side far from the "fish-bone" structure to generate a stable incident wave. The excitations along the x-axis and z-axis simulate the P-and SV-waves, respectively. Perfectly matched layers (PMLs) are used on both ends of the functional unit to avoid any reflection from the boundaries. Continuous periodic boundary conditions are applied to the upper and lower boundaries of the functional unit, as shown in Fig. 2(a). The phase shift and the transmittance are obtained by taking the averaged displacements along a line (detection line) in the far area of the right side [START_REF] Li | Tunable modulation of refracted lamb wave front facilitated by adaptive elastic metasurfaces[END_REF]. As stated above, the phase shift and transmittance properties for P-or SV-waves can be calculated in the case of the P-or SV-wave incidence, and are illustrated in Figs. 3 and4, respectively. The transmittances of P-wave ( P  ) and SV-wave ( SV  ) are defined as 0

P P P A A  = and 0 SV SV SV A A  = , (13)
where AP and ASV are the transmitted displacement amplitudes of P-and SV-waves for the functional unit, respectively; and AP0 and ASV0 are those for the uniform structure without the fish-bone structure. From the nephograms of the phase shift for the transmitted P-and SV-waves, as shown in Figs. 3(a) and 4(a), it can be found that the phase shift indeed has a full 2π coverage in a wide frequency range (28 -42 kHz). The corresponding wavelengths for the P-and SV-waves are 191 -127 mm and 111 -74 mm, respectively. The thickness of the metasurface (l = 48 mm) is less than the wavelength. The transmittance is higher than 70% below the oblique dashed line, but is very low above this line, see Figs. It is clearly seen that the phase shift covers a full 2π span and the transmittance is generally higher than 70%.

Mechanism of phase-shift tunability

When the wave propagates through the "fish-bone" functional unit, the phase shift can be obtained by

where cphase = 2πf/k is the phase velocity with k being the wave number. For a given frequency, cphase can be calculated from the dispersion curves of the infinite 1D periodic structure in which the basic element is the unit cell, as shown in Fig. 5(a). For different screw-in depths of the nuts, the dispersion curve can be calculated by using the finite element method together with the Bloch theorem. The results in the first Brillouin zone of the 1D periodic structures with the screw-in depth being 1.7, 3.0, and 4.5 mm are shown in Fig. 5(b)-(f). The color scale represents the value of

222 d d i V i x y z V uV u u u V  = ++   , ( 15
)
where i represents for x, y and z; ux, uy and uz are the displacement components in x, y, and z directions, respectively; and the integrals are taken over the volume (V) of one basic element. As shown in Figs. 5(b),(c) and (e), i → 1 (the dark red color) represents for the ui -dominated wave mode, i.e., y → 1, x → 1 and z → 1 means the flexural-, P-and SV-dominated wave modes, respectively. This also can be tested by the vibrational modes. Here, in the case of f = 35.6 kHz, Figs. 5(d) and (f) show the vibrational modes of points 1, 2 and 3, 4 which are the vibrational modes of the P-and SVdominated wave modes, respectively. It should be noticed that for the flexure-dominated wave studied in Ref. [START_REF] Yuan | Switchable multifunctional fish-bone elastic metasurface for transmitted plate wave modulation[END_REF], the operating frequency zone is 7 -13 kHz. However, for the in-plane waves, the operating frequency zone should be as high as 28 -42 kHz to ensure that the phase shifts of P-and SV-waves cover 0-2π simultaneously. Bands a and b are considered for Pand SV-dominated waves as shown in Figs. 5(c) and (e), respectively. The phase velocity can be calculated easily from the crosspoints on the dispersion bands. The phase velocity curves of the Pand SV-wave modes varying with s at f = 35.6 kHz are illustrated in Figs. 6(a) and (b), respectively. Then according to Fig. 6 and Eq. ( 14), the phase shifts can be obtained which are illustrated in Figs.

Design of the metasurface for manipulating P-and SV-waves simultaneously

In this section, we will design the metasurfaces based on the GSL for the mixed in-plane waves and the nephograms of phase shift and transmittance to realize the manipulation of the transmitted P-and SV-waves. The calculation model in a front view (xz-plane) is shown in Fig. 7. The boundaries aroud the model are set as PMLs. A surface traction along the x-or z-direction with the amplitude of

22 0 1 exp( ( ) / ) z H w    --  (N/m 2
) is applied at the yellow dashed line to excite a

Gaussian beam of the P-or SV-wave mode. Multi-functions such as directional refraction (α = 16.5 and w0 = 30H) and wave focusing (α = 0 and w0 = 40H) will be realized numerically.

Fig. 7 The model of the metasurface constructed in COMSOL Multiphysics.

Directional refraction

For the transmitted wave propagating along the designated direction, the required phase shift of each functional unit can be calculated from Eq. ( 5) and Eq. ( 6). Once the phase shifts are known, we can obtain the corresponding screw-in depths from Figs. 3(a) and 4(a). Then the metasurface with the function of directional refraction can be obtained.

a. Normal incidence of P-wave

First, we consider the normal incidence of a P-wave and design the metasurface based on Fig. 3(a) to control the refraction angle of the P-wave being 30° with f = 31.4 kHz. The phase shift curves of the P-wave in Fig. 3(a) are used to determine the screw-in depths of every functional unit of the metasurface (refer to Appendix A for the detailed values). Figs. 8(a) and (c) are the divergence and curl of displacement fields, respectively, which represent the P-and SV-wave fields; the corresponding polar directivity graphs of the transmitted waves measured at far fields are presented in Figs. 8(b) and(d). We can see that the metasurface can steer the P-wave nearly in the designed direction (red dot-dashed line). Due to the limited width of the incident beam and the influence of the near-field scattering, there is a little deviation between the numerically estimated direction (yellow arrow) and the desired direction (red dot-dashed line). And at the same time, the refraction angle of the SV-wave calculated by Eq. ( 12) is 17°; the estimated refraction angle by the curl of displacement fields is close to this value. Then, with the same target of the regulation (i.e. the refraction angle of the P-wave is 30°), the metasurface is tuned for a different operating frequency f = 35.6 kHz (the screw-in depths are presented in Appendix A). The numerical results for the divergence (P-wave) and curl (SV-wave) of displacement fields are shown in Figs. 8(e) and (g), respectively, with the corresponding polar directivity graphs shown in Figs. 8(f) and (h). It can be seen that the refraction direction is independent of the incident wave frequencies. The results verify the correctness of Eq. ( 12), i.e., for the mixed in-plane mode, if the propagation direction of one wave mode in the transmitted medium is decided, the direction of another wave mode is then determined. In other words, the transmitted P-and SV-waves are related by Eq. ( 12). Next, we try to steering the refracted P-wave in the direction of 90°. What we need to do is to tune the screw-in depth (shown in Appendix A) of every functional unit of the same metasurface to get the phase shift we need. In this case, the numerical results of the divergence and curl of displacement fields are illustrated in Figs. 9(a) and (b), respectively. It can be seen that the refracted P-wave propagates approximately along the desired direction (red dotted line) as shown in Fig. 9(a).

If the refraction angle is large, the wave front is slight curved under the influence of the near-field scattering. In addition, it is noted that the incident beam has a limited width, so the propagating-to-evanescence transformation is only effective within the beam width. The transmitted wave will spread out outside the beam width, especially when the target transmission angle is 90°. Therefore, it's very difficult to realize perfect grazing. Fig. 9(b) shows that the refracted angle of the SV-wave is also close to 35° which is calculated by Eq. ( 12) (red dot-dashed line). Generally speaking, these results show that the metasurface can steer the refracted P-wave to propagate in the targeted direction with the refraction angle of the SV-wave determined by Eq. ( 12). Eq. ( 12), which is derived from Eq. ( 10), is verified by the above numerical simulation. Next we will check Eq. ( 10) directly through comparision of P   and SV   . Fig. 10 shows the phase shifts of every functional unit when the refraction angle of the P-wave is 30° with f = 35.6 kHz. The phase shift of the P-wave (blue dots) for each functional unit is calculated from the first two items of Eq. ( 5). Then the screw-in depths of every functional unit can be obtained from the linear interpolation in Fig. 3(a). The phase shift of the SV-wave (red rhombuses) corresponding to these screw-in depths can be obtained from Fig. 4(a) directly. We can see from Fig. 10 that the slopes of the two results are approximately equal, which verifies the reliability of Eq. [START_REF] Pfeiffer | Metamaterial huygens' surfaces: tailoring wave fronts with reflectionless sheets[END_REF].

b. Normal incidence of SV-wave

In the above section, the steering of the P-wave is discussed for the P-wave incidence. In this section, the normal incidence of an SV-wave will be considered. Steering of both P-and SV-waves will be examined based on the phase shift curves of SV-wave (Fig. 4). We first consider the steering of the refracted P-wave in the angle of 30° for f = 35.6 kHz. In this case, the refraction angle of the SV-wave can be obtained by Eq. ( 12), which is 17° (red dot-dashed line). Then the phase shift curves of SV-wave (Fig. 4) are used to obtain the phase shifts and the screw-in depth of every functional unit to design the metasurface. The numerical results are shown in Fig. 11, in which (a) is the divergence of displacement field and (c) is the curl of displacement field. Figs. 11(b) and(d) are the corresponding polar directivity graphs. It is seen that our target is largely achieved, i.e., the refracted P-wave propagates at 30° with the SV-wave refracted at 17°. If the refracted P-wave is steered in the direction along the metasurface ( 90 P   = ), the refraction angle of SV-wave should be 0 =35 SV    based on Eq. ( 12). The numerical simulations are shown in Figs. 11(e) and (g) with the corresponding polar directivity graphs shown in Figs. 11(f) and(h). The evanescent behavior of the refracted P-wave is seen clearly. There are three refracted beams (yellow arrows) appearing in the transmitted field in Fig. 11(g). One is the refracted SV-wave as expected, which propagates in the theoretical direction of 35° obtained by Eq. ( 12) (red dot-dashed line). And the remaining two beams are the leaked waves caused by higher-order diffraction [START_REF] Cao | Steering elastic SH waves in an anomalous way by metasurface[END_REF]. It can be seen that the three directions shown in the polar directivity graph Fig. 11(h) agree with those in Fig. 11(g). The appearance of these undesired waves indicates that our design is not perfect yet and should be optimized. However, the optimization of a tunable metasurface is a challenging problem and deserves further research.

Fig. 11 Normalized divergence (a, e) and curl (c, g) of displacement fields for the normal incidence of SV-wave at 35.6 kHz. The desired refraction angles of P-wave are 30° (a) and 90° (e), and the refraction angles of SV-wave obtained from Eq. ( 12) are 17° (c) and 35° (g), respectively. The polar directivity graphs (b, d, f and h) of the transmitted waves measured at far fields correspond to the numerical simulation results (a, c, e and g), respectively.

It was indicated previousely that when the refraction angle of the SV-wave is bigger than the critical value 0 SV  , the refracted P-wave will becomes evanescent. For example, if we target the refraction angle of the SV-wave at 60° which is bigger than 0 =35 SV    , the refracted wave P-wave is observed to be transferred into the evanescent wave, see Figs. 12(a) and (b). When the refracted SV-wave is steered along the surface (i.e. 90 SV   = ), all transmission waves become evanescent, see Figs. 12(c) and(d). For the same reason as we explained before, the simulated directions (yellow arrows) deviate a bit from the desired directions (red dot-dashed lines), as shown in Figs. 12(a) and(c). In Fig. 12(d), the downstream wave is the leaked wave, similar to the situation in Fig. 11(g). In thses two cases, the metasurfaces are designed based on the phase shift curves of SV-wave in Fig. 4. The detailed values of the screw-in depths of the functional units are presented in Appendix A.

c. Oblique incidence of SV-wave

Consider the oblique incidence of the SV-wave at 15° with the operating frequency of 35.6 kHz. The targeted refraction angle of the SV-wave is 30° or 60° (red dot-dashed line). The results are shown in Fig. 13 which demonstrate the efficiency of the designed metasurface (the detailed screw-in depths of the units are listed in Appendix A). For the 30° refraction of the SV-wave, the P-wave propagates at the angle of 60° which is close to the value calculated from Eq. ( 12) (red dot-dashed line). Fig. 13(d) shows that the refracted P-wave becomes evanescent when the targeted refraction angle of the SV-wave is 60° which is larger than 0 SV  . This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.

PLEASE CITE THIS ARTICLE AS

Wave focusing

In this part, we design metasurface to realize the focusing of the refracted P-wave. The central focusing and arbitrary point focusing are considered. For the central focusing, the schematic diagram is illustrated in Fig. 14. The red dashed lines represent the P-wave, and the black solid lines represent the SV-wave. We suppose that the P-wave is focused at point O1 (x, 16.5H). Then according to the GSL, the required phase shift of each functional unit can be determined by [START_REF] Xia | In situ steering of shear horizontal waves in a plate by a tunable electromechanical resonant elastic metasurface[END_REF]:

( ) 2 2 2 2 1 1 2 ( 16.5 ) ( 16.5
) nn z H x z H x    -= - + - - + , (16)
where n is the serial number of the functional unit. It should be noticed that the metasurface should be symmetric about z = 16.5H for the central focusing, so we can just calculate the phase shifts of half the functional units (e.g., n = 2-17) as shown in Fig. 14. By combining Fig. 3(a), the screw-in depth for each functional unit of the metasurface can be obtained to realize the focusing of the refracted P-wave.

Then what happens to the refracted SV-wave? According to Eq. ( 12), the behavior of the SV-wave can be estimated as shown in Fig. 14. Assume that for the nth functional unit the refraction angle of the P-wave is Pn   , and that of the SV-wave is SVn   . Then 16.5 ( 1) tan

Pn H n H x  --  = . ( 17
)
From Eq. ( 12), we can get sin sin

SV Pn SVn P c c     =  . ( 18
)
As shown in Fig. 14, the coordinate of the crosspoint O2 (x′, 16.5H) of the SV-wave beams and the central axis (z = 16.5H) can be determined as

( 1) tan

SVn H n H x  --  =  . ( 19
)
Table II Similarly, the focusing of SV-waves can also be realized for the SV-wave incidence. As shown in Fig. 16, when the SV-wave with the frequency of 35.6 kHz is incident, the phase shift and screw-in depth of each functional unit are obtained by the curves of SV-wave (Fig. 4), and then the metasurface is designed to make the SV-wave focus on (250 mm, 330 mm). Next the wave focusing at an arbitrary point O3(x′′, y′′) for the refracted P-wave will be discussed. In this case, the metasurface is not symmetric. The phase shifts of the functional units can be calculated by [START_REF] Xia | In situ steering of shear horizontal waves in a plate by a tunable electromechanical resonant elastic metasurface[END_REF]:

( ) 2 2 2 2 1 1 2 ( ) ( ) nn z z x z z x        -= - + - - + , ( 20
)
where n = 2-34. The schematic diagram is shown in Fig. 17. Considering the transmitted wave beam from nth unit, we have

( 1) tan Pn n H z x   --  =  . (21)
According to Eq. ( 12), we can get sin sin

SV Pn SVn P c c     =  . ( 22
)
In that way the propagation direction of the SV-wave is determined when the P-wave is focused on (200 mm, 150 mm). On each side of the line EF there are energy dense areas for the SV-wave, as shown in Fig. 17

Conclusions

In this paper, we design a kind of tunable fish-bone elastic metasurface to manipulate transmitted P and SV-waves. We first present the generalized Snell's law for the mixed in-plane wave modes including the relation between the refraction angles of the P-and SV-waves, which allows us to predict another kind of wave's propagation behavior when we regulate one of the waves. The phase shifts of the designed tunable functional unit are obtained numerically and analytically. Then, we design the metasurfaces to achieve directional refraction and focusing for both P-and SV-wave incidence. The distinguishing feature of our design is that the metasurface is continuously tunable and broadband. The present study provides us a way to control mixed in-plane waves.

The present paper only consider the modulation of the transmitted waves. Another interesting topic is the steering of the reflected wave front, especially the simultaneous modulation of refracted and reflected waves which is very challenging.

  and SV-waves, respectively. The displacements of the incident P-or SV-wave in medium A are assumed as



  beyond which the refracted P-wave becomes evanescent.

Fig. 2

 2 Fig. 2 Schematic diagrams of the functional unit (a) and basic element (b).

Fig. 3 Fig. 4

 34 Fig. 3 The phase shift and transmittance of the transmitted P-wave. The nephograms of phase shift (a) and transmittance (b) varying with the screw-in depth and frequency for the functional unit; (c) the phase shifts obtained by numerical simulation (the black solid curve) and analytical calculation (the green dashed curve) at f = 35.6 kHz; (d) the transmittance obtained by numerical simulation at f = 35.6 kHz.

  3(b) and 4(b). Therefore, the high transmission regions below the oblique lines will be considered in the following design. As an illustration, the phase shift curves for a particular frequency of f = 35.6 kHz are plotted in Figs. 3(c) and 4(c) by the black solid curves which are obtained directly by sectioning, respectively, Figs. 3(a) and 4(a) along the dot-dashed line. And the corresponding curves of the normalized transmittance are depicted in Figs. 3(d) and 4(d).

phase 2 /

 2 fl c  = ,

  3(c) and 4(c) by the green dashed curve. It is consistent with the black solid curve obtained by numerical simulation. It is also found that the black dot-dashed lines (f = 35.6 kHz) have no crosspoint with the dispersion curve a or b for s = 1.7 mm, meaning that f = 35.6 kHz locates in the bandgap of the P-or SV-wave mode, which explains why in this case the transmissions are very low as shown in Figs. 3(b) and 4(b).

Fig. 5 (Fig. 6

 56 Fig. 5 (a) The unit with Bloch periodic boundary conditions to calculate the dispersion curves; (b), (c) and (e) the dispersion curves of the above unit for the screw-in depth s = 1.7, 3.0 and 4.5 mm; (d) the P-dominated wave modes and (f) the SV-dominated wave modes corresponding to crosspoints "x" at f = 35.6 kHz for s = 3.0 and 4.5 mm.

Fig. 8

 8 Fig. 8 Normalized divergence (a, e) and curl (c, g) of the displacement fields for the normal incidence of a P-wave with frequencies f = 31.4 kHz (a, c) and 35.6 kHz (e, g). The designed refraction angle of the refracted P-wave is 30°. The polar directivity graphs (b, d, f and h) of the transmitted waves measured at far fields correspond to the numerical simulation results (a, c, e and g), respectively.

Fig. 9

 9 Fig. 9 Normalized divergence (a) and curl (b) of the displacement fields for the normal incidence of a P-wave with frequency f = 35.6 kHz. The designed refraction angles of the P-wave is 90°.

Fig. 10

 10 Fig. 10 Phase shifts for every functional unit of the metasurface to control the refraction angle of P-wave is 30° at 35.6 kHz.

Fig. 12

 12 Fig. 12 Normalized curl (a, c) and divergence (b, d) of displacement fields for the normal incidence of SV-wave at 35.6 kHz. The manipulation targets are (a, b) the refraction angle of SV-wave is 60°, and P-wave is transfered into evanescent wave; (c, d) The transmitted SV-wave is steered along the metasurface, and P-wave has aleardy become evanescent wave.

Fig. 13

 13 Fig. 13 Normalized curl (a, c) and divergence (b, d) of displacement fields for the obliquely incident SV-wave of 35.6 kHz with the incident angle of 15°. Steer the SV-wave at 30° (a) with the refraction angle of the P-wave being 60° (b); Steer the SV-wave at 60° (c) with the P-wave becoming evanescent (d).

  illustrates the values of x′ when the P-wave with the operating frequency f = 35.6 kHz focusing on O1 (x, 16.5H) = (150 mm, 330 mm). In this case, the results are shown in Fig.14. For the sake of clearance, we have omitted some beams. As shown in Fig.14, all the beams will pass through the wedge ABC. The wave beams near AB and AC are relatively dense, indicating that the energy near lines AB and AC is relatively concentrated. The corresponding numerical results computed by Comsol Multiphysics are shown in Fig.15. We can see that the P-wave indeed focuses on the targeted position, O1, as shown in Fig.15(b). And the size and position of the focus area of the SV-wave shown in Fig.15(c) are approximately the same as the analytical results. If the operating frequency changes, the metasurface still works just by tuning the screw-in depths of the functional units. Figs.15(d), (e) and (f) show the focusing results for another choice of frequency, f = 31.4 kHz.

Fig. 14

 14 Fig. 14 Schematic diagram of central focusing for the refracted P-wave. The red dashed and black solid lines represent the refracted P-and SV-waves, respectively.

Fig. 15

 15 Fig. 15 Results of central focusing for the P-wave incidence with f = 35.6 kHz (a-c) and f = 31.4 kHz (d-f) at O 1 (150 mm, 330 mm). (a, d) Normalized total energy fields, (b, e) normalized intensity of divergence of displacement fields, and (c, f) normalized intensity of curl of displacement fields.

Fig. 16

 16 Fig. 16 Results of central focusing for the SV-wave incidence with f = 35.6 kHz at O 1 (250 mm, 330 mm). (a) Normalized total energy field, (b) normalized intensity of divergence of displacement field, and (c) normalized intensity of curl of displacement field.

  . The simulated results for the P-wave focusing at O3(200 mm, 150 mm) in case of f = 35.6 kHz and 31.4 kHz are shown in Fig.18. These results indicate that wave focusing at an arbitrary point can be achieved by the designed metasurface, and there indeed have the energy dense areas for SV-wave as shown in Figs.18(c) and (f). The numerical simulated results show a good agreement with the analysis. The screw-in depths of all functional units for central focusing and arbitrary point focusing are shown in Appendix A.

Fig. 17

 17 Fig. 17 Schematic diagram of arbitrary point focusing for the refracted P-wave. The red dashed and black solid lines represent the refracted P-and SV-waves, respectively.

  we can see that if the relation between d /d

												 	z	and d  	/d z	is known,
												P	SV
	then the transmission directions for the P-and SV-waves can be determined directly. Next, we will
	derive the relation between d  	/d z	and d  	/d z	.
					P			SV
	The phase difference between the refracted P-wave (	 ) and SV-wave (	 ) at any point (x, z)
											P		SV
	with x0 is:
			=[ ( sin k z	x	cos	c t	)	] [ k	( sin z	x	cos	c t	)	]
	P	SV	P	P		P	P	P		SV	SV		SV	SV	SV

Table I .

 I Geometric parameters (unit: mm).

	l	H	h 1	h 2	h 3	d	d 1	r	a = l/6
	48	20	7.1	2	3	5.6	2.4	3.9	8

Table II .

 II The intersection points of the SV-wave beams and the center line (z = 16.5H). This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE

	n	1	2	3	4	5	6	7	8	9
	x′ (mm)	260	263	269	277	289	302	317	334	353
	n	10	11	12	13	14	15	16		17
	x′ (mm)	373	393	415	437	460	484	508	532

CITE THIS ARTICLE AS DOI: 10.1063/5.0029045

Table IV .

 IV The screw-in depths of the nuts in the nth functional unit for directional refraction. (unit: mm).

	sn	case 1	case 2	case 3	case 4	case 5	case 6	case 7	case 8	case 9
	s1	1.22	1.94	1.94	3.81	3.81	3.81	3.81	3.81	3.81
	s2	1.29	1.95	1.98	4.13	4.48	4.83	4.97	4.07	4.52
	s3	1.40	1.98	2.05	4.48	2.17	2.26	2.31	4.36	2.18
	s4	1.52	2.01	2.14	4.83	2.26	2.53	2.75	4.65	2.28
	s5	1.66	2.05	2.27	2.17	2.40	3.05	3.47	4.93	2.45
	s6	1.88	2.10	2.49	2.21	2.68	3.79	4.56	2.17	2.77
	s7	2.30	2.14	2.85	2.26	3.06	4.81	2.24	2.21	3.18
	s8	3.68	2.19	3.86	2.32	3.52	2.25	2.55	2.25	3.71
	s9	5.06	2.27	1.95	2.40	4.11	2.52	3.18	2.30	4.39
	s10	0.96	2.36	2.01	2.53	4.82	3.04	4.14	2.36	5.07
	s11	0.98	2.49	2.10	2.68	2.21	3.77	2.19	2.44	2.26
	s12	1.00	2.65	2.19	2.86	2.32	4.78	2.40	2.55	2.41
	s13	1.03	2.85	2.36	3.06	2.53	2.25	2.94	2.68	2.71
	s14	1.05	3.16	2.64	3.27	2.86	2.51	3.76	2.83	3.11
	s15	1.08	3.86	3.15	3.52	3.27	3.02	4.91	2.99	3.62
	s16	1.11	1.94	1.94	3.81	3.79	3.75	2.30	3.16	4.27
	s17	1.15	1.95	1.97	4.11	4.46	4.76	2.72	3.35	4.99
	s18	1.21	1.97	2.05	4.47	5.08	2.25	3.43	3.56	2.24
	s19	1.29	2.01	2.14	4.82	2.26	2.50	4.51	3.81	2.37
	s20	1.39	2.05	2.26	2.16	2.39	3.01	2.24	4.06	2.65
	s21	1.51	2.10	2.48	2.21	2.67	3.73	2.53	4.35	3.04
	s22	1.65	2.14	2.83	2.26	3.04	4.73	3.15	4.64	3.52
	s23	1.87	2.19	3.79	2.32	3.50	2.24	4.08	4.92	4.15
	s24	2.27	2.27	1.95	2.40	4.09	2.49	2.18	2.17	4.88
	s25	3.54	2.36	2.00	2.53	4.80	3.00	2.38	2.21	2.22
	s26	5.05	2.48	2.09	2.68	2.21	3.71	2.91	2.25	2.34
	s27	0.96	2.64	2.19	2.86	2.31	4.71	3.71	2.30	2.60
	s28	0.98	2.84	2.35	3.05	2.52	2.24	4.86	2.35	2.97
	s29	1.00	3.15	2.63	3.27	2.85	2.48	2.29	2.44	3.43
	s30	1.03	3.82	3.13	3.51	3.25	2.98	2.69	2.55	4.03
	s31	1.05	1.94	1.94	3.80	3.77	3.69	3.39	2.68	4.77
	s32	1.08	1.95	1.97	4.10	4.44	4.68	4.45	2.83	2.21
	s33	1.11	1.97	2.05	4.46	5.07	2.24	2.23	2.98	2.32
	s34	1.15	2.01	2.14	4.81	2.25	2.47	2.50	3.15	2.55

Table V .

 V The screw-in depths of the nuts in the nth functional unit for wave focusing. (unit: mm). This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset. PLEASE

			central focusing		arbitrary point focusing
	sn	P-wave incident	SV-wave incident	P-wave incident
		35.6 kHz	31.4 kHz	35.6 kHz	35.6 kHz	31.4 kHz
	s1	1.94	1.22	3.82	1.94	1.22
	s2	1.95	1.24	3.91	3.73	1.15
	s3	1.95	1.28	4.08	3.11	1.10

CITE THIS ARTICLE AS DOI: 10.1063/5.0029045

	s4	1.97	1.35	4.36	2.84	1.08
	s5	1.99	1.46	4.73	2.69	1.06
	s6	2.04	1.59	2.17	2.58	1.04
	s7	2.09	1.83	2.24	2.53	1.03
	s8	2.15	2.46	2.33	2.49	1.03
	s9	2.24	4.97	2.54	2.49	1.03
	s10	2.38	0.97	2.87	2.53	1.03
	s11	2.61	1.00	3.32	2.58	1.04
	s12	2.97	1.04	3.95	2.69	1.06
	s13	4.06	1.08	4.79	2.84	1.08
	s14	1.95	1.15	2.23	3.11	1.10
	s15	1.99	1.26	2.39	3.73	1.15
	s16	2.07	1.45	2.79	1.94	1.22
	s17	2.15	1.70	3.35	1.96	1.31
	s18	s17	s17	s17	1.99	1.47
	s19	s16	s16	s16	2.05	1.66
	s20	s15	s15	s15	2.12	2.02
	s21	s14	s14	s14	2.19	3.74
	s22	s13	s13	s13	2.32	0.95
	s23	s12	s12	s12	2.52	0.99
	s24	s11	s11	s11	2.82	1.03
	s25	s10	s10	s10	3.46	1.07
	s26	s9	s9	s9	1.95	1.12
	s27	s8	s8	s8	1.98	1.22
	s28	s7	s7	s7	2.04	1.37
	s29	s6	s6	s6	2.13	1.59
	s30	s5	s5	s5	2.22	1.99
	s31	s4	s4	s4	2.38	4.33
	s32	s3	s3	s3	2.66	0.96
	s33	s2	s2	s2	3.12	1.00
	s34	s1	s1	s1	5.04	1.05

This is the author's peer reviewed, accepted manuscript. However, the online version of record will be different from this version once it has been copyedited and typeset.PLEASE CITE THIS ARTICLE AS DOI: 10.1063/5.0029045

Acknowledgment

This work was supported by the National Natural Science Foundation of China (Grant Nos. 11872101, 11991031, 11991032, 12021002); the Fundamental Research Funds for the Central Universities (Grant No. 2018JBM032); and the China Scholarship Council (Grant No. 201907090040).

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A. The screw-in depths of the nuts

The screw-in depths of the nuts to achieve above manipulate targets are given in the following tables. Directional refraction is more complex, so it is summarized in Table III. Tables IV andV show the screw-in depths of the nuts for directional refraction and focusing, respectively.


Edgar Tournon
email: edgar.tournon@univ-lyon1.fr

Pascal Venet

Bertrand Barbedette

Adrien Lelievre

Judicael Aubry

Ali Sari

Judicaël Aubry

Unplugged Series Hybrid Bike Study: Cost versus Performance Analysis

Keywords: Electric bike, Series Hybrid Bike, Efficiency and financial costs optimization, Powertrain modelling

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

Series Hybrid Bikes (SHBs) are chainless bicycles. The mechanical energy from the cyclist is first converted into electrical energy thanks to an electrical generator located around the crankshaft. At the same time the motion is created by a motor assembled in the wheel which is using the energy produced by the user and/or the one stored in the Energy Storage System (ESS) [START_REF] Aubry | Evaluation de l'interêt d'une hybridation serie pour un vehicule a propulsion humaine. Pre-dimensionnement d'une generatrice-pedalier[END_REF], [START_REF] Guanetti | Optimal energy management in series hybrid electric bicycles[END_REF]. The powertrain architecture of SHB is described in the figure 1.

The cyclist provides energy to the system and the power that he produces is considered to be constant during the trip. Indeed, the SHB architecture makes it possible for the cyclist to pedal in a fixed manner regardless of the road conditions thanks to the motion decoupling between wheel and crankshaft.

In order to cope with the power demand from the rear wheel motor and to keep the cyclist effort constant during all the driving phases, an ESS is needed. This project foresees using supercapacitors to store energy. Indeed, the ESS interest in SHB is not comparable to a classical e-bike battery that needs to be recharged by the power grid. On a SHB, the cyclist will manage his own energy, never recharge the ESS on the grid and will therefore need a lower storage capacity. Supercapacitors also offers a high-efficiency device, extended lifetime and a friendlier environment solution alternative to batteries that make them suitable for the application [START_REF] Raza | Recent advancements in supercapacitor technology[END_REF]- [START_REF] German | Impact of Voltage Resets on Supercapacitors Aging[END_REF].

The SHB most critical aspect is the efficiency. A previous study analysed the minimum energy efficiency criteria (𝜂 𝑆𝐻𝐵 ) and found that 𝜂 𝑆𝐻𝐵 could drop down to 60% and still travel as fast as a classical chain bicycle thanks to energy recovery and the energy management features [START_REF] Tournon | Efficiency comparison between Series Hybrid Bike and Traditional Bike[END_REF].

This paper focuses on determining design rules for SHBs. Indeed, financial cost shall stay as low as possible, performances shall be good enough to guaranty a new pedaling feeling where the power pedaling is smoother. This paper aims to propose optimal sizing solutions trying to minimize the costs and maximize the system's efficiency.

Therefore, this paper will first focus on describing the model used for the simulation, then describe the methodology used for finding optimal solutions and finally describe the various results improvements and conclusions extracted through the optimization.

I. Model description

The following sections will describe the model and the assumptions made to simulate the system. The global model is run using the VEHLIB environment which is a Simulink based tool using energy transfers principles [START_REF] Trigui | Modélisation systémique de véhicules hybrides en vue de la prédiction de leurs performances énergétiques et dynamiques[END_REF], [START_REF] Vinot | Model simulation, validation and case study of the 2004 THS of Toyota Prius[END_REF].

A. Modelling guidelines

The system modelling aims to provide performances estimation of subparts for maximizing the global efficiency and reducing the costs. Therefore, the user, both machines and their associated gearboxes, the DC-DC converter and the ESS will be described using quasistatic models. Costs will also be evaluated.

Costs are often the addition of several expenses as raw material, subsystems and manufacturing costs in addition to the commercial margin. This study expresses financial costs in a simple but realistic manner thanks to market prices analysis.

B. Energy and movement equations of bicycle

In order to track the position of the bike, the speed (𝑉 𝑘 ) is updated at every time step (𝛥𝑡) using variations of the kinetic energy (𝐸 𝑘 ):

𝐸 𝑘+1 = 𝐸 𝑘 + 𝑃 𝑘 𝛥𝑡 (1)
𝑃 𝑘 = 𝑇 𝑚𝑜𝑡 𝑉 𝑘 𝑅 𝑤ℎ𝑒𝑒𝑙 -∑ 𝐹 𝑒𝑥𝑡 𝑉 𝑘 (2)
Those variations are due to a change of power (𝑃 𝑘 ) created by the bicycle (motor torque: 𝑇 𝑚𝑜𝑡 , radius: 𝑅 𝑤ℎ𝑒𝑒𝑙 ), and external forces (𝐹 𝑒𝑥𝑡 ).

External forces are mainly due to the rolling resistance (𝐶 𝑟 ), the drag area (𝑆𝐶 𝑥 ) and the current slope (𝛼 𝑘 ) (equation (3)) [START_REF] Mariot | Mechanics of cycling[END_REF]. The bicycle and the cyclist are modelized as a unique mass (𝑀 𝑡𝑜𝑡 ) with g the standard gravity and 𝜌 the standard air mass.

𝐹 𝑒𝑥𝑡 = 𝑀 𝑡𝑜𝑡 𝑔(C r + sin(α k )) + 𝜌 2 𝑆𝐶 𝑋 𝑉 𝑘 2 (3) 𝑉 𝑘+1 = √𝑉 𝑘 2 + 2𝑃 𝑘 𝑀 𝑡𝑜𝑡 𝛥𝑡 (4)
Once the speed has been found, it is possible to extract both the bicycle position (𝑋 𝑘 ) and the current slope (𝛼 𝑘 ) by integrating the speed and updating the parameter from the ride profile (see section I.F).

𝑋 𝑘+1 = 𝑋 𝑘 + 𝑉 𝑘 𝛥𝑡

(5)

𝛼 𝑘 = 𝑆𝑙𝑜𝑝𝑒𝐿𝑜𝑜𝑘𝑈𝑝(𝑋 𝑘 ) (6)

C. Cyclist

The cyclist is more easily simulated on a SHB than on a regular bicycle. Indeed, the chainless power transmission allows the user to pedal in a regular way independent from the road profile. The user was modeled as a varying sinusoidal torque generator (𝑇 𝑢𝑠𝑒𝑟 ) depending on pedals position (𝜃 𝑢𝑠𝑒𝑟 ) (this is a fair approximation according to [START_REF] Abbiss | Optimal cadence selection during cycling[END_REF]). The pedaling speed was considered to be constant thanks to the torque control algorithm of the generator that simulates the inertia of a classic chain bicycle.

𝑇 𝑢𝑠𝑒𝑟 = 𝑇 𝑚𝑒𝑎𝑛 𝑐𝑜𝑠 (2𝜃 𝑢𝑠𝑒𝑟 - 𝜋 2 ) + 𝑇 𝑚𝑒𝑎𝑛 (7)
The torque constant 𝑇 𝑚𝑒𝑎𝑛 represents the mean torque value as well as the amplitude of the cosinus. Indeed, the user is assumed to produce no torque at both crank vertical dead positions (𝜃 𝑢𝑠𝑒𝑟 = 0 ∧ 𝜋) [START_REF] Abbiss | Optimal cadence selection during cycling[END_REF].

D. Electrical machines and gearboxes

The system has two gearboxes. One of them is in charge to increase the pedaling speed before the generator and the second multiply the motor torque (input torque is called 𝑇 𝑖𝑛 and output torque 𝑇 𝑜𝑢𝑡 ). Gearboxes may have several stages which have to fit the bicycle geometry requirement, it was considered a one stage maximum ratio of 3 (𝑘 𝑠𝑡𝑎𝑔𝑒 ≤ 3). The number of stages (𝑁𝑏 𝑠𝑡𝑎𝑔𝑒 ) has an influence on the overall gearbox ratio (𝑘 𝑔𝑏 ) but also on the costs (𝐶𝑡 𝑔𝑏 ) and on the overall gearbox efficiency (overall gearbox efficiency 𝜂 𝑔𝑏 ; one stage efficiency 𝜂 𝑠𝑡𝑎𝑔𝑒 ). The efficiency of a single stage is considered to be constant over speed and load. Moreover this efficiency doesn't take into account the manufacturing quality. [START_REF] Lundin | Efficiency Analysis of a Planetary Gearbox[END_REF], [START_REF] Henriot | Engrenages : usinage, rendements. Trains planétaires[END_REF].

𝜂 𝑠𝑡𝑎𝑔𝑒 = 0.97 (8)
𝑇 𝑜𝑢𝑡 = 𝑘 𝑔𝑏 𝑇 𝑖𝑛 𝜂 𝑔𝑏 (9)

𝜂 𝑔𝑏 = 𝜂 𝑠𝑡𝑎𝑔𝑒 𝑁𝑏 𝑠𝑡𝑎𝑔𝑒 (10)
𝑁𝑏 𝑠𝑡𝑎𝑔𝑒 = 𝑐𝑒𝑖𝑙 𝑙𝑜𝑔(𝑘 𝑔𝑏 ) 𝑙𝑜𝑔(𝑘 𝑠𝑡𝑎𝑔𝑒 ) (11)
The gearbox costs (𝐶𝑡 𝑔𝑏 ) is modeled only based on the number of stages required to achieve the desired multiplication ratio. The complexity of modeling the cost was introduced in section II.A. In this case, 𝐶𝑡 𝑔𝑏 shall include gears, bearings, manufacturing and assembly costs. In order to simplify the model, a several stages gearbox is assumed to have the same cost at each stage (𝛽 𝑔𝑏 : € 𝑠𝑡𝑎𝑔𝑒 ⁄ ). Cost per stage is determined from typical supplier's price offers (𝛽 𝑔𝑏 ≈ 30 € 𝑠𝑡𝑎𝑔𝑒 ⁄ ). Minimizing the cost is therefore equivalent as minimizing the number of stages.

𝐶𝑡 𝑔𝑏 = 𝛽 𝑔𝑏 𝑁𝑏 𝑠𝑡𝑎𝑔𝑒 [START_REF] Henriot | Engrenages : usinage, rendements. Trains planétaires[END_REF] The machines (Brushless DC motors only) and their associated converters are modeled in a simple manner thanks to a loss table (𝑃 𝑙𝑜𝑠𝑠 𝑚 ) differentiating the input power from the output power (𝑃 𝑖𝑛 𝑚 , 𝑃 𝑜𝑢𝑡 𝑚 ). Losses were assumed to fit a model including conduction losses (resistor 𝑅 𝑚 ), dry friction losses (𝑓 𝑑 ) and viscous losses (𝑓 𝑣 ). The classical torque constant was also introduced (𝐾 𝑇 : 𝑁𝑚. 𝐴 -1 ). This paper deals only with commercially available machines. The losses are then fitted (thanks to equation 14) from the available data in the datasheet. Moreover, the maximum motor torque (𝑇 𝑚𝑎𝑥 ) is also an important parameter influencing the weight, and the gearbox design. [START_REF] Gonzalez | High efficiency Electric Generator for Chainless Bicycle[END_REF].

𝑃 𝑜𝑢𝑡 𝑚 = 𝑃 𝑖𝑛 𝑚 -𝑃 𝑙𝑜𝑠𝑠 𝑚 ( 13
)
𝑃 𝑙𝑜𝑠𝑠 𝑚 = 𝑅 𝑚 𝑇 𝐾 𝑇 2 2 + 𝑓 𝑑 𝜔 + 𝑓 𝑣 𝜔 2 (14)

Table 1: Generator set

The generator will be chosen among the Maxon motors EC-FLAT series and a specific machine design realized by a manufacturer [START_REF]maxon -EC-FLAT[END_REF] (figure 2 and table 1) (𝛼 𝑔𝑒𝑛 is the choice variable: 𝛼 𝑔𝑒𝑛 = 1 for the first generator…). Moreover, a set of motor representative of the market machines diversity (direct drive (𝛼 𝑚𝑜𝑡 = 1), geared motor (𝛼 𝑔𝑒𝑛 = 2), the same special high quality motor design (𝛼 𝑚𝑜𝑡 = 3)) has been selected (see

E. Storage System and related DC-DC converter

The storage system is at the heart of the SHB, this system absorbs the cyclist's power variations and represents a major part of the costs.

ESS's performances vary along the solicitation frequencies. However, the discharge occurs within minutes and transient behaviors were then not considered. This work then focuses only on power losses (𝑃 𝑙𝑜𝑠𝑠 𝐸𝑆𝑆 ) mainly due to the pack internal resistance (𝑅 𝑠 ) and the current flowing through the storage system (𝐼 𝑠 ):

𝑃 𝑙𝑜𝑠𝑠 𝐸𝑆𝑆 = 𝑅 𝑠 𝐼 𝑠 2 (16)
Finally, the state of energy (𝑆𝑂𝐸) of an 𝐸 𝑡𝑜𝑡 energy rated ESS, is an important control parameter [START_REF] Tournon | Efficiency comparison between Series Hybrid Bike and Traditional Bike[END_REF]. It can then be directly derived from the pack open circuit voltage (𝑈 𝑝𝑎𝑐𝑘 ) assuming a constant capacity over 𝑆𝑂𝐸 and a maximum voltage of 2.7𝑉 for each cell: Low energy storage quantity lead to high variations of 𝑆𝑂𝐸 and therefore to large voltage variations. Cut off voltages are implemented to ensure the system is running safely. However, situation where the voltage is low leads to high currents in the system thus increasing the losses (see figure 7). In order to improve the system, an option (𝛼 𝐷𝐶𝐷𝐶 = 1 if DCDC activated, 0 otherwise) is to implement a DCDC converter between the ESS and machines terminals (bus) that will keep the bus voltage constant. Moreover, a low bus voltage will lead to reduced machines performing areas and therefore to deteriorated efficiency and traveling time. Yet, implementing a DCDC converter adds costs to the system and losses which is not necessary if the SOE is high, the voltage being high enough to ensure correct performances. SHBs main advantage is to smooth the developed cyclist power against ride discontinuities (hills, stops). The chosen ride used for the simulation is a virtually created trip as shown figure 3

(a).

A standard 50m hill crossing was chosen as it represents the average hill crossing in town, based on the analysis of 30 French cities elevation maps (figure 4) with Standard Elevation defined as: 𝑆𝐸 = 𝐶𝑖𝑡𝑦𝐶𝑒𝑛𝑡𝑒𝑟 𝑚𝑒𝑎𝑛 𝑒𝑙𝑒𝑣𝑎𝑡𝑖𝑜𝑛 -min 𝑒𝑙𝑒𝑣𝑎𝑡𝑖𝑜𝑛.

G. Control system

The motor control is based on torque control. An optimal motor control law has previously been found by using Dynamic Programming (DP) (see figure 3b) [START_REF] Tournon | Efficiency comparison between Series Hybrid Bike and Traditional Bike[END_REF]. The discrete torque control solution is then interpolated and implemented in the simulation.

DP's parameters have an influence on the optimal control law. All the parameters described in section II.B shall therefore remain the same in both simulations. The following parameters can be adjusted and have influence on the control law:

❖ The system energy efficiency (𝜂 𝑆𝐻𝐵 ) ❖ Power produced by the cyclist (𝑃 𝑢𝑠𝑒𝑟 )

❖ The ESS energy (𝐸 𝑡𝑜𝑡 )

The derived control law (and its associated travel time 𝑡 𝑟𝑖𝑑𝑒 ) presented in figure 3.(b) is the optimal SHB control law associated with a global efficiency of 0.7, a 8𝑊ℎ energy storage pack and a 100W constant human power generation.

H. Performance indicators

SHB performances shall be analyzed. Therefore, the overall energy efficiency (𝜂 𝑆𝐻𝐵 𝑒 ) will be computed as well as the losses for each subsystem (see figure 4) presented above (section II). The energy efficiency is the ratio between the input energy from the user (𝐸 𝑢𝑠𝑒𝑟 = 𝑃 𝑢𝑠𝑒𝑟 𝑡 𝑟𝑖𝑑𝑒 ) and the effective power given to the wheel (𝐸 𝑚𝑜𝑡𝑜𝑟 ). It also takes into account for the energy difference between the initial and the final 𝑆𝑂𝐸(𝛥𝑆𝑂𝐸).

II Methodology

A. Optimization Goal

The main purpose is to increase the system efficiency and therefore aims to reduce system's losses. On the other hand, it is necessary to reduce the costs to achieve an affordable bicycle.

The optimization problem will also be constrained to reduce the possible solutions and keep results within boundaries.

B. Optimization Method

Optimization problems have been widely addressed in the literature. Achieving the best compromises between global system costs and performances requires to solve a multiobjective problem [START_REF] Wahlroos | Optimisation of an Energy System in Finland Using NSGA-II Evolutionary Algorithm[END_REF]- [START_REF] Fonseca | Hybrid fuel cell vehicle energy management including dynamic and durability constraints[END_REF].

The optimization process will find optimal trade-off between two cost functions. Indeed, objectives are contradictory to each other. As an example, increasing the performance will also increase the costs.

Let 𝑓 1 𝑎𝑛𝑑 𝑓 2 be the functions to optimize:

𝑓 1 = 𝑚𝑖𝑛 (∑𝐶𝑡) (24) 𝑓 2 = 𝑚𝑎𝑥 (𝜂 𝑆𝐻𝐵 ) (25)
The following variables will be optimized:

The simulation also has to be constrained. Indeed, faults may occur during the trip. For example, the motor power demand can empty the ESS as a result of a to small embedded capacitance. Therefore, it is also necessary to ensure the trip as been complete at the end of the simulation.

𝑋 𝑒𝑛𝑑 ≥ 𝑚𝑎𝑥(𝐷𝑖𝑠𝑡𝑎𝑛𝑐𝑒) (29)
The NSGA II (Non dominated Sorting Genetic Algorithm) is a multiobjective optimization tool that as already be widely covered in the literature and will be used here to solve the problem described above. [START_REF] Deb | A Fast Elitist Non-dominated Sorting Genetic Algorithm for Multi-objective Optimization: NSGA-II[END_REF]

III Results

A. Reference result

The reference simulation, not yet optimized, is based on the topology used on an existing SHB prototype (generator nominal power 800𝑊, generator rated speed 1000𝑟𝑝𝑚, 𝑇 𝑚𝑎𝑥 = 7𝑁𝑚, 3 stages gearbox, 8 supercapacitors of 1200F in series). Every subsystem was designed separately from each other according to the state of the art (torque considerations, embedded energy needed...). This design will be challenged thanks to the optimization. Finally, the motor used here is a direct drive meaning there is no gearbox between the motor and the wheel (and therefore no gearbox losses).

Figure 5 shows a typical losses distribution for the standard ride defined in section I.F for the reference design bicycle configuration. The global associated bicycle energy efficiency is 𝜂 𝑆𝐻𝐵 𝑒 = 0.36. This level of performances would cause a bad user experience. Improving the efficiency up to 0.6 is then mandatory to achieve a successful bicycle [START_REF] Tournon | Efficiency comparison between Series Hybrid Bike and Traditional Bike[END_REF].

Most losses (see figure 5) occur in the machines, indeed close to 70 % of the energy is lost in the generator and the motor. Machines are then key components to optimize. Losses in the ESS represents around 15% of the losses. However, ESS's losses are strongly varying depending on its solicitation and SOE. (see figure 7).

For the hilly trip, the supercapacitors are delivering the needed power in addition to the generated power in order to follow the demand power as it is shown in figure 6. The trip profile lasts 1287 seconds. The power is oscillating due to sinusoidal user power production and represents an important losses increase (up to 13% more losses). The optimization was launched based on the NSGA II algorithm. The optimization process was involving all the parameters defined in table 3.

Figure 8 shows the results. As expected, the costs are kept between boundaries and it is possible to observe cost/efficiency variations. The transmission system's total costs (€) is plotted against the best energy efficiency achievable.

First, one motor out of the 3 versions came out of the simulation (the special design), both other don't performs well enough to reach the required performances and therefore aren't part of the optimal solution. This also point the major importance of motors for achieving a performing system. Moreover, cost are quite high especially due to the high quantity of embedded energy needed. Finally, the optimized design has better performances than the reference design as energy efficiency progressed from 0.36 up to 0.62.

The results are meeting the efficiency requirement defined in the study even if only reached by a few percents [START_REF] Tournon | Efficiency comparison between Series Hybrid Bike and Traditional Bike[END_REF] and therefore validate the possibility and the structure to build in order to create an unplugged SHB.

Nevertheless, costs are quite high for this trip and the efficiency only shortly reach the requirement. The next section will investigate the better efficiency one might reach with the best-case scenario.

C. Optimization results: the flat case

Another case was studied in order to compare the result presented in section III.B. The alternative case is a standard onekilometer flat drive including no stops and no hill. The system run on steady state point on this type of trip which enables to find an optimal maximum efficiency operating point. This configuration shows the highest efficiency level the system can reach with the available components.

The typical results are shown in figure 9. The results for this optimization goes beyond the objective. However, it is necessary to recall this efficiency is an optimum and could not be reached on a non-steady state drive cycle (on the hilly trip with those components for example).

The results confirm the major importance of the motors in the performance evaluation. A better motor automatically provides a better solution (for this drive) independent from the generator choice, supercapacitors…

D. Comparison

Both cases described in previous sections B and C are complementary. Indeed, the efficiency close to 75% confirm the possibility to create an efficient enough system in order to reach more than 60% energy efficiency. This is even possible while maintaining reasonable costs as the needed embedded energy is quite low on a flat land.

However, the steep hill crossing presented in this paper is showing some limits to those results. Indeed, the need for a bigger energy storage quantity as well as the higher power flowing through the motor combined with an eventually limited speed and therefore a limited efficiency is causing the system to have worse performances and higher cost than on the flat drive.

One challenge to be faced is to propose cost effective motors that can achieve high level performances for the application.

Conclusion

This works is a major step to the unplugged SHB technical validation. Furthermore, it has shown the limit and potential of the system pointing out it is possible to reach an efficient level system that will propose an interesting performance to cost ratio but that this achievement is strongly dependent on the drive cycle.

This study has not only validated the possibility to build a SHB reaching the energy efficiency basic conditions but it also proposes an entire range of architecture possibilities depending on the costs and performances requirement. The comparison between several cases also gives us more information about the robustness of the design choice and the capability of the system to outstand various solicitations (hills, start/stop…).

The optimization shows it is possible to reach more than 60% energy efficiency on a specific non continuous ride (the hilly ride) but it has also shown the bicycle structure could achieve up to 75% percent efficiency with optimal drive situations (flat profile). Finally, the study has also increased the previous prototype efficiency by 26% and pointed out the advantages of a global system optimization over separate standalone design. However, efficiency requirements were reached using special high-quality designs.

This study also only focuses on the structural bicycle, one major step will be to integrate the user and its physiological efficiency into the calculations. The optimization goal wouldn't be to maximize efficiency but to control user's effort and minimize user's fatigue that will be reduced thanks to steady state pedaling rate.

Moreover, DP and NSGA II algorithms are run independently and could be run together in order to get more precise and accurate results. In addition, the dynamic programming as a command can also be conflicting with the reality as it is an optimal which will never occur on real conditions drives. One further improvement could be to implement real time energy management system based on MPC (Model Predictive Control) or ECMS (Equivalent Consumption Minimization Strategy) for example.

Finally, the model and performances shall be tested on real prototype for modeling and result validations.

Figure 1 :

 1 Figure 1: Series Hybrid Bike architecture and associated efficiencies

Figure 2 :

 2 Figure 2 : Maxon Motor EC flat 160W typical efficiency map

18 )

 18 supercapacitor pack costs (𝐶𝑡 𝐸𝑆𝑆 ) are depending on the number of cells (𝑁 𝑐𝑒𝑙𝑙 ), on the supercapacitor capacitance (𝐶 𝑆𝐶 ) and on a price coefficient (𝛽 𝐸𝑆𝑆 ) expressed in €/F: 𝐶𝑡 𝐸𝑆𝑆 = 𝛽 𝐸𝑆𝑆 𝑁 𝑐𝑒𝑙𝑙 𝐶 𝑆𝐶 (

  Figure 3: (a) Trip profile and (b) motor torque command

Figure 4 :

 4 Figure 4: Nominal elevation values A second profile is also considered in this paper as a comparison case. The second trip is a one-kilometer flat profile.

  gives a direct indication about the average performance of the bicycle.

Figure 5 :

 5 Figure 5: Losses percentage overview of the reference design (non-optimized)

Figure 6 :

 6 Figure 6: power exchanges in the system for the hilly trip Finally, SOE as well as the ESS losses are plotted in the figure 7.

Figure 7 :

 7 Figure 7: SOE and ESS'losses for the hilly trip B. Optimization results: the hill case

Figure 8 :

 8 Figure 8: optimization result

Figure 9 :

 9 Figure 9: optimization result for a flat land drive

table 2

 2

		EC Flat 160 EC Flat 220	Special
				design
	Torque (𝑁𝑚)	0.457	0.729	3
	Efficiency (%)	85	85	93
	Cost (€)	121	128	700
	).

Table 1 :

 1 Motor Set

Table 2 :

 2 Parameter limits

		Description	Interval Unit
	𝑘 𝑔𝑏 𝑔𝑒𝑛 , 𝑘 𝑔𝑏 𝑚𝑜𝑡	Gearbox ratio	[1; 50]
	𝑁 𝐸𝑆𝑆	Cells number	[1; 22]
	𝐶 𝑆𝐶	Cells capacity	[60; 3500]	𝐹
	𝛼 𝐷𝐶𝐷𝐶	DC-DC	[0; 1]
		activated
	𝛼 𝑔𝑒𝑛	Generator choice [1; 2; 3]
	𝛼 𝑚𝑜𝑡	Motor choice	[1; 2; 3]
	Moreover, some constrains have been added to the system:
		400 < 𝐶𝑡 𝑆𝐻𝐵 ≤ 2000€	(26)
		𝑘 𝑔𝑏 𝑔𝑒𝑛 𝑇 𝑚𝑎𝑥 𝑔𝑒𝑛 ≥ 2𝑇 𝑚𝑒𝑎𝑛		(27)
		𝑘 𝑔𝑏 𝑚𝑜𝑡 𝑇 𝑚𝑎𝑥

𝑚𝑜𝑡

≥ 𝑚𝑎𝑥 (𝑇 𝑐𝑚𝑑 )