Pierre Romet
email: pierre.romet@utbm.fr

Romain Tabusse
email: romain.tabusse@utbm.fr

Franck Gechter
email: franck.gechter@utbm.fr

El-Hassane Aglzim
email: el-hassane.aglzim@u-bourgogne.fr

Samir Jemei
email: samir.jemei@univ-fcomte.fr

David Bouquain
email: david.bouquain@univ-fcomte.fr

Daniela Chrenko
email: daniela.chrenko@utbm.fr

Autonomous Electric Vehicle Routing Problem using Ant Colony Optimization with consideration of the Battery State-of-Health

Keywords: EVRP, ACO, autonomous vehicle, lithium-ion battery, SOE, SOH, DOD

HAL is

I. INTRODUCTION

For several years, a significant effort has been made to optimize the transport of goods in urban and peri-urban area. This improvement is mainly aimed at the reduction of the secondary effects such as road system congestion, air pollution and sound pollution. Thus, recent improvements to rationalize freight transport are mainly legislative such as the limitation of the size and weight of vehicles, the arrangement of vehicle schedules, urban distribution centers and last mile delivery. Furthermore, soft distribution circuits have emerged and new technologies as the electric car and autonomous vehicles are currently studied to resolve last mile delivery. In this context, This work has been supported by the French State through the SURATRAM project, the Eco-Campus project, the EIPHI Graduate School and the Region Bourgogne Franche-Comté. Special thanks to the Régie des Transports du Territoire de Belfort for sharing their bus fleet data.

vehicle routing problems (VRP) must be adapted to meet the specific requirements of these new technologies.

For this purpose, [START_REF] Caric | A survey on the electric vehicle routing problem: Variants and solution approaches[END_REF] presents a comprehensive overview about Electric Vehicle Routing Problem (EVRP). It begins with an introduction of EVRP issues linked to different technologies of battery. Then, the paper presents a non-exhaustive research work about EVRP possible variants, such as Energy Shortest Path Problem [START_REF] Polycarpou | Ant colony optimization for the electric vehicle routing problem[END_REF], Electric Traveling Salesman Problem, Heterogeneous or Mixed Vehicle Fleet [START_REF] Chrenko | Electric vehicles in logistics and transportation: A survey on emerging environmental, strategic,and operational challenges[END_REF], Hybrid Vehicles [START_REF] Caric | A survey on the electric vehicle routing problem: Variants and solution approaches[END_REF], Partial Recharging, Charging technologies and strategies [START_REF] Chrenko | Electric vehicles in logistics and transportation: A survey on emerging environmental, strategic,and operational challenges[END_REF], Battery Swapping, Two-Echelon Routing Problem, Charging Schedule and Dynamic Traffic Conditions.

However, the electric vehicle, due to the chemical nature of its battery, cannot be considered and used in the same way as a thermal one. Indeed, its recharge time is longer. Moreover, its performance will decrease over time in much larger proportions, especially its autonomy and power. To the authors' knowledge, this last point has never been considered in the EVRP.

In this paper, the development of a new freight transport system in urban area is proposed. It is based on autonomous electric vehicles that can move on city public transport system (dedicated bus lane). It will allow to transport goods from urban distribution center to all delivery points (bus stations) located in urban area [START_REF] Senouci | Transportation of goods in inner-city centers: Can autonomous vehicles in platoon be a suitable solution?[END_REF]. To address the issue of battery aging, the influence of battery depth-of-discharge (DOD) is investigated. For this study, real data from a bus company and a battery aging database are used.

The paper is structured as follow. Section II introduces all concepts and methods used. The section III describes experimental procedure. Then, Section IV presents the experimental results and a discussion about them. Finally, Section V summarizes the findings and states the main conclusions.

II. MATERIALS AND METHODS

This study investigates the last mile delivery of parcels within the city of Belfort using autonomous electric vehicles. For this, the city bus network is used as a support. The vehicle considered is a Tesla Model S in autonomous mode.

The objective of the study is to optimize the delivery of parcels by defining an adapted delivery strategy and an adequate sizing of the vehicle fleet. For this, the energy parameter will be particularly optimized. Indeed, for an autonomous electric vehicle, the main cost is related to the purchase of the vehicle and its energy consumption. It is also linked to the aging of the battery which will define the life span and therefore the replacement (or repurchase) of the battery.

Therefore, this work will complement the previous work on the EVRP, by taking into account the aging of the battery. For this purpose, an aging model based on the DOD, defining the vehicle's charge/charge strategies, will be used. It will allow to dimension the battery and then to calculate the financial cost of taking into account the aging. In addition, after the sizing of one vehicle, it will be possible to carry out the sizing of an entire fleet.

In the following, the main concepts for better understanding, algorithms and models used are presented. Then, the conducted experiments are described in details.

A. EVRP

This paper aims to optimize a last mile delivery logistics plan in city center, using autonomous electric vehicles. In order to determine a set of least energy cost delivery routes from an urban distribution centre (UDC) to delivery points, a Vehicle Routing Problem (VRP) algorithm adapted to electric vehicles (EVRP) is implemented. Our EVRP algorithm is based on a graph, where each edge contains the distance, the time and the energy information needed to travel through it.

The problem is solved under the following constraints:

• Each customer is visited only once by a single vehicle;

• Each vehicle must start and end its route at the urban depot; • Total volume of packages, transport of the vehicle cannot exceed its volume capacity; • Total energy consumed cannot exceed the energy capacity of the vehicle.

B. ACO

As the EVRP II-A is based on a graph, it was decided to implement an ant colony algorithm to exploit it and determine a tour that would minimize total energy consumption. For this, the ACO model implemented is the same as in [START_REF] El | Towards an environment approach to determine the ant colony termination criterion[END_REF], to which a termination criterion has been added. Thus, thanks to the maximum number of iterations calculated by the termination criterion, the optimal solution is found. The number of iterations to be performed is given by the equation 1.

