#### Abstract

In some mechanical models, the tensile armors of bent flexible pipes are treated as geodesics on a torus and, based on this hypothesis, the curvatures of these curves are calculated to obtain the acting stresses. However, a closed-form solution of the geodesic differential equations is not possible, which imposes difficulties on determining these curvatures. This work, therefore, proposes two alternative solutions to the nonlinear geodesic differential equations. The first relies on an artificial neural network (ANN) and the second is obtained by symbolic regression (SR). Both employ data from the numerical solution of the geodesic differential equations and showed good correlation with the complete dataset. Nevertheless, when tested against new data, the SR equations led to results almost equal to those obtained with the numerical solution of the differential equations and to null geodesic curvature. Despite also agreeing well with the numerical solution, the ANN indicates nonnull geodesic curvatures. Moreover, when compared to equations often employed in the design of flexible pipes, the SR equations may indicate different results, which can impact, for example, the fatigue or the instability analysis of the tensile armors of these pipes.

#### 1. Introduction

Unbonded flexible pipes (or, simply, flexible pipes), such as the one presented in Figure 1, are largely employed in offshore oil and gas exploitation. These multilayered composite structures combine high axial strength and stiffness with low bending stiffness, resulting in highly compliant pipes. Each layer of these pipes has a specific function and may be either metallic or polymeric. Whereas the metallic layers provide structural resistance, the polymeric layers are used to seal the pipe and/or to mitigate wear and friction between layers. In a typical flexible pipe, three different metallic layers are found [1]:* inner carcass*, which is made from profiled steel strips wound at angles close to 90 degrees with respect to the pipe axis that mainly resists radial inward forces;* pressure armors*, which are usually Z-shaped steel wires also wound at angles close to 90 degrees with the main function of supporting the system internal pressure and also radial inward forces; and* tensile armors*, which are constituted of several approximately rectangular steel wires laid in two or four layers, cross-wound at angles between 20 degrees and 55 degrees, which resist tension, torque, and pressure end cap effects.

The structural responses of these pipes to bending and axial loads have been extensively studied since the mid-seventies. Several analytical and numerical models, as reviewed and summarized by Tang et al. [2], are available in the public literature, but some aspects of their response remain challenging and the computation of the stresses in the tensile armors due to bending loads, as stated by Dai et al. [3] and Zhou and Vaz [4], is among them.

The bending response of a flexible pipe depends on the curvature imposed on the pipe. Several authors [5–7] describe this response as a stick-slip mechanism (Figure 2), which is activated by the contact pressures between layers due to the applied axisymmetric loads.

For small curvatures, friction between the armors and the adjacent layers prevent their relative slip. Hence, axial forces are induced in the armors and these forces are opposed by friction forces with the same magnitude. This leads to a linear bending moment versus curvature relationship with very high tangent stiffness. This tangent stiffness is called no-slip bending stiffness, .

As curvature increases, interlayer friction is overcome and progressively allows the relative movement of the layers. This slippage reduces the tension increase in the extrados of the pipe and compression decreases in its intrados, thereby reducing the tangent stiffness of the pipe. This stiffness keeps decreasing until friction forces are fully overcome and the tensile armor wires are free to slip. At this point, the tangent stiffness reaches a limit value, much lower than the no-slip one. This lower limit is called full-slip bending stiffness, , and is the value usually provided by the manufacturers of flexible pipes. Moreover, the curvature at which the interlayer friction is fully overcome is called critical curvature, , and the associated bending moment is named internal friction moment, . Figure 2 schematically shows this mechanism. If no friction is considered between layers, the bending response of the pipe is linear with bending stiffness equal to the full-slip value.

In the prediction of the stresses in the tensile armors induced by the bending of the pipe, these armors are commonly treated as naturally curved beams [6–9]. A traditional mathematical model used to describe the mechanics of a naturally curved rod was systematized by Love [10]. In there, a system of six differential equations is stated (see (1)), which holds the equilibrium of such rods concerning its local coordinate system.where and represent, respectively, the internal forces and moments of the rod in its tangential , normal , and binormal directions; the external distributed loads and moments are given by and in this order; and the natural curvature and torsion of the rod are given by and , respectively. It is important to note that the curvature in the binormal direction, for surface curves, is usually called the geodesic curvature. The comma after a symbol , in this work, indicates the derivative of the given variable with respect to the indicated one; for instance, is equivalent to .

