R Vazquez Casique

M L Viollet

M Bauerheim

V Chapin

E Gowree

FP27-AERO2022-viollet

Several studies have revealed that boxfishes, tropical fishes living in shallow waters, show exceptional aerodynamic performances despite their large volume and crosssectional area. These characteristics make the boxfish shape an interesting topic of study for drag optimisation of bluff bodies. This study focuses on the aerodynamic analysis of several mathematical shape representations of the boxfish. The effect of the edges is studied since they are responsible for the generation of vortices, which play a role in potentially manipulating the wake, resulting in an overall decrease in the drag coefficient. The effect of the sharpness or roundness of the edges on the aerodynamic performance is investigated. Aerodynamic coefficients are obtained numerically and experimentally for a range of Reynolds numbers between 3000 and 300000 and at pitch angle from -16 • to +20 • .

INTRODUCTION

In nature, animals and plants have adapted over years of evolution to their environment in order to carry out efficiently the tasks they require for survival. Bio-mimetics and bio-inspired optimisation consist of taking nature as an example of optimal performance for diverse applications. Researchers are trying to transpose some of the remarkable characteristics of natural systems into humanbuilt systems to achieve the same efficiency observed in nature. Such an example is observed in the Boxfish (Teleostei: Ostraciidae), a tropical species of fish living in shallow waters around coral reefs, recognisable by their rigid, bony, and often geometrically shaped carapaces. A very common boxfish, Ostracion cubicus, is displayed in Fig. 1.

Studies have shown that they are very stable and are able to swim on remarkably straight paths through the Figure 1: Ostracion cubicus [START_REF]Yellow boxfish -toxic and inspiring?[END_REF] highly perturbed flow conditions present in coral reefs [START_REF] Webb | Boxfishes (teleostei: ostracii-dae) as a model system for fishes swimming with many fins: kinematics[END_REF]. Other studies have also shown that they can still exhibit high manoeuvrability despite their rigid carapaces, unlike other species of fishes [START_REF] Walker | Does a rigid body limit maneuverability[END_REF]. Another important characteristic of the boxfish is its low drag coefficient despite being a bluff body with an aspect ratio close to 1/3. The boxfish's large volume and cross-sectional area coupled with their low drag coefficient would allow transport vehicles to transport payload at low drag, reducing fuel cost but more importantly emissions.

The dynamics of a swimming boxfish was studied by Gordon et al. [START_REF] Gordon | Boxfishes as unusually well-controlled autonomous underwater vehicles[END_REF]. In the study, the idea that the boxfish is handling dynamically the flow by using its fins but, more importantly, its body, to generate vortices and almost perfectly balance by continuously adapting its hydrodynamic forces and moments was presented. This work led to further studies on the effect of the body on the stability of the Boxfish. The role of the carapace of the boxfish on stability and maneuverability was investigated by Bartol et. al. [START_REF] Gharib | Flow patterns aroundthe carapaces of rigid-bodied, multi-propulsor boxfishes (teleostei: Ostraciidae)[END_REF][START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] through an analysis of the pitch and yaw control. Force balance and Digital Particle Image Velocimetry (DPIV) measurements were used to characterise the performance of the Boxfish. Results showed that the vortices generated near the anterior edges of the ventro-lateral keels grow in strength as they move to the posterior edge of the body. These regions of high vorticity generate a suction effect in the posterior part of the body, which create restoring moments for pitching motions. The vortices and the selfcorrecting moment have been shown to be directly proportional to the pitch angle of the boxfish. It has been shown that multiple parts of the carapace play important roles in the generation of vorticity [START_REF] Gharib | Flow patterns aroundthe carapaces of rigid-bodied, multi-propulsor boxfishes (teleostei: Ostraciidae)[END_REF]. The work was extended to three other species of Boxfish with varying cross-sectional shapes in 2004 [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF], for which the results show a similar tendency; the carapace of the Boxfish produces self-trimming and restoring moments through the creation and amplification of vortices in the keels. Studies, carried out with actively swimming Boxfish, have allowed to separate the contribution of the body and the fins, as well as their interactions [START_REF] Weihs | Evidence of self-correcting spiral flows in swimming boxfishes[END_REF]. The self-correcting flows observed around static carapace models were verified for actively swimming boxfishes, with a higher vorticity and irregular shapes because of the interaction between fin movement and the carapace. The formation of the vortical flows below the ventral keels was validated, and it was found that the self-correcting keel-induced leading edge vortices do not necessarily act in isolation, and can be modified with powered input from the fins [START_REF] Weihs | Evidence of self-correcting spiral flows in swimming boxfishes[END_REF]. The use of pectoral fin motions allow the Boxfish to improve stability or manoeuvrability as desired. For stability improvement, pectoral fin motions supplement the passive self-correcting control mechanism by powering it with more flow as required. To enhance manoeuvrability, the pectoral fin generates circulatory forces capable of producing vorticity of similar magnitude to that produced by the keels.

The effect of the vortices induced by the carapace of the Boxfish, defined as 'self-stabilising' by Bartol [START_REF] Weihs | Evidence of self-correcting spiral flows in swimming boxfishes[END_REF], are stated to be too small to have a significant effect by Van Wassenbergh [START_REF] Marcroft | Boxfish swimming paradox resolved: forces by the flow of water around the body promote manoeuvrability[END_REF]. It is argued that the stabilising moments created by the vortices, posterior of the centre of mass, are negligible compared with the moment caused by the frontal pressure wave on the head during pitch and yaw. This hypothesis is supported by arguing that having a destabilising flow over the boxfish's body is adequate with its swimming style and environment since it promoted manoeuvrability. Numerical results revealed that during pitching, 76% of the stabilising moments are cancelled by destabilising moments. For yawing motion, all stabilising moments are cancelled.

All previous studies found in the literature suggest that the Boxfish shape is interesting to study as it presents enhanced aerodynamic performance and, although contradicting, high manoeuvrability and stability. However, the possible low drag characteristic which potentially enables its relatively non-negligible speed, of similar range as automobiles in terms of body length displacement per second, has not been thoroughly investigated. The work presented aims to characterise the flow around a boxfish shape to understand better the flow characteristics which enable the Boxfish to swim efficiently in its environment.

The study focuses on the effect of the Reynolds number on the aerodynamic forces, for Reynolds number ranging from 3000 to 300000, in order to characterise the flow regimes and for comparison with more canonical bluff bodies. To achieve this wide range of Reynolds number both numerical and experimental studies are performed. The effect of pitching is also examined. Finally, a sensitivity of the aerodynamic forces to small but strategic geometrical characteristics of the Boxfish is carried out focusing on sharp edges on the carapace, capable of triggering in-phase vortical structures.

MATERIAL AND METHOD

Model construction

The boxfish shapes studied here are simplified mathematical geometries that have been slightly modified from a real scanned carapace of an Ostracion Cubicus. This species of boxfish is characterised by its semi-rectangular shape with four edges. The models were made using a scan fitting routine, representing the geometry by a finite collection of points in space, and then fitting a NURBS surface to the points. This allows the models to keep some concavities specific to the real boxfish shape. With this method an aspect ratio of about 1/3 is obtained. The effect of the sharp edges on the aerodynamic performance of the boxfish shape is studied by changing the sharpness or roundness of the edges. The geometries, presented in Fig. 2, are used for numerical simulations and are modelled to have no sharp edges, two sharp edges, and all four sharp edges, denominated BV0, BV2, and BV4, respectively. An experimental campaign was also carried out for the BV4 shape. The length of the boxfish is the reference length used to compute the Reynolds number. For both numerical and experimental analysis, the angle of incidence, α, refers to the angle between the free-stream flow and the centreline of the fish derived from the CAD drawing. On Fig. 2, the boxfish shapes are presented at α = 0 • . Figure 2: Simplified boxfish geometries: BV0 (left), BV2 (centre), and BV4 (right). Red edges signify edges which are modelled as sharp.

Numerical Analysis

A Lattice-Boltzmann Methods (LBM) solver is used to carry out Direct Numerical Simulations (DNS) at several Reynolds numbers. Based on the length of the boxfish shape, Reynolds number from 3000 to 105000 have been investigated. Simulation settings are set to an upstream velocity of 1 m/s, a characteristic length of 1 m, and air density of 1 kg/m 3 . The Reynolds number is varied by changing the value of the kinematic viscosity. Simulations are carried out at a Reynolds number of 3000 for BV0, BV2 and BV4 for positive α values of 0 • , 5 • , 10 • , 15 • and 20 • . Simulations on BV4 with no pitch have also been carried out for Reynolds numbers of 10000, 20000, 35000, 50000, 65000, and 105000.

The mesh is based on an octree structure, which allows cutting the domain into several cubic regions with different levels of refinement. The geometries and the wake are positioned in the highest level of refinement, to capture with accuracy the relevant flow characteristics. The rectangular domain has dimensions of 37.5 m in the x-direction, 30 m in the y-direction, and 30 m in the zdirection. A mesh independence study was carried out to determine the best compromise between computational cost and accuracy. To this effect, first it was carried out for a canonical flow around a sphere at a Reynolds number of 3000, based on the sphere diameter, for which validation data exist. Three different meshes, from M1 to M3, with different levels of refinement were simulated, and the results are shown in Table 1. Based on these results, M2 was chosen as it is able to reliably predict the expected drag coefficient at this Reynolds number and provide enough refinement at the lowest level to discretise with enough points the geometry without increasing the computational cost of the simulations. An unsteady analysis of the forces obtained by LBM show that for M2, the correct Strouhal numbers were captured in agreement with the results obtained by Yun et al [START_REF] Gowree | Vortical structures behind a sphere at subcritical reynolds numbers[END_REF]. A Strouhal number of 1, related to the vortex shedding of the sphere, and the correct mean drag coefficient of approximately 0.4 were obtained.