θ = n v=1 ω i k (1)
With v the number of vertices in the graph, ω i the number of edges starting to the current vertex and connecting it to another vertex and k the number of edges in the graph.

C. Graph Building

The graph is a replica of the public transportation system of the city of Belfort (France). Each node is a bus stop and edges are bus lanes on the road system of the city. Moreover, each edge is composed of calculating parameters: distance, time, packages volume to be dropped off and energy consumed to move between two bus stops.

These parameters were calculated thanks to data of the Régie des Transports du Territoire de Belfort (RTTB), concerning the traffic of buses on the city's public transport system. Over ten past years data, only one day was used in this work in a first approach. Their data are files storing at thirty-second intervals, successive GPS points forming the route taken by the bus during the day (KMZ file).

To calculate the parameters of each edge, all GPS points forming the path between two bus stops were extracted. Then time, distance and energy were calculated.

The distance of each edge is calculated using the WGS84 ellipsoid formula, based on the set of GPS points between two bus stops. In addition, the energy required to travel the edge is calculated with the energy model, described in section II-D.

D. Energy Model

The energy of each edge of the graph is determined using the fundamental principle of dynamics applied to the vehicle illustrated in figure 2. Knowing the characteristics of the vehicle and its speed profile on the edge covered, the mechanical power at the wheels P wheels necessary for its movement at any time can be determined by equation 2.

(2)

P wheels = v veh (M v dv veh dt + 1 2 ρ air v 2 veh SC x + M v gsin(α) + M v gC r cos(α))
where: v veh is the vehicle's speed in m.s -1 , α is the road grade in rad, g is the acceleration of gravity in m.s -2 , ρ air is the air density at sea level at 15°C, M v is the vehicle's mass in kg, Cx is the drag coefficient of the vehicle, S is the frontal area of the vehicle in m 2 , C r is the rolling resistance coefficient. The vehicle's speed v veh is determined using the RTTB's GPS database by deriving the position of the bus as a function of time. Given the low sampling rate of the GPS data used (one sample every 30 sec), only the average speed of the vehicle is considered over an edge (zero acceleration and zero energy recovery).

Similarly, the slope calculation is done over an edge from the difference in elevation between the first point of the edge and the last one by equation 3.

α edge = arctan( ∆Elevation Distance ) (3)
Thus the equation 2 simplifies and gives the average power required over an edge. After multiplying by the travel time ∆t, we find the energy required to move the vehicle on the edge. Finally, considering a constant powertrain (PWT) efficiency η P W T , the energy supplied by the vehicle's battery to cover the edge is given by equation 4.

(4) E edge = η P W T * P wheels * ∆t To calculate the corresponding energy for each edge, equation 4 is used with the Tesla Model S 90D (2017) vehicle's parameters shown in table I [START_REF]tesla model s 90d -specifications[END_REF]. The aging of a battery is commonly defined as an irreversible loss of its performance. The most commonly used indicator is the battery capacity expressed in Ah. The state of health (SOH) of a battery is defined as the ratio of its remaining capacity to its initial capacity measured when new (equation 5). It is most often expressed as a percentage. A new battery corresponds to an SOH of 100%, while a battery at end-of-life usually corresponds to an SOH of about 80% in electric cars [START_REF] Dubarry | Perspective on State-of-Health Determination in Lithium-Ion Batteries[END_REF].

SOH = currentCapacity initialCapacity (5)
The SOH of a battery can be determined by direct measurement or by modeling. Concerning the last, a large number of approaches exist [START_REF] Xiong | Towards a smarter battery management system: A critical review on battery state of health monitoring methods[END_REF]. Here an empirical data based approach is considered.

For this purpose, battery aging data from an open access database have been used [START_REF]Data for degradation of commercial lithium-ion cells as a function of chemistry and cycling conditions[END_REF]. The cells studied here are the Panasonic B in 18650 format. They are used in the Tesla Model S vehicle presented previously. Considering a linear battery aging from 100% SOH to 80%, a linear model has been established from the aging database. Thus, the life cycle of the cell considered is given for three different DOD in table II: 20%, 60% and 100%. The cycle life is given in number of cycles performed at a given DOD until reaching a SOH of 80%. Therefore, the SOH of the battery for each DOD considered is proportional to the energy that passed through it.

F. Delivery List Generator

In order to be consistent with a real last mile delivery application, it is necessary to determine the number of parcels dropped off per delivery point, their size and volume.

First of all, [START_REF] Reme-Harnay | L'emprise des sous-traitants de la messagerie vis à vis de leurs donneurs d'ordres : entre dépendance économique et reconnaissance professionnelle[END_REF], [START_REF] Ducret | Livraison de colis et logistique urbaine: quelles recompositions de la messagerie en milieu urbain? Revue franc ¸aise de gestion industrielle[END_REF], [START_REF] Erhel | La messagerie, une approche systémique pour appréhender les réseaux de transport de colis[END_REF] make it possible to identify that a set of 1 to 6 parcels is deposited on each delivery point of the tour, according to their weight and their volumes. In addition, the company Amazon provides a comprehensive list of all the packages it processes by weight and volume III. Based on the number of packages per delivery point [START_REF] Reme-Harnay | L'emprise des sous-traitants de la messagerie vis à vis de leurs donneurs d'ordres : entre dépendance économique et reconnaissance professionnelle[END_REF], [START_REF] Ducret | Livraison de colis et logistique urbaine: quelles recompositions de la messagerie en milieu urbain? Revue franc ¸aise de gestion industrielle[END_REF], [START_REF] Erhel | La messagerie, une approche systémique pour appréhender les réseaux de transport de colis[END_REF] and Amazon's package format III. The generated list of packages to be delivered is defined as being arbitrarily made up of 50% standard letter, 35% standard package and 15% oversized package.