The stresses in the tensile armors are induced by the internal forces and moments that result from the variation of the natural curvatures and torsion caused by the bending of the pipe. Therefore, the evaluation of the internal forces and moments is not possible unless these geometric components are previously known, as can be depicted from (1).

Several analytical models devoted to the bending analysis of flexible pipes consider that the curvatures and torsion induced in the tensile armors depend on the path followed by these armors during loading. Typically, two assumptions are made with respect to this path [8]: the geodesic and the loxodromic. These curves correspond to physical limit paths of the stressed tensile armors. The geodesic is the curve associated with the shortest distance between two points at a cylindrical supporting surface, while the loxodromic is the curve that makes the same angle with every meridian of this surface. The geodesic has no transverse (geodesic) curvature and both longitudinal and transverse slips relative to the loxodromic are needed to reach this path, as indicated in Figure 3.

In some analytical models proposed to predict the bending response of flexible pipes, such as those from Féret and Bournazel [9], Estrier [11], Out and Von Morgen [12], and Ramos Jr. and Pesce [13], the geodesic path is related to conditions in which no interlayer friction is considered or is disregarded, since the friction forces would restrain the slip of the armors. On the other hand, if friction forces are present, as noted by Sævik [5], the transverse slip towards the geodesic path is restrained and, therefore, the curvatures and torsion in the armor are somewhere between the loxodromic and the geodesic paths [14]. If no relative slip is assumed between the armor and the supporting surface, the loxodromic may be directly employed, as described in Witz and Tan (1992), Estrier [11], Sævik [5], or Larsen et al. [14], but, since at least the longitudinal slip is unavoidable when the critical curvature is reached, the loxodromic curvatures and torsion need to be modified even if no transverse slip occurs. In this last case, Sævik [15] and Larsen et al. [14] propose expressions to calculate these quantities by neglecting the longitudinal friction.

More general models are proposed by Féret et al. [6], Leroy and Estrier [7], and Østergaard [16]. Féret et al. [6] proposed a function dependent on an empirical term allowing the consideration of any curve in between the two limit paths. Later, Leroy and Estrier [7] followed the same approach but disregarding this empirical parameter and attributing the function to be dependent on the interlayer friction and the variation of the transverse displacement. A similar procedure was also adopted by Østergaard [16], where six differential equations based on angular coordinates are proposed to describe the mechanics of a frictionless armor sliding on a torus surface.

The use of these more generic models, however, has some difficulties regarding, for example, the stability of the solution, implementation procedures, and required approximations. Numerical models, usually based on the finite element method, are also employed in the bending analysis of flexible pipes, such as those proposed by Sævik [5, 15]. Nevertheless, these models usually demand high computational resources that may hamper, for instance, their use in the fatigue analysis of flexible pipes, as this type of analysis typically requires thousands of bending mechanical analyses [17].

In this scenario, the use of analytical models that rely on the choice of limit curves is quite attractive, but while closed-form equations for the curvatures and torsion of the loxodromic curve can be easily obtained from its differential equation, the geodesic quantities cannot be directly calculated. The geodesic is defined by stiff nonlinear differential equations of second order, which are usually solved by numerical methods [12]. Aiming at simplifying the use of this hypothesis in the analysis of the bending response of flexible pipes, this work proposes closed-form solutions of the geodesic differential equations obtained by symbolic regression (SR) over a dataset generated from the numerical solution of these equations considering different conditions. Moreover, this dataset is used to train an artificial neural network (ANN), whose predictions are then compared to the SR equations.

Therefore, initially, a brief description of the geometrical relations involved in the derivation of the geodesic differential equations is presented. After that, these equations are derived followed by a description of the numerical method (Runge-Kutta Lobatto IIIa) used to solve them. The ANN and the SR approaches are then described and the results of a case study are presented and discussed. Finally, the main conclusions of this work are stated.

#### 2. Geometric Relations