Simulations with the 3 meshes were carried out on BV4 to validate the mesh choice. The results, presented in Table 2, confirm that the choice of M2 for simulation on the BV shapes will allow a good trade-off between precision and cost of simulation.

As the flow behind the boxfish is expected to be unsteady, a time independence study was also carried out. Three different time steps were investigated, dt = dx • 10 -1 , dt = 2dx • 10 -1 and dt = 3dx • 10 -1 . This time independent study was performed using M1 to reduce the computational cost. The results, presented in Table 3, show very close results between the drag coefficient for all time steps. It was decided to keep T1 as it will allow to save computational time and it has shown that it is capable of capturing accurately both the steady and the unsteady flow over the sphere.

dt mean (C D ) T1 1.46 • 10 -3
0.4126 T2 2.92 • 10 -3 0.4077 (+1.1%) T3 4.39 • 10 -3 0.4132 (+0.1%) Table 3: Time independence study on BV4

Wind tunnel measurements

Force measurements have been conducted at ISAE-Supaero's subsonic wind tunnel on a 3D printed BV4 model in order to obtain the drag and lift. The printed BV4 shape was post-treated manually to remove unwanted surface imperfections to avoid tripping the boundary layer. The model has a body length of 160 mm, resulting in a reference cross-sectional area of A re f = 0.00337 m 2 . The wind tunnel has a 65 cm long test section with a cross-sectional area of 45 x 45 cm, resulting in a blockage ratio of 1.66% at 0 • pitch angle. The wind tunnel had stable flow range for velocities between 6 m/s (= 37 Body Lengths/s) and 28 m/s (= 175 Body Lengths/s), therefore the Reynolds number of 65000 and 300000 could be achieved.

The model is mounted on top of an in-house force balance thanks to three streamlined and tapered supports, shown in Fig. 3. The assembly allows to pitch the model from -16 • to +20 • from the geometric angle. The calibration of the balance was done at the beginning of the test campaign with weights in the range of the expected forces produced by the model in the direction of lift and drag. The contribution of the 3 stands on the measured forces are not negligible, therefore offwind and on-wind measurements have been undertaken to subtract their contribution on the force measurements during post processing. The force coefficients are then computed as such:

C l = (2 • L)/(ρ • V 2 • A re f ) and C d = (2•D)/(ρ •V 2 •A re f )
, with L and D the lift and drag forces acting on the boxfish, respectively.

Figure 3: Experimental wind tunnel setup using 3D printed and smoothed BV4 shape.

For each measurement, the output data, such as velocity, drag, and lift are recorded for 10 seconds at a frequency of 10Hz, giving 100 measurements. In the post processing the mean values are presented, with the error bar corresponding to the standard deviation taken during one realisation. The standard deviation was lower than 2% for most of the aerodynamic coefficient measurements.

RESULTS

Influence of the Reynolds number on the aerodynamic forces

Figure 4 presents numerical and experimental results on the evolution of the drag coefficient as a function of the Reynolds number for the BV4 shape at a 0 • pitch angle.

It is observed that the drag coefficient is highly dependent on the Reynolds number and that a decrease in the drag coefficient is observed with an increasing Reynolds number. We observe a mismatch in the drag coefficient estimation between the numerical and experimental results at the Reynolds numbers of 65000 and 105000. Previous investigations on the influence of the Reynolds number on the drag coefficient of a boxfish were performed by Van Wassenbergh et al. [START_REF] Marcroft | Boxfish swimming paradox resolved: forces by the flow of water around the body promote manoeuvrability[END_REF] on an Ostracion Cubicus shape. The results were obtained by measurements in a water tank for Reynolds numbers between 3400 and 68000. It is important to point out that the results obtained by Van Wassenbergh et al. [START_REF] Marcroft | Boxfish swimming paradox resolved: forces by the flow of water around the body promote manoeuvrability[END_REF] at low Reynolds numbers present high errors and are thus not completely indicative of the drag coefficient of the shape at this range. The results obtained after a Reynolds number of 30000 show a more reliable accuracy and are in a reasonable agreement with respect to simulations carried out. Figure 5 shows that the experimental results obtained in the wind tunnel follow a similar trend as the one observed by Wassenbergh et al. [START_REF] Marcroft | Boxfish swimming paradox resolved: forces by the flow of water around the body promote manoeuvrability[END_REF]. At the low Reynolds region, around 65000, both numerical and experimental estimations of the drag coefficient are in close range with the results obtained by Wassenbergh et al. [START_REF] Marcroft | Boxfish swimming paradox resolved: forces by the flow of water around the body promote manoeuvrability[END_REF] and Bartol et al. [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF]. A comparison of the average results at a Reynolds number of approximately 65000 can be made in Table 4. The differences found in the results obtained at the various Reynolds numbers compared to those found in literature can in part be attributed to the slight geometric differences of the boxfish shapes studied by Wassenbergh et al. [START_REF] Marcroft | Boxfish swimming paradox resolved: forces by the flow of water around the body promote manoeuvrability[END_REF] and Bartol et al. [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] with respect to BV4. The tendency, however, shows that the boxfish shape does not present a sudden drag crisis as other bluff bodies, like the sphere or the ellipsoid represented in Fig. 5. Instead, it shows a slower decrease in the drag coefficient with increase of the Reynolds number. The presented drag coefficient is for an ellipsoid with an aspect ratio of 1:1.8 [START_REF] Prandtl | Ergebnisse der Aerodynamischen Versuchsanstalt zu Göttingen[END_REF]. As expected, due to the more streamlined geometry of an ellipsoid at an aspect ratio of 1:1.8, the drag coefficient of the ellipsoid is lower than the drag of the boxfish. This puts the boxfish shape in a region between the drag coefficient of the sphere and of an ellipsoid.

To see if a different regime can be triggered an experiment to trip the boundary layer of the boxfish has been performed. A trip is added on the 4 faces of the body at a chord-wise length x/c = 0.38 shown in Fig. 6, in order to force the transition of the boundary layer to delay the Ostracion 65 000 0.282 [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] meleagris Table 4: Drag coefficient at Re≈65000 for multiple boxfish models.

separation. This usually leads to a decrease in the drag coefficient as the flow stays attached to the body and the major source of drag is due to the separated flow behind bluff bodies. The results of the experiment on BV4 show indeed that when the boundary layer is tripped the drag coefficient is lowered, as observed in Fig. 5. Further flow analysis is needed to better understand this behaviour.

Influence of the angle of incidence on the aerodynamic forces

The lift and drag coefficients are observed to be highly dependent on the angle of attack, as seen in Fig. 7. The lift coefficient is directly proportional to the angle of incidence and shows a similar tendency to that of lifting surfaces. It can be observed that the slope of the curves do not change greatly with varying Reynolds number. At the geometric angle of attack equal to zero, the lift coefficient is approximately zero. Results by Bartol et. al.

[3] at Reynolds number of 65000 show a similar trend. It is observed that the slope from Bartol et. al. is slightly lower than the experimental results of this study. This is due to the Reynolds number of used by Bartol et. al. No noticeable stall is observed even at very high pitch angle. This is due to the contribution from the vortices generated by the edges of the boxfish, similar to the leading-edge separation-induced vortex lift generated on delta wings. These vortices were previously reported by Bartol et. al. and were also captured by our current simulation presented in Fig. 12. The rather erratic behaviour in the lift at Re = 65000 is an artefact related to the sensitivity of the balance in detecting and capturing very small loads which were of the order of a couple of grams in this case. The drag coefficient is highly dependant on the pitch angle. The results obtained through the experimental campaign present a parabolic tendency centered around a geometric incidence angle equal to 0, which is coherent with Bartol et. al. [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] results. The plateau observed between the pitch angles of -10 • to 10 • is respected, and the drag coefficient grows again with pitching angles outside this range. The sharp rise in drag is potentially due to the increasing contribution of the form drag but also the induced drag component resulting from the non-negligible amount of lift generated with increasing angle of incidence. A difference in the angle at which the minimum drag coefficient is observed exists between the Ostracion meleagris shape [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] and the BV4 shape. It is at -4 • for the former while it is situated at +3 • for the latter. This difference can be explained by the difference in geometry between the two models (the shape, the definition of the geometric angle and the roughness) and possibly the set-up, where the mounting struts can undergo interaction with the flow. Results by [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] presented for comparison.

Effect of the body edges

After building-up some confidence from the numerical simulation capability shown in 5 we proceeded with the analysis of the effect of the edges from the boxfish's body. Numerical simulations were carried out for BV0, BV2, and BV4 at multiple pitch angles at Re = 3x10e3 as the cost of simulation could increase drastically considering the total number of cases to be simulated. The evolution of the aerodynamic coefficients with respect to angle of incidence for each shape is presented on Fig. 8. It can be observed that all shapes follow the same trend as seen for higher Reynolds number, i.e. the drag coefficient follows a near parabolic trend and the lift coefficient a linear trend with no stall behaviour was observed. It is interesting to note that the slopes of the lift coefficient evolution are different depending on the shapes, which can explain the difference in slope observed with the results from Bartol et. al. [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] as a real fish was used. In general the tendency in the lift showed strong dependence on the effect of edge modification. For a configuration which mixes both sharp and smooth edges, BV2, the lift coefficient is higher at high angles of attack compared to configurations where all edges are modelled either as smooth or sharp, BV0 and BV4 respectively. At low angles of attack BV0 and BV4 show similar increment in lift, potentially due to the symmetry in terms of edge modifications. However, the difference between these two configurations starts to increase at α > 15 • . The major difference here in terms of shape is the smoothness of the edges and with sharper edges stronger vortices will be favoured on the BV4 contributing to additional lift. The larger difference from BV2 therefore emanates from a combination of stronger suction in the presence of the edge vortices and also a larger pressure difference promoted by asymmetry in edges.