III. RESULTS GENERATION

A. From EVRP to EAHVRP

As introduced in the section II-A, an EVRP coupled with the use of an autonomous vehicle II is used on our graph II-C to optimise the tour to deliver all the packages.

The optimization phase is performed by the ant colony algorithm II-B, which at each iteration will weight the parameters distance, time, volume, energy of each edge constituting the tour generated by an ant. Each parameter is weighted differently. The total energy needed to move on an edge is weighted by a factor of 0.5. The total volume of parcels to be delivered to the next delivery point is weighted by a factor of 0.3. The time and distance to travel an edge is weighted by a factor of 0.1.

Working on last mile delivery, three lists of goods were generated according to the model defined in II-F, in order to study three typical problem sizes, as illustrated by [START_REF] Reme-Harnay | L'emprise des sous-traitants de la messagerie vis à vis de leurs donneurs d'ordres : entre dépendance économique et reconnaissance professionnelle[END_REF].

The three graphs are composed of 30, 60 and 80 delivery points randomly and homogeneously distributed on the graph.

In addition, to meet the criteria of a last mile application, the Urban Distribution Centre (UDC) was placed in the centre of the city. More precisely at the railway station, this allows the UDC to be supplied by rail so as not to impact the city's road system. Thus, after exploiting the three graphs with our AEVRP algorithm, a tour minimizing the energy needed to visit each of the delivery points was generated. Moreover, as introduced in II-E, optimising the energy needed for a tour minimizes the impact on battery aging. Finally, the phenomenon of aging of the vehicle battery is studied through the influence of different DOD strategies III-C.

Together these elements allow us to define our new variant of the EVRP, the Electric Autonomous Heatlth Vehicle Routing Problem (EAHVRP).

B. DOD strategy

Nowadays, different technologies for charging batteries are available. A slow charge made with ¡7kW chargers allows to recharge the battery between 6 and 8 hours. A fast charge, carried out in 1-2 hours, uses chargers with a power ranging from 11 to 43 kW. Finally, a super charger with a power of 50 to 250 kW can recharge Tesla's battery in a maximum of 45 minutes [START_REF] Chrenko | Electric vehicles in logistics and transportation: A survey on emerging environmental, strategic,and operational challenges[END_REF].

However, it must be taken into account that each type of recharge (slow, fast, accelerated) has a different impact on the aging of the battery [START_REF] Al | Influence de la recharge rapide sur les performances des accumulateurs lithium des véhicules électriques dans le cadre de l'utilisation postale[END_REF]. The capacity of a battery degrades at 80% of its initial value more quickly with a fast or accelerated recharge, so these contribute greatly to the acceleration of the aging processus.

Moreover, the section III-C explains the approach to sizing a battery to maximize its life. A vehicle is therefore designed to make the best use of a given capacity range, depending on the DOD strategy, to complete its tour, which varies between 5am and 18pm. Once the tour is over, the vehicle has the whole night to recharge, i.e. between 18 p.m. and 5 p.m., allowing the use of a slow charge. In this context, an electric vehicle will perform one charge/discharge cycle per day.

Furthermore, according to [START_REF]Data for degradation of commercial lithium-ion cells as a function of chemistry and cycling conditions[END_REF], three different DOD strategy has been defined to study the impact of the battery aging on the battery. The state of health of the battery defining the number of cycle and the time of use of the vehicle, the aging process impact as well as the number of packages deliver and the financial impact of the vehicle over its lifetime. The following DOD strategies are consider: 100%DOD (from 100% to 0% SOC), 60%DOD (from 80% to 20% SOC), 20%DOD (from 60% to 40% SOC).

Having defined the DOD strategies to be used in our study, it remains to size the battery capacity according to them.

C. Battery Energy Capacity Sizing

Being able to determine the energy required to complete a tour, it remains to calculate the battery energy capacity E battery according to the DOD strategy used. For this, the equation 6 is use. ( 6)

E battery = E needed /DOD(%)/20%
where the 20% coefficient corresponds to the loss of capacity of the battery at the end of its life.

IV. RESULTS

A. EAHVRP Results and Solution Set

The presentation of the research results is organised as follows. First, the solutions generated by the EAHVRP algorithm will be introduced (III-A), and then an analysis of these solutions combined with the different DOD strategies studied will be provided.

For each size of graph presented in the part II-C, the EAHVRP algorithm generated a set of possible solutions to achieve the delivery of goods IV.

Before continuing, it is necessary to consider a counterintuitive fact concerning the number of rounds. Thanks to the first line of the table IV of the set of solutions of the 80 points graph, we observe that a low number of rounds is not necessarily linked to a low energy consumption. Thus, the solution having consumed the least energy will have made one more cycle than the others to carry out its tour. The same procedure is followed for the other two data sets in order to select the best solution. These solutions can be found in the table V.

B. Battery Sizing Results

From the set of solutions selected in Table V, it is possible to calculate the minimum energy capacities required for an electric vehicle to complete its tour, depending on the DOD strategy used (II-E). Results are shown in Table VI. Thus, it is possible to define 9 vehicles operating in 9 different scenarios on three graph sizes. They provide a set that best characterizes the battery size of a vehicle. However, four vehicles are set aside for our study. Those with an energy capacity below 30 kWh and those with an energy capacity above 100 kWh.

Indeed, the energy model II-D is designed to calculate the energy necessary to move a 2199 kg vehicle with a battery of 90 kWh. Knowing that decreasing the energy capacity of the battery is equivalent to decreasing its mass, our model is no longer valid for solutions with a consumption of less than 30 kWh and conversely for consumptions greater than 100 kWh. Considering an energy density of 250 Wh/kg, reducing a battery by 1kW means reducing its mass by 4kg. As mass is a significant part of the energy model equation, it is not possible to exploit these results without distorting our study.