The geometric representation of the tensile armors as surface curves is quite common. The torus is a useful surface for this purpose, as it stands for a generic bent cylinder under constant curvature, thus representing a wide range of real structures. Figure 4 illustrates the torus and the angles necessary to its parametrization and its parametric equation is given by the following equation:

An infinite number of curves can be defined on the torus surface and an orthogonal triad must be established to orient them. In this work, the Darboux-Ribaucourt axes are employed, where is the unit tangent of the curve and and are the normal and the binormal vectors of the surface [18], as indicated in Figure 4. The variation of this triad with respect to the curve arc-length giveswhere and ( in (1)) are, respectively, the normal and geodesic curvatures and is the geodesic torsion of the curve. Generically, the curvatures and torsion can be directly determined by differential geometry formulas, named after Leonhard Euler (see (4)), Joseph Liouville (see (5)), and Sophie Germain (see (6)) [19]. where is the curve lay angle with respect to torus longitudinal direction (see in Figure 4); and are, in this order, the major and minor principal curvatures of the surface and and are the geodesic curvature for the -constant ( in Figure 4) and -constant ( in Figure 4) curves, respectively. As the angular coordinates and are orthogonal, the following formulas, which are based on the fundamental form of the torus, express these components [18]:where , , and are the first fundamental forms given by Nutbourne and Martin [18]:

Physically, and may be interpreted, for an orthogonal surface, such as the torus, as the square of the* speed* along the first and second parameters, respectively [20].

#### 3. Geodesic Equations

The evaluation of the curvature components in the tensile armors is of fundamental importance to predict the acting forces and moments on them. These curvature components, as discussed in Section 1, are usually calculated by considering two limit curves: the loxodromic and the geodesic. This work focuses on the geodesic curve.

Geodesic curves have their normal coincident with the surface normal. A more precise and mathematic accurate description of a geodesic is given by do Carmo [19], for whom geodesics are constant speed curves with vanishing geodesic curvature. They are also the shortest curves that connect two points in a surface being analogous to a straight line in a plane [18]. Many different parametrizations lead to many different differential equations of the geodesic of a surface: the first (and more widespread) was given by Leonhard Euler and Joseph Louis Lagrange [21], which consists of obtaining the minimization of the arc-length functional, , of a given surface . This minimization problem recalls the classical Euler-Lagrange equation, which is written as

This problem was also studied by Sævik [8] and Out and Von Morgen [12], stating that these differential equations cannot be solved analytically. A less straightforward, but more systematic way, to obtain the differential equations of a geodesic of a generic surface is given by the set of differential equations presented in Nutbourne and Martin [18]:

In surfaces parameterized by an orthogonal coordinate system, such as the torus, the first fundamental form is zero; thus the Christoffel symbols , which relate the first fundamental forms of a surface and their derivatives, are given by the following equations:

Substituting (10) and (12) into (11), the classical set of differential equations for the geodesic of a torus can be written as

However, do Carmo [19], using Liouville’s formula for the binormal curvature (see (14)), presents an alternative differential equation for the geodesic. The binormal curvature is initially written as

As previously stated, geodesic curves have null geodesic, or binormal, curvature, . Thus, by equating (14) to zero, a first-order differential equation written in terms of the lay angle can be written. Aiming at establishing the geodesic equations, a relation with any other torus angular parameter is still necessary. This relation may be obtained from the tangent definition, as can be seen in (15). The tangent of a unit speed curve provides the following relation:

Using this equation, it is possible to write the following equation:where is the norm of vector .

Thus, distinct differential equations of a geodesic may be obtained in terms of two functions of the arc-length , the lay angle , and the radial angular parameter , as shown in (17).

Østergaard [16] showed the equivalency between the two sets of differential equations, given by (13) and (17). In order to obtain the geodesic parameters (such as normal curvature and torsion), any of these sets of differential equations need to be solved numerically. In this work, the sixth-order Runge-Kutta (RK) Lobatto IIIa method is chosen. However, in the design of flexible pipes, a high number of local mechanical analyses are demanded and, consequently, direct approximations to the curvatures and torsion of the geodesics are desirable to define the equilibrium equations of the tensile armors (see (1)), thus minimizing the computational effort.

