Any versions for multi-criterion search of several equations system solution are based on different methods of generating multiple solutions, satisfying Pareto conditions. Choosing candidate solution out of Pareto-optimal population must be based on analysis of hydrodynamic process and is similar to identification procedure of mathematical model. In any case, procedure of multi-criterion optimization comes to solving single-criterion problems, forming multiple possible solutions. At the same time particularities of some computational approaches of fluid dynamics allows using iteration algorithms, where on

The computational procedure described below is analogous to MAC method (Fletcher, 1991), investigating possibility of NMWR application based on neural net trial functions.

Computational capabilities of the developed algorithm can be illustrated by the example of the solution of Navier-Stokes momentum equations, describing two-dimensional isothermal flows of viscous incompressible fluid. On the first stage we will be using this algorithm for

> 2 2 2 2 0 *x y*

Let us consider the flow of incompressible fluid in the channel (fig. 3).

. (12)

*E y f wb b*

*jj j j j i*

*jj j j j*

>

**xx x** ; (9)

**xx x** . (10)

*E* . (11)

,1 , *s ss s s s*

Here **w** is vector of current values of network weights and thresholds.

**3. Using NMWR for hydrodynamics equations solving**

each step solution at only one physical magnitude is generated.

*ij*

*j*

*b*

with the general formula

Statnikov & Matusov, 1995).

**3.1 Modeling flows – the first step**

**Laplace equation solution**

Laplace equation solution

*v*

Hence we have neural network approximation for given equation, which can be presented by the formula

$$y = 13.786 \frac{1}{1 + e^{1.954913 + 0.983267 \cdot x}} - 3.95569 \frac{1}{1 + e^{5.098 + 0.108345 \cdot x}} + 28.3978 \frac{1}{1 + e^{2.532 - 0.7599 \cdot x}} - 3.7751 \tag{5}$$

Using nonlinear optimization universal program products for ANN training is limited to neural networks of the simplest structure, for dimension of optimization tasks solved by data packages does not normally exceed 100; however, it frequently forms 10-20 independent variables due to the fact that efficiency of neural network optimization methods generally falls under the greater dimensions of the nonlinear programming free task. On the other hand, the same neural network training optimization methods prove efficient under much greater dimensions of vector independent variables. Within the framework of given functioning, the standard program codes of neural network models are applied, using the well-known optimization procedures, e.g. Levenberg-Markardt or conjugate gradients - and the computing block of trained neural network with those obtained by the analytical expressions for objective function of the training anti-gradient components, which in composition of the equation under investigation acts as a "teacher" is designed.

#### **2.2 Computing algorithm of minimization of neural network decision**

Let us consider perceptron operation with one hidden layer from N neuron and one output (1). As training objective function, total RMS error (4) will be considered. The objective function shall be presented as a complex function from neural network parameters; components of its gradient shall be calculated using complex function formula. Network output, therefore, is calculated by the following formula:

$$\log \left( \mathbf{x}^s \right) = \sum\_{j} w\_j \sigma\_j \left( \mathbf{x}^s \right), \tag{6}$$

where **x** - vector of inputs, *s* - number of point in training sample, *(x)* - activation function, *wj* - weights of output neuron, *j* - number of neuron in hidden layer. For activation functions, logistic sigmoid will be considered

$$\sigma\_j(\mathbf{x}) = \frac{1}{1 + e^{-t\_j \left(\mathbf{x}, b\_j\right)}} \,. \tag{7}$$

Here *bj* - threshold of *j*-number neuron of hidden layer; the function *tj(***x***,bj)* , however, has form of *j j* , *ij i j i t b vx b* **x** , where *vi* - neuron weight of hidden layer.

While training on each iterations (the epoch) we shall correct the parameters of ANN toward the anti-gradient of objective function - *E*(**v,w,b**), which components are presented in the following form:

$$\frac{\partial E^{s}}{\partial w\_{ij}} = \sigma\_{j}\left(\mathbf{x}^{s}, b\_{j}\right) \cdot \left(y^{s} - f\left(\mathbf{x}^{s}\right)\right);\tag{8}$$

206 Hydrodynamics – Optimizing Methods and Tools

Hence we have neural network approximation for given equation, which can be presented

1.954913 0.983267 5.098 0.108345 2.532 0.75393 1 11 13.786 3.95569 28.3978 3.7751

Using nonlinear optimization universal program products for ANN training is limited to neural networks of the simplest structure, for dimension of optimization tasks solved by data packages does not normally exceed 100; however, it frequently forms 10-20 independent variables due to the fact that efficiency of neural network optimization methods generally falls under the greater dimensions of the nonlinear programming free task. On the other hand, the same neural network training optimization methods prove efficient under much greater dimensions of vector independent variables. Within the framework of given functioning, the standard program codes of neural network models are applied, using the well-known optimization procedures, e.g. Levenberg-Markardt or conjugate gradients - and the computing block of trained neural network with those obtained by the analytical expressions for objective function of the training anti-gradient components, which in composition of the equation under investigation acts as a "teacher" is

Let us consider perceptron operation with one hidden layer from N neuron and one output (1). As training objective function, total RMS error (4) will be considered. The objective function shall be presented as a complex function from neural network parameters; components of its gradient shall be calculated using complex function formula. Network

> *s s j j*

where **x** - vector of inputs, *s* - number of point in training sample, *(x)* - activation function, *wj* - weights of output neuron, *j* - number of neuron in hidden layer. For activation

> , 1 1 *j j j t b e*

Here *bj* - threshold of *j*-number neuron of hidden layer; the function *tj(***x***,bj)* , however, has

While training on each iterations (the epoch) we shall correct the parameters of ANN toward the anti-gradient of objective function - *E*(**v,w,b**), which components are presented

> , *s s ss j j*

*E b yf*

*t b vx b* **x** , where *vi* - neuron weight of hidden layer.

, (6)

**x x** . (7)

**x x** ; (8)

*j y w* **x x** (5)

1 11 *x xx y e ee*

**2.2 Computing algorithm of minimization of neural network decision**

*ij*

*w*

output, therefore, is calculated by the following formula:

functions, logistic sigmoid will be considered

form of *j j* , *ij i j*

in the following form:

*i*

by the formula

designed.

$$\frac{\partial \mathbf{E}^s}{\partial \mathbf{v}\_{\dot{\eta}}} = \left(\mathbf{y}^s - f\left(\mathbf{x}^s\right)\right) \cdot \mathbf{w}\_{\dot{\jmath}} \cdot \sigma\_{\dot{\jmath}}\left(\mathbf{x}^s, b\_{\dot{\jmath}}\right) \cdot \left(\mathbf{1} - \sigma\_{\dot{\jmath}}\left(\mathbf{x}^s, b\_{\dot{\jmath}}\right)\right) \cdot \mathbf{x}\_i^s \; ; \tag{9}$$

$$\frac{\partial \hat{E}^s}{\partial b\_j} = -\left(y^s - f\left(\mathbf{x}^s\right)\right) \cdot w\_j \cdot \sigma\_j\left(\mathbf{x}^s, b\_j\right) \cdot \left(\mathbf{1} - \sigma\_j\left(\mathbf{x}^s, b\_j\right)\right). \tag{10}$$

Thereby, we have got all the components of the gradient of the objective function of minimization, comparatively which iterations will be consecutively realized in accordance with the general formula

$$
\Delta \mathbf{w} = -\varepsilon \nabla E \left( \mathbf{w} \right) \,. \tag{11}
$$

Here **w** is vector of current values of network weights and thresholds.

#### **3. Using NMWR for hydrodynamics equations solving**

Parameter optimization of neural network trial solutions is achieved by applying several optimization strategies and by subsequently choosing the maximum effective one (see Cloete & Zurada, 2000). First strategy is to apply totality of effective gradient methods "starting" from various initial points. The other strategy is to apply structural-parametrical optimization to ANN training; this method is based on indirect statistic optimization method on self-organizing basis or parameter space research (see: Egorov et al., 1998; Statnikov & Matusov, 1995).

Any versions for multi-criterion search of several equations system solution are based on different methods of generating multiple solutions, satisfying Pareto conditions. Choosing candidate solution out of Pareto-optimal population must be based on analysis of hydrodynamic process and is similar to identification procedure of mathematical model. In any case, procedure of multi-criterion optimization comes to solving single-criterion problems, forming multiple possible solutions. At the same time particularities of some computational approaches of fluid dynamics allows using iteration algorithms, where on each step solution at only one physical magnitude is generated.

#### **3.1 Modeling flows – the first step**

The computational procedure described below is analogous to MAC method (Fletcher, 1991), investigating possibility of NMWR application based on neural net trial functions.

#### **Laplace equation solution**

Computational capabilities of the developed algorithm can be illustrated by the example of the solution of Navier-Stokes momentum equations, describing two-dimensional isothermal flows of viscous incompressible fluid. On the first stage we will be using this algorithm for Laplace equation solution

$$
\frac{\partial^2 \varphi}{\partial x^2} + \frac{\partial^2 \varphi}{\partial y^2} = 0 \,\,\,\,\,\tag{12}
$$

Let us consider the flow of incompressible fluid in the channel (fig. 3).

Neural Network Modeling of Hydrodynamics Processes 209

Equation summary residual with substituted trial solutions (1) on arbitrary calculation area points with coordinates **х**s with expressions application (13) can also be calculated

2 2 *s ss s i s i*

Therefore, trial solution (1) training problem of neural network equation consists in SLP hidden layer parameter selection (weights and thresholds) at which the summary residual (14) has the minimal value limited to zero. The computer program described above, with training procedure target function being set functionally by applying analytical expressions

Efficiency of searching of neuronet learning solution parameters depends on problem dimension, i.e. weights and perceptron thresholds variable adjusted quantity. The more significant is neurons quantity in trial solution, the higher is ANN approximate capacity; however, achieving high approximation accuracy is more complicated. At the same time, neuron quantity depends not only on simulated function complexity, but also on calculation nodes quantity in which the residual equation is calculated. It is known that generally points' quantity increase in statistical set used for neural network construction is followed by increase in necessary neurons network (Galushkin, 2002; Galushkin, 2007) quantity. Consequently, the dense calculation grids application results in nonlinear programming problems; while applying rare calculation grids, it is necessary to check the solution realization between calculation nodes, i.e. there is a problem of learning solution procedure standardization. In the neuronet solution reception context on known equation, it is convenient using traditional additive parameter of training neural model quality - a control error which is calculated on the set of additional calculation nodes between calculation grid nodes. Number of these additional calculation grid nodes can be much more significant, and they should cover the whole calculation area, because the nodes number increase with control error on known network parameters does not result in essential computing expenses growth. Hence, referring to learning solution neuronet parameters reception, there exists an issue of solving twice-criterion problem of nonlinear optimization along with minimizing simultaneously both summary residual in control points, or the control error can appear as a restriction parameter, in the limited set of calculation nodes and in this case the neural network solution parameters reception is reduced to the conditional nonlinear optimization

At the first stage, residual distribution of the current equation (5) on various calculation

(5) except for calculation nodes group, for example, in the input border right point vicinity, due to a sudden change of the boundary conditions in this point. In areas with the solution insufficient exactness we will place the calculation nodes additional quantity using the

. As a whole, the received neural network solution satisfies the equation

nodes and the corresponding speed vector distribution , *T*

1 1

*E*

2 2 *x*

and

2 2 *y*

*x xy*

2 2 2 2 22

2 2

. (14)

, is used for parameter adjustment of learning model.

**v** *u v* , where speed nodes

analytically

for second derivatives

problem.

*x*

and *v*

*y*

*u*

Fig. 3. Computational area

Here's how the boundary conditions are defined: on solid walls *u=v=*0, on inflow boundary *u=*0, *v=*1, on outflow boundary 0 *u v x x* . There are no boundary conditions for pressure except for one reference point, where *p=*0 is specified (in the absolute values *p=p*0), considering which indication of incoming into the momentum equation *p x* and *p y* is

realized.

For solving flow equations by predictor method it is necessary to specify initial velocity distribution within the computational area, satisfying the equation of continuity. For this purpose, velocity potential *x*,*y* is introduced and *u x* and *v y* . As a result of Laplace equation solution, velocity distribution is generated, which can be indicated as freevortex component of the sought quantity.

If the result of learning sample neuronet calculations is defined by the following formula *s s j j j* **x x** *v f* , where , *T* **x** *x y* -input variables vector, *s* - point number in the

learning sample, *f* **х** - activation function, *j v* - output neuron weights, *j* - neuron number

in the hidden layer as activation function the logistical sigmoid is used , 1 1 *j j j t b f e* **x x** ,

where *j b* - threshold of the *j*-number neuron hidden layer, and the function *t b j j* **x**, looks like *j j* , *ij i j i t b wx b* **x** where *wij* - hidden layer neurons weights, then analytical expressions for the second speed potential derivatives can be calculated using the following formula

$$\frac{\hat{\sigma}^2 \boldsymbol{\varrho}}{\hat{\sigma} \left(\mathbf{x}\_i^s\right)^2} = \sum\_j \left( \upsilon\_j w\_{ij}^2 \Big( f\_j \left(\mathbf{x}^s, b\_j\right) - 3f\_j^2 \left(\mathbf{x}^s, b\_j\right) + 2f\_j^3 \left(\mathbf{x}^s, b\_j\right) \Big) \right). \tag{13}$$

208 Hydrodynamics – Optimizing Methods and Tools

 0 *y*

0 2 2 2 2

 *x y*

Here's how the boundary conditions are defined: on solid walls *u=v=*0, on inflow boundary