From the evolution of the drag coefficient as a function of angle of incidence shown in Fig. 8a a similar trend is depicted as that in Fig. 7a for the BV4 measured at significantly higher Reynolds number. Similar to the lift, the drag increased with increasing angle of incidence for all three cases, and once again this is due to a rapid increase in form drag and induced drag due to the rapid increase in lift. However, unlike the lift at very low angle of incidence both modified shapes, BV2 and BV4 showed closer agreement and were slightly higher than the unmodified BV0 case. The presence of sharp edges results in increases in the drag coefficient which can be significant at high incidence angles. The highest drag coefficients are associated to the configuration which combines both sharp and smooth edges, indicating that the disparity in vorticity between shapes results in a high increase of the drag of the shape, compared to shapes which have all edges modelled with the same sharpness. Between BV0 and BV4 the difference in drag only becomes significant at incidence angles around 20 • . From 8, at 0°and 5°the drag between BV2 and BV4 are closer whereas the lift are similar for all three cases. At 10 • the drag between BV0 and BV4 remained very close while BV2 experienced a higher drag which increased at a faster rate than BV0 and BV4 except at α > 15 • . The large difference arising from BV2 is due to the rapid increase in induced drag in comparison to BV0 and BV4. In this case, the increase is mainly due to the contribution of induced drag which is proportional to the lift coefficient to the square. To further understand the difference in the mean drag results, the time history of the drag coefficient at 0 • is plotted on Fig. 9a. A quasi-steady flow can be observed for all shapes, with steady oscillations of seemingly similar frequency in the drag coefficient. An average drag coefficient of approximately 0.415 is obtained for BV2 and BV4, while BV0 shows a lower value of 0.397. A spectral analysis shows that for all shapes two dominating peaks were observed corresponding to two Strouhal numbers of 0.562 and 1.185. Furthermore, it is observed in Fig. 9a that BV2 and BV4 present a higher period oscil- lation corresponding respectively to a Strouhal number of 0.042 and 0.036, which is not observed for BV0. This low Strouhal number is not present in the wake analysis so a possible hypothesis is the existence of low frequency instabilities related to the sharp edges. The existence of the higher Strouhal numbers and their relation to the wake instabilities are supported by findings in other flows around bluff bodies, such as the flow around a sphere, where a Strouhal number of 1 is found at a Reynolds number of 3700 [START_REF] Gowree | Vortical structures behind a sphere at subcritical reynolds numbers[END_REF].

The flow around the BV0 and in its wake was further analysed. It can be observed in Fig. 10 that four vortices are created at the region corresponding to the eye section (x/L=0.2) of the boxfish as found by Bartol et. al. [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] and Van Wassenbergh et. al. [START_REF] Marcroft | Boxfish swimming paradox resolved: forces by the flow of water around the body promote manoeuvrability[END_REF]. Much smaller vortices are also formed at the mid height at this section (x/L=0.2). The vortices are roughly of the same magnitude and counter-rotating between adjacent vortices. They remain attached to the body as they travel downstream. In the wake, at x/L= 1.5, the four vortices are fully developed. The interaction can be observed in the wake between the upper and lower vortices and also between vortices that are side by side. A preliminary analysis of this phenomenon seems to indicate that the absence of sharp edges on the BV0 shape, compared to BV2 and BV4, results in a slight reduction of the drag coefficient and also in the disappearance in low frequency oscillations of the flow also evident in the time series in Fig. 9a. The highest Strouhal number is related to the dynamics of the convected vortices in the wake, while the lowest Strouhal number present only on BV2 and BV4 can be directly attributed to the presence of sharp edges as BV0 is modelled entirely with smooth edges. The low frequency in BV2 and BV4 is potentially due to stronger interaction between the edge vortices which are of larger magnitude when emanating from a sharp edge as opposed to a rounded one. The propagation of these vortices in the wake and the resulting interaction would thus modify the dynamics of the wake. Further analysis of the unsteady behaviour of the wake should shed some more light on this.

This tendency in the behaviour of the drag coefficient and the dominating Strouhal numbers are maintained until a pitch angle of approximately 10 • , where a high increase of the drag coefficient for BV2 is observed. At a pitch angle of 10 • the drag coefficient evolution for BV2 seems to reach a stable value with very small oscillations while for BV0 and BV4, the drag coefficient evolution follows slightly more chaotic oscillations, as can be seen on Fig. 11a. BV0 and BV4 seem to result in similar behaviour between the upper and lower regions and thus lower drag coefficients. It is likely that the similarity between the upper and lower edges in the geometries cause vortices which are of similar magnitude and undergo weaker interactions. At high incidence angles, the upper vortices grow in strength to create a selfcorrecting moment due to the low pressure region associated to them, this is supported by findings in [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] and [START_REF] Marcroft | Boxfish swimming paradox resolved: forces by the flow of water around the body promote manoeuvrability[END_REF].

At a pitch angle of 10 • , the low frequencies are favoured as they increase in amplitude while the peaks observed at 0 • at high frequencies are much reduced for all shapes. The peaks from the spectra in Fig. 9a where still present for the case of BV0 and BV4 despite being at significantly lower magnitudes as seen in Fig. 11a, unlike the BV2. Since these peaks are associated with the dynamics of the wake, their gradual elimination with increasing angle of incidence would suggest that they undergo stronger interaction resulting in faster break-down. From Fig. 12 this effect is seen to be more pronounced in BV2 due to the asymmetry in the magnitude of the vortices between the lower and upper surface and thus accelerating the break-down process. Traditionally on bluff bodies a break-down in the wake could result in the reduction of form drag, but here the rapid increase in lift resulting in induced drag makes it difficult to postulate on the form drag.

SUMMARY

Numerical and experimental wind tunnel studies have been carried out for different boxfish shapes created from mathematical representations and optimisations based on a real Ostracion cubicus specimen. Three shapes with different sharpness or curvature of the edges have been modelled and analysed to evaluate the effect of the edges in the aerodynamic performance of the boxfish shape. Experiments and numerical simulations done for BV4, the shape with all edges modelled as sharp, for a range of Reynolds number showed that the drag coefficient evolution does not show a pronounced drag crisis like observed in other shapes such as a sphere or an ellipsoid. Instead, a smooth decrease in the drag coefficient is observed with increasing Reynolds number. Comparison with the results obtained by Bartol et al. [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] and Van Wassenbergh et al. [START_REF] Marcroft | Boxfish swimming paradox resolved: forces by the flow of water around the body promote manoeuvrability[END_REF] show that this tendency is expected for similar boxfish shapes.

An analysis of the evolution of the aerodynamic coefficients as a function of positive pitching angle shows that for the angles, up to 20 • , there is no stall behaviour on the BV4 shape. This finding is supported by study carried out by Bartol et al. [START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF]. The drag coefficient evolution with pitch angle follows a near parabolic behaviour. A plateau, where the drag coefficient does not change greatly, is observed for low pitching angles experimentally for BV4.

A numerical analysis on the aerodynamic performance for all shapes at a Reynolds number of 3000 shows that for low incidence angles the shapes present similar drag and lift coefficients regardless of the sharpness or curvature of the edges. At low incidence angles the flow is characterised by the presence of four vortices being formed at the edges of the shapes which grow as they travel downstream. In the wake, the interaction of the system of vortices drives the oscillation of the wake at a Strouhal number of 0.562. The vortex shedding between the vortices of the right and the left of the shape occur at two Strouhal numbers of 0.562 and 1.185. At higher pitch angle, starting at 10 • , a high increase in the drag coefficient of BV2 is observed. An analysis of the vorticity shows that the upper and the lower being of significantly different magnitudes, while the sharp edge ones being more dominant. This asymmetry in vorticity magnitude affects the dynamics of the wake significantly and potentially provokes an earlier breakdown. Further analysis of the propagation of the edge vortices and their interaction in the wake needs to be conducted both from a steady and an unsteady framework.

Figure 4 :

 4 Figure 4: BV4 drag coefficient function of the Reynolds number, from simulations (△) and experiments (•)

Figure 6 :

 6 Figure 6: Location of the trip position. Tripping was implemented by employing a sufficiently thin tape wrapped around the shape.

  Experimental drag coefficient polar curves for BV4 at varying Reynolds number. Experimental lift coefficient polar curves for BV4 at varying Reynolds number.

Figure 7 :

 7 Figure 7: Experimental aerodynamic coefficient polar curves for the BV4 shape at varying Reynolds numbers.Results by[START_REF] Weihs | Body-induced vortical flows:a common mechanism for self-corrective trimming control in boxfishes[END_REF] presented for comparison.

Figure 8 :

 8 Figure 8: Numerical predictions of the aerodynamic coefficients for the BV0, BV2, and BV4 shapes at positive pitch angles.

  Time history of drag coefficient for BV0 (blue), BV2 (red), and BV4 (yellow) at 0 • positive pitch angle and Re=3000.

  Spectral analysis of the drag coefficient variation for BV0 (blue), BV2 (red), and BV4 (yellow) at 0 • positive pitch angle.

