PREDICTION OF EQUILIBRIUM SCOUR HOLES USING OPTIMIZATION PROCEDURES

This paper presents an optimization procedure that finds the equilibrium scour depth under a pipeline. Even though much knowledge on scour is available for the most typical marine structures such as a vertical circular monopile of a horizontal pipelines the calculation of the scour depth complex and time-consuming as the transient solution is often modelled as well. In this paper we present a optimization procedure that combined with a computational fluid dynamics, and a model of the bed load finds the equilibrium shape of a scour hole. This can potentially speed up the calculation of the shape of the equilibrium scour hole with a factor of 100. However, it comes with a coast as we will not model the transition and the time scale of the scour hole development.Recorded Presentation from the vICCE (YouTube Link): https://youtu.be/LpKq9Twj7zo


INTRODUCTION
Scour around and adjacent to marine and coastal structures is an important design issue. The study of scour and design of scour protection is often based on theoretical considerations and hydraulic experiments. An introduction to the field is given in (Sumer and Fredsøe, 2002). During the last couple of decades computational fluid dynamics (CFD) has been used to model the scour development, see for instance (Roulund et al., 2005), (Baykal et al., 2015), and (Fuhrman et al., 2014). The use of CFDmodels has shown promising results. They can predict transition from for instance a flatbed to a fully developed scour hole with good accuracy. A major drawback of such models is the long computational time needed for modelling the whole transition from flatbed to fully developed scour whole. These calculations can take weeks or in some cases even months.
As an alternative, only the depth of the scour whole could be modelled. In (Dixen et al., 2012) CFD was used to calculate scour/deposition rates for prescribed frustum shaped scour holes around a monopile. The comparison with experimental results was encouraging. The main objective is to model the scour depth over long time span as suggested in (Nielsen and Hansen, 2007). In this method, scour and back-filling rates are needed. These can be found from either experiments or computational fluid dynamics as for instance in (Fuhrman et al., 2014). To use the method outlined in (Dixen et al., 2012) an approximate shape of the scour hole is needed. This could be based on the equilibrium scour hole. In this study we used an optimization procedure to assess the shape of the equilibrium scour hole shape for horizontal pipelines above a sandy seabed.
The remaining part of this paper will described the methodology, which three parts the; the optimization procedure, the CFD-model and the bed-load model. This is followed by examples af finding the equilibrium scour depth under a pipe line and under tandem pipeline arrangement. Finally a chapter summarizes and concludes the study.