except for one reference point, where *p=*0 is specified (in the absolute values *p=p*0),

For solving flow equations by predictor method it is necessary to specify initial velocity distribution within the computational area, satisfying the equation of continuity. For this

Laplace equation solution, velocity distribution is generated, which can be indicated as free-

If the result of learning sample neuronet calculations is defined by the following formula

learning sample, *f* **х** - activation function, *j v* - output neuron weights, *j* - neuron number

in the hidden layer as activation function the logistical sigmoid is used ,

where *j b* - threshold of the *j*-number neuron hidden layer, and the function *t b j j* **x**, looks

expressions for the second speed potential derivatives can be calculated using the following

 2 2 23 2 ,3 ,2 ,

*j ij j j j j j j s j i vw f b f b f b*

*sss*

**xxx** . (13)

*t b wx b* **x** where *wij* - hidden layer neurons weights, then analytical

*x*,*y* is introduced and *u*

 0 *y*

. There are no boundary conditions for pressure

*x*

**x** *x y* -input variables vector, *s* - point number in the

and *v*

*x*

*x*

*y*

 and *p y* is

. As a result of

1 1 *j j j t b f e* **x x** ,

1 *y*

*x x*

considering which indication of incoming into the momentum equation *p*

 0 *x*

*u=*0, *v=*1, on outflow boundary 0 *u v*

vortex component of the sought quantity.

*x*

**x x** *v f* , where , *T*

Fig. 3. Computational area

purpose, velocity potential

like *j j* , *ij i j i*

 *s s j j j*

realized.

formula

*y*

Equation summary residual with substituted trial solutions (1) on arbitrary calculation area points with coordinates **х**s with expressions application (13) can also be calculated analytically

$$E = \frac{1}{2} \sum\_{s} \sum\_{i} \left( \frac{\hat{\sigma}^2 \rho}{\hat{\sigma} \left( \mathbf{x}\_i^s \right)^2} \right)^2 = \frac{1}{2} \sum\_{s} \left( \frac{\hat{\sigma}^2 \rho}{\hat{\sigma} \left( \mathbf{x}^s \right)^2} + \frac{\hat{\sigma}^2 \rho}{\hat{\sigma} \left( \mathbf{y}^s \right)^2} \right)^2 \,. \tag{14}$$

Therefore, trial solution (1) training problem of neural network equation consists in SLP hidden layer parameter selection (weights and thresholds) at which the summary residual (14) has the minimal value limited to zero. The computer program described above, with training procedure target function being set functionally by applying analytical expressions

for second derivatives 2 2 *x* and 2 2 *y* , is used for parameter adjustment of learning model.

Efficiency of searching of neuronet learning solution parameters depends on problem dimension, i.e. weights and perceptron thresholds variable adjusted quantity. The more significant is neurons quantity in trial solution, the higher is ANN approximate capacity; however, achieving high approximation accuracy is more complicated. At the same time, neuron quantity depends not only on simulated function complexity, but also on calculation nodes quantity in which the residual equation is calculated. It is known that generally points' quantity increase in statistical set used for neural network construction is followed by increase in necessary neurons network (Galushkin, 2002; Galushkin, 2007) quantity. Consequently, the dense calculation grids application results in nonlinear programming problems; while applying rare calculation grids, it is necessary to check the solution realization between calculation nodes, i.e. there is a problem of learning solution procedure standardization. In the neuronet solution reception context on known equation, it is convenient using traditional additive parameter of training neural model quality - a control error which is calculated on the set of additional calculation nodes between calculation grid nodes. Number of these additional calculation grid nodes can be much more significant, and they should cover the whole calculation area, because the nodes number increase with control error on known network parameters does not result in essential computing expenses growth. Hence, referring to learning solution neuronet parameters reception, there exists an issue of solving twice-criterion problem of nonlinear optimization along with minimizing simultaneously both summary residual in control points, or the control error can appear as a restriction parameter, in the limited set of calculation nodes and in this case the neural network solution parameters reception is reduced to the conditional nonlinear optimization problem.

At the first stage, residual distribution of the current equation (5) on various calculation nodes and the corresponding speed vector distribution , *T* **v** *u v* , where speed nodes

*u x* and *v y* . As a whole, the received neural network solution satisfies the equation

(5) except for calculation nodes group, for example, in the input border right point vicinity, due to a sudden change of the boundary conditions in this point. In areas with the solution insufficient exactness we will place the calculation nodes additional quantity using the

Neural Network Modeling of Hydrodynamics Processes 211

solution (12). Calculation grid points are formed as centers of Cohonen network clusters constructed on units coordinates of the uniform rectangular scale and on the right part of the equation (16) corresponding to these units values Δ. Fig. 6 (а) presents formation re sults of the calculation grid and the speed distribution on the pseudo-non-stationary algorithm first iterative step of the Navier-Stokes equation solution. Here it was possible to receive an exact neural network solution for the whole calculation area without using additional set of

Let us now study am incompressible fluid internal flow within a channel with a stream turning (fig. 3). Navier-Stokes equation system describing two-dimensional isothermal flows

> 0 *u v x y*

Re *p u u uu*

Re *p v v vv*

Here *u, v* – nodes speed, Re - Reynolds number. Hydrodynamics equations system is written

can be chosen in the current field, for example an input fluid speed value in the channel and

Boundary conditions are stated as follows: on solid walls *u=v=0*, on the input border *u=0*,

[a,b][c,d] within the plane *XY*, and there is a rectangular analytical grid, specified by

We will understand neural net functions *uvp f xy* ,, ( ,,) *NET* **w** as the (17)-(19) system solution giving minimum of the total squared residual in the knot set of computational grid. The trial

*xxy x y*

*yxy x y*

in the non-dimensional view; i.e. it includes non-dimensional values \* *u*

*u v*

*u v*

2 2 2 2 1 0

2 2 2 2 1 0

. Quality of *u* and any speed and linear size values

. Let us consider that there is rectangular region

**w** ; (20)

**w** ; (21)

1 2

1 2

*i ii i v*

*i ii i u*

are received as a result of the Laplace equation

; (17)

; (18)

. (19)

*u* , \* *v v*

*u* ,

*u*

*y*

where speed nodes *u*

calculation nodes.

\* *r r*

*D* , \*

2

*v=1*, on the output border 0 *u v*

 *u D*

*x x*

Cartesian product of two one-dimensional grids *{xk}, k=l,…,n* and *{yl}, l*=l,…,*m*.

solution (fig. 5) of the system (17)-(19) *u, v, p* can be presented in the form equation (1):

( ,,) ( )

*u xy vf b w x w y b*

*v xy vf b w x w y b*

1

1 ( ,,) ( )

*i*

2

*i q*

*q*

*q*

*p p u* , Re

the channel width *h* .

*x*

of the viscous incompressible fluid (Fletcher, 1991):

also *v*

following algorithm. Let us formulate the Cohonen neural network with three inlet variables presented by the coordinates of available computation nodes *x* and *y*, and also the equations (5) residual value in these nodes, along with the required cluster centers quantity equal to the additional nodes quantity. The cluster center coordinates which will generally be placed in areas with the learning solution low precision (Prokhorov et al., 2001) we will consider additional computation nodes coordinates. The number of these additional nodes in each case is different and defined by iterations, until the decision error does not accept comprehensible value. As a result of the additional formation of received neural network learning solution using additional computation nodes, it turned out to be possible to increase the solution local accuracy in the point *B* vicinity while maintaining the accuracy high in all other points.

Fig. 4. Formation of additional computation nodes for Laplace equation solution

Therefore, not only has the computing experiment proven reception opportunity of the general neural network solution in the calculation area, but also defined coordinate calculation logic of computational nodes for increasing the accuracy of neural network initial equation solution. Let us study a reception opportunity of the Poisson equation solution using an irregular computational grid, i.e. equation total residual with solutions (1) will be calculated in nodes located in the casual image or certain algorithm, which use has not been connected with the necessity of computational grid coordination and computational area borders.

#### **Poisson equation**

Let us study a neural network solution precision on irregular calculation scales for Poisson equation

$$\frac{\partial^2 p}{\partial \mathbf{x}^2} + \frac{\partial^2 p}{\partial y^2} = \boldsymbol{\Delta} \neq \mathbf{0} \,. \tag{15}$$

This equation is particularly used for calculating the pressure distribution as well as for time iterations organization at the Navier-Stokes equations solution by pseudo-non-stationary algorithms (Fletcher, 1991). For the solution we shall use an irregular calculation grid, because, in contrast to fluid dynamics classical numerical methods, it does not result in the neutral network learning functions algorithm complication. Meanwhile, advantages using calculation nodes located in calculation area for the complex geometry study current are obvious. The decision is defined by the equations (15) with the right side as follows

$$\frac{\partial^2 p}{\partial \mathbf{x}^2} + \frac{\partial^2 p}{\partial y^2} = \Delta = 2 \left( \frac{\partial u}{\partial \mathbf{x}} \frac{\partial v}{\partial y} - \frac{\partial v}{\partial \mathbf{x}} \frac{\partial u}{\partial y} \right), \tag{16}$$

210 Hydrodynamics – Optimizing Methods and Tools

following algorithm. Let us formulate the Cohonen neural network with three inlet variables presented by the coordinates of available computation nodes *x* and *y*, and also the equations (5) residual value in these nodes, along with the required cluster centers quantity equal to the additional nodes quantity. The cluster center coordinates which will generally be placed in areas with the learning solution low precision (Prokhorov et al., 2001) we will consider additional computation nodes coordinates. The number of these additional nodes in each case is different and defined by iterations, until the decision error does not accept comprehensible value. As a result of the additional formation of received neural network learning solution using additional computation nodes, it turned out to be possible to increase the solution local accuracy in the point *B* vicinity while maintaining the accuracy

Fig. 4. Formation of additional computation nodes for Laplace equation solution

Therefore, not only has the computing experiment proven reception opportunity of the general neural network solution in the calculation area, but also defined coordinate calculation logic of computational nodes for increasing the accuracy of neural network initial equation solution. Let us study a reception opportunity of the Poisson equation solution using an irregular computational grid, i.e. equation total residual with solutions (1) will be calculated in nodes located in the casual image or certain algorithm, which use has not been connected with the necessity of computational grid coordination and

Let us study a neural network solution precision on irregular calculation scales for Poisson

2 2 0 *p p*

This equation is particularly used for calculating the pressure distribution as well as for time iterations organization at the Navier-Stokes equations solution by pseudo-non-stationary algorithms (Fletcher, 1991). For the solution we shall use an irregular calculation grid, because, in contrast to fluid dynamics classical numerical methods, it does not result in the neutral network learning functions algorithm complication. Meanwhile, advantages using calculation nodes located in calculation area for the complex geometry study current are

> 2 2 2 *p p uv vu x y x y x y*

. (15)

, (16)

2 2

*x y*

obvious. The decision is defined by the equations (15) with the right side as follows

2 2

high in all other points.

computational area borders.

**Poisson equation**

equation

where speed nodes *u x* also *v y* are received as a result of the Laplace equation

solution (12). Calculation grid points are formed as centers of Cohonen network clusters constructed on units coordinates of the uniform rectangular scale and on the right part of the equation (16) corresponding to these units values Δ. Fig. 6 (а) presents formation re sults of the calculation grid and the speed distribution on the pseudo-non-stationary algorithm first iterative step of the Navier-Stokes equation solution. Here it was possible to receive an exact neural network solution for the whole calculation area without using additional set of calculation nodes.

Let us now study am incompressible fluid internal flow within a channel with a stream turning (fig. 3). Navier-Stokes equation system describing two-dimensional isothermal flows of the viscous incompressible fluid (Fletcher, 1991):

$$\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 \; ; \tag{17}$$

$$
\left\| \frac{\partial p}{\partial \mathbf{x}} + u \frac{\partial u}{\partial \mathbf{x}} + v \frac{\partial u}{\partial y} - \frac{1}{\text{Re}} \left| \frac{\partial^2 u}{\partial \mathbf{x}^2} + \frac{\partial^2 u}{\partial y^2} \right| = 0 \; ; \tag{18}
$$

$$
\left\{\frac{\partial p}{\partial y} + \mu \frac{\partial v}{\partial \mathbf{x}} + v \frac{\partial v}{\partial y} - \frac{1}{\text{Re}} \left| \frac{\partial^2 v}{\partial \mathbf{x}^2} + \frac{\partial^2 v}{\partial y^2} \right| = \mathbf{0} \right\}.\tag{19}
$$

Here *u, v* – nodes speed, Re - Reynolds number. Hydrodynamics equations system is written in the non-dimensional view; i.e. it includes non-dimensional values \* *u u u* , \* *v v u* ,

$$\text{For } ^\circ = \frac{r}{D}, \text{ } \notp = \frac{p}{\rho \cdot u\_w^2}, \text{ } \text{Re} = \frac{\rho \cdot u\_w \cdot D}{\mu}. \text{ Quality of } \ u\_w \text{ and any speed and linear size values}$$

can be chosen in the current field, for example an input fluid speed value in the channel and the channel width *h* .

Boundary conditions are stated as follows: on solid walls *u=v=0*, on the input border *u=0*, *v=1*, on the output border 0 *u v x x* . Let us consider that there is rectangular region [a,b][c,d] within the plane *XY*, and there is a rectangular analytical grid, specified by Cartesian product of two one-dimensional grids *{xk}, k=l,…,n* and *{yl}, l*=l,…,*m*.