From this, it must be understood that the extremes are amplified. Thus, for batteries with an energy capacity of less than 30 kWh, their real energy capacities will be revised downwards and for batteries with an energy capacity of more than 100 kWh, their real energy capacities will be revised upwards.

Finally, concerning the analysis of the impact of aging on the battery, and thanks to the exploitation of open access data [START_REF]Data for degradation of commercial lithium-ion cells as a function of chemistry and cycling conditions[END_REF], it was possible to calculate the number of cycles achievable given the degradation factor, due to the DOD strategy used VII. Thanks to the table VII, we can observe that, depending on the DOD strategy chosen, it has a major impact on the life of the battery. We can see that, for an equivalent energy capacity, between a DOD strategy of 100% and a strategy of 20%, the life span in terms of number of cycles (III-B) is multiplied by 18. Thus, due to the DOD strategy implemented III-B, the battery's lifetime in number of years varies between 1 year and over 18 years.

C. Financial Cost

Having been able to determine the impact of aging on the life of a battery, we now seek to determine its impact on the cost of purchase. Considering the price of 1 kWh battery of 100$, it is possible to put a price on each of the previously calculated energy capacities VI. In the table VIII, it can be seen that the price of a battery using a DOD strategy of 100% is 4994 $. For a battery using a 20% DOD strategy, its price is 8484$. Furthermore, based on the table VII, we can see that using a DOD strategy of 20% makes it possible to minimize the aging phenomenon of a battery and thus multiply its life by 18. Thus, it appears that implementing a DOD strategy of 20% certainly implies multiplying the price of a battery by 1.7, but the lifetime of the latter is multiplied by 18.

We are now interested in the impact of aging on the financial cost per delivery point. It can be concluded that, although the "oversizing" linked to the DOD strategy of 20% multiplies the price of the battery by 1.7, the financial impact of a delivery point is divided by 4 (0.174/0.043=4).

D. Fleet sizing

Being able to size a vehicle, it is possible to continue this work to size a fleet of electric vehicles. The table X provides information on the energy cost of a graph according to the number of delivery points. As can be seen in the table X, the graph size that requires the least amount of energy to explore is 60 delivery points.

This result may seem counter-intuitive, but under-utilisation of the transport system (30 delivery points) implies long distances between each point, increasing consumption. Furthermore, over-utilisation (80 delivery points) leads to more detours between each point, thus increasing consumption.

This observation highlights the need for a detailed study on the sizing of a fleet of autonomous electric vehicles with heterogeneous energy capacity.

V. CONCLUSION

This article has studied the consideration of the problem of the aging of the battery of an electric vehicle through its use and its dimensioning.

Exploiting here a VRP problem with autonomous electric vehicles, we propose an Autonomous Electric Health VRP algorithm to address this problem. Having generated a set of solutions for three sizes of graphs II-C, traversed by vehicles exploiting different DOD strategies VI. We sought to define the impact of aging on the battery.

Through the results IV, we can conclude that the aging of lithium batteries is a key point for a good use of this technology in the transport field. Indeed, we have been able to show the importance of aging through its impact on the life of the vehicle IV-B, as well as on the financial aspect IV-C.

Thus, it can be concluded that the choice of the DOD strategy is the major element in the aging of a battery VII. In addition, although the use of a low percentage DOD strategy implies a larger battery size. Its lifetime is, however, multiplied by 18. It can therefore be concluded that "over-consumption" of raw material is in this case preferable for sustainable management of natural resources. It should also be noted that since the battery degrades slowly, by minimizing the aging factor, this favours the second life of the battery.

Finally, taking account of aging is not a financial problem, as the factor for increasing battery life is more than 10 times higher than the factor for increasing battery price.

Finally, by analysing the energy cost per delivery point X, we notice that the energy required depends on the size of the studied graph. This element opens the way to the continuation of our work, concerning the dimensioning of a fleet of autonomous electric vehicles of heterogeneous energy capacity.

Fig. 1 .

 1 Fig. 1. Public transportation system of Belfort city

Fig. 2 .

 2 Fig. 2. Balance of forces applied to the vehicle. The inertia of the rotating masses is not taken into account.

TABLE II BATTERY

 II CYCLE LIFE AS A FUNCTION OF DOD. END-OF-LIFE IS REACH BY 80% SOH.

	DOD(%)	100	60	20
	Cycle life 363	1220 6595

TABLE VI BATTERY

 VI ENERGY CAPACITY ACCORDING TO DOD STRATEGIES

	Delivery points 100% DOD 60% DOD	20% DOD
	30	16.96 kWh	28.28 kWh	84.84 kWh
	60	30.50 kWh	50.84 kWh	152.50 kWh
	80	49.93 kWh	83.23 kWh	249.67 kWh


Xinyue Wu

Yabin Jin

Abdelkrim Khelif

Xiaoying Zhuang
email: zhuang@iop.uni-hannover.de

Timon Rabczuk

Bahram Djafari-Rouhani

Topological surface wave metamaterials for robust vibration attenuation and energy harvesting

Keywords: Surface wave metamaterial, topological insulator, vibration attenuation, energy harvesting, robustness

Metamaterials are extensively utilized to manipulate ground surface waves for vibration isolation within the bandgap frequency ranges whereas topological crystals allow the creation of robust edge states immune to scattering by defects. In this work, we propose a topological surface wave metamaterial working in the Hertz frequency range, constituted of triangle shape concrete pillars arranged in a honeycomb lattice and deposited on the soil ground. Based on the analogue of quantum valley Hall effect, a non-trivial bandgap is formed from the degeneracy lifting of the Dirac cone at the K point of the Brillouin zone by breaking the inversion symmetry of the two pillars in each unit cell. A topological interface is created between two different crystal phases and a topological edge state based on surface acoustic wave propagation is demonstrated. The robustness of the topologically protected edge state is quantitatively analyzed in presence of various defects and disorders. Finally, we take advantage of the robust and compact topological edge state for designing a harvesting energy device. The results demonstrate the functionality of the proposed structure for both robust 2 surface vibration reduction as well as energy harvesting by designing proper topological waveguides.