Figure 9 :

 9 Figure 9: Time history of drag coefficient and spectral analysis for all shapes at 0 • .

Figure 10 :

 10 Figure 10: Vorticity along BV0 shape and wake at 0 • pitch angle. Cuts (1, 2, 3, 4, 5) at (x/L=0.2, 0.5, 1.5, 2.0, 2.5) from leading-edge.

  Time history of drag coefficient for BV0 (blue), BV2 (red), and BV4 (yellow) at 10 • positive pitch angle and Re=3000.

  Spectral analysis of the drag coefficient variation for BV0 (blue), BV2 (red), and BV4 (yellow) at 10 • positive pitch angle.

Figure 11 :

 11 Figure 11: Time history of drag coefficient and spectral analysis for all shapes at 10 • .

Figure 12 :

 12 Figure 12: Vorticity along BV2 shape at 10 • pitch angle. Cuts (1, 2, 3, 4, 5) at (x/L=0.2, 0.5, 1.5, 2.0, 2.5) from leading-edge.

Table 2 :

 2 Mesh independence study on BV4

	Mesh	mean (C D )
	M1	0.4126 (+0.56%)
	M2	0.4152 (+1.1%)
	M3	0.4107

Table 1 :

 1 Mesh independence study

	Mesh Points on the chord dx at finest level dt at finest level lattice nodes mean (C D )
	M1	68	0.0146	0.00146	4551680	0.409
	M2	137	0.0073	0.00073	36413440	0.401
	M3	204	0.0049	0.00049	86676480	0.400

  Figure 5: Drag coefficient of several geometries considered as bluff bodies as a function of Reynolds numbers.

				BV4 experiment
				BV4 tripped experiment
				BV4 simulation
				O.meleagris, Bartol 2004
				O.cubicus, Van Wassenberg 2015
				Sphere, NASA 2010
	10 0			Sphere simulation Elipsoid, Prandtl 2009
	D
	C
	Author	boxfish	Re	Cd
	Numerical	BV4	65 000 0.2829
	Experiment	BV4	65 000 0.2300
	Van Wassenbergh 2014 Ostracion 68 000 0.2814
	[8]	cubicus
	Bartol 2004

ACKNOWLEDGMENTS

We thank Marc Grellet (ISAE-Supaero) for providing assistance and guidance during the experimental wind tunnel testing campaign carried out. We also thank Sylvain Belliot for his precious help on the coating of the boxfish model use in the experiment.


Arnaud Heibig
email: arnaud.heibig@insa-lyon.fr

Pavel Krejčí
email: pavel.krejci@cvut.cz

Adrien Petrov
email: apetrov@math.univ-lyon1.fr

Solvability for a dynamical beam problem on a frictionally damped foundation under a moving load

Keywords: Euler-Bernoulli beam equation, Coulomb friction law, existence and uniqueness results 2008 MSC: 35B45, 35D30, 49J40, 74H20, 74M10

This paper deals with a dynamic Euler-Bernoulli beam of infinite length subjected to a moving concentrated Dirac mass. The beam relies on a foundation composed of a continuous distribution of linear elastic springs associated in parallel with a uniform distribution of Coulomb friction elements and viscous dampers. The problem is stated in distributional form, and the existence and uniqueness results are established by means of a combination of L ∞ -L 2 -L 1 estimates together with a monotonicity argument. Traveling wave solutions are studied in detail in the case without Coulomb friction, and they are shown to be globally exponentially stable under positive viscous damping.

Introduction

The increasing demand for sustainable development allowed a comeback of railway transport. This is a clear trend especially in Asia and Europe. The share in passenger and freight traffic increases during the last decades and it is accompanied with the development of high-speed railways. This has in parallel stimulated mathematical models trying to describe the interactions between the rails and foundation which are in persistent frictional contact. The physical properties of the foundation, the so-called ballast composed of stones having different shapes and sizes, play a substantial role in the problem. The ballast track structures are successfully used nowadays all around the world in high-speed railway tracks. However, the track settlement is one of main causes of track deformations leading to a considerable maintenance effort. Thus a good understanding of this complex deterioration track mechanism is crucial to preserve the quality required by traffic safety and passenger comfort and to reduce track maintenance costs. The reliability of a railway track depends essentially on the limitation of the dynamic amplifications of the foundation behavior [START_REF] Madshus | High-speed railway lines on soft ground: dynamic behaviour at critical train speed[END_REF][START_REF] Kaynia | Ground vibration from high-speed trains: prediction and countermeasure[END_REF]. It depends strongly on the ability of this structure to dissipate the energy transmitted by the transit of the moving loads. The train velocity, traffic loads, number of load cycles and other factors increase the energy that should be dissipated. For all these reasons, many studies are devoted to the foundation behavior (see [START_REF] Suiker | The Mechanical Behaviour of Ballasted Railway Tracks[END_REF][START_REF] Kaynia | Ground vibration from high-speed trains: prediction and countermeasure[END_REF] and the references therein). It is also highlighted in [START_REF] Suiker | The Mechanical Behaviour of Ballasted Railway Tracks[END_REF][START_REF] Drożdziel | Simulation of railway track deterioration influenced by ballast stiffness and dry friction[END_REF] that the properties of the foundation contribute substantially to possible deformations of the rail infrastructure, while the reduced scale experiments in [START_REF] Shaer | Experimental settlement and dynamic behavior of a portion of ballasted railway track under high speed trains[END_REF] established that the sleepers subjected to sufficiently high accelerations increase strongly the settlement of the ballast.

Most mathematical models treat the rails as Euler-Bernoulli beams subjected to moving loads, and their properties are intensively studied. In the linear case, dynamic stresses are identified for beams under a moving load in [START_REF] Krylov | Über die erzwungenen Schwingungen von gleichförmigen elastischen Stäben (On forced vibrations of prismatic elastic beams)[END_REF][START_REF] Timoshenko | Erzwungene Schwingungen prismatischer Stäbe (Forced vibrations of prismatic beams)[END_REF] and the transverse vibrations are investigated in [START_REF] Frýba | Vibration of Solids and Structures under Moving Loads[END_REF]. Note that an analytical solution is exhibited for the free vibrations problem of a beam lying on an elastic foundation in [START_REF] Timoshenko | Vibration Problems in Engineering[END_REF], and [START_REF] Dumont | Vibrations of a beam between obstacles. Convergence of a fully discretized approximation[END_REF][START_REF] Dumont | Dynamic contact of a beam against rigid obstacles: convergence of a velocity-based approximation and numerical results[END_REF] contain some discrete approximations and numerical results for the problem of a vibrating beam between two rigid obstacles. However, the observed irreversible ballast accumulation in railway tracks experiments indicates that the main mechanism governing the interaction between the infrastructure constituents is driven by frictional contact. For this reason, it makes sense to study the dynamic rail behavior on a foundation including both frictional and visco-elastic dissipative mechanisms. Questions related to contact problems for Euler-Bernoulli beams have been studied in [START_REF] Ahn | An Euler-Bernoulli beam with dynamic frictionless contact: penalty approximation and existence[END_REF][START_REF] Andrews | On the dynamic vibrations of an elastic beam in frictional contact with a rigid obstacle[END_REF]. More precisely, frictional contact at each end of the beam is treated in [START_REF] Andrews | On the dynamic vibrations of an elastic beam in frictional contact with a rigid obstacle[END_REF], while the distributed contact along the beam is investigated in [START_REF] Ahn | An Euler-Bernoulli beam with dynamic frictionless contact: penalty approximation and existence[END_REF].

We consider in this paper a straight infinite linear elastic Euler-Bernoulli beam with a constant cross sectional area A, central moment of inertia I, mass density , and Young's modulus E, in persistent visco-elasto-frictional contact with a foundation. The mechanical properties of the foundation are represented by a continuous distribution of linear elastic springs with stiffness per unit length k associated in parallel with frictional and viscous dampers with yield limit m ≥ 0 and damping coefficient c ≥ 0, respectively. The present work is an extension of the results obtained in [START_REF] Heibig | A dynamic Euler-Bernoulli beam equation frictionally damped on an elastic foundation[END_REF], where the force F depends on a time variable t ∈ (0, T ) and space variable x ∈ R, and acts on the beam moving from left to right. Here the beam is activated by a concentrated force F moving at a constant velocity v ∈ R (see Figure 1). by a rather standard approach. Finally, some travelling wave solutions are exhibited which fit in our mathematical analysis of the problem.

The mathematical model

We focus here on a horizontal infinite linear elastic Euler-Bernoulli beam, which is characterized by a cross sectional area A, a central moment of inertia I, a mass density and a Young's modulus E, lying on a foundation. The beam is connected to a fixed foundation by system of linear elastic springs, with stiffness per unit length k associated in parallel with frictional and viscous dampers, with a maximum frictional force per unit length m and with damping coefficient c, respectively. The beam is acted by a concentrated force F moving from left to right at a constant velocity v (see Figure 1). The transverse displacement of the beam u(t, x) is governed by a partial differential equation that can be formally written as

+∞

-∞ EI, µ k, c, m F u(t, x) x vt
µu tt + EIu xxxx + m sign (u t ) + cu t + ku = Fδ(x -vt), (2.1)
with Cauchy initial data

u(0, •) = u 0 and u t (0, •) = u 1 , (2.2)
where (•)