We will understand neural net functions *uvp f xy* ,, ( ,,) *NET* **w** as the (17)-(19) system solution giving minimum of the total squared residual in the knot set of computational grid. The trial solution (fig. 5) of the system (17)-(19) *u, v, p* can be presented in the form equation (1):

$$u(\mathbf{w}, \mathbf{x}, y) = \sum\_{i=1}^{q} v\_i f\_{\sigma}(b\_i + w\_{i1}\mathbf{x} + w\_{i2}y) + b\_u;\tag{20}$$

$$w(\mathbf{w}, x, y) = \sum\_{i=q+1}^{2q} v\_i f\_\sigma(b\_i + w\_{i1}\mathbf{x} + w\_{i2}y) + b\_v \; ; \tag{21}$$

Neural Network Modeling of Hydrodynamics Processes 213

2 1 1 (,) *n n NET p f xy t*

Vector ( , )*T* **F** *F G* introduced to this algorithm can be calculated by the momentum equations (18-19); or Poisson equation (16) can be used for pressure calculation. Thus, the solution for pressure received from the equation (16) results in the continuity equation realization at the moment of time 1 *n* . Once *n* 1 *p* is calculated, substitution of these values in the formula *nn n* 1 1 *t p* **u F** allows to determine *n* 1 *v* . The iterative process

Let us briefly generalize the above mentioned results of calculation experiment in methodic form to set dynamic calculation scale at the Navier-Stokes equations (17)-(19) solution by the establishment method. First, Laplace equation (12) solution is calculated with the received scale or with rectangular scale, or by means of random numbers generator, or by using the

equation residual in the grid nodes; third, the additional multitude of calculation nodes is generated with using the Cohonen network; then, if the precision is not achieved, points 2 and 3 are realized and the additional calculation components quantity grows until exact neural network solution for the whole calculation area is found; further, pseudo-nonstationary algorithm iterations are organized by using the equations (23)-(24) where equation (16) solution is found on each time step on the multitude of calculation components, the coordinates of which vary depending on the Poisson equation right part distribution for each iterative step; finally, at steps 3 and 5 realization of structural optimization algorithms and learning neural network solutions standardization formation

a)

b)

c)

. (24)

(Statnikov & Matusov, 2002); then, distribution solving

 **F**

goes on until speed distribution stops varying.

Sobol - Statnikov generator of *LP*

Fig. 6. Net velocity distribution

$$p(\mathbf{w}, \mathbf{x}, y) = \sum\_{i=2q+1}^{3q} v\_i f\_{\sigma}(b\_i + w\_{i1}\mathbf{x} + w\_{i2}y) + b\_p \,. \tag{22}$$

Here again, **w** is the vector of all the weights and thresholds of the net. In this case the amount of *q* neurons in the trial solutions remains the same for each decision variable set. This is the parameter on which depend approximated capabilities of neural net trial solution. Result of computational algorithm functioning should be achievement of necessary accuracy level of solution at *q* minimum value*.*

Fig. 5. Neural net trial solution

Let us name the residuals of equations (17)-(19) *R*1, *R*2 and *R*3 correspondingly, then for the vase of NMWR realization for parameters setup of the trial solution it is necessary to minimize three objective functions 222 123 *RRR* ; ; min . In the simplest case, the only solution of the multi-criterion problem of minimization can be generated substituting of three criterions by one, presented in compression form; for example, 2222 *RRRR* 123 min . Presenting the trial solution in the form of continuous functions (20)-(22) allows to define analytically the first and the second differential coefficient in the equations (17)-(19), knowing which one can generate analytic expressions of the function of residuals *Rs(***w**,*x,y)* and further for antigradient component of the total residual in the *s*-reference point at ANN

inner parameters *j R v* , *ij R w* and *j R b* , being later used in the minimization algorithm in

accordance with anti-gradient direction.

For the momentum equations solution by MAC method, it is necessary to specify an initial speeds distribution in the calculation area satisfying to the continuity equation. For this purpose, the speed potential *x y*, is introduced, *u x* and *v y* . As a Laplace equation solution result, we obtain speed distribution which can be called non-vortex required value.Final speeds and pressure distribution are results of the momentum equations solution according to the following algorithm.

The speed distribution on the following time layer is calculated according to the formula

$$\mathbf{u}^{n+1} = \mathbf{F}^n - \Delta t \cdot \nabla p^{n+1} \,, \tag{23}$$

where pressure distribution to each iterative step is defined upon the Poisson equation solution

212 Hydrodynamics – Optimizing Methods and Tools

( ,,) ( )

*p xy vf b w x w y b*

Here again, **w** is the vector of all the weights and thresholds of the net. In this case the amount of *q* neurons in the trial solutions remains the same for each decision variable set. This is the parameter on which depend approximated capabilities of neural net trial solution. Result of computational algorithm functioning should be achievement of necessary

Let us name the residuals of equations (17)-(19) *R*1, *R*2 and *R*3 correspondingly, then for the vase of NMWR realization for parameters setup of the trial solution it is necessary to

of the multi-criterion problem of minimization can be generated substituting of three criterions by one, presented in compression form; for example, 2222 *RRRR* 123 min . Presenting the trial solution in the form of continuous functions (20)-(22) allows to define analytically the first and the second differential coefficient in the equations (17)-(19), knowing which one can generate analytic expressions of the function of residuals *Rs(***w**,*x,y)* and further for antigradient component of the total residual in the *s*-reference point at ANN

For the momentum equations solution by MAC method, it is necessary to specify an initial speeds distribution in the calculation area satisfying to the continuity equation. For this

equation solution result, we obtain speed distribution which can be called non-vortex required value.Final speeds and pressure distribution are results of the momentum

where pressure distribution to each iterative step is defined upon the Poisson equation

*x y*, is introduced, *u*

The speed distribution on the following time layer is calculated according to the formula

1 2

**w** . (22)

123 *RRR* ; ; min . In the simplest case, the only solution

, being later used in the minimization algorithm in

and *v*

*nn n* 1 1 *t p* **u F** , (23)

*y*

. As a Laplace

*x*

*i ii i p*

3

*i q*

accuracy level of solution at *q* minimum value*.*

Fig. 5. Neural net trial solution

inner parameters

solution

minimize three objective functions 222

*j R v* ,

accordance with anti-gradient direction.

purpose, the speed potential

*ij R w*

equations solution according to the following algorithm.

and

*j R b*

*q*

2 1

$$
\nabla^2 p^{n+1} = \frac{1}{\Delta t} \nabla \mathbf{F}^n = f\_{\text{NET}}(\mathbf{x}, y) \tag{24}
$$

Vector ( , )*T* **F** *F G* introduced to this algorithm can be calculated by the momentum equations (18-19); or Poisson equation (16) can be used for pressure calculation. Thus, the solution for pressure received from the equation (16) results in the continuity equation realization at the moment of time 1 *n* . Once *n* 1 *p* is calculated, substitution of these values in the formula *nn n* 1 1 *t p* **u F** allows to determine *n* 1 *v* . The iterative process goes on until speed distribution stops varying.

Let us briefly generalize the above mentioned results of calculation experiment in methodic form to set dynamic calculation scale at the Navier-Stokes equations (17)-(19) solution by the establishment method. First, Laplace equation (12) solution is calculated with the received scale or with rectangular scale, or by means of random numbers generator, or by using the Sobol - Statnikov generator of *LP* (Statnikov & Matusov, 2002); then, distribution solving equation residual in the grid nodes; third, the additional multitude of calculation nodes is generated with using the Cohonen network; then, if the precision is not achieved, points 2 and 3 are realized and the additional calculation components quantity grows until exact neural network solution for the whole calculation area is found; further, pseudo-nonstationary algorithm iterations are organized by using the equations (23)-(24) where equation (16) solution is found on each time step on the multitude of calculation components, the coordinates of which vary depending on the Poisson equation right part distribution for each iterative step; finally, at steps 3 and 5 realization of structural optimization algorithms and learning neural network solutions standardization formation

Fig. 6. Net velocity distribution

Neural Network Modeling of Hydrodynamics Processes 215

At the same time, neuron quantity depends not only on simulated function complexity, but also on calculation nodes quantity in which the residual equation is calculated. It is known that generally points' quantity increase in statistical set used for neural network construction is followed by increase in necessary neurons network (Galushkin, 2002; Galushkin, 2007) quantity. Consequently, the dense calculation grids application results in nonlinear programming problems; while applying rare calculation grids, it is necessary to check the solution realization between calculation nodes, i.e. there is a problem of learning solution procedure standardization. In the neuronet solution reception context on known equation, it is convenient using traditional additive parameter of training neural model quality - a control error which is calculated on the set of additional calculation nodes between calculation grid nodes. The number of these additional calculation grid nodes can be much more important and they should cover all calculation area because the nodes number increase in which the control error on known network parameters does not result in essential computing expenses growth. Hence, at the learning solution neuronet parameters reception there is a problem of solving the multi-criterion problem of the nonlinear optimization and to minimize simultaneously both the summary residual in control points, or the control error can appear as a restriction parameter, in the limited set of calculation nodes and in this case the neural network solution parameters reception is reduced to the

On figure 8, velocity distribution vector within the current random computation nodes of revolving ring on one of neural network training iteration of equation (25)-(26) decisions is presented. Contours of stream function within the computational zone obtained by using

Fig. 7. Rotating Domain

conditional nonlinear optimization problem.

NMWR are presented on fig. 9.

stated in (Kretinin et al., 2010). On fig. 6 (а-c), changing dynamics of calculation components during realizing various moments of time of equations (23)-(24) algorithm is shown, the speed vector distribution reorganization during the transition from the equation (12) solution to the equations (17)-(19) solution is also illustrated.

Finally, analytical solution (as neural network function) for Navier-Stokes equations (17)- (19) systems within the channel with stream turning at Re 100 which is expressed by the formula (23) on last iterative layer and neural network dependence of the pressure distribution on this layer

$$p(\mathbf{w}, \mathbf{x}, y) = \sum\_{i=2q+1}^{3q} v\_i f\_{\sigma}(b\_i + w\_{i1}\mathbf{x} + w\_{i2}y) + b\_p$$

with weights calculated array and the boundary network **w** *wvb* , , .

#### **3.2 Modeling flows in rotating ring zone - the equations that are applicable to rotating reference frame**

For this flow, NMWR solves conservation equations for mass and momentum describing incompressible flows of viscous Newtonian fluid.

Continuity equation

$$\frac{\partial \overline{u}\_j}{\partial \mathbf{x}\_j} = \mathbf{0} \tag{25}$$

Momentum equations

$$
\frac{\partial}{\partial \mathbf{x}\_j} \left( \overline{\mathbf{u}\_i \boldsymbol{\mu}\_j} \right) + \frac{\partial}{\partial \mathbf{x}\_j} \left( \overline{\mathbf{u}\_i' \boldsymbol{\mu}\_j'} \right) = -\frac{\partial p}{\partial \mathbf{x}\_i} + \frac{\partial}{\partial \mathbf{x}\_j} \left[ \mu \left( \overline{\frac{\partial \mathbf{u}\_i}{\partial \mathbf{x}\_j}} + \overline{\frac{\partial \mathbf{u}\_j}{\partial \mathbf{x}\_i}} \right) \right] + f\_i \,\tag{26}
$$

For flows in rotating domain (fig. 7), the equations for conservation of mass and momentum are written for the relative velocity formulation, where *fi* in right hand side is given by

$$\vec{f}\_i = -\rho \left( 2\vec{\alpha} \times \vec{u} + \vec{\alpha} \times \left( \vec{\alpha} \times \vec{r} \right) \right). \tag{27}$$

The absolute velocity formulation is used in the non-rotating domain, and *fi=*0. Thus, the standard *k* turbulence model is used.

Let us consider the neural network functions ,,, ( ,,) *u pk f xy i NN* **w** as equations solution, where **w**– all weights and network thresholds vector that assure summary quadratic residuals minimums for each equation in optional totality of computation nodes coordinating each neural network solution learning iteration, are generated using the random number generator.

Neuronet learning solution parameter search process efficiency depends on problem dimension, i.e. weights and perceptron thresholds varying adjusted quantity. The more significant is neurons quantity in trial solution, the higher is ANN approximate capacity; however, achieving high approximation accuracy is more complicated.

214 Hydrodynamics – Optimizing Methods and Tools

stated in (Kretinin et al., 2010). On fig. 6 (а-c), changing dynamics of calculation components during realizing various moments of time of equations (23)-(24) algorithm is shown, the speed vector distribution reorganization during the transition from the equation (12)

Finally, analytical solution (as neural network function) for Navier-Stokes equations (17)- (19) systems within the channel with stream turning at Re 100 which is expressed by the formula (23) on last iterative layer and neural network dependence of the pressure

( ,,) ( )

*p xy vf b w x w y b*

**3.2 Modeling flows in rotating ring zone - the equations that are applicable to rotating**

For this flow, NMWR solves conservation equations for mass and momentum describing

0 *j j u x*

 *j i i j i j i*

*uu uu f*

For flows in rotating domain (fig. 7), the equations for conservation of mass and momentum are written for the relative velocity formulation, where *fi* in right hand side is given by

The absolute velocity formulation is used in the non-rotating domain, and *fi=*0. Thus, the

where **w**– all weights and network thresholds vector that assure summary quadratic residuals minimums for each equation in optional totality of computation nodes coordinating each neural network solution learning iteration, are generated using the

Neuronet learning solution parameter search process efficiency depends on problem dimension, i.e. weights and perceptron thresholds varying adjusted quantity. The more significant is neurons quantity in trial solution, the higher is ANN approximate capacity;

