Optimum Shape Design of Geometrically Nonlinear Submerged Arches Using the Coral Reefs Optimization with Substrate Layers Algorithm

: In this paper, a novel procedure for optimal design of geometrically nonlinear submerged arches is proposed. It is based on the Coral Reefs Optimization with Substrate Layers algorithm, a multi-method ensemble evolutionary approach for solving optimization problems. A novel arch shape parameterization is combined with the Coral Reefs Optimization with Substrate Layers algo-rithm. This new parameterization allows considering geometrical parameters in the design process, in addition to the reduction of the bending moment carried out by the classical design approach. The importance of considering the second-order behaviour of the arch structure is shown by different numerical experiments. Moreover, it is shown that the use of Coral Reefs Optimization with Substrate Layers algorithm leads to nearly-optimal solutions, ensuring the stability of the structure, reducing the maximum absolute bending moment value, and complying with the serviceability structural restrictions.


Introduction
Arches are structures that provide strength through their geometry. In particular, funicular arches can span long distances without significant bending stresses [1]. However, these stresses may be large when real design conditions are considered, such as irregularities of the arch shape, support configuration, and load conditions, which implies a high effort in both design and construction, with a consequent increasing in costs [2].
The term submerged arch is used to denote an arch immersed in a fluid (usually water). These arches are very important in a number of applications in marine and civil structures [2] or aquarium tunnels. The design of these submerged structures has been previously approached from the funicular point of view. However, the funicular shape does not exist under certain design parameters [3].
The nonlinearity involved when the buckling is considered in the design of arches has been also widely analyzed in the literature. In this context, it stands out especially in [4], where the nonlinear buckling and post-buckling of elastic shallow arches via a nonlinear curved finite element model is studied. Recently, the buckling phenomenon has been studied for different scenarios, for example, with the aim of showing the elastic and elasto-plastic flexural-torsional buckling and post-buckling analysis in combined bending and compression actions [5], the validity of approximations on analytical solutions [6], the influence of non-homogeneous arches [7], or the effect of considering the shear deformation [8].
The optimal design in the context of elasticity and structural systems usually leads to very hard optimization problems [9][10][11][12][13][14][15], with extremely large search spaces and nonlinear or black-boxes objective functions [12,16], which in many occasions come from computational simulations of the potential solutions. Thus, the use of meta-heuristic approaches is completely justified in the particular case of submerged arches design, where traditional optimization procedures do not provide results in all the cases, as pointed out above. In addition, considering nonlinear effects in the problem leads to even harder optimization problems in the structural design, as specific parameterization and adaptations in the objective function must be carried out.
Different optimization techniques have been used in the optimal design of submerged arches problem, such as the reduced gradient method, which was applied in a design problem that aims to reduce the arch weight with respect to the axial force and different geometrical arch parameters [17]. It was shown in [18] that many of the arch shapes proposed in the literature may reach their critical buckling load for service conditions. Thus, a buckling design of submerged arches was proposed by using a genetic algorithm. The arch design may be approached as a multi-objective optimization problem [19]. The optimization of the arch design by reducing the bending moment along the arch curve and increasing the airspace of the arch is shown in [20]. For that purpose, an evolutionary algorithm that uses artificial neural networks with extreme learning machine was used.
In this work, the optimization problem is approached by means of the Coral Reefs Optimization with Substrate Layers (CRO-SL) method [21]. The CRO-SL is a multi-method ensemble [22,23] based on evolutionary computation, where different search procedures are applied within a single population of solutions to the problem. The CRO-SL evolves then possible solutions to a given optimization problem by applying different search operators, depending on the location of the solutions within the population. This procedure results in a very effective approach to solve hard optimization problems [24][25][26][27].
The optimal design of arch shapes by using the CRO-SL has been previously carried out in [28], obtaining excellent results for different case studies. However, a complete linear analysis was considered in that work, thus missing important nonlinear phenomena which may condition the optimal arch shape. In this paper, the problem of optimal design of submerged arches considering nonlinear effects is studied. Specifically, the design cases in [18] (where a parametric shape function is formulated using two conical shapes involved in the funicular design of submerged arches) have been extended to the context of the geometrical nonlinear analysis, through the use of a nonlinear finite element model. Thus, the optimal shape of the arch can be obtained by an iterative strategy [10], based on minimizing a functional by means of an optimization algorithm or strategy. In this work, the previous procedure is extended to the simultaneous inclusion of design criteria such as bending moment, vertical deflection, or the amount of area inside the arch curve, checking the structural stability for all the scenarios considered.
In this work, a set of search operators are defined, which mix well-known traditional operators, such as two-point crossover or Differential evolution, with other operators based on nonlinear processes [29], such as Firefly or Water Wave optimization. The proposed novel design procedure is integrated with CRO-SL in order to find the optimum arch geometry in terms of minimum bending moment and maximum airspace in the submerged arch subjected to maximum deflection and bucking restrictions. Two submerged arches designs are addressed, in order to illustrate the good performance of the proposed procedure optimized by CRO-SL.
The rest of this article is structured as follows. Section 2 introduces the problem of nonlinear submerged arches optimization and design. Section 3 presents the CRO-SL algorithm, its most important characteristics, and a description of the search procedures implemented in the algorithm. Section 4 presents some simulation results on the optimiza-tion of nonlinear submerged arches with the CRO-SL algorithm. Finally, Section 5 outlines some conclusions, remarks on the research carried out, and some future lines to continue this work.