t def = ∂(•) ∂t and (•) x def = ∂(•)
∂x denote the partial derivatives with respect to t and x, respectively, µ def = A is the mass of the beam per unit length of the beam, δ is the Dirac delta function. and µ, E, I, k are positive, and 3 The transverse displacement of the beam u(t, x) for (t, x) ∈ (0, T ) × R is governed by a partial differential equation that can be formally written as

µu tt + EIu xxxx + mξ + cu t + ku = Fδ(x -vt), ξ ∈ sign (u t ), (1.1)
with positive constants µ, E, I, k, F. The constant µ is defined as the product µ def = A of the mass density and the cross section area A and represents the mass of the beam per unit length. We denote by (•)

t def = ∂(•)
∂t and (•)

x def = ∂(•)
∂x the partial derivatives with respect to t and x, respectively. The symbol δ represents the Dirac measure on real line centered at 0, so that δ(xvt) should be interpreted in distributional sense as

T 0 R δ(x -vt)φ(t, x) dx dt def = T 0 φ(t, vt) dt
for all φ ∈ L 1 (0, T ; C 0 (R)). Last, as customary, "sign" denotes the multival-ued sign function given by the formula

sign(z) =    z |z| for z = 0, [-1, 1] for z = 0,
and we prescribe initial data

u(0, •) = u 0 and u t (0, •) = u 1 . (1.
2

)
The paper is organized as follows. Section 2 is devoted to a rigorous formulation of (1.1) in distributional form in Problem 2.1 below, and to the statement of the main existence and uniqueness theorem. The first step in the existence argument consists in choosing in Section 3 an approximation parameter ε > 0, and replacing the singular terms in (1.1) by suitable regularizations. More specifically, we regularize the sign function by its Yosida approximation f ε , and the Dirac measure by a smooth probability density ψ ε with compact support on R, with the intention to prove sufficient estimates independent of ε which will allow us to let ε tend to 0 in Section 4.

The identification of the limit in Section 5 represents the most delicate step of the proof which consists in proving, using finer properties of the Yosida approximation, that the limit velocity u t belongs to L 1 ((0, T ) × R) although the sequence of its approximations u ε t is not bounded in this space. Uniqueness is proved in Section 6 by a standard monotonicity argument. Finally, Section 7 is devoted to an explicit computation of travelling wave solutions in the case without friction. We obtain not only complete agreement with previous numerical experiments (see [START_REF] Castro | Dynamics of beams on non-uniform nonlinear foundations subjected to moving loads[END_REF], [START_REF] Dimitrovová | Critical velocity of a uniformly moving load[END_REF] and the references therein), but we prove in addition a global exponential stability result.

The mathematical model

We present in this section a rigorous meaning to (1.1). To this aim, we fix a final time T > 0 and introduce the following spaces

• W def = H 2 (R) ∩ L 1 (
R) endowed with the following norm:

w W def = R |w(x)| 2 + |w xx (x)| 2 dx 1/2 + R |w(x)| dx, • V def = L 1 (0, T ; W), • K def = u ∈ C 0 ([0, T ] × R) : u t ∈ L 1 ((0, T ) × R) ∩ L ∞ (0, T ; L 2 (R)), u, u xx ∈ L ∞ (0, T ; L 2 (R)), u tt ∈ L ∞ (0, T ; W ) ,
where W is the dual space to W. This non-standard functional framework allows us to reformulate (1.1) rigorously as follows:

Problem 2.1. For given initial data u 0 ∈ H 2 (R) and u 1 ∈ L 2 (R) find u ∈ K and ξ ∈ L ∞ ((0, T ) × R) such that (1.2) holds, and for all τ ∈ (0, T ) and for all φ ∈ V we have

• τ 0 µu tt , φ + R (EIu xx φ xx + mξφ + cu t φ + kuφ) dx dt (2.1a) = F τ 0 φ(t, vt) dt = F Φ(τ, vτ ) -v τ 0 Φ x (t, vt) dt , • ξ(t, x) ∈ sign (u t (t, x)) a. e., (2.1b)
where we put Φ(t, x) def = t 0 φ(s, x) ds, and •, • denotes the duality between W and W .

Notice that Φ and Φ x are continuous functions, so that (2.1b) is meaningful. Indeed, by using Sobolev embeddings, we have, for any x ∈ R and almost every s ∈ (0, T ):

|φ x (s, x)| ≤ C φ(s, •) H 2 (R) . We thus have sup (t,x)∈(0,T )×R |Φ x (t, x)| ≤ C T 0 φ(t, •) W dt. ( 2.2)
Furthermore, for τ > t and x ∈ R we have

|Φ x (τ, x) -Φ x (t, x)| = τ t |φ x (s, x)| ds → 0
as t → τ uniformly with respect to x ∈ R. Similarly, we have also that

|Φ x (t, y) -Φ x (t, x)| ≤ t 0 y x |φ xx (s, z)| dz ds ≤ |x -y| 1/2 T 0 φ(t, •) W dt uniformly with respect to t ∈ [0, T ].
The positive constants µ, E, I, F, k in (1.1) play no role in the mathematical analysis carried out in the next sections. Thus without loss of generality, we set them equal to 1.

The main result of this paper can be stated as follows.

Theorem 2.2. Let v ∈ R, m ≥ 0, and c ≥ 0 be given constants and let u 0 ∈ H 2 (R) and u 1 ∈ L 2 (R) be given. Then there exists a unique u ∈ K and a not necessary unique

ξ ∈ L ∞ ((0, T ) × R) solving Problem 2.1 with µ = E = I = F = k = 1.
The proof of Theorem 2.2 is given in the next sections. We first choose a small regularization parameter ε > 0 and solve a regularized problem. A priori estimates independent of ε and a compactness and monotonicity technique allow us to pass to the limit as ε → 0 and prove that the limit is the desired solution of the problem. Uniqueness is obtained by a standard monotonicity argument.

The approximate problem

The starting point in the proof of existence result is to approximate (1.1)-(1.2) by using the Yosida regularization for the sign function and a convolution-type mollifier for the delta function. In other words, for ε ∈ (0, 1), the sign function in (1.1) is replaced by a stiff response represented by the function

f ε (u) def =      1 for u ≥ ε, -1 for u ≤ -ε, u ε for u ∈ (-ε, ε). (3.1)
The Dirac measure is approximated by a probability density ψ ε of the form

ψ ε (x) = 1 ε ψ(x/ε), where ψ : R → [0, 1] is an even function of class C 2 such that supp ψ ⊂ [-1, 1] and R ψ(x) dx = 1. We then replace (1.1)-(1.2) with the equation for a new unknown function u ε τ 0 R (u ε tt φ + u ε xx φ xx + mf ε (u ε t )φ + cu ε t φ + u ε φ) dx dt = 1 ε τ 0 R ψ x -vt ε φ(t, x) dx dt (3.2)
for every φ ∈ L 2 (0, T ; H 2 (R)), coupled with initial conditions

u ε (0, x) = u 0,ε (x) and u ε t (0, x) = u 1,ε (x), (3.3)
where

u 0,ε ∈ H 4 (R) and u 1,ε ∈ H 2 (R) satisfy the conditions lim ε→0 R |u 0,ε -u 0 | 2 + |u 0,ε xx -u 0 xx | 2 (x) dx = 0, (3.4a) lim ε→0 R |u 1,ε -u 1 | 2 (x) dx = 0. (3.4b)
The existence and uniqueness of solutions to (3.2)-(3.3) have been established in [START_REF] Heibig | A dynamic Euler-Bernoulli beam equation frictionally damped on an elastic foundation[END_REF]Thm 3.2] in a more general setting using the classical semigroup theory. For our purposes, we recall the result in the following special form.

Theorem 3.1. Let T > 0. Let α ∈ C 1 ([0, T ]; L 2 (R)), u 0 ∈ H 4 (R)
, and u 1 ∈ H 2 (R) be given, and let β : R → R be a bounded nondecreasing Lipschitz continuous function such that β(0) = 0. Then there exists a unique solution u ∈ ∩ 2 i=0 C 2-i ([0, T ]; H 2i (R)) to the system

u tt + u xxxx + β(u t ) + cu t + u = α, u(0, •) = u 0 and u t (0, •) = u 1 . (3.5a) (3.5b) Remark 3.2.
We can deduce from Theorem 3.1 that for the solution u of

(3.5) we have u ∈ C 0 ([0, T ] × R) and u x ∈ C 0 ([0, T ] × R).
As an immediate consequence of Theorem 3.1 we see that Problem (3.2)-(3.3) has for every ε ∈ (0, 1) a unique solution u ε with the regularity as in Theorem 3.1. The convergence as ε → 0 and the identification of the limit will be established in the next sections.

4 Limit as ε → 0

We first derive some estimates independent of ε for solutions to (3.2)-(3.3) which will be subsequently used for passing to the limit in the approximate system and for proving that the limit is a solution to (2.1). Let us summarize the estimates in the following statement.

• T 0 R mf ε (u ε t )u ε t dx dt (4.1) + sup τ ∈(0,T ) R |u ε t | 2 + |u ε xx | 2 + |u ε | 2 (τ, x) dx ≤ C, • T 0 R u ε tt (t, x)φ(t, x) dx dt ≤ C T 0 φ(t, •) W dt (4.2)
with a constant C > 0 independent of ε.

Proof. We test (3.2) with φ(x, t) = u ε t (x, t). For every τ ∈ (0, T ) we get 1 2 R |u ε t | 2 +|u ε xx | 2 +|u ε | 2 (τ, x) dx - 1 2 R |u 1,ε | 2 +|u 0,ε xx | 2 +|u 0,ε | 2 (x) dx + τ 0 R mf ε (u ε t )u ε t + c|u ε t | 2 dx dt = 1 ε τ 0 R ψ x -vt ε u ε t dx dt. (4.3)
We rewrite the right-hand side of (4.3) as