*j j ij ji*

*fi*

*x x xx x x*

**w**

1 2

. (25)

. (26)

*i ii i p*

*p u u*

2 *u r* . (27)

*NN* **w** as equations solution,

3

*i q*

with weights calculated array and the boundary network **w** *wvb* , , .

incompressible flows of viscous Newtonian fluid.

turbulence model is used.

Let us consider the neural network functions ,,, ( ,,) *u pk f xy i*

however, achieving high approximation accuracy is more complicated.

*q*

2 1

solution to the equations (17)-(19) solution is also illustrated.

distribution on this layer

**reference frame**

Continuity equation

Momentum equations

standard *k*

random number generator.

#### Fig. 7. Rotating Domain

At the same time, neuron quantity depends not only on simulated function complexity, but also on calculation nodes quantity in which the residual equation is calculated. It is known that generally points' quantity increase in statistical set used for neural network construction is followed by increase in necessary neurons network (Galushkin, 2002; Galushkin, 2007) quantity. Consequently, the dense calculation grids application results in nonlinear programming problems; while applying rare calculation grids, it is necessary to check the solution realization between calculation nodes, i.e. there is a problem of learning solution procedure standardization. In the neuronet solution reception context on known equation, it is convenient using traditional additive parameter of training neural model quality - a control error which is calculated on the set of additional calculation nodes between calculation grid nodes. The number of these additional calculation grid nodes can be much more important and they should cover all calculation area because the nodes number increase in which the control error on known network parameters does not result in essential computing expenses growth. Hence, at the learning solution neuronet parameters reception there is a problem of solving the multi-criterion problem of the nonlinear optimization and to minimize simultaneously both the summary residual in control points, or the control error can appear as a restriction parameter, in the limited set of calculation nodes and in this case the neural network solution parameters reception is reduced to the conditional nonlinear optimization problem.

On figure 8, velocity distribution vector within the current random computation nodes of revolving ring on one of neural network training iteration of equation (25)-(26) decisions is presented. Contours of stream function within the computational zone obtained by using NMWR are presented on fig. 9.

Neural Network Modeling of Hydrodynamics Processes 217

Method of leakage zonal location (Zverev & Lurie, 2009), which generalizes the known method of leakage detection by hydraulic gradient line salient point, can be used to define of pipeline leakage position. When using the "base" zone location method variant, it is supposed that in case of stationary process, hydraulic gradients 1*i* and 2*i* before and after the leakage are constant and can be designed by the known values of liquid charge at the ends of controllable site. I.e., the distribution indignation of hydrodynamic parameters is not taken into account in comparison with the established flow which arises because of the environment outflow influence through the leakage aperture. It results in speed profile deformation within the cross sections of the pipeline near the leakage and the nonlinear dependence of hydraulic gradient in the function of distance from position of leakage downwards the stream. In order for pipe cross section speed distribution of the pipe to become appropriate for the established flow, the distance up to 40-60 calibres from leakage coordinate *у x* can be required. Defining nonlinear dependence of hydraulic gradient function near the leakage is possible either based on special experimental researches, or based on the numerical decision of hydrodynamics equations in three-dimensional statement. The given function can be added to the algorithm of

zone location methode to reduce the leakage position coordinate definition error.

8

*u*

*R* - design section average speed,

One pipeline section of 300 m length and diameter *D* 1 m is used for the nonlinear dependence of the total pressure drop determination along the pipeline length in the leakage neighborhood through the wall hole. The leakage position coordinate is fixed 155.5 m *у x* . On the segment *x* 140,200 m computational faces are formed with an interval *h m* 1 . In segments *x* 100,140 m and *x* 200,250 m computational faces are disposed with an interval *h* 10 m . Such a leakage location (approximately in the middle of the pipeline sector into a question) and computational faces disposition have been chosen for the influence elimination of the boundary conditions setting in the computational zone inlet and outlet on the flow distribution near leakage position. The leakage hole *D Dl* relative diameter during the computational planed experiment is constant. The mathematical model includes conservation equations for mass and momentum enclosed by

turbulence model. The local hydraulic friction coefficient in any section

, (27)

, (28)

2

*u*

*r R*

*u r*

Turbulent flow speed distribution within the pipe cross section is described by the universal

log *u u y A B*

**4. Modeling leakage in a fuel transfer pipeline**

the standard *k*

where 2

where

viscosity.

0 1 2 *R u urdr*

8

logarithmic dependence

*u u*

of the pipeline is defined by the ratio

Fig. 8. Velocity distribution on one of neural network training iteration

Fig. 9. The contours of stream function

216 Hydrodynamics – Optimizing Methods and Tools

Fig. 8. Velocity distribution on one of neural network training iteration

Fig. 9. The contours of stream function

#### **4. Modeling leakage in a fuel transfer pipeline**

Method of leakage zonal location (Zverev & Lurie, 2009), which generalizes the known method of leakage detection by hydraulic gradient line salient point, can be used to define of pipeline leakage position. When using the "base" zone location method variant, it is supposed that in case of stationary process, hydraulic gradients 1*i* and 2*i* before and after the leakage are constant and can be designed by the known values of liquid charge at the ends of controllable site. I.e., the distribution indignation of hydrodynamic parameters is not taken into account in comparison with the established flow which arises because of the environment outflow influence through the leakage aperture. It results in speed profile deformation within the cross sections of the pipeline near the leakage and the nonlinear dependence of hydraulic gradient in the function of distance from position of leakage downwards the stream. In order for pipe cross section speed distribution of the pipe to become appropriate for the established flow, the distance up to 40-60 calibres from leakage coordinate *у x* can be required. Defining nonlinear dependence of hydraulic gradient function near the leakage is possible either based on special experimental researches, or based on the numerical decision of hydrodynamics equations in three-dimensional statement. The given function can be added to the algorithm of zone location methode to reduce the leakage position coordinate definition error.

One pipeline section of 300 m length and diameter *D* 1 m is used for the nonlinear dependence of the total pressure drop determination along the pipeline length in the leakage neighborhood through the wall hole. The leakage position coordinate is fixed 155.5 m *у x* . On the segment *x* 140,200 m computational faces are formed with an interval *h m* 1 . In segments *x* 100,140 m and *x* 200,250 m computational faces are disposed with an interval *h* 10 m . Such a leakage location (approximately in the middle of the pipeline sector into a question) and computational faces disposition have been chosen for the influence elimination of the boundary conditions setting in the computational zone inlet and outlet on the flow distribution near leakage position. The leakage hole *D Dl* relative diameter during the computational planed experiment is constant. The mathematical model includes conservation equations for mass and momentum enclosed by the standard *k* turbulence model. The local hydraulic friction coefficient in any section of the pipeline is defined by the ratio

$$\lambda = -\left| \begin{matrix} 8\mu \frac{\partial u}{\partial r} \Big|\_{r=R} \\\\ \rho \overline{u}^2 \end{matrix} \right| \tag{27}$$

where 2 0 1 2 *R u urdr R* - design section average speed, - dynamic factor of liquid viscosity.

Turbulent flow speed distribution within the pipe cross section is described by the universal logarithmic dependence

$$\frac{\mu}{\mu\_{\star}} = A \log \frac{\mu\_{\ast} y}{\nu} + B \text{ .} \tag{28}$$

where 8 *u u* - the dynamic speed, *y* - the distance from the wall, - the kinematic viscosity.

Neural Network Modeling of Hydrodynamics Processes 219

*p g ix ix*

*gi i*

Let us present the total pressure drop nonlinear dependence for its calculation in leakage neighbourhood in nonlinear function form 1 2 , , *i i f iix x NN NN* on the pipeline sector *xx x l NN NN NN* , , where for our computational model the left area limit coordinate of the nonlinear leakage function determination *x NN* 150 m and the determination length fragment area *l NN* 60 m (leakage is situated in the point with coordinate *x у* 155,5 m ). In this case, total pressure drop in the controlled sector *x x* 1 2 , is calculated by the

> 0 1 1 22 1 2

In the known function 1 2 , , *i f iix x NN NN* , the present equation is nonlinear with one unknown value *xNN* . After solving this equation, we can define the pipeline leakage

Starting from above-mentioned, computational experiment is effectuated with the purpose of function determination 1 2 , , *i f iix x NN NN* , i.e. the total pressure drop dependence in the leakage neighbourhood *i* (criterion) of three variables 1 1

2 21 / *NN i ii* and *x xx NN NN* (factors). Earlier, the interval forfactor variation *xNN* 0,60 was determined. Let us determine the intervals 1 10,20 *NN i* , 2 0,5,0,95 *NN i* . The plan of computing experiment is made with the points received with the help of Sobol-Statnikov generator of quasiuniform number sequences (Statnikov & Matusov, 2002). The working hypercube of space **R**3 is filled with the points

algorithm for the experiment plan formation is made by high efficiency of the research method of parameter space based on the sounding of computational area by points of the uniform distributed sequence.For each variable, parameter combinations 1 2 , *NN NN i i* from the experiment plan limitation boundary conditions in the computational model are selected by determining values corresponding tothe inlet mass flow *min* and leakage intensity *ml* , where the hydraulic friction coefficient is calculated by the implicit formula

1 2,82 2,04log

As a result of computational hydrodynamics equations decision, we obtain the value distribution of total pressure drop on the length of the controlled pipeline sector, from which we are exterminating the criteria value *i* for the corresponding value *xNN* from the experiment plan. Let us now register the obtained vector 1 2 ,, , *NN NN ii xi NN* in database for

For the illustration of numerical calculations, the distribution of speed in the leakage neighbourhood is presented on Fig. 11, and on fig. 12 the distribution of the hydraulic

.

Re

the posterior generation of neural network dependence 1 2 , , *i NN NN if i i x NN NN* .

gradient on the pipeline sector *x* 100,250 m for factor values 1

2 21 / 0,5; 0,57; 0,64; 0,71; 0,8; 0,9 *NN i ii* is presented.

*NN NN NN NN NN p ix x ix x l f i i x x dx*

0

*lNN i*

. (30)

*у*

*x*

following formula

coordinate *l x* .

*g*

1 2 , , *NN NN ii x NN* according to *LP*

of Altschul for "smooth" pipes

 0 11 22 1 2

. (29)

, ,

algorithm (Statnikov & Matusov, 2002). Choice of this

*NN i gi* ,

Pa 14 m *NN i* ,

At the statement of boundary conditions on the input, mass charge *m*1 or the speed *u*1 appropriate to this charge is set. In the grid units belonging to the leakage aperture, liquid outflow speed *uу* is appointed. The given parameter at numerical researches varies

for the modeling leakages of various intensity. For any point of the calculating area, any operational value of pressure *x* of all the hydrodynamical parameters are equaled to zero, for the concerning of which differences of pressure will be calculated, is set. On the output from the calculating area the conditions of the established current are set, i.e. the derivatives modeling of the boundary layer the standard wall functions for parameters of turbulence *k* and ε (Fletcher, 1991) are used.

The equations discretization is effectuated on the base of finite volumes method in the combination with hexagonal grid. On fig. 10 the formation of the grid in the calculating sections of the pipeline is represented.

Fig. 10. Calculation grid in the pipeline sections

#### **4.1 The plan of computing experiment**

If 1 2 *x x*, are the coordinates of the beginning and the end of the controllable pipeline sector, then 1 2 *i i* , are hydraulic gradients before and after the leakage respectively, 0 *p* is total pressure drop of the sector, so in case of 1*i const* from the beginning of the sector to the leakage and 2*i const* from the leakage up to the end of the sector, in the stationary hydrodynamic mode, leakage coordinate shall be defined by the formula

218 Hydrodynamics – Optimizing Methods and Tools

At the statement of boundary conditions on the input, mass charge *m*1 or the speed *u*1 appropriate to this charge is set. In the grid units belonging to the leakage aperture, liquid outflow speed *uу* is appointed. The given parameter at numerical researches varies for the modeling leakages of various intensity. For any point of the calculating area, any operational value of pressure *x* of all the hydrodynamical parameters are equaled to zero, for the concerning of which differences of pressure will be calculated, is set. On the output from the calculating area the conditions of the established current are set, i.e. the derivatives modeling of the boundary layer the standard wall functions for parameters of

The equations discretization is effectuated on the base of finite volumes method in the combination with hexagonal grid. On fig. 10 the formation of the grid in the calculating

If 1 2 *x x*, are the coordinates of the beginning and the end of the controllable pipeline sector, then 1 2 *i i* , are hydraulic gradients before and after the leakage respectively, 0 *p* is total pressure drop of the sector, so in case of 1*i const* from the beginning of the sector to the leakage and 2*i const* from the leakage up to the end of the sector, in the stationary

hydrodynamic mode, leakage coordinate shall be defined by the formula

turbulence *k* and ε (Fletcher, 1991) are used.

Fig. 10. Calculation grid in the pipeline sections

**4.1 The plan of computing experiment**

sections of the pipeline is represented.

$$\propto\_{y} = \frac{\Delta p\_0 + \rho g \left(i\_1 x\_1 - i\_2 x\_2\right)}{\rho g \left(i\_1 - i\_2\right)}.\tag{29}$$

Let us present the total pressure drop nonlinear dependence for its calculation in leakage neighbourhood in nonlinear function form 1 2 , , *i i f iix x NN NN* on the pipeline sector *xx x l NN NN NN* , , where for our computational model the left area limit coordinate of the nonlinear leakage function determination *x NN* 150 m and the determination length fragment area *l NN* 60 m (leakage is situated in the point with coordinate *x у* 155,5 m ). In this case, total pressure drop in the controlled sector *x x* 1 2 , is calculated by the following formula