Therefore, here, two machine learning techniques are employed to obtain these approximations, that is, ANNs and SRs, based on Genetic Programming (GP) model. The results of the accurate numerical solution serve as a database to train the ANNs and as input for the SR approach. Next, these approaches are discussed.

#### 4. Runge-Kutta (RK) Lobatto IIIa Method

Lobatto IIIa method derives from the family of implicit RK methods. The method was proposed to overcome some implementation difficulties, stability issues, and convergence properties of traditional explicit and implicit Runge-Kutta methods. Later, other families of RK methods, such as Radau IIa, were developed. They present better convergence and stability than Lobatto IIIa method, which is A-stable but not algebraically or L-stable [22]. Despite the apparent drawbacks, Lobatto IIIa gives quite accurate solutions for stiff differential equations with global error lower than other methods, even the Radau IIa [22].

Consider the stiff ordinary differential equations and their boundary condition; for example,

The formulas that define general RK methods of stages are given bywhere is the number of stages , is an approximation of the function that represents the integral of the analyzed differential equation (main numerical RK approximation), and are the so-called internal stages (an approximation of the solution at ). The coefficients , , and define the chosen RK method. These coefficients are usually presented in a* Butcher tableau* (see (20)), which relates the coefficients , , and to any of the RK family methods, as follows:

It is interesting to note that when Lobatto IIIa method is considered in its simpler form , the well-known implicit trapezoidal rule is obtained. Lobatto RK family of methods (IIIa, IIIb, IIIc, and ) are characterized by having the same coefficients and , known as* weight* and* node* coefficients, respectively. Moreover, for all these family of methods, and . The coefficients for the Lobatto IIIa method are given by [23]where is the Legendre polynomial [23].

The sixth-order Lobatto IIIa method is obtained considering = 4, , and for . The* Butcher tableau* for this method is, therefore, given as follows:

In this work, the differential equations given by (17) were solved according to the algorithm proposed by Pinto et al. [24], where the RK Lobatto IIIa method is generalized for systems with any number of first-order differential equations. The applied boundary conditions are equal to and .

#### 5. Artificial Neural Networks (ANNs)

Artificial neural networks (ANNs) constitute a form of mathematical processing of information with architecture and operation similar to the biological neural networks. The learning ability is the remarkable feature of the ANNs, and, in practice, this ability applies pattern recognition, classification, clustering, optimization, and function approximation, among others [25].

In this work, ANNs are employed as function approximators. These ANNs define relationships between the bent flexible pipe tensile armor geometrical parameters (network inputs) and the geodesic curve parameters, that is, the radial angle and the lay angle given as functions of the armor arc-length (network outputs). However, the ANN needs sample data of the inputs and outputs to define the relation between them (ANN training) and, here, this data is provided by the numerical integration of the geodesic curve differential equations (see (17)) for some training cases.

The main goal is training an ANN and checking whether its extrapolated results are valid or not. This is one of the most critical issues when it comes to ANN applications. In order to do so, the radial and the lay angles along the armor arc-length for a given pipe configuration, which is not present in the training examples, are compared to the numerical solution.

The most common ANN architecture employs a single hidden layer between the input and output layers. This configuration is outlined in Figure 5. The first layer reads the input values of the network, while the output layer delivers the responses of the network. Both hidden and output layers are composed of neurons (numerical functions). Each connection between elements in neighboring layers contains a weight (synapse weighting). Linked to the hidden and output layers, there are also bias parameters with unitary values.

In a typical ANN, firstly, a neuron in the hidden layer receives an input given bywhere is the total number of network inputs, is the bias weight of the input layer, is the weight between the input and the respective neuron in the hidden layer (, is the total number of hidden layer neurons), and is the network input. Each neuron returns an output written aswhere is the so-called activation or transfer function. The activation functions are usually sigmoid functions, like the hyperbolic tangent or the logistic functions, as shown in Figure 6. Both functions are monotonic, differentiable over their domain, and limited. These features are required to introduce a nonlinear mapping in the network [26]. Due to the nature of the limit values, the logistic function is usually employed in pattern recognition applications, while the hyperbolic tangent is employed to function approximation and time series prediction. The weights control the steepness and the offset function. As the use of ANN in this work is related to function fitting, the hyperbolic tangent is chosen for the hidden neurons of the networks.