1 ε τ 0 R ψ x -vt ε u ε t (t, x) dx dt = 1 ε τ 0 d dt R ψ x -vt ε u ε (t, x) dx dt + v ε 2 τ 0 R ψ x -vt ε u ε (t, x) dx dt = 1 ε R ψ x -vτ ε u ε (τ, x) -ψ x ε u 0,ε (x) dx - v ε τ 0 R ψ x -vt ε u ε x (t, x) dx dt, (4.4)
which implies that there exists a constant C > 0 independent of ε such that

1 ε τ 0 R ψ x -vt ε u ε t (t, x) dx dt ≤ C sup x∈R |u 0,ε (x)| + sup x∈R |u ε (τ, x)| + τ 0 sup x∈R |u ε x (t, x)| dt . (4.5)
Using the L 1 -L ∞ convolution inequality and Sobolev embedding H

1 (R) → L ∞ (R), we get sup x∈R |u 0,ε (x)| ≤ C u 0 H 1 (R) (4.6)
The same embedding provides

sup x∈R |u ε (t, x)| + |u ε x (t, x)| ≤ C u ε (t, •) H 2 (R) (4.7)
for almost every t ∈ (0, T ) and a constant C > 0 independent of ε. According to (4.4)-(4.7), we may deduce from (4.3) that

1 2 R |u ε t | 2 + |u ε xx | 2 + |u ε | 2 (τ, x) dx + τ 0 R mf ε (u ε t )u ε t + c|u ε t | 2 dx dt ≤ C 1 + τ 0 R |u ε xx | 2 + |u ε | 2 dx dt 1/2
with a constant C > 0 independent of ε. This is an inequality of the form

ẇ(τ ) + a(τ ) ≤ C(1 + w(τ )) 1/2 , (4.8)
where (˙) denotes the derivative with respect to τ . Note that (4.8) can easily be integrated and we get sup τ ∈[0,T ] (w(τ ) + ẇ(τ ) + a(τ )) ≤ C, which yields in particular that

τ 0 R (mf ε (u ε t )u ε t + c|u ε t | 2 ) dx dt + sup τ ∈[0,T ] R |u ε t | 2 + |u ε xx | 2 + |u ε | 2 (τ, x) dx ≤ C, (4.9)
and (4.1) follows.

To prove the inequality (4.2), we first notice that (3.2) holds for all φ ∈ V as a consequence of (4.1) and of the identity Due to the fact that the beam is assumed to be infinite, we cannot expect that the estimates in Lemma 4.1 would imply any compactness on the whole real axis. Instead, for any R > 0, we denote I R def = [-R, R] and consider the restrictions of u ε to (0, T ) × I R . We first notice that (u ε ) ε>0 are bounded in L ∞ (0, T ; H 2 (I R )) and (u ε t ) ε>0 are bounded in L ∞ (0, T ; L 2 (I R )). Using the Aubin-Lions lemma ([16, Theorem 5.1]), we can extract a subsequence (ε i ) i∈N such that lim i→∞ ε i = 0 and