Optimisation methodology
In order to find a find an optimal value of a continuous function f(x), one can use a numerical technique such as Newton's method to iterate towards a local minimum. The method expresses a change in the function Δx towards a local extrema, with a second order expression. Considering the following Taylor series expansion of the function: Taylors formula gives and approximate solution to a function: A function of as ( )= ( + ) is defined, and extrema at ( ( )/ ( ) = ) One can equate the left-hand side to a function of as ( )= ( + ), then by taking the derivative of this function with respect to and setting it to zero to find a extrema ( ( ) ( ) = ) yields: where is the Hessian matrix, and by rearranging, results in the following: where is the directional change towards the optima. New iterate = + , where is a selected step size.
The following section presents an implementation of the Boyden-Fletcher-Goldfarb-Shanno (BFGS) method (Nocedal and Wright, 2006), (Arora, 2004), which is classified as a 'quasi-Newton' method, as the Hessian matrix ( ) is approximated at each iteration rather than directly calculated. The implementation here was performed in MATLAB  , where the algorithm set-up, ran and analyzed OpenFOAM simulations on a high performance cluster to find a minimized cost function pertaining to an equilibrium scour bed.
The optimization problem was formulated such that a cost function or also termed an 'objective function' had to be minimized. In the case of the live scour bed, we proposed using the following: Where ∞ is the free-stream bed-load, ( ) is the bed load along the flow direction, is an integral gain coefficient, and [ − ] is the extents of the bed. Other functions could also be considered such as minimizing the Exner equation along the bed, i.e. find a bed shape, which minimizes the integrated change of the bed height. The optimal value of the cost function was found with the Boyden-Fletcher-Goldfarb-Shanno (BFGS) method, (Bonnans et al., 2006).he following outlines the implemented algorithm where: 1. The procedure is initialized with an iteration counter = 0, and an initial shape of the bed with design variables: = [ 1 2 3 … . . ]. These design variables are input to a shape function that describe the shape of the seabed, which is outlined further later. At the initialization, the Hessian matrix is equated to an identity matrix ( = ). The gradient vector = ∇f( ), is then computed linearly by calculating a new value of the cost function by perturbing each of the design variable by ' ' (a user defined small step size of each design variable) as demonstrated for the design variable 1 by: Calculating gradients requires ′ ′ OpenFOAM calculations to evaluate the gradients for each design variable. On a high performance cluster, this can be achieved in parallel. 2. Evaluate the norm of the gradient function, as a stopping criterion of the algorithm. If it is larger than a small tolerance continue (‖ ‖ > ), else terminate. 3. Calculate the direction of descent by the following = −( ) − , where is the iteration number. 4. Evaluate the step size ( ). This is achieved from a user defined array of 10 or so alpha values. After computation for each case, the step size value, which minimizes the cost function, is selected. Armijo's rule (Arora, 2004) is one of many more formal approaches to finding the step size. 5. Update the design variables + = + 6. Calculate the new gradient function + = ∇f( ) based on a linear evaluation of the gradients as shown in eq. (5), now for the updated design variables. 7. Update the Hessian matrix using the following: Increment the iteration number ( = + 1) and Go to Step 2)

CFD-method
The CFD model was based on the open source CFD library/code OpenFOAM, ("OpenFOAM, The open source CFD toolbox," 2018). The model solved the Navier-Stokes equations with a RANS turbulence model, and a sediment transport model.
The following outlines the setup for each OpenFOAM simulation, where the hydrodynamics and sediment transport is calculated for each alteration of the seabed. Fig. 1 shows the general geometry for the computational runs.

Fig. 1 Computational model domain for analyzing 2D pipeline scour
The boundary conditions applied to the 'Inlet' in the figure for the velocity and turbulence quantities was mapped from a separate precursor solution of a fully developed boundary layer, driven by a body force similar to the methodology described in (Petersen et al., 2020), while a zero-gradient condition was applied for the pressure ( = 0).
Pressure was set to a reference value of zero and zero gradient conditions were applied to all other variables in the 'Outlet'. The 'Sides' were specified as 'empty' to enforce a two dimensional solution. The 'Bed' and 'Cylinder' walls were modelled with no-slip conditions and used the standard wall functions in OpenFOAM 2.2.2 for the turbulent quantities, with the k-ω SST model (Menter, 2003). A setup using the same turbulence model and wall conditions for forces on a cylinder were presented in (Petersen et al., 2020).
The 'Top' boundary was modelled as a friction-less lid with a slip condition. The dimensions of the domain are 20D from the pipeline center to the 'Inlet' and 30D to the 'Outlet', while the height was 10D tall. The diameter of the pipeline was = 0.1 .
The SIMPLE (Semi-Implicit Method for Pressure Linked Equations) algorithm (Ferziger and Peric´, 2002) was used to calculate the hydrodynamics. This iterates to a quasi-steady state, as vortex shedding behind the cylinder is expected, a true steady state cannot be attained. A model simulation is deemed 'steady' once an equilibrium of systematic vortex shedding pattern behind the cylinder is achieved. Using a steady-state solver is a practical measure, as model runs were comparably faster than the transient equivalents.

Bed load
The non-linear vectorial bed-load model (Kovacs and Parker, 1994) was implemented in OpenFOAM and used in the following analysis. Properties of the sediment modelled were as follows, sediment diameter, = 0.26 , specific weight, = 2.65, and coulomb friction factor = 0.73. This model will estimate the sediment transport.

