Damage detection in a fiber reinforced polymers based tower of a floating offshore wind and tidal power platform using structural dynamic parameters

<!-- metadata commented in wiki content
==DAMAGE DETECTION IN A FIBER REINFORCED POLYMERS BASED TOWER OF A FLOATING OFFSHORE WIND AND TIDAL POWER PLATFORM USING STRUCTURAL DYNAMIC PARAMETERS==

'''P. Sierra<math>^{1*}</math>, X. Martinez<math>^{1,2}</math> y R. Chacón<math>^{3}</math> 



























'''
-->

==Abstract==

The development of offshore renewable energy structures is growing at a fast pace, and their maintenance will become critical in further years. Because steel corrosion on offshore structures is one of their highest maintenance cost, researchers are working to replace steel by Fibre Reinforced Polymers (FRP), due to their immunity to corrosion. This is done in the EU H2020 funded project FIBREGY, which also seeks improving maintenance costs by incorporating Structural Health Monitoring (SHM) strategies in the structure. A critical issue in SHM is the damage detection process, which is addressed in this paper. This study aims to characterize the dynamic behaviour of an offshore windmill FRP tower, and its evolution due the presence of structural damage. This analysis will be carried out using a FEM model of the structure, and the Serial-Parallel mixing theory to characterize the material performance. A sensitivity analysis of the changes in the dynamic response of the tower (frequencies and modal shapes) produced by structural damage is carried out. This analysis will provide the optimum design variables and objective functions required by a damage detection algorithm. The algorithm is optimization based, since it converts the damage detection problem into a mathematical one, and finds the best design variables to minimize the objective function. The algorithm is also model based, as it uses the structural model to obtain the objective function for the proposed design variables.Finally, this work presents the analysis of potential damage scenario, in order to evaluate if the damage is effectively detected.

==1 INTRODUCTION==

Offshore renewable energy is a cornerstone of the clean energy transition in the EU. It consists of many different sources that are abundant, natural and clean, obtaining energy from wind, wave and tidal. This can be harnessed by modern technologies without emitting any greenhouse gases, (European Commission 2020). There is no doubt that the offshore renewable energy exploitation has a large potential to grow. But, the open sea is a very aggressive environment, which may largely affect the maintenance costs of the installations and therefore the overall cost of offshore energy generation. The owners of offshore assets are aware of that and are paying a steep price. A massive amount of steel goes into those assets, and all this metal is subject to degradation, which explains why corrosion accounts for approximately 60% of offshore maintenance costs. In order to address this problem, the main objective of the EU H2020 952966 project FIBREGY, (EU 2023), is enabling the extensive use of FRP materials in the structure of the next generation of large Offshore Wind and Tidal Power (OWTP) platforms, due to the convenient immunity to corrosion and superior fatigue performance of these materials. Another innovative aspect that will be considered is the use of multifunctional FRP materials in the OWTP structure, by embedding sensors in the material for continuous Structural Health Monitoring (SHM).

During FIBREGY project, a scale (1:50) of an FRP tower with an embedded SHM system will be tested. In this paper, a dynamic numerical analysis of this scaled structure is carried out. A Finite Element Model (FEM) is made and analysed to obtain its dynamic behaviour, characterized by frequencies and modal shapes. Also, a model updating and optimization based damage detection method is proposed for this structure and tested using a numerical damaged scenario.

==2 Numerical Model==

<div id='img-1'>
===2.1 Geometry===

The structure under study is a cone with the following geometric characteristics, showed in Fig. 1. The bottom face of the tower is fully fixed (displacement in x, y and z).

* Height: <math display="inline">h = 1200 \  mm</math>
* Bottom external diameter: <math display="inline">\Phi _b = 117 \ mm</math>
* Top external diameter: <math display="inline">\Phi _t = 84 \ mm</math>
* Thickness: <math display="inline">t = 4 \ mm</math>

The tower is divided into <math display="inline">10</math> sectors in its height (x axis) and in <math display="inline">6</math> layers evenly distributed in the thickness, in correspondence with the laminate that forms it, resulting in <math display="inline">60</math> volumes.
</div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Sierra_755961050-cono.png|500px|Geometry of the analysed tower, dimensions in mm.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Geometry of the analysed tower, dimensions in mm.
|}