Problem Statement: Nonlinear Submerged Arches Design
The objective of this section is to explain how geometrical nonlinear analysis is included in the optimal arch design methodology, used afterwards for encoding arches in the CRO-SL algorithm. This section starts with funicular equilibrium of submerged arches. The buckling phenomenon, as a particular case of the structural stability in submerged arches, is then explained. The adopted shape parameterization used in the geometrical optimal design is also presented. Finally, the geometrical nonlinear analysis and the optimal design methodology are detailed.

Funicular Equilibrium of Submerged Arches
The funicular equilibrium of a submerged arch is governed by the following set of differential equations in Cartesian coordinates (x, y) [1]: where D p is the depth (i.e., the distance from the water surface to the supporting line, see Figure 1); N is the axial internal compression force; σ adm is the compression stress in the arch cross section under a fully stressed condition; h is the arch thickness at the apex; and, finally, γ h and γ s are the unit weights of the water and the arch's material, respectively. Solutions to Equation (1) may be obtained by numerical integration until x = B (see Figure 1), at which y(x) = 0. The funicular arch shape can be represented in Cartesian coordinates, simplifying the traditional geometric expressions based on the arch length and the angle of the tangent [30]. From these analytical expressions, the maximum span and height of the funicular submerged arch, as well as the minimum compression force (and therefore, the minimum amount of airspace), can be determined. The hydrostatic pressure and self-weight load were considered in [31] to show the influence of the self-weight in the funicular shape, resulting in a sufficiently accurate solution for practical applications.

Structural Stability: Buckling Analysis of Submerged Arches
Structural stability represents a type of structural analysis which explores the behaviour of the structures subjected to compression loads. In essence, an equilibrium state achieved under compression load is considered unstable when the original equilibrium state is not achieved after removing the compression load.
In case of plastic behaviour, equilibrium state is assumed to be stable, if a structure even tends to return to the initial equilibrium state [32]. Nevertheless, elastic members subjected to static loads are only considered in this work.
The minimum force at which the structure no longer returns to the initial state, if all disturbing factors are removed, is the so-called critical force. If we note as N the internal axial force due to an external load, then the linear buckling load factor (α crit ) can be expressed as in which N crit is the critical axial force. Then, for α crit < 1, the buckling phenomenon occurs in the arch structure, while for α crit > 1, the stability is guaranteed, as the axial force along the structure is lower than the critical axial force. As a general rule, a loss of stability produces a change in the initial form of equilibrium or buckling, that usually leads to the structural collapse. The critical load of the arches may be determined in an analytical way only in the simplest cases, such as uniform circular arches subjected to uniform pressure normal to the axis. Thus, Timoshenko [33] only provided analytical solution of the buckling problem for circular two-hinged arches subjected to a uniform distributed radial pressure. Approximate methods of solving stability problems for arches consist of approximating the arch by a framed structure, and proceeding with analysis by the displacements method (also called the eigenvalue buckling test [34]).
The importance of considering the nonlinear effects in the submerged arch structure has been shown in [18]. Particularly, the arches corresponding to shallow waters usually present linear buckling factors lower than 1.0, thus becoming unstable under design loads. This leads to modify the slenderness of the arch (i.e., augmenting its thickness), which implies large volume and, consequently, large spans due to self-weight effect in the funicular equilibrium.
This last problem may be approached in an alternative way through a shape parameterization of the submerged arch. This strategy allows searching the optimal shape and, simultaneously, controlling its second-order behaviour during the optimization process, as explained below.