$$\frac{\Delta p\_0}{\rho \mathcal{R}} = \mathbf{i}\_1 \left( \mathbf{x}\_{\text{NN}} - \mathbf{x}\_1 \right) + \mathbf{i}\_2 \left( \mathbf{x}\_2 - \left( \mathbf{x}\_{\text{NN}} + l\_{\text{NN}} \right) \right) + \int\_0^{l\_{\text{NN}}} f\_{\text{NN}}^i \left( \mathbf{i}\_1, \mathbf{i}\_2, \mathbf{x} - \mathbf{x}\_{\text{NN}} \right) d\mathbf{x} \, \tag{30}$$

In the known function 1 2 , , *i f iix x NN NN* , the present equation is nonlinear with one unknown value *xNN* . After solving this equation, we can define the pipeline leakage coordinate *l x* .

Starting from above-mentioned, computational experiment is effectuated with the purpose of function determination 1 2 , , *i f iix x NN NN* , i.e. the total pressure drop dependence in the leakage neighbourhood *i* (criterion) of three variables 1 1 *NN i gi* , 2 21 / *NN i ii* and *x xx NN NN* (factors). Earlier, the interval forfactor variation *xNN* 0,60 was determined. Let us determine the intervals 1 10,20 *NN i* , 2 0,5,0,95 *NN i* . The plan of computing experiment is made with the points received with the help of Sobol-Statnikov generator of quasiuniform number sequences (Statnikov & Matusov, 2002). The working hypercube of space **R**3 is filled with the points 1 2 , , *NN NN ii x NN* according to *LP* algorithm (Statnikov & Matusov, 2002). Choice of this algorithm for the experiment plan formation is made by high efficiency of the research method of parameter space based on the sounding of computational area by points of the uniform distributed sequence.For each variable, parameter combinations 1 2 , *NN NN i i* from the experiment plan limitation boundary conditions in the computational model are selected by determining values corresponding tothe inlet mass flow *min* and leakage intensity *ml* , where the hydraulic friction coefficient is calculated by the implicit formula of Altschul for "smooth" pipes

$$\frac{1}{\sqrt{\lambda}} = -2,04\log\left(\frac{2,82}{\text{Re}\sqrt{\lambda}}\right)\cdot 1$$

As a result of computational hydrodynamics equations decision, we obtain the value distribution of total pressure drop on the length of the controlled pipeline sector, from which we are exterminating the criteria value *i* for the corresponding value *xNN* from the experiment plan. Let us now register the obtained vector 1 2 ,, , *NN NN ii xi NN* in database for the posterior generation of neural network dependence 1 2 , , *i NN NN if i i x NN NN* .

For the illustration of numerical calculations, the distribution of speed in the leakage neighbourhood is presented on Fig. 11, and on fig. 12 the distribution of the hydraulic gradient on the pipeline sector *x* 100,250 m for factor values 1 Pa 14 m *NN i* , 2 21 / 0,5; 0,57; 0,64; 0,71; 0,8; 0,9 *NN i ii* is presented.

Neural Network Modeling of Hydrodynamics Processes 221

For the regressive dependence construction, 1 2 , , *i NN NN if i i x NN NN* artificial neural

networks device is used. For the formation of display 1 2 , , *i NN NN if i i x NN NN* , the standard structure of the multilayer perceptron (MLP) with 3 inputs, one output and two latent layers with 7 and 5 neurons accordingly was used. At training, MLP algorithm of Levenberg-Markardt was used. The total RMS error on 512 points of statistical sample was *E=0.001*. The comparative analysis of dependences calculation results obtained for inlet parameters

*i* is presented on Fig. 13 and 1

 *xNN* 0,60 m from the differential equation computational solution– continuous lines and determined with neural network dependence – markers. It is also necessary to note the

Fig. 13. Comparative analysis of neural network data and numerical solution results

depending on learning solution local accuracy is offered.

Results presented in this chapter concern one of many applications of neural network learning functions to the mathematical physics problems solution, which is solution of hydrodynamics equations. Using learning functions of neural network allows to exclude the errors in solving differential equations caused by derivatives discretization and the borders low precision representation. Numerical algorithm of Navier-Stokes equations system solution differing from known fluid dynamics methods by the possibility of computational grid arbitrary formationis developed. The formation technique of the additional nodes file

Pa 10 m

1 0.8 *i i* for

*xNN*

*gi* , 2

**4.2 The neural network regressive model**

Pa 16 m

1 0.9 *i*

*gi* , 2

high approximation precision.

values 1

1 *i i*

**5. Conclusion**

Fig. 11. Speed distribution in the leakage neighbourhood

Fig. 12. Total pressure drop distribution for 1 meter of pipeline length

#### **4.2 The neural network regressive model**

220 Hydrodynamics – Optimizing Methods and Tools

Fig. 11. Speed distribution in the leakage neighbourhood

*gi*

Fig. 12. Total pressure drop distribution for 1 meter of pipeline length

*x, m*

For the regressive dependence construction, 1 2 , , *i NN NN if i i x NN NN* artificial neural networks device is used. For the formation of display 1 2 , , *i NN NN if i i x NN NN* , the standard structure of the multilayer perceptron (MLP) with 3 inputs, one output and two latent layers with 7 and 5 neurons accordingly was used. At training, MLP algorithm of Levenberg-Markardt was used. The total RMS error on 512 points of statistical sample was *E=0.001*. The comparative analysis of dependences calculation results obtained for inlet parameters values 1 Pa 16 m *gi* , 2 1 0.9 *i i* is presented on Fig. 13 and 1 Pa 10 m *gi* , 2 1 0.8 *i i* for

 *xNN* 0,60 m from the differential equation computational solution– continuous lines and determined with neural network dependence – markers. It is also necessary to note the high approximation precision.

Fig. 13. Comparative analysis of neural network data and numerical solution results

#### **5. Conclusion**

Results presented in this chapter concern one of many applications of neural network learning functions to the mathematical physics problems solution, which is solution of hydrodynamics equations. Using learning functions of neural network allows to exclude the errors in solving differential equations caused by derivatives discretization and the borders low precision representation. Numerical algorithm of Navier-Stokes equations system solution differing from known fluid dynamics methods by the possibility of computational grid arbitrary formationis developed. The formation technique of the additional nodes file depending on learning solution local accuracy is offered.

**Part 3**

**Complex Hydraulic Engineering Applications**

The formation of the finite-elemental models for all extended pipelines and the numerical decision of equations of liquid movement for modeling the hydrodynamics processes now is limited to computer resources as by the quantity of the final elements of model calculated grid, so as by the time of finding of the decision. Even with such model, its use for system functioning operative analysis is quite doubtful. However, the numerical models of hydrodynamics processes can be used for the formation of information databases based on the neural network computing architecture which, after "training", have high speed of calculation. Neural networks are the universal approximation tool of multivariate nonlinear dependences, capable "to be arranged" under appearing of the new information of the researched process, i.e. they can serve as an intellectual tool for monitoring, which is constantly filled up and clarified. Thus, the introduction of neural algorithms in the approved and well-recommending methods of pipeline leakage detection is expedient for increased accuracy and the efficiency of accepted decisions.

#### **6. References**

Bishop, C.M. (1995) *Neural networks for pattern recognition*, Oxford University Press

## **Part 3**

## **Complex Hydraulic Engineering Applications**

222 Hydrodynamics – Optimizing Methods and Tools

The formation of the finite-elemental models for all extended pipelines and the numerical decision of equations of liquid movement for modeling the hydrodynamics processes now is limited to computer resources as by the quantity of the final elements of model calculated grid, so as by the time of finding of the decision. Even with such model, its use for system functioning operative analysis is quite doubtful. However, the numerical models of hydrodynamics processes can be used for the formation of information databases based on the neural network computing architecture which, after "training", have high speed of calculation. Neural networks are the universal approximation tool of multivariate nonlinear dependences, capable "to be arranged" under appearing of the new information of the researched process, i.e. they can serve as an intellectual tool for monitoring, which is constantly filled up and clarified. Thus, the introduction of neural algorithms in the approved and well-recommending methods of pipeline leakage detection is expedient for

increased accuracy and the efficiency of accepted decisions.

Galushkin, A.I. (2002). *Neuromathematics,* IPRJR, Moscow Galushkin, A.I. (2007). *Neural networks theory,* Springer

1-4244-6916-1, Barcelona, Spain, 2010.

[*Russian Aeronautics* (Engl.Transl.)]

*Oil industry*, No 8, 2009, pp. 85-87

Academic Publishers, Dordrecht/ Boston / London

Bishop, C.M. (1995) *Neural networks for pattern recognition*, Oxford University Press Cloete, L. & Zurada, J.M. (2000) *Knowledge-Based Neurocomputing,* MIT Press

*Sciences,* Vol. 359, No. 3, 1998, pp. 330-333. ISSN 0869-5652

Egorov, I. N.; Kretinin, G.V.; Matusov I.B. & Statnikov, R.B. Problems of design and

Fletcher, K. (1991). *Computational techniques for fluid dynamics,* Mir, Moscow, ISBN 5-03-

Kretinin, A.V. The weighted residuals method based on neuronet approximations for

Kretinin, A.V.; Bulygin, Yu. A. & Kirpichev, M.I. Method of Weighted Residuals on the Base

Kretinin, A.V.; Bulygin, Yu. A.; Kirpichev, M.I. & Darneva, T.I. Neuronet's method of

Prokhorov, D.; Feldkamp, L. & Feldkamp, T. A new approach to cluster weighted modeling.

Statnikov, R.B. & Matusov, J.B. (1995). *Multicriteria Optimization and Engineering,* Chapman &

Statnikov, R.B. & Matusov, J.B. (2002). *Multicriteria Analysis in Engineering,* Kluwer

Stogney, V.G. & Kretinin*, A.V.* Simulation of Flows in a Channel with Permeable Wall Using

Zverev, F.C. & Lurie, M.V. Generalized zonal location method of pipeline leakage detection,

*Cybernetics ICCC 2008*. pp. 237-240, Stara Lesná, Slovakia, 2008.

multicriteria control of controlled engineering systems. *J. Reports of RF Academy of*

simulation of hydrodynamics problems, *Siberian J. Num. Math*, 2006, Vol. 9, No.1,

of Neuronet's Approximations for Computer Simulation of Hydrodynamics Problems, *Proceedings of IEEE 6th International Conference on Computational*

weighted residuals for computer simulation of hydrodynamics problems, Proceedings of 2010 International Joint Conference on Neural Networks. ISBN 978-

*Proceedings of International Joint Conference on Neural Networks (IJCNN),* Washington

Artificial Neural Networks. *Izv.Vuz. Av. Tekhnika*, Vol. 48, No. 1, 2005, pp. 34–38

**6. References**

001881-6

pp. 23-35.

DC, July 2001

Hall, N.Y.

**11**

*Argentina*

**Interaction Between Hydraulic and Numerical**

**Models for the Design of Hydraulic Structures**

The design of complex hydraulic structures requires its testing through hydraulic models (i.e., reduced scale physical representations). The main practical limitation of hydraulic models are the so called 'scale effects', i.e., the fact that only the primary physical mechanisms can be correctly represented, while the secondary ones are distorted. In particular, for free surface flows the gravitational driving forces – primary mechanism – must be correctly scaled in relation to inertia (Froude scaling), leading to an incorrect representation of viscous forces (no Reynolds scaling) – usually the leading secondary mechanism – as the fluid in the hydraulic model is the same as in the prototype (water). Though for most applications Reynolds number effects introduce only small quantitative deviations, which can be readily absorbed within the margin of safety assumed for design, this is not always the case. In fact, they can for example accumulate, in such a way that the effects compete with those arising from the primary mechanism. In those cases, being the Reynolds effects distorted in the hydraulic model, the observed response deviates from the

Numerical modeling is the appropriate tool to help solving in a rigorous way this type of difficulty. A sound numerical model should be able to correctly represent both the primary and secondary mechanisms, i.e., it is not subject to 'scale effects'. Its main limitations might arise from insufficient resolution, or from inaccurate representation of turbulence effects. The first limitation could be overcome by reducing the spatial step of the numerical grid; the

Based on these observations, the following strategy is proposed: (i) the flow in the hydraulic model is numerically simulated, i.e., the dimensions of the hydraulic model are used (thus accounting for the 'spurious' scale effects); this constitutes a way of validating the theoretical model; eventually, adjustments in the representation (higher resolution, more elaborated theoretical approaches) are introduced in order to improve the comparison; (ii) the flow in the prototype is numerically simulated, by introducing the dimensions of the prototype in the validated numerical model (i.e., distortion of secondary mechanisms is now avoided); this

Two problems (with quite different levels of complexity) are presented as case studies in order to illustrate the proposed approach, both of them associated to the design of the Third Set of Locks of the Panama Canal (communicating the Atlantic and Pacific Oceans), for which the present authors were responsible: (a) the determination of the time for water level

constitutes the adequate way of extrapolating the results to the prototype dimensions.

one corresponding to the prototype, thus needing some empirical correction.

second one, by resorting to more elaborated theoretical approaches.

**1. Introduction**

Angel N. Menéndez and Nicolás D. Badano *INA (National Institute for Water), Hydraulics Laboratory*

## **Interaction Between Hydraulic and Numerical Models for the Design of Hydraulic Structures**

Angel N. Menéndez and Nicolás D. Badano *INA (National Institute for Water), Hydraulics Laboratory Argentina*