For an ANN with a single output , this parameter is determined bywhere is the activation function of the neuron in the output layer, usually a linear function. The remaining parameters are defined similarly to those in the hidden layer expressed by (23).

For an ANN, considering the input as a group of parameters values, the relation between these values and the network output is expressed bywhere is a vector containing the parameters values of the input and is a matrix containing all ANN connecting weights.

The so-called network training consists in the search of a set of weights which leads to a good fitting of the approximated function. In this manner, a sample set of inputs and outputs is needed to adjust the network weights by an optimization method. The optimization problem is expressed by the minimization of an error function, usually the mean square error , defined bywhere is the length of the training sample and is the desired output sample.

An iterative optimization procedure is employed to obtain optimum values for the weights . In this work, the Levenberg-Marquardt algorithm is employed in the networking training. This algorithm is a variation of Newton’s method designed to minimize functions that are sums of squares of other nonlinear functions [27] as the error function (see (27)). Through the variation of a parameter called learning rate, this algorithm provides a good equilibrium between the speed of Newton’s method and the convergence guaranteed by the steepest descent procedure. Each iteration in the optimization procedure is named* epoch*.

All the optimization methods (also named training algorithms) begin with setting the initial weights (a point in the error surface; see (27)). Usually, the weights are initialized using small random values [26]. Hence, the outputs of the activation functions fall out of the saturation region (low magnitude derivatives region), turning the training procedure faster.

Since the input and output parameters can present different order of magnitude values, they must be scaled in some way. As in this work the parameters values will be in a defined range, the training data is scaled to fall in the range. This is done using the mathematical expression:where the minor subscripts and refer to the maximum and minimum values present in the training sample, respectively. For some applications (as time series), the data is normalized with respect to its mean and deviations values. The data normalization also aids in turning the training process faster. For a network trained with scaled data, the network application for new inputs will result in scaled output values, which must be reversed for an adequate use.

In general, two types of errors may occur during the network training. The first is related to the stopping of training in the beginning of the iterative procedure, far from a good set of weights, resulting in high approximation errors (even for trained data). The second is related to the overfitting of the training sample: the network adapts too much to the particular training set. Both cases lead to poor estimation for new inputs (nontrained data) [26, 27]. Therefore, there is an optimum level of training to be achieved. A common and simple procedure applied to avoid both errors is to separate part of the training set for the training itself (weights optimization) and the other part for validation, with separated error function evaluations. The usual behavior of these errors is depicted in Figure 7. The training process must be stopped when the error of the validation set begins to increase, avoiding the overfitting and ensuring good network generalization.

#### 6. Symbolic Regression (SR)

Genetic Programming (GP) consists of searching and optimizing executable expressions modeled as natural evolution [28]. The GP is inspired by Darwin’s theory of evolution, which states that the fittest species are the ones to survive. The three concepts stated by Darwin (*reproduction*,* crossover*, and* mutation*) are the main operations followed in GP. These operations control the evolution of the* population* or the set of solutions to the problems. Each of these solutions is named* chromosome* and, at each* generation*, it will undergo a* fitness* process. This process determines if each of the chromosomes is more likely to be maintained in the next generation (the higher the fitness rate, the higher the odds) [29] or not. Crossover is a reproduction process that accounts for the combination of genes of two successful chromosomes. Mutations occur in the same generation and alter the genes of its chromosomes. These processes are controlled by the crossover and mutation rates, respectively. Figure 8 presents a flowchart that exemplifies the main steps of a generic genetic algorithm.

Symbolic regression (SR) is a data-driven approach to extract an appropriate model from the space of mathematical expressions, , which is defined by a set of binary operations and mathematical functions [30]. This regression technique, unlike traditional regression methods, does not require a mathematical model of a given form being a kind of nondeterministic polynomial problem, which simultaneously optimizes the structure and coefficients of a target model [30].