Bed Shape and Shape Functions
Ideally the shape of the bed is modelled in a simplistic fashion, in order to minimise the number of design variables. Employment of two cubic Hermite style spline functions provides a good general form for a hole in 2D, applied in the following piece-wise fashion (Pedersen, 2016):  The bed could also more generally comprise of piecewise polynomials evaluated as basis splines or 'B-splines'. Fig. 3 demonstrates the general approximation of a B-spline to the experimental scour hole from (Mao, 1986). The red dots in the figure show the control points used to manipulate the shape of the B-spline. Calculation of the B-spline was done via de Boors algorithm (De Boor, 1977), implemented in MATLAB  by (Hunyadi, 2012). For the optimisation procedure, the design variables are the control points, where the shape of B-spline can be manipulated. The five control points in the figure, will result in 10 design variables as each point has a vertical and horizontal coordinate.

Outline of the final algorithm
The BFGS algorithm as outlined in the previous section was code as a MATLAB  -script, for gradient evaluations of the cost function (Step 6) and step size evaluations (Step 4), new OpenFOAM runs needed to be set up, executed and analyzed automatically for the algorithm to proceed. OpenFOAM runs were set up on the high performance computing (HPC) cluster at DTU, such that the runs could be executed in parallel over many computational cores. The folder structure on the cluster was organized in terms of the iteration number from the algorithm (i.e. folder 'Iter1' will contain all the relevant simulations for that iteration). A seed simulation was set up, this contained the correct boundary conditions, solver options, numerical schemes used for each simulation. When the algorithm needed to create a new OpenFOAM run, it simply copied the seed folder and replaced the geometry files (to calculate a new shape of the bed). The shape of bed was interpolated using a selected shape function and converted to a Stereolithographic format (.stl), using a script available from the MathWorks user File Exchange (Holcombe, n.d.). The .stl file was then set in the OpenFOAM run directory and used to create a new computational mesh using the automated meshing tool snappyHexMesh. snappyHexMesh is standard in the OpenFOAM release and creates a predominately hexahedral element mesh, with cut cells used to conform the mesh to the geometry surface described by the .stl file. Layered elements at the boundary were added to the surface for better modelling of the boundary layers. Fig. 4 shows an example of a constructed computational mesh for a bed and pipeline. The .stl surfaces used for the bed and pipeline are shown in red. The graded boundary elements start with a cell center roughly situated at + = 50 in the logarithmic boundary layer.  The new .stl geometries are then put into their respective simulation folders. The snappyHexMesh procedure executes on the cluster for each of the runs in order to create new computational meshes. On completion of the mesh creation, the runs are executed in parallel. The MATLAB  script continuously queries the cluster to check for completion of the simulations. Once finished, the bed load profiles from each run are transferred back to MATLAB  for the algorithm to calculate the gradients as presented in eq. (5) and proceed with the subsequent steps of the BFGS algorithm.