The FEM design is done with GiD,(CIMNE 2020), a pre and post processor for numerical simulations in science and engineering. The structure is discretized in <math display="inline">54\,096</math> nodes and <math display="inline">46\,080</math> hexahedra linear elements and is shown in Fig. 2.
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Sierra_755961050-modelGid.png|191px|Numerical model discretization]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 2:''' Numerical model discretization
|}

===2.2 Material===

The FRP used in the scaled tower is composed of a vinylester based resin with glass fiber reinforcement, with a <math display="inline">20 %</math> on volume participation. The main characteristics are listed in Tab 1. The composite has <math display="inline">6</math> layers of unidirectional fibres in the following orientations, angles from the axis of the tower, see Fig. 1. From external to internal layer: <math display="inline">0^{\circ }   /  55^{\circ }   /  0^{\circ }   /  {-55}^{\circ }   /  0^{\circ }   /  90^{\circ } </math>.
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Materials properties
|- style="border-top: 2px solid;"
| colspan='1' | Material
| colspan='1' | Young's Modulus
| colspan='1' | Poisson
| colspan='1' | Density
|-
|      
| colspan='1' | E (MPa)
| colspan='1' | <math>\nu </math>
| colspan='1' | <math>\rho  \left(\frac{t}{m^{3}} \right)</math> 
|- style="border-top: 2px solid;"
|           Resin 
| <math>3\,350</math> 
| <math>0.260</math> 
| <math>1.12</math>  
|- style="border-bottom: 2px solid;"
|      Fibre 
| <math>73\,100</math> 
| <math>0.238</math> 
| <math>2.62</math>  

|}

===2.3 Numerical Analysis===

A dynamic analysis of the model is done using an academic multi-physics FEM code, FEMUSS, developed by the International Centre for Numerical Methods in Engineering (CIMNE) and the Polytechnical University of Catalunya (UPC). The different orthotropic composite layers, are treated using the serial/parallel mixing theory, (Martinez, Oller, Rastellini, & Barbat 2008).

===2.4 Damaged Case===

A reference damaged tower is proposed to observe changes in the dynamic behaviour of the model and to use its modal frequencies as input data in the damage detection method. Numerically, the damaged is represented by a reduction of the Young's Modulus of both materials, fibre and resin, in layers <math display="inline">1</math>, <math display="inline">3</math>, and <math display="inline">5</math> (zero degree fibre orientation layers of the laminate), to <math display="inline">100  MPa</math> in the sector <math display="inline">5</math> of the tower, as shown in Fig. 1.

==3 Dynamic behaviour==

The embedded SHM system that is going to be used in tests allows up to <math display="inline">1\,000 \  Hz</math> of sampling rate. In consequence, the analysis is limited to <math display="inline">500 \ Hz</math>, the half of the sampling rate value, corresponding to the Nyquist frequency.

In this range, <math display="inline">7</math> modes were found. The frequency values obtained in the analysis of the undamaged scenario and the described damaged reference case, are listed and compared in Tab. 2. Due to the rotational symmetry of the structure, modes 1-2, 3-4 and 6-7 have almost the same values between them, but not equal because the asymmetry in the fibre's orientation. These are flexural modes, while mode 5 is expansive, as is shown in Figs. 3 and 4. There, dark blue colour represents no displacement, red biggest displacement.
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-2'></span>Table. 2 Frequencies, in <math>[Hz]</math>
|- style="border-top: 2px solid;"
|         Mode 
|  Main 
|  Undamaged 
|  Damaged 
| <math>\Delta _{rel}</math> 
|-
|     
|  Direction 
|  Case 
|  Case 
| <math>[%]</math> 
|- style="border-top: 2px solid;"
|         1 
|  Z 
|  43.16 
|  37.97 
|  12.02 
|-
|     2 
|  Y 
|  43.34 
|  38.07 
|  12.14 
|-
|     3 
|  Z 
|  209.12 
|  177.76 
|  14.99 
|-
|     4 
|  Y 
|  209.82 
|  178.49 
|  14.93 
|-
|     5 
|  - 
|  292.47 
|  279.26 
|  4.51 
|-
|     6 
|  Z 
|  488.72 
|  470.05 
|  3.81 
|- style="border-bottom: 2px solid;"
|     7 
|  Y 
|  489.78 
|  470.93 
|  3.84 

|}

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Sierra_755961050-SD.png|600px|Modal shapes of selected modes in undamaged tower.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Modal shapes of selected modes in undamaged tower.
|}

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Sierra_755961050-D1.png|600px|Modal shapes of selected modes in reference damaged tower.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Modal shapes of selected modes in reference damaged tower.
|}

The first <math display="inline">4</math> modes are more affected for the damage, with changes in frequencies above <math display="inline">10 %</math>, than the last 3 modes analysed. That is because the damage is located close to a nodal point of the modal shape in modes <math display="inline">6</math> and <math display="inline">7</math> and minimally affects the pressure stiffness, observed in mode <math display="inline">5</math>.

From this analysis is concluded that the frequencies are sensitive to the damage and its location is related to the more affected modes. So, this parameter can be used in a damage detection method.

==4 Damage detection==

In this paper, structural damage is detected using changes in the natural frequencies. The method is a FEM updating based, that consist in modifying the mass, stiffness and damping parameters of the numerical model to match the behaviour in the model as in the test. To find it, an optimization algorithm is used, described in section 4.1. The problem is converted to a minimization of an objective function <math display="inline">f\mathbf{(x)}</math>, where <math display="inline">\mathbf{(x)}</math> are the design variables, in this case the Young modulus of different volumes of the FEM. The basic scheme of the method is illustrated in Fig. 5, where <math display="inline">\left\lbrace \mathbf{f} \right\rbrace </math> is the vector of natural frequencies.
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Sierra_755961050-DDSchemehalf.png|452px|Damage detection method scheme]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Damage detection method scheme
|}

Several interfaces have been developed to couple the optimization algorithm with the structural analysis software. For example, to convert the proposed numerical design variables <math display="inline">\mathbf{x}</math> by the optimizer in an valid analysable FEM by FEMUSS, other to read the results of the dynamic analysis and evaluate the objective function.

===4.1 Optimization===

The general optimization problem can be formulated as solving the function:

<span id="eq-1"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>minimize \quad f\mathbf{(x)}, \quad \mathbf{x} = [x_1, x_2, \dots , x_n]  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

Where <math display="inline">f\mathbf{(x)}</math> is the objective function, <math display="inline">n</math> is the dimension of the problem and <math display="inline">(x_1, x_2, \dots , x_n)</math> are the design variables. These may be subject to constraints in order to preserve physical sense.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>lb_i \leq x_i \leq ub_i \quad for \quad  i = 1,2, \dots , n </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

There <math display="inline">lb</math> and <math display="inline">ub</math> represents the lower and upper bound limit, respectively. The global optimization is referred to find the set of <math display="inline">\mathbf{x}</math> that minimize <math display="inline">f\mathbf{(x)}</math> over the entire feasible region. The problem is harder as the dimension <math display="inline">n</math> increases and, in cases when the behaviour of the objective function is unknown, is not possible to be certain whether one has found the true global optimum.

Several global optimization algorithms were developed that are more or less efficient, depending on the behaviour of the problem. In this paper, the Multi Level Single Linkage (MLSL) algorithm for global optimization from the NLopt library, (Johnson 2020), is used for the damage detection problem. It consist of a sequence of local optimizations starting from random points, proposed by Rinnooy Kan & Timmer (1987a), (1987b). An alternative implementation is applied. This uses a Sobol low discrepancy sequence as starting points, instead of random numbers, providing more accurate results, practically and theoretically guaranteed. This application was presented by Kucherenko & Sytsko (2005). The local optimization is carried out by the derivative - free and bound - constrained BOBYQUA local optimizer. This algorithm was presented by Powell (2009). It consist of a quadratic approximation of the objective function, constructed iteratively.

===4.2 Design variables===

In this numerical case, the design variables represent changes in the Young's modulus of the materials that are in the composite, fibre and resin. These variables can take any value between <math display="inline">0</math> and <math display="inline">1</math>.

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x} = [x_1, x_2, \dots , x_n], \quad x_i \in \left[0, 1 \right] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|}

The zero value represents no damage, so the <math display="inline">\mathbf{E}</math> remains unalterable, with the values listed in Tab. 1. A value of one, represents fully damaged, to avoid convergence problems in the numerical analysis, the <math display="inline">\mathbf{E}</math> is reduced to <math display="inline">100  MPa</math> for both, the resin and the fibre. Intermediate values are linearly interpolated, as shown in Tab. 3, with <math display="inline">x=0.5</math>.
{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-3'></span>Table. 3 Young's Modulus in function of design variable
|- style="border-top: 2px solid;"
| <math>x_i</math> 
|  Fibre 
|  Resin  
|-
|     
| <math>\mathbf{E}</math> [MPa] 
| <math>\mathbf{E}</math> [MPa] 
|- style="border-top: 2px solid;"
|         0 
|  73\,100 
|  3\,350  
|-
|     0.5 
|  36\,600 
|  1\,725 
|- style="border-bottom: 2px solid;"
|     1 
|  100 
|  100 

|}

At first, a simpler scenario is evaluated, where are only <math display="inline">10</math> design variables. Each one of these represents the <math display="inline">3</math> layers with <math display="inline">0^{\circ }</math> fibre orientation of <math display="inline">1</math> zone in the structure. Divided as shown in Fig 1. If, for example, the design variable <math display="inline">x_7 = 1</math>, it implies that the <math display="inline">3</math> layers <math display="inline">0^{\circ }</math> in zone <math display="inline">7</math> are damaged, so it's Young's modulus is equal to <math display="inline">100  MPa</math>. The damaged case, described in Section 2.4, is represented by the vector of Eq. 4.<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x} = [0.0,  0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  0.0,  0.0,  0.0] </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}

Another scenario is considered, with <math display="inline">60</math> design variables. These are correlated with the material properties of each one of the volumes that the structure is divided. So, first <math display="inline">10</math> variables represent the external layer, of <math display="inline">0^{\circ }</math> fibre orientation, in the different zones. From <math display="inline">11</math> to <math display="inline">20</math>, the second layer, of <math display="inline">55^{\circ }</math> fibre orientation, and so on for the following variables. The damaged case, described in Section 2.4, is represented by all values equal to <math display="inline">0.0</math>, except the listed in Eq. 5.<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>x_5 = 1.0, \quad x_{25} = 1.0, \quad x_{45} = 1.0 </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

===4.3 Objective function===

The aim of the optimization process is to obtain the set of variables which minimize the objective function, Eq. 1. The success of the process depends on the correct behaviour of this function. This function must compare measurable values of a damaged structure (reference - ref), with the results of a proposed by the optimizer model of the structure (evaluation - ev).
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Sierra_755961050-freq.res.png|340px|Results optimization with 10 design variables]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Results optimization with <math>10</math> design variables
|}

In this paper Eigenfrequencies Method is used. The natural frequencies of the structure are very sensitive to the damage and it is easy to implement a monitoring system, thus is widely used. The variations in the frequencies and which mode is more or less affected depends on the position and severity of the damage. But it represents the global behavior of the structure, since doesn't use spatial information. The proposed objective function is Eq. 6.<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> f_{(\mathbf{x})} = \left\|\mathbf{f}_{ref} - \mathbf{f}_{ev} \right\|  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{f}_{ref} = [f_{ref,1}, f_{ref,2}, \dots , f_{ref,7}] = [37.97, 38.07, \dots , 470.93]  Hz </math>
|-
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{f}_{ev} = [f_{ev,1}, f_{ev,2}, \dots , f_{ev,7}]  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

Where <math display="inline">\mathbf{f}_{ref}</math> and <math display="inline">\mathbf{f}_{ev}</math> are vectors of the first <math display="inline">7</math> natural frequencies of the reference damaged structure and the results of the evaluation, respectively. This function is <math display="inline">0</math> and is the global minimum when the vectors are equal, so the structure proposed and the damaged case matches.

===4.4 Results===

In the first scenario, with <math display="inline">10</math> design variables as explained before, the optimization process is set up with a limit of <math display="inline">1\,000</math> evaluations. In that range, the algorithm found a minimum in the evaluation <math display="inline">255</math>, with an objective function value of <math display="inline">0.088</math>. The evolution of the objective function value is shown in Fig. 6.

The design variables of that minima are listed in Eq. 9 and plotted in the histogram of Fig. 6. Except for minimal differences in zones <math display="inline">2</math> and <math display="inline">10</math> (less than <math display="inline">2  %</math>), the obtained structure and the damaged case are equal, so the damage is found.<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{x}_{res}  = [0.000, 0.020, 0.000, 0.000, 0.998, 0.000, 0.000, 0.000, 0.000, 0,019]</math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|-
|}
|}

In the <math display="inline">60</math> design variables scenario, the optimization is limited to <math display="inline">3\,500</math> evaluations. In this scenario the minimum value of the objective function was found in the evaluation <math display="inline">2\,801</math>, with an objective function value of <math display="inline">1.148</math>. The evolution of this function and the histogram representing the best sample found are shown in Fig. 7.
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_Sierra_755961050-freq60.res.png|340px|Results optimization with 60 design variables]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 7:''' Results optimization with <math>60</math> design variables
|}