SR is used in a wide branch of disciplines, varying from physical to psychological and from chemistry to astronomy. In engineering, it also has a wide branch of applications, such as deriving implicit or explicit parameters relations [31], simplifying complex mathematical relations into concise expressions [32], fitting numerical data to an equation to avoid dealing with unstable and hard to implement mathematical models [33], distilling natural laws from experimental data [34], and automating of reverse engineering problems in dynamics of nonlinear systems [35].

Symbolic regression initially works forming expressions by randomly combining what is called* mathematical building blocks*. These blocks comprise algebraic operators , analytical functions (i.e., sine, cosine, and exponential), constants, and state variables [34]. Symbolic equations, in each step, compete to model the experimental data retaining equations that best fit the given dataset while abandoning the unpromising solutions. This procedure can be written as [30]where and represent the sampling data; is the result equation and is the target model.

New expressions are formed by a combination of previously accepted equations with probabilistic adjustment of their subexpressions [34]. After a considerable amount of computational time, the algorithm returns a set of equations, their size and fitness corresponding to the given data series. These equations vary from the most accurate and sophisticated to the simplest and less accurate. Paraphrasing the suggestion given by Schmidt and Lipson [34], in order to obtain a physical significant equation, the adopted solution should pass Occam’s razor principle, which states that, given more than one solution that satisfactorily fits the same problem, that is, solutions with desirable level of accuracy, the simplest one shall be adopted.

In this work, the commercial software Eureqa® [34] was employed for the symbolic regression of the data generated by the numerical solution of the geodesic equation. The software allows the user to choose between different error metrics [34]. In this work, the mean squared error (MSE) (see (30)) was chosen, because the noise in the data would, by hypothesis, follow a normal distribution. The same error metric is used in the ANN training.

#### 7. Case Study

This work proposal consists in training ANNs and obtaining a physically meaningful equation by SR, both capable of predicting the path of the geodesic curve followed by a tensile armor when a flexible pipe is under bending. Aiming at predicting this path, geometrical parameters such as the torus minor (mean radius of the armor, ) and major (bend radius of the pipe, ) radii, initial lay angle of the armor, length of the armor , and the normalized arc-length are taken as input parameters, while the angles that define the geodesic differential equations (see (17)), and , are the aimed outputs. This can be represented by the following expressions:

When dealing with machine learning algorithms, a typical approach is to feed the algorithm with inputs and target outputs for training. Thus, multiple cases are selected to compose the training dataset for the ANN and the SR approaches. These cases comprise variations in the armor geometrical parameters and the imposed bending radius. By considering these parameters and the geodesic differential equations, (see (17)), the radial angle and the lay angle responses are obtained. The training dataset can thus be represented in the following form:where is the total number of selected cases and is the total number of points in the arc-length in each case. The first five columns are the inputs of the network/genetic algorithm, while the sixth and seventh columns are the outputs (see (31)). In this manner, the inputs can be seen like vectors. The outputs for training the network and evolving the population of the genetic algorithm are two vectors of values of and , respectively*.*

The selection of cases must be done with caution. The limit values of the parameters must be present in the training set to cover their entire domain. Furthermore, specifically for the ANN, cases with a combination of values inside the domain are required to enable the network learning of the cross relationship between the parameters. Additionally, when selecting the input data, the identification of possible outliers is necessary. These outliers, which may arise from numerical errors and/or chosen unreal armor parameters, may conceive bias to correlation and, therefore, hamper the modeling of the data.

It is important to highlight that the geodesic lay angle is a periodic cosine-like shape function, while the parameter is close to a linear function, both with respect to the arc-length [16]. Therefore, only one pitch of the armor is considered as representative of the entire length, which significantly reduces the computational efforts of the finite difference problem (RK solution). As the shapes of the functions and are approximately known, initial forms of these functions, known as base equations (or* seeds*), may be stated for the search of the SR program. This procedure consists in giving a short-cut in the trajectory followed by the algorithm, which speeds up the convergence of the model. These base equations should be used with caution, because misleading base equations may result in greater computational time effort and even unreliable responses of the SR program. In this work, the base equations for and are written, respectively, as