Introduction

Mechanical vibrations in urban areas resulting from engineering construction and ground transportation may damage buildings and bridges, impact the operation of sophisticated instruments in high-tech laboratories and intrude residents' slumber. Apart from the airborne acoustic wave, the energy of these vibrations is generally concentrated in the depth of one wavelength underneath the ground surface, that is, primarily propagating in the form of surface waves, and the frequency range is usually below 20.0 Hz [START_REF] Miniaci | Large scale mechanical metamaterials as seismic shields[END_REF][START_REF] Cai | Propagation and attenuation of Rayleigh and pseudo surface waves in viscoelastic metamaterials[END_REF].

Mechanical metamaterials [START_REF] Jin | Physics of surface vibrational resonances: Pillared phononic crystals, metamaterials, and metasurfaces[END_REF] whose configurations can be artificially fabricated and modulated have been proven to manipulate the propagation of mechanical waves with new effects such as wave isolation [START_REF] Khelif | Locally resonant surface acoustic wave band gaps in a two-dimensional phononic crystal of pillars on a surface[END_REF][START_REF] Liu | Locally Resonant Sonic Materials[END_REF][START_REF] Jin | Tunable waveguide and cavity in a phononic crystal plate by controlling whispering-gallery modes in hollow pillars[END_REF][START_REF] Oudich | Rayleigh Waves in Phononic Crystal Made of Multilayered Pillars: Confined Modes, Fano Resonances, and Acoustically Induced Transparency[END_REF][START_REF] Chen | Isolation of Surface Wave-Induced Vibration Using Periodically Modulated Piles[END_REF][START_REF] Achaoui | Seismic waves damping with arrays of inertial resonators[END_REF], topological insulator [START_REF] Cai | Symmetry-enforced three-dimensional Dirac phononic crystals[END_REF][START_REF] Chaunsali | Subwavelength and directional control of flexural waves in zonefolding induced topological plates[END_REF][START_REF] Ma | Topological phases in acoustic and mechanical systems[END_REF][START_REF] Nguyen | Tunable topological bandgaps and frequencies in a pre-stressed soft phononic crystal[END_REF][START_REF] Chaunsali | Experimental demonstration of topological waveguiding in elastic plates with local resonators[END_REF][START_REF] Fan | Elastic Higher-Order Topological Insulator with Topologically Protected Corner States[END_REF][START_REF] Jin | Topological states in twisted pillared phononic plates[END_REF][START_REF] Chen | Self-ordering induces multiple topological transitions for in-plane bulk waves in solid phononic crystals[END_REF][START_REF] Huang | Topologically protected zero refraction of elastic waves in pseudospin-Hall phononic crystals[END_REF][START_REF] Huang | Subwavelength elastic topological negative refraction in ternary locally resonant phononic crystals[END_REF], Fano resonance [START_REF] Wang | Robust Fano resonance in a topological mechanical beam[END_REF][START_REF] Goffaux | Evidence of Fano-Like Interference Phenomena in Locally Resonant Materials[END_REF][START_REF] Jin | Acoustic analogue of electromagnetically induced transparency and Autler-Townes splitting in pillared metasurfaces[END_REF][START_REF] Jin | Tunable Fano resonances of Lamb modes in a pillared metasurface[END_REF], chirality [START_REF] Wen | Robust edge states of subwavelength chiral phononic plates[END_REF], focusing and imaging [START_REF] Jin | Elastic Metasurfaces for Deep and Robust Subwavelength Focusing and Imaging[END_REF], among others. Especially, surface wave metamaterials [START_REF] Miniaci | Large scale mechanical metamaterials as seismic shields[END_REF][START_REF] Brûlé | Emergence of seismic metamaterials: Current state and future perspectives[END_REF][START_REF] Colombi | A seismic metamaterial: The resonant metawedge[END_REF][START_REF] Colombi | Forests as a natural seismic metamaterial: Rayleigh wave bandgaps induced by local resonances[END_REF][START_REF] Li | Seismic metamaterial barriers for ground vibration mitigation in railways considering the train-track-soil dynamic interactions[END_REF][START_REF] Mu | A Review of Research on Seismic Metamaterials[END_REF][START_REF] Muhammad | Built-up structural steel sections as seismic metamaterials for surface wave attenuation with low frequency wide bandgap in layered soil medium[END_REF][START_REF] Pu | Seismic metasurfaces on porous layered media: Surface resonators and fluid-solid interaction effects on the propagation of Rayleigh waves[END_REF][START_REF] Wootton | An asymptotic hyperbolic-elliptic model for flexuralseismic metasurfaces[END_REF][START_REF] Zeng | Low-frequency broadband seismic metamaterial using I-shaped pillars in a half-space[END_REF][START_REF] Liu | Ultra-wide band gap metasurfaces for controlling seismic surface waves[END_REF][START_REF] Zeng | Subwavelength seismic metamaterial with an ultra-low frequency bandgap[END_REF] are designed mainly for transmission suppression of low-frequency surface wave in certain frequency ranges based on the conception of hybridized bandgaps. Nonetheless, energy harvesting and signal detection of the subwavelength ground surface wave can be additional functions for surface wave metamaterials, remaining a meaningful and challenging issue that has hitherto received few attentions. It is possible to transform the surface wave into available electric energy and detect ground dynamics for potential catastrophe early warning.