1 ε τ 0 R ψ x -vt ε φ(t, x) dx dt = 1 ε R ψ x -vτ ε Φ(τ, x) dx - v ε τ 0 R ψ x -vt ε Φ x (t,
lim i→∞ u ε i -u L ∞ (0,T ;L 2 (I R )) = 0 (4.11)
for some u ∈ L ∞ (0, T ; H 2 (R)) and for every R > 0. We need in the sequel the following stronger convergence result.

Lemma 4.2. Let u and ε i be as in (4.11), and let R > 0 be arbitrarily chosen. Then (u ε i ) i∈N and (u ε i x ) i∈N converge strongly in C 0 ([0, T ] × I R ) to u and u x , respectively.

Proof. Using the embedding H 5/3 (I R ) → W 1,∞ (I R ) and the interpolation inequality, we deduce that there exists a constant C > 0 such that for all i, j ∈ N we have

u ε i -u ε j L ∞ ((0,T )×I R ) + u ε i x -u ε j x L ∞ ((0,T )×I R ) ≤ C u ε i -u ε j L ∞ (0,T ;H 5/3 (I R )) ≤ C u ε i -u ε j 1/6 L ∞ (0,T ;L 2 (I R )) u ε i -u ε j 5/6 L ∞ (0,T ;H 2 (I R )) .
(4.12)

Since (u ε iu ε j ) (i,j)∈N 2 are bounded in L ∞ (0, T ; H 2 (I R )), we conclude using (4.11) that (u ε i ) i∈N , (u ε i x ) i∈N are Cauchy sequences in C 0 ([0, T ] × I R ), and the assertion follows.

As a direct consequence of Lemma 4.1 and Lemma 4.2, we have the following result.

Lemma 4.3. There exist u ∈ K, ξ ∈ L ∞ ((0, T )×R), and a sequence

ε i → 0 as i → ∞ such that u ε i , u ε i
x converge for i → ∞ locally uniformly to u, u x , respectively, and

• u ε i t * u t in L ∞ (0, T ; L 2 (R)), • u ε i xx * u xx in L ∞ (0, T ; L 2 (R)), • f ε i (u ε i t ) * ξ in L ∞ ((0, T ) × R),
where the symbol * denotes the weak-star convergence.

Indeed, the identification of ξ represents a major problem. Note that the argument of [START_REF] Heibig | A dynamic Euler-Bernoulli beam equation frictionally damped on an elastic foundation[END_REF] cannot be applied for lack of pointwise convergence, and we have to proceed instead in several consecutive steps. The first one is the L 1 bound for u t established in Lemma 4.4 below. Lemma 4.4. Under the hypotheses of Lemma 4.3 we have

• u ε i tt * u tt in V , (4.13)
• T 0 R |u t | dx dt ≤ lim inf i→∞ T 0 R f ε i (u ε i t )u ε i t dx dt ≤ C, (4.14)
• 1 2 R |u t | 2 + |u xx | 2 + |u| 2 (τ, x) dx + c τ 0 R |u t | 2 dx dt (4.15) + lim inf ε→0 τ 0 R mf ε (u ε t )u ε t dx dt ≤ u(τ, vτ ) -u(0, 0) -v τ 0 u x (t, vt) dt + 1 2 R |u 1 | 2 + |u 0 xx | 2 + |u 0 | 2 (x) dx.
Proof. The convergence (4.13) follows from (4.2). To prove (4.14), we notice that for each η ∈ (0, 1) and R > 0 we have

T 0 R -R f η (u t )u t dx dt = lim i→∞ T 0 R -R f η (u t )u ε i t dx dt ≤ lim inf i→∞ T 0 R -R |u ε i t | dx dt.
We have

T 0 R -R |u ε i t | dx dt ≤ T 0 R f ε i (u ε i t )u ε i t dx dt + 2ε i RT.
Hence, we find

T 0 R -R f η (u t )u t dx dt ≤ lim inf i→∞ T 0 R f ε i (u ε i t )u ε i t dx dt
for every R > 0 and η > 0. Passing to the limit as η → 0, and by monotone convergence for R → ∞ we obtain (4.14).

To prove (4.15), we use the weak-star convergences in Lemma 4.3 to get for all τ ∈ (0, T ) that

1 2 R |u t | 2 + |u xx | 2 + |u| 2 (τ, x) dx + c τ 0 R |u t (t, x)| 2 dx dt ≤ lim inf i→∞ 1 2 R |u ε i t | 2 + |u ε i xx | 2 + |u ε i | 2 (τ, x) dx + c τ 0 R |u ε i t (t, x)| 2 dx dt . (4.16)
The next goal is to pass to the limit in (4.3) as ε = ε i → 0. It remains to prove convergence of the last term on the right-hand side of (4.3) rewritten in the form (4.4). Note that the support of ψ is contained in [-1, 1]. Hence, assuming that ε i < 1 for all i ∈ N we have

v ε i τ 0 R ψ x -vt ε i u ε i x (t, x) dx dt = v ε i τ 0 |v|T +1 -|v|T -1 ψ x -vt ε i u ε i x (t, x) dx dt = v ε i τ 0 |v|T +1 -|v|T -1 ψ x -vt ε i (u ε i x -u x )(t, x) dx dt + v ε i τ 0 |v|T +1 -|v|T -1 ψ x -vt ε i u x (t, x) dx dt.
(4.17)

Let us estimate the right-hand side of (4.17). First, by using Lemma 4.2 we observe that

v ε i τ 0 |v|T +1 -|v|T -1 ψ x -vt ε i (u ε i x -u x )(t, x) dx dt ≤ C u ε i x -u x L ∞ ((0,T )×(-|v|T -1,|v|T +1)) → 0 (4.18)
as i → ∞. The fact that u ε i xx are bounded in L ∞ (0, T ; L 2 (R)) yields that there exists a constant C > 0 such that

|u x (t, x) -u x (t, vt)| = x vt u xx (t, y) dy ≤ C|x -vt| 1/2 , (4.19)
which implies that

lim i→+∞ v ε i τ 0 |v|T +1 -|v|T -1 ψ x -vt ε i u x (t, x) dx dt = v τ 0 u x (t, vt) dt.
(4.20)

By using (4.18) and (4.20), we may conclude that lim i→+∞

1 ε i τ 0 R ψ x -vt ε i u ε i t (t, x) dx dt = u(τ, vτ ) -u(0, 0) -v τ 0 u x (t, vt) dt.

Identification of the limit

There is nothing to prove if m = 0. Indeed, using the convergences in Lemmas 4.3 and 4.4 we can directly pass to the limit as ε → 0 in (3.2) and (4.10) and check that the limit is a solution to (2.1a). Hence, assume in the sequel that m > 0. We can still pass to the limit as ε → 0 in (3.2) and (4.10) and obtain (2.1a). However, (2.1b) is not a direct consequence of the weak convergence results, and a monotonicity argument will be necessary.

A natural idea would be to choose φ = u t in (2.1a). This is, however, not possible, as u t is not sufficiently regular. Instead, we take the function ψ as in (3.2) and for α > 0 put

u α (t, x) = 1 α 2 R 2 ψ x -y α ψ z -y α u(t, z) dz dy, (5.1a) u α (t, x) = 1 α R ψ x -y α u(t, y) dy, (5.1b)
and choose φ = u α t in (2.1a). We claim that we obtain

1 2 R | u α t | 2 + | u α xx | 2 + | u α | 2 (τ, x) dx + τ 0 R mξu α t + c| u α t | 2 (t, x) dx dt = 1 2 R | u α t | 2 + | u α xx | 2 + | u α | 2 (0, x) dx + τ 0 u α t (t, vt) dt.
(5.2)

The only term which may perhaps need some explanation is the identity

τ 0 u tt , u α t (t) = τ 0 R ( u α tt u α t )(t, x) dx dt.
This follows from the general identity

w, ψ * (ψ * φ) = R (ψ * w)(x)(ψ * φ)(x) dx, (5.3)
where ψ * φ is the convolution product (ψ * φ)(x) = R ψ(xy)φ(y) dy. Formula (5.3) holds for all w, φ ∈ L 2 (R) ∩ L 1 (R) and by density of L 2 (R) ∩ L 1 (R) in W it can be extended to all elements w ∈ W .

In the identity (5.2), the convergences

u α t (t, •) → u t (t, •), u α xx (t, •) → u xx (t, •), u α (t, •) → u(t, •), u α t (t, •) → u t (t,
•) as α tends to 0 are strong in L 2 (R) for almost all t ∈ (0, T ), and u α t → u t strongly in L 1 ((0, T ) × R). The right-hand side of (5.2) can be rewritten similarly as in (4.4), more specifically,

τ 0 u α t (t, vt) dt = u α (τ, vτ ) -u α (0, 0) - v α 3 τ 0 R 2 ψ vt -y α ψ z -y α u(t, z) dz dy dt = u α (τ, vτ ) -u α (0, 0) + v α 3 τ 0 R 2 ψ vt -y α ψ z -y α u(t, z) dz dy dt = u α (τ, vτ ) -u α (0, 0) - v α 2 τ 0 R 2 ψ vt -y α ψ z -y α u x (t, z) dz dy dt.
We now argue as in Section 4 and pass to the limit in (5.2) as α tends to 0 to get the identity

1 2 R |u t | 2 + |u xx | 2 + |u| 2 (τ, x) dx + τ 0 R mξu t + c|u t | 2 dx dt = u(τ, vτ ) -u(0, 0) -v τ 0 u x (t, vt) dt + 1 2 R |u 1 | 2 + |u 0 xx | 2 + |u 0 | 2 (x) dx.
(5.4)

From (4.15), (5.4) it follows that lim inf

ε→0 τ 0 R f ε (u ε t )u ε t dx dt ≤ τ 0 R ξu t dx dt (5.5)
for all τ ∈ (0, T ). Hence, by (4.14) and (5.5) we have the inequality

T 0 R |u t | dx dt ≤ lim inf i→∞ T 0 R f ε i (u ε i t )u ε i t dx dt ≤ T 0 R ξu t dx dt,
which implies that ξ(x, t) ∈ sign (u t (x, t)) almost everywhere. We have thus proved the following existence statement of Theorem 2.2.

Proposition 5.1. Under the hypotheses of Theorem 2.2, there exists a solution to Problem 2.1.

We conclude this section by proving the following higher regularity result.

Proposition 5.2. Under the additional hypothesis

u 0 ∈ H 4 (R) and u 1 ∈ H 2 (R)
the solution to Problem 2.1 has the regularity u tt ∈ L ∞ (0, T, L 2 (R)), and (2.1a) can be stated in the form

τ 0 R (u tt φ + u xx φ xx + mξφ + cu t φ + uφ) dx dt = τ 0 φ(t, vt) dt (5.6)
for all φ ∈ V and τ ∈ (0, T ).

Proof of Proposition 5.2. We consider (3.2) in the form (3.5a), that is,

u ε tt -u ε xxxx + mf ε (u ε t ) + cu ε t + u ε = 1 ε ψ x -vt ε (5.7)
with initial conditions u ε (0, x) = u 0 (x), u ε t (0, x) = u 1 (x). By virtue of Theorem 3.1, the identity (5.7) is satisfied in the sense of L 2 ((0, T ) × R). For (t, x) ∈ (0, T )×R and σ > 0 sufficiently small we introduce the functions

u ε,σ (t, x) def = 1 σ u ε (t + σ, x + vσ) -u ε (t, x) . (5.8)
It follows from (5.7) that the functions u ε,σ satisfy the equation

u ε,σ tt -u ε,σ xxxx + m σ (f ε (u ε t (t + σ, x + vσ))-f ε (u ε t (t, x)))+cu ε,σ t +u ε,σ = 0 (5.9)
in the sense of L 2 ((0, T ) × R). We now multiply (5.9) by u ε,σ t , integrate over (0, t) × R for an arbitrary t ∈ (0, T ), and using integration by parts and the monotonicity of f ε we obtain

R (|u ε,σ t | 2 + |u ε,σ xx | 2 + |u ε,σ | 2 )(t, x) dx ≤ R (|u ε,σ t | 2 + |u ε,σ xx | 2 + |u ε,σ | 2 )( 0, x) dx.
(5.10)

Still keeping ε > 0 fixed, we pass to the limit as σ → 0 using the regularity of u ε , which yields

R (|u ε tt + vu ε xt | 2 + |u ε xxt + vu ε xxx | 2 + |u ε t + vu ε x | 2 )(t, x) dx ≤ R (|u ε tt + vu ε xt | 2 + |u ε xxt + vu ε xxx | 2 + |u ε t + vu ε x | 2 )(0, x) dx ≤ C u 0 H 4 (R) + u 1 H 2 (R) 2 (5.11)
with a constant C > 0 independent of ε. Taking into account the fact that

R |u ε xt + vu ε xx | 2 (t, x) dx ≤ R (|u ε xxt + vu ε xxx | 2 + |u ε t + vu ε x | 2 )(t, x) dx and (4.1) we obtain that sup t∈(0,T ) R |u ε tt + vu ε xt | 2 + |u ε xt + vu ε xx | 2 + |u ε xx | 2 (t, x) dx ≤ C
with a constant C > 0 depending only on the initial data, hence,

sup t∈(0,T ) R |u ε tt | 2 + |u ε xt | 2 + |u ε xx | 2 (t, x) dx ≤ C. (5.12)
The regularity (5.12) is preserved in the limit as ε → 0, which proves that u tt ∈ L ∞ (0, T, L 2 (R)). Notice also that this implies compactness of u ε t in C 0 ([0, T ] × [-R, R]) for every R > 0, which would thus simplify the convergence argument.

Uniqueness

We establish in this section the uniqueness result by using a classical approach. Namely, let u 1 and u 2 be two solutions of Problem 2.1, and let ū def = u 1u 2 and ξ def = ξ 1ξ 2 . Then for every test function φ as in Theorem 2.2 the function ū satisfies the identity

τ 0 ūtt , φ + R ūxx φ xx + m ξφ + cū t φ + ūφ dx dt = 0. ( 6.1)
We define for α > 0 similarly as in (5.1)

ūα (t, x) def = 1 α 2 R 2 ψ x -y α ψ z -y α ū(t, z) dz dy, (6.2a
)

u α (t, x) def = 1 α R ψ x -y α ū(t, y) dy, (6.2b)
and choose φ = ūα t in (6.1). Repeating the argument of (5.2)-( 5.4) we pass to the limit as α tends to 0 and get

1 2 R |ū t | 2 + |ū xx | 2 + |ū| 2 (τ, x) dx + τ 0 R m ξ ūt + c|ū t | 2 dx dt = 0. (6.3)
We have indeed ξ(t, x)ū t (t, x) ≥ 0 a. e., hence ū = 0 and u 1 = u 2 .

More generally, we can consider different initial conditions u 0 i , u 1 i for i = 1, 2. Putting ū0 = u 0 1u 0 2 , ū1 = u 1 1u 1 2 we obtain by the same argument the inequality

R |ū t | 2 +|ū xx | 2 +|ū| 2 (τ, x) dx ≤ R |ū 1 | 2 +|ū 0 xx | 2 +|ū 0 | 2 (x) dx, (6.4)
which implies continuous dependence of the solution on the initial data.

Traveling wave solutions

In order to make a link between the theory developed in the previous sections and the classical approach to (2.1a), we first notice that if the nonlinearity ξ does not appear in (2.1a), then for some special initial conditions there may exist traveling wave solutions of the form u(t, x) = w(xvt). Let us first consider one special case of (1.1) with m = c = 0, that is,

µu tt + EIu xxxx + ku = δ(x -vt). (7.1)
In terms of the new unknown function w, this equation reads

µv 2 w ii (z) + EIw iv (z) + kw(z) = δ(z), (7.2)
where w ii and w iv denote the second and fourth derivatives with respect to z. The solution can be constructed in such a way that we solve the linear autonomous equations

µv 2 w ii -(z) + EIw iv -(z) + kw -(z) = 0 for z < 0, (7.3
)

µv 2 w ii + (z) + EIw iv + (z) + kw + (z) = 0 for z > 0, (7.4)
and prescribe at z = 0 the transition conditions w + (0+) = w -(0-),

w i + (0+) = w i -(0-), w ii + (0+) = w ii -(0-), w iii + (0+) = w iii -(0-) + 1, (7.5)
where w i and w iii denote the first and third derivatives with respect to z.

The characteristic equation of both (7.3) and (7.4) reads EIλ 4 + µv 2 λ 2 + k = 0, (

hence

λ 2 = -µv 2 ± µ 2 v 4 -4EIk 2EI . (7.7)
The condition that u(t, •) belongs to L 2 (R) requires that w -(x) → 0 as x → -∞ and w + (x) → 0 as x → +∞. If µ 2 v 4 ≥ 4EIk, then the right-hand side of (7.7) is real and negative for both choices of the sign, so that the characteristic equation (7.6) has four roots on the imaginary axis, and no solution of (7.3)-(7.4) can asymptotically vanish at ±∞. The situation is different if µ 2 v 4 < 4EIk. Then (7.7) is of the form λ 2 = -p ± ri with p > 0, r > 0, so that λ = ±a ± bi with a = 1 2 p 2 + r 2p and b = r 2a .

The condition that w -has to vanish at -∞ and w + has to vanish at +∞ implies that we search for solutions of the form

(-a 3 + 3ab 2 ) α + (3a 2 b -b 3 ) β = (a 3 -3ab 2 )α + (3a 2 b -b 3 )β + 1, (7.9d)
which is equivalent to

α = α = 1 4a(a 2 + b 2 ) , β = -β = 1 4b(a 2 + b 2 )
.

(7.10)

The above computations show that a traveling wave solution in L 2 (R) exists if and only if |v| stays below the critical velocity

v c def = (4EIk/µ 2 ) 1/4 . (7.11)
A similar formula to (7.11) for the critical velocity of load travelling on finite and infinite beams was exhibited in [START_REF] Dimitrovová | Critical velocity of a uniformly moving load[END_REF] and it is in accordance with experimental and numerical results (see, for instance, [START_REF] Shaer | Experimental settlement and dynamic behavior of a portion of ballasted railway track under high speed trains[END_REF][START_REF] Castro | Dynamics of beams on non-uniform nonlinear foundations subjected to moving loads[END_REF]). The critical velocity is defined as the load velocity inducing the highest upward and/or downward deflections of the rail. In particular, upward displacements of the rail may have a negative impact on the stability and the lifetime of the system. In practice, however, the theoretical critical speed is still far from standard train velocities. Note also that the nonexistence of traveling wave solutions for supercritical velocities is not in contradiction with the existence statement in Theorem 2.2. The point is that initial conditions of the traveling wave solution in the "bad case" µ 2 v 4 ≥ 4EIk do not belong to L 2 (R), so that we can only conclude that the solution constructed in Theorem 2.2 does not have the form of traveling wave. On the other hand, in the "good case" µ 2 v 4 < 4EIk, a traveling wave solution exists in L 2 (R) only for one particular choice of the initial condition, namely

u 0 (x) def = w -(x) for x < 0 w + (x) for x > 0 and u 1 (x) def = w i -(x) for x < 0 w i + (x) for x > 0 ,
which satisfies the hypotheses of Theorem 2.2.

If the damping coefficient c > 0 is present in the equation, the counterpart of the equation (7.2) for the traveling wave solution w is of the form µv 2 w ii (z) + EIw iv (z)cvw i (z) + kw(z) = δ(z).

(7.12)

As a first observation, we see that only the sign of the product cv has an influence in the solution, that is, we get the same solution for negative damping (which physically makes no sense) when the train goes backward in space.

The characteristic equation has the form

EIλ 4 + µv 2 λ 2 -cvλ + k = 0. (7.13)
Consider for example the case cv > 0. The intersection of the straight line ω = cvλ with the "parabola" ω = EIλ 4 + µv 2 λ 2 + k for λ ∈ R is either empty, or consists of one or two points with positive values of λ. Moreover, the sum of the four complex roots is zero. Hence, for the roots of (7.13), only one of the following three cases can occur:

(1) complex roots a ± bi, p ± ri, p > 0, a + p = 0;

(2) complex roots a ± bi, two real roots 0 < p 1 < p 2 , 2a + p 1 + p 2 = 0;

(3) complex roots a ± bi, one double real root p > 0, a + p = 0.

Then (note that a < 0) the counterpart of (7.8) reads w + (z) = e az (α cos(bz) + β sin(bz)),

w -(z) =     
e pz ( α cos(rz) + β sin(rz)) in case (1), αe p 1 z + βe p 2 z in case ( 2), e pz ( α + βz) in case (3).

(7.14) System (7.5) is a linear system of four equations for four unknowns α, β, α, β, and the matrix of the system is nonsingular, hence it has a unique solution.

The same conclusion holds if cv < 0 with the difference that p and p i are negative and a is positive, so that we can only simply interchange w + and w -. We conclude that even negative dissipation may give a "good" traveling wave solution. The difference is that while for c ≥ 0 the traveling wave solution is stable with respect to perturbations of the initial condition by virtue of (6.4), for c < 0 the traveling wave solution is unstable. In fact, in the case c > 0, the traveling wave solution is globally exponentially stable, and the result can be stated as follows.

Proposition 7.1. There exist constants κ > 0 and M > 0 such that if u 1 , u 2 ∈ K are two solutions to (2.1) with µ = E = I = F = k = 1, m = 0, and c > 0, corresponding to initial conditions u 0 i , u 1 i for i = 1, 2, then the function We conclude that the inequality E 0 (T ) ≤ e -(cη/2)T E 0 (0) holds for all T > 0, and the assertion follows from (7.21).

E(t) = R (|ū t | 2 + |ū xx | 2 + |ū| 2 )(t,

Figure 1 :

 1 Figure 1: An infinite beam on a frictionally and viscously damped continuous elastic foundation under a moving load.

Figure 1 :

 1 Figure 1: An infinite beam on a visco-elastic foundation subjected to frictional damping driven by a moving load.

Lemma 4 . 1 .

 41 The solution u ε to (3.2)-(3.3) satisfies the estimate

  with (3.4), (4.16), and (4.21), we obtain (4.15).

w

  -(z) def = e az (α cos(bz) + β sin(bz)), (7.8a) w + (z) def = e -az ( α cos(bz) + β sin(bz)), (7.8b) where (α, α, β, β) ∈ R 4 . Conditions (7.5) read α = α, (7.9a) -a α + b β = aα + bβ, (7.9b) (a 2b 2 ) α -2ab β = (a 2b 2 )α + 2abβ, (7.9c)

def= 2 t 0 E 0 0 E 0 2 T 0 E 0 0 E 0

 2000020000 E 0 (t) + cη (s) ds is therefore nonincreasing in (0, +∞), hence the Stieltjes integral inequality 0 ≥ T 0 q(t) dG(t) = q(t)G(t) positive absolutely continuous function q. Choosing in particular q(t) = e (cη/2)t (s) ds , and we obtain from(7.23) that e (cη/2)T E 0 (T ) + cη (s) ds -E 0 (0) ≤ cη 2 e (cη/2)T T (s) ds.

  x) dx,(7.15) where ū = u 1u 2 , satisfies the estimateE(t) ≤ M e -κt .(7.16)Proof. The argument of Section 6 yields for every τ > t ≥ 0 the following counterpart of (6.3) with m = 01 2 R |ū t | 2 + |ū xx | 2 + |ū| 2 (τ, x) dx + c

	for all τ > t ≥ 0. There exists a constant η ∈ (0, 1) such that
					ηE 0 (t) ≤ E(t) ≤	1 η	E 0 (t)	(7.21)
	for all t ≥ 0. Hence, it follows from (7.20) that
					E 0 (τ ) -E 0 (t) +	cη 2	t	τ	E 0 (s) ds ≤ 0	(7.22)
	for all τ > t ≥ 0. The function
					G(t)
							τ t R	|ū t | 2 (s, x) dx ds
	=	1 2 R	|ū R		2 c	|ū| 2 + ūt ū (t, x) dx.	(7.19)
	We now multiply (7.19) by c/2, add the result to (7.17), denote
	E 0 (t)	def =	1 2 R	|ū t | 2 + |ū xx | 2 + 1 +	c 2 2	|ū| 2 + cū t ū (τ, x) dx,
	and obtain the inequality
					E 0 (τ ) -E 0 (t) +	c 2	t	τ	E(s) ds ≤ 0	(7.20)

t | 2 + |ū xx | 2 + |ū| 2 (t, x) dx. (7.17) Furthermore, we can test (6.1) with m = 0 by φ = ū, which is an admissible test function, and obtain d ds R c 2 |ū| 2 + ūt ū (s, x) dx + R |ū xx | 2 + |ū| 2 (s, x) dx = R |ū t | 2 (s, x) dx. (7.18) Integrating (7.18) from t to τ yields that R c 2 |ū| 2 + ūt ū (τ, x) dx + τ t R |ū xx | 2 + |ū| 2 (s, x) dx = τ t R |ū t | 2 (s, x) dx +

and by the European Regional Development Fund, Project No. CZ.02.1.01/0.0/0.0/16 019/0000778.