A total of 106 cases were selected with discretization in the arc-length of 51 points (from 0 to 1, normalized) totalizing 5406 training (inputs and outputs) data points. The range values of the inputs are presented in Table 1 and comprise the common usage scenarios of flexible pipes.

Moreover, for the ANN to work properly, it is crucial to randomly sort the input dataset for its training. This is necessary because discontinuities on the parameter arise from case to case and they may end up governing the fitting of the network. On the other hand, in the SR approach, the data should be organized case by case in an ascendant manner to obtain better accuracy.

As discussed in Section 5, the network architecture consists in the number of inputs, outputs, and nodes in the hidden layer. While the numbers of inputs and outputs are defined by the problem, the quantity of hidden nodes must be chosen. The capacity of the network to represent functions increases with the number of nodes employed. However, an excess of hidden nodes turns the network prone to overfitting and weak generalization [26]. One way to estimate the sufficient quantity of nodes is to perform a sensitivity analysis by iteratively varying this parameter and checking the validation fitting. As the network weights are initialized randomly, the training procedure stops at a different point in the error surface (a local minimum) every time the training is done. Hence, in the sensitivity analysis, the training was repeated 15 times for each number of nodes to compose a statistical view of the errors (mean and variation). From this study and some initial test trainings, a hidden layer with 20 nodes was found to be enough to reach good fitting for both geodesic angles.

Therefore, in this work, two ANNs, one for the radial angle and the other for the lay angle , were trained with 20 nodes each using 80% of the data for the training itself and the remaining 20% is separated for validation. The selection was done randomly in the shuffled dataset generated by numerically solving (17). All networks were trained using the Neural Network Toolbox® from MATLAB® [36].

On the other hand, the SR approach results in a series of equations as viable solutions. Each of these equations carries a fit measure and a complexity size. The best fit is also the most complex equation and, therefore, it is wise to choose the less complex equation with enough accuracy among all solutions provided. Moreover, as the geodesic has no geodesic (or binormal) curvature, only equations that nullify (14) are considered. In the performed SR analyses, the equations that satisfy all these requirements are

With the trained ANNs and the SR equations, new cases are selected to test their performance. These cases have parameters in ranges between the values applied in the training of the networks and the evolution of the genetic algorithm. Their values are presented in Table 2. Altogether, 49 cases were generated, which constitute 2499 testing data points.

In Figure 9, the results of the ANN approximation and the SR equations (see (34)) are plotted for a case in which = 0.17 m, = 8 m, = 31 degrees, and = 2.074 m, while, in Figure 10, the results obtained in a case with = 0.35 m, = 13 m, = 22 degrees, and = 5.871 m are presented. These figures indicate that the ANN and SR results showed agreement with the numerical integration results of (17) (named RK).

**(a)**

**(b)**

**(a)**

**(b)**

Figures 11 and 12 present the correlation coefficient and the MSE for the radial and lay angles in all test cases, respectively. These figures indicate that the correlation of the radial angle is very high for both the SR and ANN results, as expected from its almost linear behavior. For the lay angle, the ANN approach led to higher correlation coefficients than those from the SR approach. The MSEs for the radial angle are higher than the MSEs for the laying angle, despite its high fitting.

Furthermore, in what concerns the mechanical analysis of flexible pipes, the proposed models (ANNs and SR equations) may be employed to evaluate the curvatures (geodesic and normal) and torsion of the geodesics considering (4) and (5). Moreover, these results can also be compared to those obtained using the approximate expressions proposed by Sævik [8]:

Figures 13 and 14 present the variation of the curvatures and torsion for the cases illustrated in Figures 9 and 10. These quantities were predicted by the ANNs, the SR equations, the approximate expressions proposed by Sævik [8], and the numerical solution of (17). It is important to highlight, however, that the values presented in Figures 13 and 14 are variations of the curvatures and torsion with respect to the initial configuration of the considered pipes. In these cases, the pipes are initially straight and the armors can be represented by cylindrical helices. The variations of these quantities are, thus, obtained by subtracting the initial curvature parameters (see (36)) from their final values (see (4)–(6)).

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