For practical applications, defects and disturbances are ubiquitous in surface wave metamaterials, hence the surface wave control with high robustness is a crucial challenge. Recently, the rapid development of metamaterial topological insulator provides an efficient way to realize robust surface wave metamaterials. Topology provides a method to describe global wave properties over a band structure which is able to be conserved under certain local perturbations. In recent years, the development of topological insulators extended from condensed matter physics to classic wave systems such as acoustic and elastic waves. In general, Dirac cone dispersion results from the lattice symmetry which is protected by the space inversion and the timereversal symmetry. Different mechanisms like the quantum Hall, quantum spin Hall and valley Hall effects are proposed to open a bandgap from the Dirac cone [START_REF] Ma | Topological phases in acoustic and mechanical systems[END_REF].

According to the elastic analog of the quantum valley Hall effect, topological phases associated to opened bandgaps from the degeneracy lifting of the Dirac cone can be achieved via the inversion symmetry breaking [START_REF] Torrent | Elastic analog of graphene: Dirac cones and edge states for flexural waves in thin plates[END_REF]. When combining two bulk metamaterials with the opposite topological phases, topologically protected edge states can be induced showing strong energy localization at the interface within certain frequency ranges. The effect has been employed successfully in Lamb wave manipulation with significant robustness [START_REF] Torrent | Elastic analog of graphene: Dirac cones and edge states for flexural waves in thin plates[END_REF][START_REF] Chaplain | Topological Rainbow Trapping for Elastic Energy Harvesting in Graded Su-Schrieffer-Heeger Systems[END_REF][START_REF] Chen | Topological valley transport of platemode waves in a homogenous thin plate with periodic stubbed surface[END_REF][START_REF] Jin | Robustness of conventional and topologically protected edge states in phononic crystal plates[END_REF][START_REF] Pal | Edge waves in plates with resonators: an elastic analogue of the quantum valley Hall effect[END_REF][START_REF] Vila | Observation of topological valley modes in an elastic hexagonal lattice[END_REF][START_REF] Yan | On-chip valley topological materials for elastic wave manipulation[END_REF], revealing a great potential in the ground surface waves manipulation. For the application to surface acoustic waves, Wang et al. reported robust guiding valley-dependent edge states for Rayleigh waves excited by chiral sources at around 30 MHz [START_REF] Wang | Guiding robust valleydependent edge states by surface acoustic waves[END_REF]; Topological chiral edge state is also realized on a periodically corrugated surface [START_REF] Inoue | Topological band structure of surface acoustic waves on a periodically corrugated surface[END_REF].

For ground surface waves in macro scale, industrial and traffic activities induced vibrations are different from chiral excitations. In this work, we study the topological protected edge modes based on surface acoustic waves within the non-trivial bandgap of a pillared metamaterial. The edge states are excited by a force point source with either out-of-plane or in plane polarization. Indeed, these are the main polarizations of the typical surface waves such as Rayleigh or Love waves. In Sec. 2, we describe the models and the implemented methods for designing non-trivial bandgaps and topological edge modes of the pillared metamaterial based on surface acoustic waves.

In Sec. 3, we quantitatively study the robustness of the designed topological edge mode against three kinds of perturbations. In Sec. 4, we adopt the robust and compact topological edge mode whose path can be designed, according to the need, to study energy harvesting properties by using piezoelectric materials. Finally, a summary of the main results is presented in Sec. 5.

Models and Methods

We propose a design constituted by a honeycomb lattice of equilateral triangular section pillars made of concrete on the surface of soil. We are interested by ground surface waves with frequency range under 20.0 Hz, so we chose a lattice constant of a = 2.0 m and the parameters of the pair of identical pillars in the unit cell are fixed as edge length d = 0.8 m and height h = 2.0 m. The design and the unit cell for the 3D simulation model are shown in FIG. 1(a). To prevent the physical surface modes being perturbed by the unphysical Lamb type waves of the thick substrate plate in the simulation model, the depth of soil in the unit cell is chosen to be H = 30.0 m and a numerical low-reflecting boundary is applied at the bottom surface to avoid wave reflection. The soil (Young's modulus Es = 10 MPa, Possion's ratio νs = 0.3, mass density ρs = 1800 kg/m 3 ) and concrete pillars (Young's modulus Ec = 40 GPa, Possion's ratio νc = 0.25, mass density ρc = 2500 kg/m 3 ) are assumed to be homogeneous, isotropic and linearly elastic materials [START_REF] Miniaci | Large scale mechanical metamaterials as seismic shields[END_REF]. Therefore, the equation of motions in solids can be written with harmonic time dependence assumption as

ω 2 ρu i = -σ ij, j (1)
σ ij = c ijkl u k, l (2)
where ω is the angular frequency, u i the displacement component, σ ij the stress tensor and c ijkl the stiffness tensors that can be expressed as a function of the Young's modulus E and Poisson ratio v. To obtain the dispersion relations, Bloch periodic conditions are applied on the four vertical surfaces of the unit cells based on Bloch theorem [START_REF] Hofstadter | Energy levels and wave functions of Bloch electrons in rational and irrational magnetic fields[END_REF] as