#### **1. Introduction**

The design of complex hydraulic structures requires its testing through hydraulic models (i.e., reduced scale physical representations). The main practical limitation of hydraulic models are the so called 'scale effects', i.e., the fact that only the primary physical mechanisms can be correctly represented, while the secondary ones are distorted. In particular, for free surface flows the gravitational driving forces – primary mechanism – must be correctly scaled in relation to inertia (Froude scaling), leading to an incorrect representation of viscous forces (no Reynolds scaling) – usually the leading secondary mechanism – as the fluid in the hydraulic model is the same as in the prototype (water). Though for most applications Reynolds number effects introduce only small quantitative deviations, which can be readily absorbed within the margin of safety assumed for design, this is not always the case. In fact, they can for example accumulate, in such a way that the effects compete with those arising from the primary mechanism. In those cases, being the Reynolds effects distorted in the hydraulic model, the observed response deviates from the one corresponding to the prototype, thus needing some empirical correction.

Numerical modeling is the appropriate tool to help solving in a rigorous way this type of difficulty. A sound numerical model should be able to correctly represent both the primary and secondary mechanisms, i.e., it is not subject to 'scale effects'. Its main limitations might arise from insufficient resolution, or from inaccurate representation of turbulence effects. The first limitation could be overcome by reducing the spatial step of the numerical grid; the second one, by resorting to more elaborated theoretical approaches.

Based on these observations, the following strategy is proposed: (i) the flow in the hydraulic model is numerically simulated, i.e., the dimensions of the hydraulic model are used (thus accounting for the 'spurious' scale effects); this constitutes a way of validating the theoretical model; eventually, adjustments in the representation (higher resolution, more elaborated theoretical approaches) are introduced in order to improve the comparison; (ii) the flow in the prototype is numerically simulated, by introducing the dimensions of the prototype in the validated numerical model (i.e., distortion of secondary mechanisms is now avoided); this constitutes the adequate way of extrapolating the results to the prototype dimensions.

Two problems (with quite different levels of complexity) are presented as case studies in order to illustrate the proposed approach, both of them associated to the design of the Third Set of Locks of the Panama Canal (communicating the Atlantic and Pacific Oceans), for which the present authors were responsible: (a) the determination of the time for water level

Interaction Between Hydraulic and Numerical Models for the Design of Hydraulic Structures 227

In order to allow an even filling or emptying of the chambers, thus minimizing longitudinal water surface slopes and consequently hawser forces, water enters or exits each lock chamber through 20 ports located on each lateral wall. They are connected to secondary culverts which, in turn, connect to the main culvert at the midpoint of each chamber, through a carefully designed hydraulic component called Central Connection (CC) (Fig. 2.2). The CC was designed so as the flow coming from the main culvert is evenly split between both secondary culverts, before filling the chamber through the ports. Additionally, during the opposite operation, i.e., chamber emptying, equal discharges should flow through both secondary culverts, before its confluence towards the main culvert. The symmetry of the flow with respect to the chamber midpoint, for both the filling and emptying operations, is what ultimately guarantees low longitudinal water slopes and, consequently, low

The two branches of the CC are called 'U' and 'S' branches, in reference to the trajectory followed by the incoming flow (Fig. 2.2). Note that the role as a U or a S branch depends on

A conduit arising from the midpoint of each secondary culvert, and running below the chambers, connects them with the WSBs, after traversing a 'trifurcation' and a set of valves

The hydrodynamic model of the F/E system was built using Flowmaster V7 (http://www.flowmaster.com/), a commercial code which solves one-dimensional transient flow over a network of conduits. Incompressible fluid and rigid pipe hypothesis were made, without compromising accuracy. For convenience, the whole F/E system was divided into sub models, one for each type of operation, such as equalization of the Upper Chamber with the Lake, equalization of two contiguous chambers, equalization of chamber

longitudinal hawser forces.

Fig. 2.2. Detail of the F/E system and Central Connection

the hydraulic operation (filling or emptying).

**2.2 Description of numerical model**

with its Water Saving Basins, and so on.

(Fig. 2.2).

equalization between chambers, for which a one-dimensional numerical model was used; (b) the calculation of the amplitude of free surface oscillations in the lock chambers (which leads to increments in the hawser forces) due to close-to-resonance conditions under interaction with an oscillation in a flow partition component of the filling/emptying system (triggered by large turbulent eddies), for which a full three-dimensional numerical model – i.e., a CFD approach – was applied.

#### **2. Equalization times**

The equalization time, named as Filling/Emptying (F/E) time during the study, is a key parameter in establishing the system performance of a Lock Complex, as it has a direct impact on the vessel throughput, measured as the number of vessels passing through the system per day. Contractual requirements existed for the design of the Third Set of Locks of the Panama Canal, imposing maximum allowable F/E times for different scenarios. Minimization of these F/E times, through the reduction of local head losses, was the main strategy used during the design optimization process. Consequently, scale effects affecting these F/E times were carefully studied.

#### **2.1 Description of the physical system**

The Third Set of Locks of the Panama Canal, presently under construction, comprises twin lock complexes located near each ocean. Each complex has three lock chambers in series (Fig. 2.1). These lock complexes allow vessels to be transported up or down between Gatun Lake and the oceans, spanning a 27 m water level difference. Each lock chamber has three side pools, called Water Saving Basins (WSB). This WSBs store part of the water used during the equalization operations, that otherwise would be flushed downstream towards the ocean. This stored water is then utilized to refill part of the lock, allowing a reduction of freshwater consumption during dry hydrological seasons.

The lock chambers are connected through two longitudinal main culverts located within the lateral walls, running from the lake to the ocean. Four sets of valves are located along these culverts, which can isolate each chamber (Fig. 2.1). The operation of these valves allows the successive equalization of the water level between chambers by gravity flow.

226 Hydrodynamics – Optimizing Methods and Tools

equalization between chambers, for which a one-dimensional numerical model was used; (b) the calculation of the amplitude of free surface oscillations in the lock chambers (which leads to increments in the hawser forces) due to close-to-resonance conditions under interaction with an oscillation in a flow partition component of the filling/emptying system (triggered by large turbulent eddies), for which a full three-dimensional numerical model –

The equalization time, named as Filling/Emptying (F/E) time during the study, is a key parameter in establishing the system performance of a Lock Complex, as it has a direct impact on the vessel throughput, measured as the number of vessels passing through the system per day. Contractual requirements existed for the design of the Third Set of Locks of the Panama Canal, imposing maximum allowable F/E times for different scenarios. Minimization of these F/E times, through the reduction of local head losses, was the main strategy used during the design optimization process. Consequently, scale effects affecting

The Third Set of Locks of the Panama Canal, presently under construction, comprises twin lock complexes located near each ocean. Each complex has three lock chambers in series (Fig. 2.1). These lock complexes allow vessels to be transported up or down between Gatun Lake and the oceans, spanning a 27 m water level difference. Each lock chamber has three side pools, called Water Saving Basins (WSB). This WSBs store part of the water used during the equalization operations, that otherwise would be flushed downstream towards the ocean. This stored water is then utilized to refill part of the lock, allowing a reduction of

The lock chambers are connected through two longitudinal main culverts located within the lateral walls, running from the lake to the ocean. Four sets of valves are located along these culverts, which can isolate each chamber (Fig. 2.1). The operation of these valves allows the

successive equalization of the water level between chambers by gravity flow.

Fig. 2.1. Panama Canal Third Set of Locks F/E System – Bottom up view

i.e., a CFD approach – was applied.

these F/E times were carefully studied.

**2.1 Description of the physical system**

freshwater consumption during dry hydrological seasons.

**2. Equalization times**

In order to allow an even filling or emptying of the chambers, thus minimizing longitudinal water surface slopes and consequently hawser forces, water enters or exits each lock chamber through 20 ports located on each lateral wall. They are connected to secondary culverts which, in turn, connect to the main culvert at the midpoint of each chamber, through a carefully designed hydraulic component called Central Connection (CC) (Fig. 2.2). The CC was designed so as the flow coming from the main culvert is evenly split between both secondary culverts, before filling the chamber through the ports. Additionally, during the opposite operation, i.e., chamber emptying, equal discharges should flow through both secondary culverts, before its confluence towards the main culvert. The symmetry of the flow with respect to the chamber midpoint, for both the filling and emptying operations, is what ultimately guarantees low longitudinal water slopes and, consequently, low longitudinal hawser forces.

Fig. 2.2. Detail of the F/E system and Central Connection

The two branches of the CC are called 'U' and 'S' branches, in reference to the trajectory followed by the incoming flow (Fig. 2.2). Note that the role as a U or a S branch depends on the hydraulic operation (filling or emptying).

A conduit arising from the midpoint of each secondary culvert, and running below the chambers, connects them with the WSBs, after traversing a 'trifurcation' and a set of valves (Fig. 2.2).

#### **2.2 Description of numerical model**

The hydrodynamic model of the F/E system was built using Flowmaster V7 (http://www.flowmaster.com/), a commercial code which solves one-dimensional transient flow over a network of conduits. Incompressible fluid and rigid pipe hypothesis were made, without compromising accuracy. For convenience, the whole F/E system was divided into sub models, one for each type of operation, such as equalization of the Upper Chamber with the Lake, equalization of two contiguous chambers, equalization of chamber with its Water Saving Basins, and so on.

Interaction Between Hydraulic and Numerical Models for the Design of Hydraulic Structures 229

0 60 120 180 240 300 360 420 480 540 600 660 720

**Time (secs)** Physical Model - Far Side Mathematical Model - Far Side Valve Aperture

Fig. 2.3. Comparison of physical and numerical models: Discharge in the Main Culvert, for a

0 60 120 180 240 300 360 420 480 540 600 660 720

**Time (secs)**

Fig. 2.4. Comparison of physical and numerical models: Water levels in the chambers, for a

Physical Model Mathematical Model Valve Aperture

Lock to Lock operation with 21 m of initial head difference.

Lock to Lock operation with 21 m of initial head difference.

**Water Level (mPLD)**

**Main Culvert Dischage (m3/s)**

0%

0%

20%

40%

60%

**Valve Aperture (%)**

80%

100%

25%

50%

**Valve Aperture (%)**

75%

100%

The F/E system was represented as a network of interconnected component elements, namely:

#### **2.3 Numerical modeling of physical model**

The Third Set of Locks has been subject to physical modeling, both during the development of the conceptual design, and later during the design for the final project. Both physical models where commissioned to the Compagnie Nationale du Rhône (CNR), Lyon, France.

The physical models were built at a 1/30 scale, comprising 2 chambers and one set of three WSBs. Extensive tests were made for various normal and special operations, measuring water levels, discharges, pressures and water slopes in the chambers. Some tests included the presence of a design vessel model, measuring hydraulic longitudinal and transversal forces over its hull. Based on these tests, a correlation between forces on the ship, and water surface slopes in the chamber in the absence of the ship (easier to measure and allegedly more repeatable), was established. This correlation was used to impose maximum values to the longitudinal and lateral water surface slopes, as contractual requirements.

The flow in the hydraulic model was numerically simulated. Real physical dimensions of the physical model components (culverts, conduits, chambers) were used.

Local head loss coefficients for the special hydraulic components were obtained through steady CFD modeling (see Section 3 for more details on CFD modeling), by calculating the difference between upstream and downstream mechanical energy, and subtracting energy losses due to wall friction. Most parts of the physical model were made out of acrylic (with a 0.025 mm roughness height), which behaves as a hydraulically smooth surface, for which the roughness height is completely submerged within the viscous sublayer (White, 1974). Some of the special hydraulic components, though, were built with Styrofoam (enclosed inside of acrylic boxes), as the initial expectations were that many alternative geometries would have to be tested, so this system would allow swapping with relative ease (very few alternatives were finally tested, due the great success of the optimization process carried out with CFD models). As it was later demonstrated that Styrofoam behaves as hydraulically rough at the physical model scale, most of it had to be coated with a low roughness layer of paint in order to avoid a spurious response (a scale effect in itself).

The results obtained with the numerical model (water levels, discharges, pressures) showed a very good agreement with physical model measurements, for different operations and conditions. As an illustration, Figs 2.3 and 2.4 show comparisons for a typical Lock to Lock operation, with maximum initial head difference. All comparisons are presented with results scaled up to prototype dimensions.

228 Hydrodynamics – Optimizing Methods and Tools

The F/E system was represented as a network of interconnected component elements, namely: Reservoirs, representing the lock chambers, WSBs, the lake, and the oceans. Level-Area relations were specified for each one of them. The lake and oceans were considered as

 Rigid rectangular pipes, representing primary and secondary culverts, WSB conduits, etc. The calculation of friction losses was made using Darcy-Weisbach and Colebrook-White equations, as a function of the flow Reynolds number and the effective roughness

 Local energy losses parameterized with a cross-section area and a head loss coefficient, representing most of the special hydraulic components, such as bends, bifurcations,

 Local energy losses expressed as laws for cross-section area and head loss coefficient in terms of a control parameter, representing valves for which the control parameter is the

The Third Set of Locks has been subject to physical modeling, both during the development of the conceptual design, and later during the design for the final project. Both physical models where commissioned to the Compagnie Nationale du Rhône (CNR), Lyon, France. The physical models were built at a 1/30 scale, comprising 2 chambers and one set of three WSBs. Extensive tests were made for various normal and special operations, measuring water levels, discharges, pressures and water slopes in the chambers. Some tests included the presence of a design vessel model, measuring hydraulic longitudinal and transversal forces over its hull. Based on these tests, a correlation between forces on the ship, and water surface slopes in the chamber in the absence of the ship (easier to measure and allegedly more repeatable), was established. This correlation was used to impose maximum values to

The flow in the hydraulic model was numerically simulated. Real physical dimensions of

Local head loss coefficients for the special hydraulic components were obtained through steady CFD modeling (see Section 3 for more details on CFD modeling), by calculating the difference between upstream and downstream mechanical energy, and subtracting energy losses due to wall friction. Most parts of the physical model were made out of acrylic (with a 0.025 mm roughness height), which behaves as a hydraulically smooth surface, for which the roughness height is completely submerged within the viscous sublayer (White, 1974). Some of the special hydraulic components, though, were built with Styrofoam (enclosed inside of acrylic boxes), as the initial expectations were that many alternative geometries would have to be tested, so this system would allow swapping with relative ease (very few alternatives were finally tested, due the great success of the optimization process carried out with CFD models). As it was later demonstrated that Styrofoam behaves as hydraulically rough at the physical model scale, most of it had to be coated with a low roughness layer of

The results obtained with the numerical model (water levels, discharges, pressures) showed a very good agreement with physical model measurements, for different operations and conditions. As an illustration, Figs 2.3 and 2.4 show comparisons for a typical Lock to Lock operation, with maximum initial head difference. All comparisons are presented with

the longitudinal and lateral water surface slopes, as contractual requirements.

the physical model components (culverts, conduits, chambers) were used.

paint in order to avoid a spurious response (a scale effect in itself).

results scaled up to prototype dimensions.

infinite area constant level reservoirs.

**2.3 Numerical modeling of physical model**

height of the conduit walls.

transitions, etc.

aperture.

Fig. 2.3. Comparison of physical and numerical models: Discharge in the Main Culvert, for a Lock to Lock operation with 21 m of initial head difference.

Fig. 2.4. Comparison of physical and numerical models: Water levels in the chambers, for a Lock to Lock operation with 21 m of initial head difference.

Interaction Between Hydraulic and Numerical Models for the Design of Hydraulic Structures 231

For example, the Reynolds number in the primary culvert (in which most of the friction losses are produced) changes in time following the flow hydrograph, from zero to the peak discharge, and back to zero again. The peak discharge for 21 m initial head difference in a Lock to Lock operation is around 425 m3/s (the corresponding flow velocity is 7.87 m/s). This leads to a Reynolds number of around 6.5 107 for the prototype. When scaled to the physical model, the Reynolds number is only 3.9 105, i.e., a drop of more than two orders of magnitude. The associated friction coefficients are then below 0.008 for the prototype, and about 0.014 for the physical model. The consequently higher friction losses produced in the physical model, exclusively due to scale effects, reduce the flow velocities, then increasing the F/E times. The numerical model contemplates the variation of frictional losses with the Reynolds number. Hence, it allows to be used in order to extrapolate the physical model results to those expected for the prototype, overcoming the distortion introduced by scale effects in

For the Panama Canal Third Set of Lock, the validated 1D model was scaled up to prototype dimensions. Variations in local head loss coefficients, indicated by 3D models, were also introduced. Relatively little effects were observed in the simulations because of the change in local head loss coefficients. On the contrary, friction losses decreased significantly, as already explained. Consequently, for a typical Lock to Lock operation with maximum initial

0 60 120 180 240 300 360 420 480 540 600 660 720

**Time (secs)**

Fig. 2.6. Comparison of physical model scale and prototype scale numerical models: Water levels in the chambers, for a Lock to Lock operation with 21 m of initial head difference.

Additionally, a 5% increase in the peak discharge of the main culverts was also observed (Fig. 2.7). This has an effect over the pressures on the vena contracta, downstream of the main culvert valves (Fig. 2.8), which had to be contemplated during the design stage, as air intrusion had to be avoided (for contractual reasons), and because piezometric levels downstream of the valves were close to the roof level of the culvert for various special operating conditions. So avoiding scale effects was also significant to correctly deal with

Prototype Scale Physical Model Scale Valve Aperture Series4

0%

20%

40%

60%

**Valve Aperture (%)**

80%

100%

head difference, F/E times showed a 10% decrease (61 seconds) (Fig. 2.6).

the physical model results.

these two limitations.

**Water Level (mPLD)**

#### **2.4 Numerical modeling of the prototype**

Practical knowledge exists about the discrepancies between F/E times as measured in a physical model and those effectively occurring at the prototype. For instance, USACE manual on hydraulic design of navigation locks (2006) states:

..."A prototype lock filling-and-emptying system is normally more efficient than predicted by its model"..."The difference in efficiency is acceptable as far as most of the modeled quantities are concerned (hawser forces, for example) and can be accommodated empirically for others (filling time and over travel, specifically)."...

In the specific commentaries about F/E times, it suggests quantitative corrections:

..."General guidance is that the operation time with rapid valving should be reduced from the model values by about 10 percent for small locks (600 ft or less) with short culverts; about 15 percent for small locks with longer, more complex culvert systems; and about 20 percent for small locks (Lower Granite, for example) or large locks having extremely long culvert systems."...

The alternative, rigorous strategy proposed in the present paper is to numerically simulate the flow in the prototype. This means using the physical dimensions of the prototype, the corresponding local head loss coefficients for the special hydraulic components, and the roughness height for concrete. Though the concrete wall also behaves as hydraulically smooth, the friction coefficient for smooth pipes is a function of the flow Reynolds number, as indicated by the "smooth pipe" curve in the Moody chart (Fig. 2.5).

Fig. 2.5. Friction coefficient as a function of Reynolds number (Moody chart)

230 Hydrodynamics – Optimizing Methods and Tools

Practical knowledge exists about the discrepancies between F/E times as measured in a physical model and those effectively occurring at the prototype. For instance, USACE

In the specific commentaries about F/E times, it suggests quantitative corrections:

as indicated by the "smooth pipe" curve in the Moody chart (Fig. 2.5).

Fig. 2.5. Friction coefficient as a function of Reynolds number (Moody chart)

..."A prototype lock filling-and-emptying system is normally more efficient than predicted by its model"..."The difference in efficiency is acceptable as far as most of the modeled quantities are concerned (hawser forces, for example) and can be accommodated empirically for others (filling time and over travel, specifically)."...

..."General guidance is that the operation time with rapid valving should be reduced from the model values by about 10 percent for small locks (600 ft or less) with short culverts; about 15 percent for small locks with longer, more complex culvert systems; and about 20 percent for small locks (Lower Granite, for example) or large locks having

The alternative, rigorous strategy proposed in the present paper is to numerically simulate the flow in the prototype. This means using the physical dimensions of the prototype, the corresponding local head loss coefficients for the special hydraulic components, and the roughness height for concrete. Though the concrete wall also behaves as hydraulically smooth, the friction coefficient for smooth pipes is a function of the flow Reynolds number,

**2.4 Numerical modeling of the prototype**

extremely long culvert systems."...

manual on hydraulic design of navigation locks (2006) states:

For example, the Reynolds number in the primary culvert (in which most of the friction losses are produced) changes in time following the flow hydrograph, from zero to the peak discharge, and back to zero again. The peak discharge for 21 m initial head difference in a Lock to Lock operation is around 425 m3/s (the corresponding flow velocity is 7.87 m/s). This leads to a Reynolds number of around 6.5 107 for the prototype. When scaled to the physical model, the Reynolds number is only 3.9 105, i.e., a drop of more than two orders of magnitude. The associated friction coefficients are then below 0.008 for the prototype, and about 0.014 for the physical model. The consequently higher friction losses produced in the physical model, exclusively due to scale effects, reduce the flow velocities, then increasing the F/E times.

The numerical model contemplates the variation of frictional losses with the Reynolds number. Hence, it allows to be used in order to extrapolate the physical model results to those expected for the prototype, overcoming the distortion introduced by scale effects in the physical model results.

For the Panama Canal Third Set of Lock, the validated 1D model was scaled up to prototype dimensions. Variations in local head loss coefficients, indicated by 3D models, were also introduced. Relatively little effects were observed in the simulations because of the change in local head loss coefficients. On the contrary, friction losses decreased significantly, as already explained. Consequently, for a typical Lock to Lock operation with maximum initial head difference, F/E times showed a 10% decrease (61 seconds) (Fig. 2.6).

Fig. 2.6. Comparison of physical model scale and prototype scale numerical models: Water levels in the chambers, for a Lock to Lock operation with 21 m of initial head difference.

Additionally, a 5% increase in the peak discharge of the main culverts was also observed (Fig. 2.7). This has an effect over the pressures on the vena contracta, downstream of the main culvert valves (Fig. 2.8), which had to be contemplated during the design stage, as air intrusion had to be avoided (for contractual reasons), and because piezometric levels downstream of the valves were close to the roof level of the culvert for various special operating conditions. So avoiding scale effects was also significant to correctly deal with these two limitations.

Interaction Between Hydraulic and Numerical Models for the Design of Hydraulic Structures 233

Free surface oscillations in the lock chambers leads to forces in the hawsers. Based on results from the physical model constructed during the development of the conceptual design, a correlation was found between these forces and the free surface slope in the absence of the vessel, as already mentioned in Section 2. Hence, the free surface slope was used as an indicator for the hawser forces. As a design restriction, a maximum value of 0.14 ‰ was

Free surface oscillations in the lock chambers are triggered by asymmetries both in the flow

a) Flow distribution according to 1D model

b) Plan view of chamber