Alternative Approach: Shape Parameterization
In the case of a nonuniform arch of arbitrary shape, it becomes impossible to present the buckling solution in analytical form [32]. For these cases, the solution may be obtained using only approximate procedures, such as the variational methods. In the case of arches in shallow waters, the solution of differential equation of stability may be obtained in closed form, as previously indicated in the first epigraph of this Section. Many important solutions related to stability of arches have been obtained by Dinnik [35] or Tadjbakhsh [1]. Likewise, buckling of funicular arches under gravity loading and overburden (the so-called Inglis problem [36,37]) has been solved numerically by several authors [1,38], under different support conditions, reaching in all the cases a good adjustment. The buckling loads were calculated for an arch with a relatively high slenderness greater than 100. Moreover, an approximated solution to the buckling load has been also developed for the hinged support problem [39].
On the other hand, the optimal design of the submerged arches in terms of their secondorder response may be approached using a parameterization of the arch's centerline. In this sense, note that, in practice, the submerged arches are not strictly funicular due to aspects like geometric imperfections, or design load variations (e.g., tidal variation), among other reasons [31]. For these cases, bending moments cannot be completely eliminated, but they can be reduced resulting in a well-dimensioned arch with a negligible bending stress. In this sense, alternative design procedures, based on the arch shape optimization using parametric functions [40], are feasible. This is the strategy adopted in this work.
The solution of Equation (1) without considering self-weight load lays halfway between parabolic and elliptical types for intermediate values of ratios D p /H and/or H/B [41]. A geometrical parameterization over these two conical forms is proposed in [20]. Thus, the equations of the parabola and the ellipse intersecting at the apex and the support of the arch (see Figure 2), and having the same pole, may be expressed, respectively, as where r and θ are the radial coordinate and the angular coordinate, respectively, and H and B are the lengths of the semi-minor axis and the semi-major axis of the ellipse, respectively. From the above considerations, the funicular solution is easily approximated using a linear combination of the two conic types: where t a ∈ [0, 1] and r pe is the radial coordinate of the intermediate curve corresponding to the parameter t a and the angular coordinate θ. Equation (5) represents a parametric family of curves between parabolic (t a = 0) and elliptical (t a = 1) forms, where the shape parameterization is posed as a function of four parameters: t a , H , B and B . However, only three of them are independent, as the arch support is assumed to be contained in the ellipse defined in Equation (4), and in this case, H and B must verify Appl. Sci. 2021, 11, 5862