u(r + a) = u(r)e ik•a ( 3
)
where a is the lattice vector, r the position vector, and k = [ kx , ky ] the wave vector. In addition, stress-free conditions are applied on the other top boundaries of the simulation model. We sweep the wave vector along the high-symmetry axes Γ -K -M -Γ in the first irreducible Brillouin zone of the honeycomb lattice (see FIG. To demonstrate the robustness of the topological edge modes, we design a zigzag interface separating the two bulk metamaterials 1 and 2 with the L-L and S-S interfaces respectively. Since typical surface waves such as Rayleigh and Love waves have different polarizations, we study the surface wave-based topological edge modes by using both out-of-plane and in-plane force excitations.

A point source represented by an out-of-plane force (marked as the green star in for the two structures. The mode exhibits a high confinement near the surface, with most of the energy localized within one wavelength depth from the ground surface. This behavior supports the conclusion that the designed topological edge states are based on nonleaky surface waves that do not radiate into the bulk. In addition, the energy density at the dot 2 near the exit is a bit smaller than that at the dot 1 near the entrance, remaining at relatively the same level, which also supports high transmission and backscattering immunity of the topological edge mode.

We now perform similar calculations for an excitation by an in-plane point source force with the same amplitude as above. The results are presented in FIGs.

3(e)-(h). The generated surface wave displays a higher weight of the in-plane displacement over total displacement than the case of the out-of-plane force excitation. One can see in FIGs.

3(e)-(f) similar high transmission and backscattering immunity of the topological edge state as above. However, the kinetic energy density (FIGs. 3(g)-(h)) is now smaller than before.

From the dispersion diagram shown in FIG. 2, one can note that the edge mode of the S-S type interface covers a certain frequency range, which makes it a good candidate for energy harvesting. In the following section, we focus on the edge mode of the S-S type excited by the out-of-plane force source for robustness analysis and energy harvesting. We investigate the robustness of the edge state to propagate along the waveguide, by calculating the loss in the transmission coefficient in presence of a perturbation.

Robustness of the topological edge mode

Therefore, we define the ratio Δ1/Δ2, where Δ1 and Δ2 represent the integrated displacements along the green arcs at the exit of the waveguide with and without perturbations, respectively. We consider three kinds of perturbations to quantitatively study the robustness of the topological edge state.

Topological edge mode for energy harvesting

The robust and compact topological edge mode has a great potential for energy harvesting from a practical point of view. We adopt piezoelectric PZT-5H patches and attach them to the top surface of soil. With the electro-mechanical coupling effect, the vibrating energy of soil can be transformed to electric power. We choose five identical square PZT patches with sides l = 0.1 m and thicknesses t = 0.001 m, and arrange them alternately as shown by the five yellow squares in FIG. 5(b). A simple circuit powered by each PZT patch with a resistance R = 1000 Ω is devised, whose average output power Pi is calculated by

P i = V i √ 2 I i √ 2 = RI i 2 2 ( 5
)
where Vi is the output voltage, Ii is the loop current to be measured. The detailed information can be found in Appendix B. We sum the total power 4)) in FIG. 5(a). From FIG. 5(a), the topological edge mode with bright color is successfully captured within the frequency range [10.33, 10.97] Hz and the peak is measured at 10.40 Hz which coincides to the design. The electric power at the peak is about two orders of magnitude higher than that within the frequency ra [9.50, 10.00] Hz outside the non-trivial bandgap. Besides, it is more than one order of magnitude higher than that of a bare surface, showing a big advantage of the energy harvesting by topological edge mode. We also present the displacement field at 10.40 Hz in FIG. 5(b), almost identical to the field without patches (see FIG. 3(b)).

In fact, all the vibrating pillars at the interface are available for the energy conversion in practice. Given that the topological edge mode is highly robust to the studied three kinds of perturbations, the corresponding induced energy harvesting also has an efficient robustness.

At the macro scale, the basic functionality of the structure is to robustly isolate surface waves with the non-trivial low-frequency bandgap. Energy harvesting is an additional functionality by designing the proper interface path inside the metamaterial.

With maximum vibrating displacements at the interface of the order of 0.5mm, and the small size of the designed PZT patch in FIG. 5, an electric power of about 10 -5 W is achieved. However, if the surface vibrations increase in practice, the PZT patch would be able to generate much higher electric power, making the metamaterial as a potential ground electric generator. FIG. 5

Conclusion

A topological surface wave metamaterial composed of triangular section concrete pillars arranged in a honeycomb lattice on the soil is proposed for robust vibration attenuation and energy harvesting at low frequency ranges, below 20 Hz. By breaking the inversion symmetry of the unit cell while preserving C3v symmetry, the degeneracy of the Dirac cone at K point is lifted to open a non-trivial bandgap. By combining two metamaterials with different topological phases, two types of topological edge modes (L-L and S-S types) with backscattering suppression and compact properties can be designed around 10 Hz. We further quantitatively analyzed the robustness of the S-S type topological edge mode against three kinds of perturbations, namely, (i) stripe defects in the vicinity of the interface towards the bulk of the crystal, (ii) unit defects at the interface and (iii) disorders in the height of pillars at the interface. For perturbation (i), the transmission of the topological edge mode keeps stable at relative high level no matter how large is the defect in the bulk media, showing a highly compact property of the topological edge mode. For perturbation (ii), the transmission first decreases very fast before becoming stable when the unit-defect number reaches 4. For perturbation (iii), the transmission slowly decreases in quasi-linear shape. For the latter two cases, the topological edge modes are still preserved at high considered perturbation degree, exhibiting a good robustness. Finally, piezoelectric patches are designed and attached on top of the pillars at the interface. The robust and compact topological edge mode makes it possible for energy harvesting within the frequency range [10.33, 10.97] Hz of the topological edge mode. The entire size of the pillared metamaterial can be further minimized according to the potential practical demands. For instance, a thickness such as 4 or 5 rows of units is sufficient to effectively attenuate surface waves. The unit cells' geometric parameters can also be optimized for different working frequency range. The employed concrete and PZT are commercial materials with low price, making the metamaterial realizable in practical engineering. The proposed new surface wave metamaterial has great potential in robust surface vibration attenuation with the nontrivial low frequency bandgap and robust energy harvesting with the topological edge mode.