The damage in zones <math display="inline">5</math>, <math display="inline">25</math> and <math display="inline">45</math> of the reference damaged case are well captured, but, there are also significant damage values in zones <math display="inline">60</math> and <math display="inline">20</math>, <math display="inline">0.956</math> and <math display="inline">0.175</math> respectively. The error in the remaining values is negligible, less than <math display="inline">0.051</math>. To obtain better results, the optimization algorithm needs more evaluations, but the obtained results can be considered acceptable, considering the difficulty of the problem.

==5 Conclusions==

In this paper a numerical model of offshore windmill FRP tower is presented and its dynamic behaviour is characterized by its natural frequencies and modal shapes. A reference damaged case of the tower is analysed. The sensitivity to damage of the natural vibration frequencies of the structure has been verified.

Also, a damage detection method is presented, based on optimization algorithms and model updating techniques. This method includes detection, localization and extension of the damage.  Based on the comparison of the natural frequencies, <math display="inline">2</math> cases are evaluated with different space discretization. At first one, with the structure divided into <math display="inline">10</math> zones to apply damage, while in the second scenario into <math display="inline">60</math> zones. In both, satisfactory results were obtained. It was demonstrated that the proposed damage detection method could be used in SHM for this kind of structures.

==6 Acknowledgments==

The founding of FIBREGY, “Development, engineering, production and life-cycle management of improved FIBRE-based material solutions for structure and functional components of large offshore wind enerGY and tidal power platform” an H2020 project under agreement 952966 is gratefully acknowledged.