Geometrical Nonlinear Analysis and Optimal Design
In practice, the arch bending moment cannot be completely suppressed. Therefore, the real buckling factors (α crit ) may be lower than the ones obtained by the eigenvalue method.
A more accurate stability analysis may be performed taking into account the geometrically nonlinear behaviour of the submerged arch. For this purpose, the updated Lagrangian formulation [42][43][44] may be applied, at which all variables are referred to the current (or updated) configuration of the arch. In this context, the equilibrium of the arch may be expressed through the principle of virtual displacements in the following terms: where S is the 2nd Piola-Kirchhoff stress tensor and is the Green-Lagrange strain tensor from the configuration at time t to the configuration at time t + δt and referred to the configuration at time t, and R is the external virtual work. Equation (7) is a nonlinear one in the unknown incremental displacements. Nevertheless, higher-order terms may be neglected and the previous weak form may be linearized through the corresponding incremental decomposition of stresses and strains. Moreover, the resulting linearized expression of Equation (7) may be solved using the iso-parametric finite element method, where the coordinates of the nodal points and the corresponding nodal displacements are interpolated in order to evaluate the displacement derivatives, required in the integrals of the equilibrium equation. Thus, the linearized expression of Equation (7) may be written in terms of finite element matrices as [43] matrix), F is the vector of external nodal loads at time t + ∆t and t t F r is the vector of nodal point forces equivalent to the current element stresses (or the vector of restoring forces) at time t. U is the vector of increments in nodal displacements, T is the Cauchy's stress tensor and C is the incremental material property matrix.
The linearization of the equilibrium equations may introduce errors which ultimately result into solution instability. Due to this reason, an iterative strategy is necessary. To this aim, the design load may be divided into a set of steps. At each load step the out-of-balance load vector (i.e., the difference between the loads corresponding to the element stresses and the applied loads) is evaluated in order to check the convergence criteria. If these criteria are not satisfied, a new seed is evaluated and a new solution is obtained until the problem converges. Consequently, obtaining the second-order load-displacement history of the submerged arch it is a path-dependent problem at which loads must be applied through several steps in order to get a high accuracy. In the context of the Modified Newton's method (where the stiffness relation is only evaluated at the start of the load step: this method is implemented in the finite element software ANSYS [45]), the previous iterative process may be formulated as [46] where k is the number of iterations until verifying the convergence criteria.
As was pointed out in the first epigraph of this section, two types of loads have been considered in this work: self-weight and hydrostatic pressure. The first has been introduced using the gravitational acceleration, whereas the hydrostatic action over the submerged arch has been introduced as a pressure load. In a large-deflection analysis, pressure loads present a singular behaviour (with regard to other load types as accelerations or nodal forces), as they change their original orientation during their application, acting always normal to the deflected element surface and "following" the finite element as it undergoes large rotations. In this case, the load step must be small enough, so that the external virtual work can be approximated to sufficient accuracy using the load corresponding to the step t + ∆t and integrating over the last area calculated in the corresponding iteration [47]. This is also the scheme implemented for geometrically nonlinear analysis in the software ANSYS.
A symmetric plane simply-supported arch has been modeled using the BEAM188 element, which is well suited for large displacements analysis [48]. As the discretization into finite elements of each shape is performed by constant increments of the polar angle, the finite element size varies depending on the values of the shape parameterization, oscillating between 0.05 m and 0.3 m. An isotropic and linear elastic stress-strain relationship has been considered as constitutive model for the submerged arch, where the material properties have been characterized by the Young's modulus and the Poisson's ratio. Moreover, in order to clarify the proposed methodology, the unfavorable effects of possible deviations in the geometry of the structure and the position of loads, such imperfections, are neglected.
The weighted functional used in this work can be expressed as in which (max |M s |) is the maximum absolute bending moment value along the arch structure, while A s the area under the arch curve. This area has been also chosen as design parameter since it is directly related to the serviceability of the underwater installation [20]. The parameters p m and p a are the bending moment and area weights, respectively, such that p m + p a = 1. Thus, a dual objective can be achieved: minimization of the maximum bending moment (which determines the cross section of the arch ) and maximization of the area under the arch. One of the advantage of this formulation is the easy way to implement different design criteria. The maximum deflection allowed has been chosen to be max d y < B/150 (11) in which d y represents the deflection of the curved beam along the guideline of the arch. The calculation process is represented in Figure 3. It starts with random parameter values, which are contained in a specific range. These arch parameters are used to calculate the area under the arch (A s ), the maximum absolute bending moment (max |M s |), and the maximum vertical deflection (max d y ). The parameters (max |M s |) and (max d y ) are calculated using ANSYS software. The arch parameters are changed by the CRO-SL in the next interactions in order to minimize the functional defined into Equation (10). Thus, the numerical model is updated with the new proposed geometrical parameters. In addition, a penalty value is applied to the functional value when this iterative strategy does not converge or max (d y ) > B/150. Thus, the quality of each arch is evaluated through the functional value, so the CRO-SL evolves the population towards high quality solutions, i.e., arches which minimize the objective functional value, given by Equation (10). The CRO-SL algorithm and its substrates (search procedures) will be fully described in the next section.

The CRO-SL: A Low-Level Ensemble Evolutionary Algorithm
The concept of ensemble in optimization refers to a combination of algorithms, search strategies, or operators to better deal with optimization problems [23]. The basic idea is that the ensemble strategy can obtain better results than any single strategy on its own in the same problem. Following the work in [23], it is possible to classify optimization ensembles attending to different characteristics of the technique, mainly the type of the constituent elements and its implementation: when the ensemble technique combines different types of search strategies, operators, algorithms, etc., it is classified as a low-level ensemble. On the other hand, high-level ensembles refer to methods which try to select the best optimization algorithm for a given problem.
The CRO-SL approach is a recently proposed low-level multi-method ensemble, based on the previously defined Coral Reefs Optimization algorithm [49], a meta-heuristic approach of the evolutionary-computation family. Specifically, in the CRO-SL different search operators are jointly applied and self-adapted within the same population, by means of defining a set of "substrates" in the CRO reef [21]. In this case, each substrate represents the application of a different search operator in the same population.
The CRO-SL algorithm was proposed as an advanced version of a basic original version of the CRO [49]. The CRO is an evolutionary-type algorithm in which the search operators are based on the processes occurring in a coral reefs, including reproduction, fight for space, or depredation, see in [50] for details on the algorithm. The CRO-SL, in turn, was first introduced in [51], and further developed in [21]. The CRO-SL is a general ensemble approach which promotes competitive co-evolution, where each substrate layer represents different processes, in this case, different search operators. The CRO-SL has been successfully applied to a large number of hard optimization problems in science and engineering, including different problems related to structural design [24,25].