FIG. 2

 2 FIGs. 1(e)-(f). The elastic energy flux travels anticlockwise and clockwise around the pillars for D1 and D2 valleys, respectively. If d1 = 1.0 m and d2 = 0.7 m, the chirality of the valley states for the lower and higher bands associated to the opened bandgap becomes inverted. Therefore, the valley Chern number associated with this lower and upper bands have opposite signs[START_REF] Jin | Topological states in twisted pillared phononic plates[END_REF][START_REF] Pal | Edge waves in plates with resonators: an elastic analogue of the quantum valley Hall effect[END_REF], which further supports an elastic analog of quantum valley Hall effect. As a consequence, the opened bandgap becomes non-trivial.The wavelength of referenced Rayleigh wave propagating along the free surface without pillars at the center frequency of the non-trivial bandgap is about 4 m, also demonstrating a subwavelength property of the pillared metamaterial. The induced nontrivial bandgap is able to provide robust surface wave attenuation.

FIGs. 3

 3 FIGs. 3(a)-(b)) is excited at the entrance of the zig-zag interface. Low-reflecting boundaries are applied to the surrounding boundaries to eliminate the reflection effects. The edge modes are excited at 10.75 Hz and 10.40 Hz and their propagation presented for the L-L (FIG. 3(a)) and S-S (FIG. 3(b)) types interfaces, respectively. In both structures, the wave entering the waveguide travels towards the end and passes the sharp bending corners without a significant backscattering. To further demonstrate that the edge mode is essentially based on surface acoustic waves, we plot the kinetic energy density at the green dots 1 and 2 as a function of the depth in the soil in FIGs. 3(c)-(d)

FIG. 4

 4 FIG. 4 The robustness of the edge state at the S-S type interface in presence of defects and disorder in height of the pillars at 10.40 Hz. The upper panels show the out-of-plane displacement along the waveguide for a particular choice of disorder whereas the lower panels give the attenuation in the transmission for various degrees of disorder. (a) & (d) The effect of removing m strips on each side of the interface in the middle part of the waveguide. The blue strip supercell constituting the interface area contains 6 unit cells. The out-of-plane displacement is detected at the black arcs at the exit of the zigzag path. (b) & (e) The effect of introducing n defects in the unit cells contained in the blue area close to the interface. (c) & (f) The effect of height disorder by a factor η for the 24 units in the blue area. The light blue domain in (f) demonstrates the standard deviation of the transmission for different realizations of the disorder at each η. In the lower panels, the calculations are performed for discrete values of the abscissa and the continuous lines are a guide for the eyes.

  The first perturbation consists of removing m strips on each side of the interface in the middle part of the waveguide as shown in FIG.4(a). The area delimited by two blue strips on each side of the waveguide is constituted by 6 unit cells. When one strip of defect is introduced (m = 1), FIG.4(d) shows that the transmission decreases suddenly to about 75%; then, by increasing further m, it slightly oscillates around the level (70% -75%), exhibiting a good robustness to stronger disorder. The displacement field for m = 7 is calculated and presented in FIG.4(a), showing a preserved zigzag edge state. Such high robustness is mainly attributed to the strong localization of the topological edge mode in the vicinity of the interface. The second perturbation (FIGs. 4(b) and (e)) consists of randomly removing n units at the left side of the middle part of the zigzag waveguide. From FIG. 4(e) one can observe that when the defect number n increases from 0 to 4, the transmission decreases by a large amount to about 20% of the unperturbed situation, due to the appearance of backscattering effect; then the level of transmission keeps at about 20% when the defect number continues to increase until 7. The transmission field at n = 7 is shown in FIG. 4(b) which indicates that although the major part of the wave is reflected, there still remain weak monopolar edge state after the 7-unit-defect along the last horizontal part of the zigzag path, showing a good robustness in edge mode shape.The third perturbation is dealing with the effect of disorder in the height of a set of 24 pillars located in the blue shaded region of FIG.4(c). These heights are randomly distributed within the interval [h -Δh, h + Δh], thus defining a disorder degree η = Δh/h. For a given disorder degree, the average transmission coefficient is evaluated over 30 random samplings and the results are shown as the solid line in FIG.4(f). We also evaluate the standard deviation of transmissions and plot it as the shaded blue domain around the solid line in FIG.4(f). The stable level of the standard deviation supports that the random sampling number of 30 should be sufficient. Generally, when the disorder degree η increases from 0 to 60%, the average transmission decreases quasilinearly, exhibiting that the backscattering effect also becomes stronger. We show the transmission field at η = 60% in FIG.4(c). Despite the backscattering effect, the relative weak propagating wave along the last horizontal part of the zigzag path still conserves the monopolar edge state property, showing a strong robustness of the topological edge mode.

  all the six PZT patches, and present the transmission spectra together with the color information of localization ratio β (Eq. (

  FIG. 5 (a) Piezoelectric energy harvesting for the edge modes with (b) the displacement fields at 10.40 Hz. The red dashed curve in (a) represents the electric power measured without pillars. The 5

Acknowledgment

This work is supported by the National Natural Science Foundation of China (11902223), the Shanghai Pujiang Program (19PJ1410100), the program for professor of special appointment (Eastern Scholar) at Shanghai Institutions of Higher Learning, the Fundamental Research Funds for the Central Universities, high-level of foreign expert program of Tongji University and Shanghai municipal peak discipline program (2019010106). This work is also partially supported by the French EIPHI Graduate School (contract "ANR-17-EURE-0002").

Data Availability

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

Appendix A. The strip's dispersion calculation

The strip's dispersion calculation is modeled via the finite element method. In FIG.

A1, a stripe consisting of lattices red and blue with an interface at the center is deposited on top of the substrate. The valley Chern number of the two bands associated with the non-trivial bandgap are opposite for the red and blue lattice bulk media. Low reflection boundary is applied to the bottom surface. Periodic conditions are applied to the twoblue side-faces of the substrate along x direction. Stress free conditions are applied to the two-grey side-faces of the substrate along y direction. The built model contains about 350 thousand domain element meshes with about 2.9 million degrees of freedom.

Appendix B. Simulating information for the energy harvesting model

In this section, we detail the simulation process of the transmission with piezoelectric energy harvesting along the zig-zag interface in the COMSOL