OPTIMISATION OF SCOUR UNDER A PIPELINE
The following presents test cases of the optimization routine. Fig. 6 shows the iteration of the BFGS algorithm for the bed shape represented by the Hermite interpolation as outlined in eqs. (6) to (9). An initial shape of the bed assumes a small hole under the pipeline, this is visible in the figure with the line corresponding to 'Iter 1'. An initial bed scour bed shape allows for a gap between the cylinder and the bed, which makes the meshing of the problem easier. Making an initial guess of a hole around a bottom-mounted structure could be made from viewing the extent of the shear stress amplification around the object on the bed. From Fig. 6 it can be seen that the algorithm finds a direction of decent by bringing the shape of the hole downwards. This would be expected due to the constriction between the pipeline and the bed, resulting in high shear stress and consequently enhanced bed load transport. Proceeding iterations push the lee side of the boundary away from the pipeline, this is to reduce the enhanced bedload from the shear stresses of the shed vortices from the pipeline. Progressive reduction of the cost function with each iteration is shown in Fig. 7.
Comparison of the equilibrium scour depth shows that the optimized depth is larger than the experiment of (Mao, 1986). However, it can be argued that it is close to the higher end of the scatter of equilibrium scour depths as presented in (Sumer and Fredsøe, 1990), where the mean scour depth and standard deviation are derived to be S/D=0.6 ,σ/D=0.2 (Sumer and Fredsøe, 2002). Generally, the shape of the hole was well captured. Fig. 8 presents the optimization using B-splines, with five control points as demonstrated in Fig. 3. Here the shape of the final iteration also provides a reasonable comparison against the experiments. As there are more degrees of freedom than the optimization shown in Fig. 6, the front end of the scour hole has a slight convex bend, which would not be physical. Improvement of the algorithm could incorporate a constraint to the bed surface with regard to the angle of repose, by potentially implementing a 'sand-slide' routine as used in morphological simulations (Jacobsen, 2011). The progression of the cost function is shown in Fig. 9.  Fig. 10 plots the optimization run for two pipelines in tandem, exposed to current in the live bed regime. The pipelines are positioned such that they are placed on the flatbed ( = 0) and are ( = 3), where is the distance between the pipelines from their respective edges. The diameter of the pipelines was = 0.1 , as used in the previous optimisation runs. The figure compares the optimisation run with experiments and numerical results for the problem presented in (Zhao et al., 2015). The numerical results are from a morphological CFD model, where both suspended and bed load were modelled. It is apparent from the figure that the optimization routine captures the scours depths well in comparison to the experiments. The small rise of the bed between the pipelines within hole observed from the experiments was not resolved by the optimization routine, or very well by the numerical results from the paper. It should be noted that such a bump in the bed could be the result of a transient development of the hole, which may not be possible to resolve at all by an optimization routine which seeks an equilibrium shape of the scour hole, essentially bypassing the realistic development of the hole. Fig. 11 shows the development of the cost function against iterations for the tandem problem. The optimization problems presented in this section took approximately 1-1.5 days to compute on the HPC cluster at DTU. To get an idea of the number of runs required, for example, if we were solving a problem with 10 design variables, each iteration would require 10 evaluations (OpenFOAM simulations) for determining of the gradients of the cost function. Furthermore, runs for determining the size of the step size( ), could mean an additional 10-15 simulations per iteration. As the computation ran on a cluster, the runs per iteration can be set off in parallel, but it does become apparent that a practicle balance between the number of design variables and representation of the bed geometry must be kept in order to keep the routine efficient.
The routine at present has at times issues with robustness, where it can have difficulty finding search directions for extending the sides of the hole. However, the scour depth of the hole was found relatively effectively by the routine. It was found by the author that changes of the gain coefficient of the cost function ( = 2, eq. (4)) and refined step size evaluations appeared to improve the prediction of the hole width.

SUMMARY
To fast estimate the shape of a scour hole, scour development was formulated in terms of an optimization problem. The technique was applied in 2D for the scour around pipelines. Results showed that the method has potential to predict an equilibrium scour hole, in a relatively fast manner in comparison to a fully morphological CFD calculations. A series of B-splines described the bed form. We only incorporated bed-load in the calculations. Future work could include suspended load as well. Furthermore, the analysis could be extended to a transient solution using the PISO algorithm instead of the steady-state SIMPLE algorithm, to model the equilibrium hole shape under waves. This might require a phase averaging of the cost function over the wave periods, in order to find an equilibrium depth.
The shape functions such as the B-splines needs to be extended to 3D to form panels in order to describe a 3D scour problem, such as the scour around a monopile or 3D scour around a pipeline. The use of the BFGS method in 3D a balance must be attained between a reasonable number of design variables and adequate representation of the bed surface.