Substrate Layers Defined in the CRO-SL
Very different search strategies can be defined in the CRO-SL as part of the multimethod approach. They are usually defined at the practitioner's discretion. In previous approaches, it has been shown that combinations of well-known traditional operators with novel search techniques provide the best set of techniques (substrates) in the CRO-SL. Some of the most used traditional operators in previous applications of the CRO-SL are Harmony Search (HS), Differential Evolution (DE), classical two-points crossover (2Px), classical multi-points crossover (MPx), and Gaussian-based mutation (GM). In this paper, we combine these traditional search strategies with two new operators for search and optimization, recently proposed: the Firefly algorithm (Fa) and the Water Wave Optimization (WWo) approach. Previous simulations were developed considering each substrate individually. All the functional values for both cases analyzed in the paper were above the functional value obtained using a combination of all the substrates. Thus, demonstrating a worst performance when they work independently than when working assembled in the CRO-SL. We provide here a short definition of these substrates for the CRO-SL ensemble. DE: Mutation from Differential Evolution algorithm (with F = 0.6). This operator is based on the evolutionary algorithm with the same name [53], a method with powerful global search capabilities. DE introduces a differential mechanism for exploring the search space. Therefore, new larvae are generated by perturbing the population members using vector differences of individuals. Perturbations are introduced by applying where F determines the evolution factor weighting the perturbation amplitude) for each encoded parameter on a random basis. After this perturbation, final perturbed vector x is combined with an alternative coral in the reef following a classical 2-points crossover, as defined next. 3.
2Px: Classical 2-points crossover. The crossover operator is the most classical exploration mechanism in genetic and evolutionary algorithms [54]. It consists of coupling to individuals at random, choosing two points for the crossover, and interchanging the genetic material in-between both points. In the CRO-SL, one individual to be crossed is from the 2Px substrate, whereas the couple can be chosen from any part of the reef.

4.
MPx: Multi-points crossover. Similar to the 2-points crossover, but in this case a number k of crossover points are selected, and a binary template decides whether parts of the individuals are interchanged.
The reason of adapting the value of σ along the generations is to provide a stronger mutation in the beginning of the optimization, while fine-tuning with smaller displacements nearing the end. The mutated larva is thus calculated as where N i (0, 1) is a random number following the Gaussian distribution. 6.
Firefly Optimization (Fa): The Fa is a kind of swarm intelligence algorithm, first introduced by Yang in [55]. The Fa is based on the flashing patterns and behaviour of fireflies in nature. The pattern movement of a firefly i attracted to another (brighter) firefly j is calculated as follows: where β 0 stands for the attractiveness at distance r = 0. The specific Fa mutation implemented in the CRO-SL is a modified version of the algorithm known as Neighborhood Attraction Firefly Algorithm (NaFa) [56]. It has been implemented as follows. When a coral (solution) in the reef belongs to the NaFa substrate, it is updated by following Equation (12). All the parameters of the equation are tuned during the CRO-SL evolution. The corals in the NaFa substrate consider as swarm a neighborhood among all the other corals in the reef (not only the NaFa substrate). Thus, the corals in the NaFa substrate are updated taking into account some solutions from other substrates, as all the corals in the reef share the same objective function (brightness for the NaFa substrate). Note that the Fa algorithms has been applied to different structural problems in previous works [57,58], showing a good behaviour. 7.
Water Wave Optimization (WWo): The WWo [59] is a recently proposed meta-heuristic based on the phenomena of water waves, such as propagation, refraction, and breaking. Three different procedures are then defined in this algorithm: Wave propagation: at each generation of the algorithm, each wave x in the population is propagated, to generate another wave x in the following way: where rand(−1, 1) is a uniformly distributed random number within the range [−1,1], and L(d) is the length of the dth dimension of the search space (see in [59] for reference). Then, Wave refraction is simulated as where x * (d) is the best solution found so far, and N(·) is a Gaussian random number with mean µ and standard deviation σ. Finally, a wave breaking process is also simulated as where β w is the breaking coefficient.
The implementation of the algorithm as CRO-SL is straightforward, as any solution within the substrate is just applied the set of operators described above, with the algorithm's parameters described in [59].