A 2D (vertically averaged) hydrodynamic model, based on code HIDROBID II developed at INA (Menéndez, 1990), was used to simulate the surface waves. It was driven by the inflow from the ports, specified as boundary conditions through time series for each one of them,

Fig. 3.2 shows the comparison between the calculated longitudinal free surface slope (using the dimensions of the physical model) and the recorded one at the physical model, for a case

Fig. 3.1. Asymmetries which trigger free surface oscillations.

that were obtained with the 1D model described in the Section 2.

contractually established for the longitudinal water surface slope.

distribution among ports, and in the geometry of the chambers (Figure 3.1).

**3.1 Description and modeling of phenomenon**

**3. Free surface oscillations**

Fig. 2.7. Comparison of physical model scale and prototype scale numerical models: Discharge in the Far Main Culvert, for a Lock to Lock operation with 21 m of initial head difference.

Fig. 2.8. Comparison of physical model scale and prototype scale numerical models: Piezometric level at the vena contracta, for a Lock to Lock operation with 21 m of initial head difference.

#### **3. Free surface oscillations**

232 Hydrodynamics – Optimizing Methods and Tools

0 60 120 180 240 300 360 420 480 540 600 660 720

**Time (secs)** Culvert Discharge - Physical Model Scale Culvert Discharge - Prototype Scale Valve Aperture - Physical Model Scale Valve Aperture - Prototype Scale

0 60 120 180 240 300 360 420

**Time (sec)**

Fig. 2.8. Comparison of physical model scale and prototype scale numerical models: Piezometric level at the vena contracta, for a Lock to Lock operation with 21 m of initial

Prototype Scale - Piezometric Level Physical Model Scale - Piezometric Level Valves

Fig. 2.7. Comparison of physical model scale and prototype scale numerical models: Discharge in the Far Main Culvert, for a Lock to Lock operation with 21 m of initial head

difference.

head difference.

**Piezometric Level (mPLD)**

**Main Culvert Dischage (m3/s)**

0%

0%

20%

40%

*Culvert Roof Elevation -7.13 mPLD*

*Culvert Sill Elevation -13.63 mPLD*

60%

**Valve Aperture (%)**

80%

100%

120%

25%

50%

**Valve Aperture (%)**

75%

100%

Free surface oscillations in the lock chambers leads to forces in the hawsers. Based on results from the physical model constructed during the development of the conceptual design, a correlation was found between these forces and the free surface slope in the absence of the vessel, as already mentioned in Section 2. Hence, the free surface slope was used as an indicator for the hawser forces. As a design restriction, a maximum value of 0.14 ‰ was contractually established for the longitudinal water surface slope.

#### **3.1 Description and modeling of phenomenon**

Free surface oscillations in the lock chambers are triggered by asymmetries both in the flow distribution among ports, and in the geometry of the chambers (Figure 3.1).

b) Plan view of chamber

Fig. 3.1. Asymmetries which trigger free surface oscillations.

A 2D (vertically averaged) hydrodynamic model, based on code HIDROBID II developed at INA (Menéndez, 1990), was used to simulate the surface waves. It was driven by the inflow from the ports, specified as boundary conditions through time series for each one of them, that were obtained with the 1D model described in the Section 2.

Fig. 3.2 shows the comparison between the calculated longitudinal free surface slope (using the dimensions of the physical model) and the recorded one at the physical model, for a case

Interaction Between Hydraulic and Numerical Models for the Design of Hydraulic Structures 235

numerical model. This indicates that turbulence scales are exerting a significant influence, so a more elaborated theoretical approach is needed. Hence, 3D modeling of the combination Central Connection + Secondary Culvert + Ports + Lock Chamber (actually, only half of the chamber, assuming that the flow is symmetrical with respect to the longitudinal axis) was