Figures 13 and 14 indicate that the quantities calculated with the proposed ANNs and SR equations are very accurate in comparison to the numerical values. On the other hand, the proposed ANNs, in some cases, indicate nonnull geodesic curvatures. The ANNs consider that the dataset employed in their training is a continuous sign, which is approximated by a combination of weights and hyperbolic tangential functions. The overall fitting to the training data is thus approximately the same along the whole sign. However, the employed dataset (this sign) is composed of several cases, each one corresponding to a tensile armor in a different flexible pipe, resulting in unmeaningful discontinuities from one case to another that the ANNs try to represent. Hence, the ANNs overfit these discontinuities and may lose accuracy in the overall behavior of the meaningful part of the response. On the other hand, the SR searches for closed form functions that appear to satisfy the training data, thus guaranteeing the well representation of the behavior of the geodesic among the meaningful discrete data points. Therefore, the SR equations, in this case, seem to be more suitable for the mechanical analysis of flexible pipes.

In the bending mechanical analysis of flexible pipes, equations similar to those proposed by Sævik [8] are usually employed. Figures 13 and 14 indicate that these equations present differences with respect to the SR, ANN, and numerical results. Figure 14, for instance, indicates that the SR and the numerical results are almost coincident, while the maximum normal curvature predicted by (35) is about 12% higher than the value estimated by the SR equation and the geodesic torsion is 8% lower. Moreover, differences occur not only in the mid-section of the armor but also along all its length.

The bending stresses in the tensile armors of flexible pipes are proportional to the calculated normal curvature and the torsion in these armors [11]. Furthermore, bending stresses are important in the computation of their fatigue lives and slight changes in these stresses may lead to significant modifications in the armor’s service life [17]. Moreover, the normal curvatures are usually considered as initial imperfections in the buckling analysis of the tensile armors and, again, the buckling loads may be affected by changes in the values of these curvatures [16]. Altogether, the use of the proposed SR equations may be adequate in both types of analyses.

#### 8. Conclusions

The stiff nonlinear differential equations of the geodesic curve of a torus do not have exact solutions and, in this work, two different approaches are proposed and compared to approximate this curve, that is, artificial neural networks (ANNs) and equations obtained by symbolic regression (SR). These approximated approaches rely on a dataset generated by numerically solving the geodesic differential equations, in terms of the radial and lay angles, using the Runge-Kutta (RK) Lobatto IIIa method.

Both the ANN and the SR approaches agreed quite well with the dataset results. However, when tested with new data, the SR equations indicate null geodesic curvatures, while the ANN calculates nonnull values. The difference between the ANN and SR predictions may be probably due to discontinuities in the dataset which are caused by the different case conditions considered in solving the geodesic equations. These unmeaningful discontinuities affect the overall behavior of the ANN in the meaningful part of the response. The SR approach, nonetheless, indicated several possible equations that have the same correlation coefficient with the RK results, but the admissible equation must have null geodesic curvature. The solution proposed in this work is thus the one that has the maximum correlation coefficient with null geodesic curvature.

These approximated approaches are of fundamental importance in the local mechanical bending analysis of flexible pipes. In this analysis, the geodesic of a torus is assumed as a limit path for the tensile armors of flexible pipes and the curvatures and torsion of the geodesic need to be determined to calculate the stresses in these armors. Moreover, the definition of a limit curve is also important in the stability analysis of these armors when the pipe is under axial compression. The results obtained with the proposed approaches were compared to those from equations usually employed in the mechanical bending analysis of flexible pipes. The curvatures and torsion predicted with the SR approach agree with the numerical (RK) calculation, but not all cases coincide with the usual equations. Deviations between the analytical and SR normal curvature and the geodesic torsion up to 12% and 8% were observed, respectively.

To sum up, it is authors’ belief that the ANNs and, mainly, the SR equations proposed here can be used as efficient approximations to the geodesics of a torus, which is of importance not only in the analysis of flexible pipes but also in other engineering applications.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this article.

#### Acknowledgments

The authors acknowledge the National Council for Scientific and Technological Development (CNPq) for the research grants PQ 309278/2013-9 and PQ 308403/2016-9, which are essential for the accomplishment of this work.