Numerical Examples
Different design restrictions may be applied depending on the application type for which the arch is used. To approach a general perspective of the problem, two general cases have been analyzed: In the first one, the height of the arch H is fixed, while the parameters H , B, h, and t a change within a range. For the second scenario, the span of the arch B is constant, and the parameters H, H , h, and t a are set for respective ranges of values. The search for the best arch shape is based on the optimization of the functional given by Equation (10). Thus, the maximum absolute bending moment is reduced whilst the airspace of the arch is increased. Note that in this work, the reduction of the bending moment has been set as the most important objective as it implies some benefits, such as in terms of the construction cost (the amount of material used for its construction is reduced). Moreover, the thickness of the arch may be reduced, thus increasing its slenderness, which is an important aspect in many cases, as in leisure applications. If the importance of the airspace (or another parameter which can be included in the functional expression) is comparable with the bending moment, a normalization of the expressions should be addressed to guarantee that the airspace weight is clearly represented in the functional value.
There are two main restrictions in the problem. The first one is inherent to the numerical calculation. As nonlinear analysis is developed, some of the solutions may not converge, thus discarding the proposed solution. The second one is based on Equation (11). If that condition is not achieved either for linear and nonlinear analysis, the arch shape is not considered.
The results illustrate how the combined use of the curve parameterization and the CRO-SL leads to optimal solutions. In addition, the importance of considering the nonlinear effects in the arch design is also evaluated. The CRO-SL parameters are based on trial and error tests. They are shown in Table 1. The Young's modulus of the arch material is E = 27,664.49 MPa and the Poisson's ratio has been chosen to be ν = 0.2 (uncracked concrete) [60], with a density ρ h = 2400 kg/m 3 . The self-weight of the arch structure is considered in the analysis.

Case in Which the Height H Is Fixed
In this scenario, the height H and the depth D p are fixed, while the parameters B, H , h, and t a are tuned within a range (see Table 2). The functional weights used for the calculation are p m = 0.8 (N·m) −1 , p a = 0.2 m −2 . Thus, the reduction of the bending moment is more important than increasing the airspace. The results have been obtained for a single run simulation.