As sub-grid scale (SGS) model for the LES approach, a sub-grid kinetic energy equation eddy viscosity model was used (Sagaut, 2001). Deardorff's method was selected to define the filter cutoff length (Sagaut, 2001). A wall model was considered to treat the boundary conditions at solid borders; Spalding law-of-the-wall – which encompasses the logarithmic law (overlap region), but it holds deeper into the inner layer – was selected for the velocity (White, 1974), while a zero normal gradient condition was taken for the remaining variables. At the inflow boundary, in addition to the ensemble-averaged velocity (which arises from the 1D model), the amplitude of the stochastic components were provided (Sagaut, 2001): 4% for the longitudinal component, and 1.3% for the transversal one, values associated to a fully developed flow, very appropriate for the present problem; additionally, a weighted average of the previous and present generated stochastic components was imposed in order to add some temporal correlation; for the turbulent kinetic energy, a zero normal gradient was taken. For the free surface at the Chamber, the rigid-lid approximation was used, where uniform pressure was imposed, together with zero normal gradient conditions for the remaining quantities. The model was implemented using OpenFOAM (Open Field Operation And Manipulation), an open source toolbox for the development of customizable numerical solvers and utilities for the solution of continuum mechanics problems (Weller et al., 1998). The model solves the integral form of the conservation equations using a finite volume, cell centered approach in the spirit of Rhie and Chow (1983). PISO (Pressure Implicit with Splitting of Operators) algorithm is

Fig. 3.4 presents a view of the model domain. The computational mesh was composed by 1.5 million elements. Special considerations were made for the mesh near the wall, as the center of the first cell has to lie within a distance range to the wall – 30 *y+* 300 – to rigurously apply the logaritmic velocity profile as boundary condition (Sagaut, 2001). Typical computing times for stabilization with a steady discharge, in a Core i7 PC running 8 parallel processes, were 3 to 8 days. When complete hydroghaphs were simulated (of approximately 550 secs), 15 to 30 days of computing time were required. By parallelizing the simulation

using more than one PC, computing times were reduced, though non-linearly.

undertaken using a Large-Eddy Simulation (LES) approach (Sagaut, 2001).

**3.2 Improved theoretical approach**

used for time marching (Ferziger & Peric, 2001).

Fig. 3.4. Model domain for 3D model.

with a relatively low initial head difference (9 m in prototype units) between the Lower Chamber and the Ocean. The agreement is considered as very good, taking into account that the numerical model does not include the resolution of turbulent scales (which introduce a smaller-amplitude, higher-frequency oscillation riding on the basic oscillation).

Fig. 3.2. Longitudinal water surface slope using 1D model input. Low initial head difference.

However, the 2D model completely fails to correctly predict the longitudinal free surface slope for higher initial head differences, as observed in Fig. 3.3 for a Lock to Lock operation with an initial head difference of 21 m. More specifically, the recorded oscillation indicates a quite more irregular response, with a much higher amplitude than the one calculated with the

Fig. 3.3. Longitudinal water surface slope using 1D model input. High initial head difference.

numerical model. This indicates that turbulence scales are exerting a significant influence, so a more elaborated theoretical approach is needed. Hence, 3D modeling of the combination Central Connection + Secondary Culvert + Ports + Lock Chamber (actually, only half of the chamber, assuming that the flow is symmetrical with respect to the longitudinal axis) was undertaken using a Large-Eddy Simulation (LES) approach (Sagaut, 2001).

#### **3.2 Improved theoretical approach**

234 Hydrodynamics – Optimizing Methods and Tools

with a relatively low initial head difference (9 m in prototype units) between the Lower Chamber and the Ocean. The agreement is considered as very good, taking into account that the numerical model does not include the resolution of turbulent scales (which introduce a

0 60 120 180 240 300 360 420 480 540 600

**Time (seconds)** Physical Model Numerical Model Valve Opening

Fig. 3.2. Longitudinal water surface slope using 1D model input. Low initial head difference. However, the 2D model completely fails to correctly predict the longitudinal free surface slope for higher initial head differences, as observed in Fig. 3.3 for a Lock to Lock operation with an initial head difference of 21 m. More specifically, the recorded oscillation indicates a quite more irregular response, with a much higher amplitude than the one calculated with the

0 60 120 180 240 300 360 420 480 540 600

**Time (seconds)** Physical Model Numerical Model Valve Opening

Fig. 3.3. Longitudinal water surface slope using 1D model input. High initial head

0.00

0.00

0.25

0.50

**Valve Opening ratio**

0.75

1.00

0.25

0.50

**Valve Opening ratio**

0.75

1.00

smaller-amplitude, higher-frequency oscillation riding on the basic oscillation).

**Longitudinal Water Surface Slope (o/oo)**

difference.

**Longitudinal Water Surface Slope (o/oo)**

As sub-grid scale (SGS) model for the LES approach, a sub-grid kinetic energy equation eddy viscosity model was used (Sagaut, 2001). Deardorff's method was selected to define the filter cutoff length (Sagaut, 2001). A wall model was considered to treat the boundary conditions at solid borders; Spalding law-of-the-wall – which encompasses the logarithmic law (overlap region), but it holds deeper into the inner layer – was selected for the velocity (White, 1974), while a zero normal gradient condition was taken for the remaining variables. At the inflow boundary, in addition to the ensemble-averaged velocity (which arises from the 1D model), the amplitude of the stochastic components were provided (Sagaut, 2001): 4% for the longitudinal component, and 1.3% for the transversal one, values associated to a fully developed flow, very appropriate for the present problem; additionally, a weighted average of the previous and present generated stochastic components was imposed in order to add some temporal correlation; for the turbulent kinetic energy, a zero normal gradient was taken. For the free surface at the Chamber, the rigid-lid approximation was used, where uniform pressure was imposed, together with zero normal gradient conditions for the remaining quantities. The model was implemented using OpenFOAM (Open Field Operation And Manipulation), an open source toolbox for the development of customizable numerical solvers and utilities for the solution of continuum mechanics problems (Weller et al., 1998). The model solves the integral form of the conservation equations using a finite volume, cell centered approach in the spirit of Rhie and Chow (1983). PISO (Pressure Implicit with Splitting of Operators) algorithm is used for time marching (Ferziger & Peric, 2001).

Fig. 3.4 presents a view of the model domain. The computational mesh was composed by 1.5 million elements. Special considerations were made for the mesh near the wall, as the center of the first cell has to lie within a distance range to the wall – 30 *y+* 300 – to rigurously apply the logaritmic velocity profile as boundary condition (Sagaut, 2001). Typical computing times for stabilization with a steady discharge, in a Core i7 PC running 8 parallel processes, were 3 to 8 days. When complete hydroghaphs were simulated (of approximately 550 secs), 15 to 30 days of computing time were required. By parallelizing the simulation using more than one PC, computing times were reduced, though non-linearly.

Fig. 3.4. Model domain for 3D model.

Interaction Between Hydraulic and Numerical Models for the Design of Hydraulic Structures 237

0 60 120 180 240 300 360 420 480 540 600

**Time (seconds)** Test 1 Test 2 Test 3 a) Experimental realizations

0 60 120 180 240 300 360 420 480 540 600

**Time (seconds)** Physical Model Simulation 1 Simulation 2 b) Numerical realizations

Fig. 3.6. Variability of longitudinal water surface slope. High initial head difference.

should be fully understood. This is performed in the following.

Before proceeding to simulate prototype conditions, it is relevant to analyze the response provided by the numerical model, in order to be confident about using this tool to make such a prediction. Specifically, the physical mechanisms involved in the present problem

Fig. 3.7a shows the time series of the discharges through the U and S branches of the Central Connection (in prototype units), according to the 3D numerical model. It is observed that, for the higher discharges, they present oscillations, which seem coherently out-of-phase. The difference between those discharges is shown in Fig. 3.7b (together with the total discharge, i.e., the one through the Main Culvert). It is effectively observed that this difference oscillates, and that during the time window of higher discharges (above about 250 m3/s) there is a dominant period of oscillation which spans from 40 to 80 seconds, approximately.

0 0.035 0.07 0.105 0.14

**Longitudinal Water Surface Slope (o/oo)**

**Longitudinal Water Surface Slope (o/oo)**

Note that the rigid-lid approximation implies that the free surface oscillations are not solved by the 3D model; this was done in order to avoid extremely high computing times. Instead, the 3D model provided the time series of the flow discharge for each port, which were used to drive the 2D model of the chamber. Alternatively (and less costly in post-processing), the time series of the discharges at the U and S branches of the Central Connection, provided by the 3D model, were used to feed the 1D model, from which the discharge distribution among ports was obtained, and used to feed the 2D model.

Fig. 3.5 shows the longitudinal water surface slope obtained with the two approaches (using the dimensions of the physical model), and their comparison with the results from the physical model, for the high initial head difference case. It is observed that both numerical simulations are now able to capture the high amplitude oscillations, indicating that large eddies must be responsible for this amplification phenomenon. Note that the numerical results with input straight from the 3D model show oscillations, associated to large eddies, which are not present in the ones with input through the 1D model (which filters out those oscillations), but they are quite compatible between them.

Fig. 3.5. Longitudinal water surface slope using 3D-LES model input. High initial head difference.

The differences between the numerical results and the measurements at the physical model are due essentially to the variability of the system reponse (variations in amplitude and phase of the oscillations), under the same driving conditions, due to the stochastic nature of turbulence. This was verified both experimentally (Fig. 3.6a) and numerically (Fig. 3.6b) by repeating the same test (in the case of the numerical model, using the 'through 1D model' approach, and different initializations for the stochastic number generator). This behavior puts a limit to the degree of agreement that can be attained between the results from the numerical and physical models. In any case, the maximum amplitudes for any of the experimental or numerical realizations are relatively consistent among them.

236 Hydrodynamics – Optimizing Methods and Tools

Note that the rigid-lid approximation implies that the free surface oscillations are not solved by the 3D model; this was done in order to avoid extremely high computing times. Instead, the 3D model provided the time series of the flow discharge for each port, which were used to drive the 2D model of the chamber. Alternatively (and less costly in post-processing), the time series of the discharges at the U and S branches of the Central Connection, provided by the 3D model, were used to feed the 1D model, from which the discharge distribution

Fig. 3.5 shows the longitudinal water surface slope obtained with the two approaches (using the dimensions of the physical model), and their comparison with the results from the physical model, for the high initial head difference case. It is observed that both numerical simulations are now able to capture the high amplitude oscillations, indicating that large eddies must be responsible for this amplification phenomenon. Note that the numerical results with input straight from the 3D model show oscillations, associated to large eddies, which are not present in the ones with input through the 1D model (which filters out those

0 60 120 180 240 300 360 420 480 540 600

**Time (seconds)**

Physical Model From 3D Numerical Model Through 1D Numerical Model

The differences between the numerical results and the measurements at the physical model are due essentially to the variability of the system reponse (variations in amplitude and phase of the oscillations), under the same driving conditions, due to the stochastic nature of turbulence. This was verified both experimentally (Fig. 3.6a) and numerically (Fig. 3.6b) by repeating the same test (in the case of the numerical model, using the 'through 1D model' approach, and different initializations for the stochastic number generator). This behavior puts a limit to the degree of agreement that can be attained between the results from the numerical and physical models. In any case, the maximum amplitudes for any of the

Fig. 3.5. Longitudinal water surface slope using 3D-LES model input. High initial head

experimental or numerical realizations are relatively consistent among them.

among ports was obtained, and used to feed the 2D model.

oscillations), but they are quite compatible between them.

difference.

**Longitudinal Water Surface Slope (o/oo)**

0.035

0.07

0.105

0.14

0

Fig. 3.6. Variability of longitudinal water surface slope. High initial head difference.

Before proceeding to simulate prototype conditions, it is relevant to analyze the response provided by the numerical model, in order to be confident about using this tool to make such a prediction. Specifically, the physical mechanisms involved in the present problem should be fully understood. This is performed in the following.

Fig. 3.7a shows the time series of the discharges through the U and S branches of the Central Connection (in prototype units), according to the 3D numerical model. It is observed that, for the higher discharges, they present oscillations, which seem coherently out-of-phase. The difference between those discharges is shown in Fig. 3.7b (together with the total discharge, i.e., the one through the Main Culvert). It is effectively observed that this difference oscillates, and that during the time window of higher discharges (above about 250 m3/s) there is a dominant period of oscillation which spans from 40 to 80 seconds, approximately.

Interaction Between Hydraulic and Numerical Models for the Design of Hydraulic Structures 239

Before pursuing with the analysis, it is worth to confirm that the close-to-resonance conditions are responsible for the amplification of the free surface oscillations. Hence synthetic hydrographs for the U and S branches of the Central Connection were built, introducing a purely sinusoidal oscillation to their difference during the higher-discharges time window, as indicated in Fig. 3.8a for 2 m3/s amplitude of oscillation, and two different

0 60 120 180 240 300 360 420 480 540 600

**Time (seconds)** 60 secs 120 secs a) Discharge difference (driving force)

0 60 120 180 240 300 360 420 480 540 600

**Time (seconds)** 60 secs 120 secs

b) Longitudinal water surface slope (system response)

Fig. 3.8. Synthetic discharge difference and system response for different periods of

0.000

**Longitudinal Water Surface Slope (o/oo)**

oscillation.

0.035

0.070

0.105

0.140

**Discharge difference (m3/s)**

Now, these periods are close to, and include, the period of free surface oscillations in the Chamber (around 70 seconds), indicating that conditions close to resonance are achieved, thus resulting in an amplification of the free surface oscillation, which is the observed effect on the water surface slope. As in the numerical simulation the free surface was represented like a rigid lid, the oscillation in the discharge difference between the two branches of the Central Connection is not influenced at all by free surface oscillations themselves, i.e., the dominant period arises from the flow properties in the Central Connection. This dominant period must then be associated to the largest, energy-containing eddies (the ones resolved with the LES approach).