===BIBLIOGRAPHY===
<nowiki>CIMNE (2020, Sep. 30, 2020). GiD. https://www.gidhome.com.</nowiki>

<nowiki>EU (2021-2023). H2020 952966 FIBREGY - Development, engineering, production and life-cycle management of improved FIBRE-based material solutions for structure and functional components of large offshore wind enerGY and tidal power platform. https://fibregy.eu/.</nowiki>

European Commission (2020). An eu strategy to harness the potential of offshore renewable energy for a climate neutral future. Technical report, Commision to the European Parliament, the Council, the European Economic and Social Committee and the Committee of the regiions. https://eur-lex.europa.eu/legalcontent/EN/TXT/?uri=COM:2020:741:FIN&qid=1605792629666.

<nowiki>Johnson, S. G. (2020, Sep. 30, 2020). NLopt - the NonLinear - optimization package. http://github.com/stevengj/nlopt.</nowiki>

Kucherenko, S. & Y. Sytsko (2005). Application of Deterministic Low-Discrepancy Sequences in Global Optimization. Computational Optimization and Applications 30, 297–318. 

Martinez, X., S. Oller, F. Rastellini, & A. Barbat (2008). A numerical procedure simulating RC structures reinforced with FRP using the serial/parallel mixing theory. Computers and

Structures 86, 1604–1618.

Powell, M. (2009). The BOBYQA Algorithm for Bound Constrained Optimization without Derivatives. Technical report, Department of Applied Mathematics and Theoretical Physics.

Rinnooy Kan, A. H. G. & G. T. Timmer (1987a). Stochastic global optimization methods part I: Clustering methods. Mathematical Programming 39, 27–56.

Rinnooy Kan, A. H. G. & G. T. Timmer (1987b). Stochastic global optimization methods part II: Multi level methods. Mathematical Programming 39, 57–78.