Type of Parameter Parameter Value [m]
Fixed parameters D p 10 H 6 Variable parameters The best solutions for both cases, linear and nonlinear, are shown in Table 3. The arch shape found by the CRO-SL is basically the same for the linear and the nonlinear cases (see Figure 4). However, the thickness of the structure (h) for the nonlinear based design is almost the double than for the linear design. This implies that, although the shape of the arch is the same for both scenarios, the increment of the bending stiffness of the arch has lead to a greater thickness of the structure to avoid instabilities issues. The optimal curve shape has been found between the parabola (t a = 0) and the ellipse (t a = 1), for both cases.   The parameters used to calculate the functional can be found in Table 3. As expected, the area for both arches is the same. However, the maximum bending moment for the nonlinear case is greater than in the linear case. Then, its functional value reaches a higher value. As it was aforementioned [18], instabilities issues can occur under service loads for many of the arch shapes proposed so far, and this example shows the importance of considering the nonlinear effects in the arch design.
The bending moment and the vertical deflections per unit arch length are shown in Figure 5a,b, respectively. The nonlinear analysis implies that instabilities issues may appear for small displacement increments. The CRO-SL has found a stable solution with low vertical deflection. sense, such increase not only augments the bending stiffness as well as bending stresses (as predictable, due to the self-weight increasing), but causes the characteristic redistribution of the bending moment already shown in previous works [23]. However, in this particular case, such redistribution is more relevant for B/4 ≤ x ≤ 3B/4, where the non-linear bending stresses decreases respect to the linear one, that produces a generalized reduction of the vertical displacements in the arch; this is due to the different stress redistribution that the arch experiments depending on the geometrical behaviour hypotheses: linear or nonlinear. Therefore, under certain design conditions (as the corresponding to this scenario), the second order optimization of submerged arches is not only feasible (as shown previous works, using other evolutionary strategies, such as the genetic algorithms [23,25]), but it is necessary in order to obtain a significant reduction in the structural vertical displacements. Regarding the computational performance of the CRO-SL in this problem, Figure  6a,b show the number of larvae (potential solutions to the problem) getting into the reef per substrate, over the iterations. Note that a new individual can enter in the population on a random place, if its associated fitness value is better than the occupied one, or if the place is empty. Due to the initial hollows in the reef, in the first iterations of the algorithm's run there is a higher probability to enter, independently of the quality of solutions, so this number is usually greater at the beginning of the algorithm's run. In Figure 6c,d the accumulative ratio of times that each substrate is producing the best larva is shown. It is remarkable that the best larva might not be better than the fittest coral of the reef. Regarding both figures, it can be pointed out that DE and NaFa substrates are good in In this first scenario, the coincidence of both shapes (linear and nonlinear solution) allows clarifying the effectiveness of the optimization procedure adopted in this work, as the structural stability is achieved mainly through the arch thickness increase. In this sense, such increase not only augments the bending stiffness as well as bending stresses (as predictable, due to the self-weight increasing), but causes the characteristic redistribution of the bending moment already shown in previous works [18]. However, in this particular case, such redistribution is more relevant for B/4 ≤ x ≤ 3B/4, where the nonlinear bending stresses decreases respect to the linear one, that produces a generalized reduction of the vertical displacements in the arch; this is due to the different stress redistribution that the arch experiments depending on the geometrical behaviour hypotheses: linear or nonlinear. Therefore, under certain design conditions (as the corresponding to this scenario), the second-order optimization of submerged arches is not only feasible (as shown previous works, using other evolutionary strategies, such as the genetic algorithms [18,20]), but it is necessary in order to obtain a significant reduction in the structural vertical displacements.
Regarding the computational performance of the CRO-SL in this problem, Figure 6a,b shows the number of larvae (potential solutions to the problem) getting into the reef per substrate, over the iterations. Note that a new individual can enter in the population on a random place, if its associated fitness value is better than the occupied one, or if the place is empty. Due to the initial hollows in the reef, in the first iterations of the algorithm's run there is a higher probability to enter, independently of the quality of solutions, so this number is usually greater at the beginning of the algorithm's run. In Figure 6c,d, the accumulative ratio of times that each substrate is producing the best larva is shown. It is remarkable that the best larva might not be better than the fittest coral of the reef. Regarding both figures, it can be pointed out that DE and NaFa substrates are good in search space exploration and exploitation. On the one side, they get into the reef a higher amount of larvae, which means, that they usually produce fitter individuals that the existed ones in the reef. It can be observed that GM does not show a good performance, as the solutions proposed by this method never fit in the reef. In addition, note that around 50% and 30% of the iterations these substrates, respectively, are producing the best larva in the linear model. Furthermore, in the nonlinear model DE substrate reaches the 80%. Other substrates such as HS, GM, MPx, or 2Px, are not producing fitted individuals because they are collaborating in a further exploration of the search space, instead of its exploitation. (d) Figure 6. Number of larvae entering to the reef per substrate for the linear (a) and nonlinear (b) cases, and ratio of times that each substrate generates the best larva in the linear (c) and nonlinear (d) cases when the height H is fixed. Figure 7 shows the fitness evolution for the linear and nonlinear cases. Note that the linear solution produces the lowest value of the functional, though it needs 37 iterations to be close of the final functional value.

