LOW COMPLEXITY SHAPE OPTIMIZATION & A POSTERIORI HIGH FIDELITY VALIDATION

. In this paper, shape optimization tools are applied to coastal engineering. A barred beach protection study is described, using a new type of water wave attenuator device, namely geotextile tube, to reduce the mobilization of sediments. A global optimization algorithm, able to pursue beyond local minima, is used to search for the optimum properties of the geotextile tube.The optimal conﬁguration is then post-validated for generated currents.


1.
Introduction. Despite beach protection becomes a major problem, shape optimization techniques have never been used yet for coastal engineering applications. Groynes, breakwaters or many other structures are used to settle down water waves impact or to control sediment flows but their shapes are usually determined using simple hydrodynamical assumptions, structural strength laws and empirical considerations. In this study, we expose a general modelling for a new type of defense structures, namely geotextile tubes, and apply a semi-deterministic algorithm to the problem of beach protection between Sète and Marseillan (NW Mediterranean sea, Gulf of Lion, France). We also discuss the circulation of water waves, ie unidirectional currents generated by waves and velocity, induced by external forcing. This investigation appears to be crucial, since sand mobilization is directly depending on the intensity of the orbital velocity, corresponding to the motion of the fluid driven by waves, on immersed beach.
The paper is structured as follows. In section 2, we recall oceanographer observations of erosion and introduce geotextile tubes. Section 3 presents the flow model used for the water waves propagation and the quasi-3D hydrodynamic model used for currents circulation simulations. Section 4 is dedicated to the description of the minimization problem. Section 5 displays optimization results for beach protection. Finally, in section 6, numerical results are discussed.
2. Geotextile tubes . Water waves propagating toward the coast are characterized by their height H, period T , wavelength L and direction θ. The period does not change during the propagation but the direction, the wavelength and the height may vary by scattering on devices.
Among these characteristics, the height H is the more crucial for the study of an erosion problem. Indeed, the mobilization of sediments produced by water waves action is the main cause of erosion process. This is essentially linked to the associated mechanical energy E = 1 8 ρgH 2 by area unit where ρ is the water density and g the gravity acceleration [1,2].
We split water waves into two categories depending if the height H of the wave is higher or lower than a given critical value H lim . Former waves are destructive because they generate big mechanical energy and accelerate the erosion process. Latter waves are called "constructive waves", since they generate a suitable dynamic through the natural beach evolution. These two observations are the basis of the principle of optimization for the protection of immersed beach. Destructive storm water waves must be mitigated and constructive weather water waves propagation must be favored. From this principle, we establish, for the four dominant propagation direction on the site of Sete (South, South South East, East South East, East), data tables for constructive and destructive water waves (See Table 1 In order to reduce water waves impact along the shoreline many devices have been tested mostly based on emerged break-waters or groynes, built with rocks or concrete. However, these techniques are expensive and are not a long-term solution for beach protection, as they usually only shift the erosion process to the close neighborhood.
Currently, a new type of attenuator device is considered acting more softly than traditional emerged break-waters. These devices are geotextile tubes, also called geotubes, and are mainly used to act on destructive water waves and prevent so the extracting of sediments. These geotubes are long cylinders made of synthetic textile and filled with sand.
3. State equations. In this section, we describe the set of state equations considered during optimization and the a posteriori validation of other aspects, mainly flow current pattern, not accounted for in the design process.
3.1. Wave propagation by Ref/Dif model. The water wave propagation and the transformation of a forward scattered wave field along an irregular mild slope is computed by using the refraction-diffraction model REF/DIF developed at the Center for Applied Coastal Research (University of Delaware, US) [3,4]. It must be stated that this model does not account for the reflection phenomenon, which appears when water waves reach an obstacle like a harbor with vertical emergent structures. However, applications introduced here only concern water waves propagation toward sandy beaches.
This model is based on a combination between the following elliptic mild-slope equation (so-called Berkhoff's equation) [5,6,7] ∇ · (CC g ∇ξ) where h is the depth, k the wave number, g the gravity acceleration, C = g k tanh kh the velocity of the wave and C g = C (1+ 2kh sinh 2kh ) 2 the group velocity, and the following parabolic model for the diffraction [8,9], where A the complex amplitude related to the water surface displacement by ξ(x, y) = A(x, y)e ikx and Moreover, the wave number satisfies the dispersion equation where ω is the angular frequency.

3.2.
Current analysis by Shorecirc model. Shorecirc is a so-called quasi-3D hydrodynamic model that computes the nearshore circulation on a rectangular grid, under the combined action of wind and waves. It was developed by the University of Delaware in the framework of the Nearcom project [10]. Now, several modified versions can be found [11,12], including the one developed at University Montpellier II [13], which can accurately deal with irregular shoreline, emerged structures, Coriolis forces and Neumann boundary conditions. A wave field on the whole simulation grid is mandatory and can be imported from any wave propagation tool that can compute monochromatic waves and provides wave forcing (like REFDIF for instance). Non-linear wave interactions, including generation of infra-gravity waves, are simulated [14].
Like many other nearshore circulation tools, SHORECIRC solves Euler equations integrated over the water column and over the period of a monochromatic wave, assuming first that the total velocity U is splitted in U = U circ + U orb + U turb where U circ refers to large scale and low frequency motion of the fluid (driven by wind), U turb is the turbulent motion (the small scale phenomena that may not be considered at nearshore scale) and U orb refers to the motion of the fluid driven by waves (which is general varying over depth).
This velocity is also splitted into two parts U orb = V + U w where U w corresponds to the fluid velocity induced by particules movement. This quantity is computed using the linear theory of water waves [2]: U w = −∇φ where φ is a potential defined by φ(x, y, z, t) = Hg 2w where k = (k x , k y ) is the vector wave, H is the height and w angular frequency. These quantities are computed by REF-DIF. V represents the vertical variation of the current and satisfies: where i = 1, 2 denote horizontal coordinates, ζ is the instantaneous water surface elevation, h 0 is still water depth and ζ t corresponds to the surface elevation of the wave trough level. This equality conveys that any varialibity is according to the vertical is provided by water waves forcing. Next short-wave averaged, we finally obtain the following equations, with Einstein notation: where S ij are the radiation stresses induced by the wave. Again, this quantity is computed using REF-DIF. Indeed, in (8) the radiation stress is defined by p being the total pressure and δ ij is the Kronecker delta function. T ij , τ B i et τ S i are the depth-integrated Reynolds'stresses, the bottom shear stresses, and the surface shear stresses, respectively. The term L ij arises from the depth-varying currents and is computed using a perturbation method, we refer to [15,16].
Finally, equations (7) and (8) provide information over the mean water surface and U circ . 4. Optimization problem. As already stated above, the longshore sediment drift is mainly due to water waves breaking on the barred beach with a large amount of mechanical energy.
The aim of the device is to prematurely release this energy by breaking the water waves sufficiently far away from the coastline.
An area between 100 and 250 meters far from the coastline is considered for the cost function computation, corresponding to the trench between the first and the second natural sand bar.
As mentioned in section 2 and Table 1, we consider two categories of water waves for the study, constructive and destructive, defined by comparison with the critical water height H lim = 2m. For a given direction of propagation θ, the following cost function for minimization is considered where, for a given point (x i , y i ) in D = [100, 250] × [0, 2400], the energy is: where ρ is the water density and H(x i , y i ) = 2A(x i , y i ), E H>H lim and E H<H lim are respectively energies corresponding to destructive (large) and constructive (small) waves.
Numerical results shown in the newt section are obtained using a multi-directional optimization. In this case, we set as new cost function where p θ is the observation frequency for a given wave direction. From a physical point of view, this cost function aims at decreasing the energy for the destructive water waves and to be non active for the constructive water waves.
The absolute minimun of this cost function J is the best configuration, following the optimization principle stated in 2. We can note that this principle corresponds to a minimization to the potential energy of the water waves, which is equivalent to a minimization of the orbital velocity in the theory for linear water waves [2]. Consequently, the aim of the shape optimization is to find the optimal configuration to minimize the orbital velocity, ad therefore minimize the percentage of mobilized sediments [17].  The studied area is 2.4 km long. The propagation is computed for water waves data at 1.2 km far from the shoreline. We would like to avoid the tube to be located at the top of the second natural sandy bar. This is imposed as a constraint for the optimization. Also the height and width of the tube are imposed. Starting from data given by CANDHIS (National Center Archive for In Situ Wave Data, http://www.cetmef.equipement.gouv.fr/donnees/candhis/), we compute, for the two categories of water waves, the following mean height H s , mean period T s and mean frequency of observation p for four significant directions of propagation. Corresponding data are displayed in Table 2 We can observe, in Table 2, that the optimized geotube is non-active for the constructive water waves and also equally efficient for the four propagation directions considered.  Figure 2 shows the comparison between the height H resulting from the propagation over the optimized geotube and the height for a water wave propagation over the initial barred beach bathymetry. The water wave considered is the destructive SSE water wave (See Table 1). We notice that H is indeed lower between 100 and 250 meters far from the shoreline when we locate a geotube to 353 meters of the coast. This shows that a geotube laid immediately downstream from the second natural sandy bar enables to break the water wave close to the shore. Figure 3 shows that this optimized configuration does not increase the bottom orbital velocity, i-e the sea bottom fluid particle velocity, compared to the initial configuration. We therefore carry out a numerical study to the circulation for weather and storm water waves propagation, with variable angle of incidence. We realize this study on differents configurations. See figure 4.

5.2.1.
Presentation of the simulation context. We use SHORECIRC on a rectangular domain 1210m×1340m and a rectangular grid of 121 × 134 cells. The wave module REFDIF is used to compute the wave field over the domain with a rectangular grid of 1210 × 1340. The refinement of the REFDIF grid is mandatory to get relevant wave numbers and wave directions when fair weather waves are propagated. Neumann condition (with derivative of velocity U and free surface ζ set to zero) were systematically used on the lateral boundaries [13]. A no-flux condition was used at the shoreline.
Forcing used for the simulations include both fair-weather and storms conditions for each realistic wave incidence (see Tables 1 and 5.2.1). Simulations are performed respectively for 10000 and 30000 seconds with a constant wave forcing. Results at 10000 are very similar to those at 30000: wave angle does not shift, velocities increase very slightly. This suggests that a quasi-stationary state is reached for the circulation before 30000 seconds. Results presented in this study are that at 30000 seconds and correspond to a quasi-stationary flow. Computations with three distinct bottom conditions are performed. The bottom condition 1 corresponds to the unprotected scenario, with a water depth profile that fits the shape of the mean seabottom profile during fair weather periods. Bottom condition 2 corresponds to the same water depth profile, with a 12 meters large geotube deployed at 550 meters from the coastline all along the domain. Bottom condition 3 corresponds to the same water depth profile with a 12 meters large geotube deployed at 350 meters from the coastline. In addition, two alternative bottom conditions are tested. Bottom condition 2bis corresponds to bottom condition 2 with a limited geotube extension, with no protection from the lateral boundary up to 300 meters to the centre of the  Figure 4. Bottom conditions for the simulations of current velocity with SHORECIRC. Bottom condition 1 corresponds to a natural immersed beach with no protection device. Bottom condition 2 corresponds to a geotube deployment on the outer sand bar, relatively far from the crest of the bar (in this case at 550 m far from the coast). Bottom condition 3 corresponds to the same type of deployment at 350 meters from the coastline. This last bottom condition is the one that results from the process of shape optimization on the basis of wave energy. Alternative extruded geotubes were also considered to test the numerical effect of geotube boundaries on the circulation. Dash line on the colored maps sign the location of the cross shore water depth profile plotted above.
domain. Bottom condition 3bis corresponds to bottom condition 3 with a geotubed limited in the same way. Bottom conditions 2bis and 3bis are introduced to analyze the circulation patterns at the boundaries of the geotube, while bottom conditions 2 and 3 highlight circulation features for an infinitely long geotube. Figure 4 displays the bottom conditions used in this study.

5.2.2.
Results and analysis. Figure 5 displays typical circulation for a wave angle normal to the coast,H = 2.5m and T = 8s. Figure 6 displays typical circulation for θ = 15 degrees oblique to the coast, the same wave amplitude and same period. All maps are plotted with an unique colour scale. Such results are representative of all simulations performed. Modify the wave amplitude or the wave period do not change the interpretation and the inferred conclusions. A preliminary analyze of the whole set of computations highlighted two distinct types of wave climate: normal wave incidence and oblique wave incidence (angle varying between a few degrees up to 30). Normal wave incidence. Barotropic circulation under normal wave incidence is characterized by a concentration of hydrodynamic features in the uppermost immersed beach. Maximal velocities observed are about 0.15m.s −1 for fair weather waves and reached 1.6m.s −1 locally for storm waves. In both cases, peculiar patterns developed: 20 to 30 cells large gyres occurred in less than 1 meter of water depth. These patterns repeated several times along the profile. They occurred since the first steps of the simulations until the quasi-stationary state at 30000 seconds. They do not vanish if the resolution of the cells is changed. For bottom condition 3bis, some stronger and small-scale circulation patterns with respect to bottom condition 1 occur in the lee of the geotube, close to its lateral boundaries. Such patterns are limited along shore to 30 to 40 cells. On the simulations for bottom condition 3, such patterns disappear. These evidences limited but significant perturbations of the circulation at the lateral terminations of the geotube if it deployed in the optimized configuration. Having this exception in mind, numerical results for normal wave incidence demonstrate that the wave-driven circulation is not significantly altered in presence of geotube, whatever may be their location on the cross-shore profile.
Oblique wave incidence. Barotropic circulation for 15 degrees incident waves is characterized by a well-expressed along shore component for any computation.  Last, for bottom conditions 3 and 3bis, the length along any cross-shore profile where the velocity is greater than e.g. 2m.s −1 is smaller than this one observed for bottom condition 1. This is not true for bottom conditions 2 and 2bis. This result points out a significant reduction of the along shore current at the lee of the geotube if it deployed close to the crest of the sand bar. Results presented herein concern the barotropic circulation. SHORECIRC being a quasi-3D circulation tool, it may have been possible to discuss on the basis of surface of sea bottom currents. However, a preliminary analyze of such quasi-3D velocities highlighted that for a wave incidence greater than 3 or 4 degrees, surface current, bottom current and barotropic current display the same orientation with decreasing norm of the velocity from the surface to the bottom. As a consequence, results presented here may have not changed if one have considered surface current or bottom current, and the quasi-3D development is not mandatory.
Simulations for normal wave incidence demonstrate that the optimized configuration (bottom condition 3 and 3bis) may introduced some significant perturbations of the circulation in the lee of the lateral boundaries of the geotube. These alterations could lead to near field scour. However, a normal wave incidence is very rare, and such results may be considered with caution. Such peculiar features disappear as soon as the wave angle is about 3-4 degrees, that is a minimum in real case conditions usually. As a consequence, simulations for 15 degrees incident wave, representative of any realistic storm forcing, demonstrate that optimized configuration does not introduce stronger cross shore or long shore current with respect to the unprotected configuration. On the contrary, it may reduce the longshore current in the lee of the coastal defense structure. In other words, the geotube configuration resulting from the optimization with a cost function derived from wave energy (wave orbital velocity) is relevant from a circulation point of view. This circulation analyze validates a posteriori the choice derived from the numerical optimization study.
6. Conclusion. Shape optimization procedure has been applied to a beach protection device based on geotextile tubes. This device permits to reduce the wave energy, orbital velocity and free surface elevation. The optimal configuration is post-validated for generated currents. It is highlighted that secondary longshore currents are not dominant when considering the protected configuration. This numerical design has then been experimentally validated in wave channel and wave tank. We show that the device acts as a low-pass filter, highly efficient for beach protection and that shape optimization tools can be accurately and efficiently applied to reduce beach erosion effect. 7. Appendix: Global optimization algorithm. In this appendix, we briefly present our global optimization method [18,19,21,20,22,23,24].
Consider the minimization of a functional J(ψ) ≥ 0, ψ ∈ O ad , where O ad is an admissible domain.We suppose the problem admissible (i.e. there exist at least one solution ψ m to the problem: J(ψ m ) = J m ). Most minimization algorithms can be seen as discretizations of: M is aimed to be positive definite and M −1 d is built to be an admissible direction. Global solution of (12) means, for instance, finding ψ(1) verifying This is an overdetermined boundary value problem and it tells us why one should not solve global optimization problems with methods which are discrete form of Cauchy problems for first order differential systems. Except if one can provide an initial condition in the attraction basin of the global optimum. Hence, in the context of global optimization the initial condition (ψ(0) = ψ 0 ,) is misleading. On the other hand, one would like to realize the optimality condition (J (ψ(1)) = 0). To remove the overdetermination we consider the following second order system with two final conditions: This can be solved using solution techniques for BVPs with free surface. An analogy can be given with the problem of finding the interface between water and ice which is only implicitly known through the iso-value of zero temperature. 7.1. Recursive minimization algorithm. Global optimization of J has a solution if, for a given precision ε in the functional, one can build at least one trajectory (ψ(t), 0 ≤ t ≤ T ε ) passing in finite time T ε enough close to ψ m . In what follows, these trajectories will be generated by first or second order Cauchy problems starting from an initial value v. These will be denoted (ψ v (t), 0 ≤ t ≤ T ε , v ∈ O ad ). Existing minimization methods build such trajectories in discrete form. This can be summarized as: If J m is unknown T ε defines the maximum calculation complexity wanted. In these cases, setting J m = −∞, one retains the best solution obtained over [0, T ε ]. In other words, one solves: T ε ] such that J(ψ v (τ )) − J m ≤ ε (16) Below, we propose, a recursive algorithm to solve (15) or (16): Given ε > 0, J m and 0 ≤ T ε < ∞, we minimize the functional h ε,Tε,Jm : O ad → IR + below: h ε,Tε,Jm (v) = min ψv(τ )∈O ad , τ ∈[0,Tε] Hence, global minimization becomes a nested minimization problem where one looks for v to provide a better initial value to generate the trajectory (ψ v (τ ), 0 ≤ τ ≤ T ε ) (i.e. passing closer to ψ m ). Below, one shows a recursive solution algorithm to this problem (denote h ε,Tε,Jm by h). Consider the following algorithm A 1 (v 1 , v 2 ): -(v 1 , v 2 ) ∈ O ad × O ad given -Find v ∈ argmin w∈O(v1,v2) h(w) where O(v 1 , v 2 ) = {v 1 + t(v 2 − v 1 ), t ∈ IR} ∩ O ad using a line search method -return v The line search minimization in A 1 is defined by the user. It may fail. For instance, a secant method degenerates on critical points. In this case, in order to have a multidimensional search, we add an external layer to the algorithm A 1 minimizing h (v ) = h(A 1 (v , w )) with w chosen randomly in O ad . This is the only stochastic feature of the algorithm.
This leads to the following two-layer algorithm A 2 (v 1 , v 2 ): , t ∈ IR} ∩ O ad using a line search method -return v Again, the line search minimization in A 2 is user-defined. This construction can be pursued by building recursively h i (v i 1 ) = h i−1 (A i−1 (v i 1 , v i 2 )), with h 1 (v) = h(v) and h 2 (v) = h (v) where i = 1, 2, 3, ... denotes external layers. Due to the stochastic feature mentioned, this is a semi-deterministic algorithm.