Case in Which the Span B Is Fixed
In this scenario, the span of the arch B and the depth D p are fixed. The ranges of values for the parameters in which the CRO-SL seeks the optimal solution are defined in Table 4. The functional weights used for this experiment are p m = 0.7 (N·m) −1 and p a = 0.3 m −2 , then the importance of the arch airspace is slightly greater than for the last scenario.
The obtained functional values can be found in Table 5. As it can be observed, the area of the arch for the linear analysis based design is lower than for the nonlinear cases. However, the maximum bending moment for the nonlinear case is more than twice as great as the bending moment of the linear case.
The parameters of the optimal solution for both cases, linear and nonlinear, are shown in Table 5. The difference in the contour ellipse height H is very small, 0.50 m. However, the difference between the heights of the arches is more than 1.00 m, see Figure 8. Furthermore, the aim of reducing the vertical deflection in the nonlinear solution to avoid buckling phenomenon has resulted in the optimal solution having a greater thickness than in the linear solution. The shape parameter t a is very similar in both cases. Both solutions tend to be a parabola curve since t a are equal to 0.22 and 0.18 for the linear and nonlinear cases, respectively.  The weighted functional expression leads to an increment of the airspace of the arch and a reduction of the maximum bending moment. In this scenario, as the spam is fixed, an increment in the height implies an increment in area. Furthermore, this implies a reduction of the water column above the arch structure. However, this does not imply that the maximum bending moment is reduced, especially if the nonlinear behaviour is considered. In fact, the maximum bending moment of the arch when the nonlinear behaviour is considered is greater than in the linear case. This evidences the importance of the self-weight effect in the optimal design of a submerged arch in shallow waters [31], particularly in their second-order behaviour [18]. Moreover, if the optimal linear solution is simulated under nonlinear conditions, the stability is not reached. The bending moment produces an irreversible increment of the vertical deflection. Therefore, it is shown the importance of considering the nonlinear effects in the arch design.
As in the previous scenario, the bending moment and the vertical deflections per unit arch length are shown in Figure 9a,b, respectively. The nonlinear analysis also implies that instabilities issues may appear for small displacement increments. Thus, the CRO-SL has found a stable solution, with a low vertical deflection. In addition, the bending moment along the structure is greater for the nonlinear case. Nevertheless, once more, the redistribution effect in the nonlinear bending stresses produces a significant decrease in nonlinear vertical displacements (see Figure 9b); in this last case, the redistribution effect on the vertical displacements is more significance than in the first scenario due to the greater magnitude of the nonlinear bending moments. In addition, the results obtained by CRO-SL conjugates such significant reduction of the second-order vertical displacements with the increasing of the airspace enclosed by the arch, which allows a generalized improvement of the serviceability conditions in the underwater installation.
In a similar way than the previous case, the computational performance of the CRO-SL in this problem is discussed by means of the number of larvae entering to the reef per substrate and the ratio of times that each substrate generates the best larva, displayed in Figure 10. It is possible to see that the behaviour of the substrates is also very similar to the previous cases. DE and NaFa substrates have a more important weight in the search space exploitation and algorithm's evolution, while the rest of them are useful to further explore it. Figure 11 shows the best fitness evolution in the CRO-SL, for linear and nonlinear cases. In this case, the submerged arch design only considering linear behaviour in the calculation of the funicular arch always obtains a better value of the fitness functional (Equation (10)) than the nonlinear design.   (d) Figure 10. Number of larvae entering to the reef per substrate for the linear (a) and nonlinear (b) cases, and ratio of times that each substrate generates the best larva in the linear (c) and nonlinear (d) cases when the span B is fixed.

Conclusions
This work considers a recent parameterization of the funicular solution, in order to propose a novel approach for the optimal design of submerged arches in which geometrical nonlinear behaviour is considered. As novelty, such parameterization has been combined with a recently proposed ensemble multi-method meta-heuristic algorithm, the Coral Reefs Optimization with Substrate Layer algorithm (CRO-SL). A weighted functional has been proposed in order to guide the CRO-SL search, optimizing the arch serviceability conditions according to the following design criteria: maximum bending moment reduction and area under the arch increment. In addition, a vertical deflection restriction and the buckling load factor have been calculated to discard unrealistic solutions.
To evidence the effectiveness of CRO-SL algorithm in the optimal nonlinear design of submerged arches, two scenarios have been studied. In the first one, its height has been fixed and the span has been modified. While in the second scenario, the span of the arch has been fixed and its height has been tuned. Both cases have been analyzed considering only linear and then nonlinear effects. From the direct comparison between the optimized linear and nonlinear solutions, the main conclusion obtained is that the characteristic redistribution of the bending moments along the arch length may diminish in a significant way the nonlinear vertical displacements under certain design conditions. This last fact represents a great improvement of CRO-SL algorithm with respect to previous experiments performed using genetic algorithms. In addition, the results also show how the contribution of the CRO-SL substrates (search procedures) depends on the problem, illustrating that the proposed multi-method approach is able to adapt to the problem at hand.