Mathematical Problems in Engineering
Volume 2017 (2017), Article ID 6482527, 14 pages
https://doi.org/10.1155/2017/6482527

Research Article

Analysis of the Hydroelastic Performance of Very Large Floating Structures Based on Multimodules Beam Theory

Jin Xu,1,2 Yonggang Sun,² Zhifu Li,³ Xiantao Zhang,⁴ and Da Lu⁵

¹PLA University of Science and Technology, Nanjing, China
²Wuxi First Scientific Research Institute, Wuxi, China
³College of Shipbuilding Engineering, Harbin Engineering University, Heilongjiang, China
⁴Center for Offshore Foundation Systems, School of Civil, Environmental and Mining Engineering, University of Western Australia, Perth, WA, Australia
⁵State Key Laboratory of Ocean Engineering, Shanghai Jiao Tong University, Shanghai, China

Correspondence should be addressed to Da Lu; 1045733188@qq.com

Received 30 December 2016; Revised 30 March 2017; Accepted 2 May 2017; Published 24 May 2017

Academic Editor: Yuri Vladimirovich Mikhlin

Copyright © 2017 Jin Xu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The hydroelastic behavior of very large floating structures (VLFSs) is investigated based on the proposed multimodules beam theory (MBT). To carry out the analysis, the VLFS is first divided into multiple submodules that are connected through their gravity center by a spatial beam with specific stiffness. The external force exerted on the submodules includes the wave hydrodynamic force as well as the beam bending force due to the relative displacements of different submodules. The wave hydrodynamic force is computed based on three-dimensional potential theory. The beam bending force is expressed in the form of a stiffness matrix. The motion response defined at the gravity center of the submodules is solved by the multibody hydrodynamic control equations; then both the displacement and the structure bending moment of the VLFS are determined from the stiffness matrix equations. To account for the moving point mass effects, the proposed method is extended to the time domain based on impulse response function (IRF) theory. The method is verified by comparison with existing results. Detailed results through the displacement and bending moment of the VLFS are provided to show the influence of the number of the submodules and the influence of the moving point mass.

1. Introduction

A very large floating structure (VFLS) is a unique type of oceanic structure that embraces a range of unprecedented parameters. VLFSs are designed primarily for floating airports and as a dock for calm waters on open seas.

Because it is larger than existing floating structures, a VLFS’s flexibility must be taken into account. Also, for a VLFS to be used as an airport or bridge, time-varying moving loads must be considered, for example, the analysis of the displacement and bending moment of a very large floating bridge in calm water with some moving point masses [1, 2] or the displacement response of a VLFS in waves with a moving load point [3, 4].

Simulations of the displacement and bending moment of a VLFS are mostly based on hydroelastic theory to date. In the frequency domain, hydroelastic theory was initially developed in two dimensions [5, 6]. Wu [7] extended the work of Bishop and Price to three dimensions, and Chen et al. (2003a) developed three-dimensional nonlinear hydroelastic theory. Hydroelastic theory can be also extended to the time domain [8, 9].

To account for the effects of moving point masses, simulations must be carried out in the time domain. Chakrabarti [10] and Jacobsen and Clauss [11] analyzed the hydrodynamic interactions of a multibody system based on the impulse response function (IRF) method, but the bodies are not connected to each other. Shen et al. [12] presented procedures for numerical solutions of system motion responses of a multi-rigid-body in the time domain. A moving load on a rigid body has also been studied [13–15].

An approximation theory for the analysis of the displacement and bending moment of a VLFS, based on multimodules beam theory (MBT), was recently developed by Lu et al. [16]. In this method, the VLFS is divided into submodules that are connected by a beam with specific stiffness through the gravity center of the submodules. The external force exerted on the submodules includes the wave hydrodynamic force as well as the beam bending force due to the relative displacements of different submodules. Then, the final hydroelastic response equations can be established based on multibody hydrodynamic theory and beam theory. The displacement along the VLFS can be obtained by solving these equations, and then the structural deformation as well as the section loads can be computed.

In this paper, Lu’s methodology for the computation of the displacement along the VLFS has been improved, and some new methods have been developed for the computation of the section loads. We discard the beam bending method of Lu et al. [16] and calculate the displacement response of the structure by the method of matrix transformation. To obtain the complete bending moment distribution, the high-order difference is obtained for a few accurate moment values. This method is more rigorous, and the displacement and bending moment along the whole VLFS can be obtained. This paper avoids the resonance problem caused by small gap reasonably when calculating the hydrodynamic coefficient [17].

The methodology has been extended to time domain analysis based on IRF theory, and this paper attempts to simulate the hydroelastic response of the VLFS in waves with some moving point mass. Study on the moving load of an elastic floating body under wave action is the foundation of research on landing pontoons, but there is little literature. This paper combines MBT (Lu, 2015), IRF theory [18, 19], and the analysis method of moving load on rigid floating body. MBT and IRF are used in the time domain analysis of a VLFS; the floating body in the moving load is treated as a rigid body, and the inertia force introduced by the moving load is added to the wave force of this floating body. Although the method is very approximate, the results can reflect the response characteristics of the structure.

2. Mathematical Formulation

Lu’s method is based on potential theory and multibody theory. For the hydrodynamic aspect, the ideal fluid assumption is adopted; that is, the fluid is inviscid, irrotational, and incompressible. The incident wave amplitude is assumed to be small relative to a characteristic wavelength and body dimension. For the structural aspect, a stiffness matrix is introduced between two adjacent modules to consider the flexible structure.

2.1. Hydrodynamic Model

Based on the assumptions of an ideal fluid and linearity, the velocity potential can be decomposed into three parts as follows:

where , , and denote, respectively, the incident wave potential, diffraction wave potential, and radiation wave potential. The incident, diffraction potential, and radiation potential satisfy the following boundary conditions:

As shown in Figure 1, is the fluid domain, and , , , and are the free surface, bottom surface, wetted body surface of the kth module, and boundary surface at infinity of the fluid, respectively. represents the outward-directed unit vector normal to the wetted surface of the kth module. is the speed on the wetted surface of the kth module. is the velocity potential. After the velocity potential is obtained, the added mass and the radiation damping of the bodies, as well as the wave excitation force, can be calculated.

Figure 1: Definition of the fluid and structure boundaries.

In the frequency domain, the excitation forces are related to the incident and diffracted potential as follows:

where , , , , and are the imaginary unit, wave excitation forces, fluid density, wave frequency, and average wetted surfaces, respectively.

The added mass and damping coefficient is given by

2.2. Motion Equation

According to Lu’s method, a VLFS is divided into many modules, and the equation of six-degree-of-freedom (DOF) motion in the frequency domain for a given wave amplitude and wave frequency can be written as follows:

where is the number of floating bodies, is the mass (or inertia moment) matrix of body k, is the added mass (or added mass moment) matrix of body caused by the motion of the body itself, is the radiation damping matrix of body caused by the motion of the body itself, is the hydrostatic stiffness matrix of body k, is the added mass matrix of body k induced by the motion of body j, is the radiation damping matrix of body caused by the motion of body , is the wave excitation force (or moment) of body , is the displacement matrix of bodies, is the stiffness matrix that shows the relationship between displacement and force, and is the six-degree-of-freedom (DOF) displacement of body expressed as . The number of equations is the same as the number of floating bodies.

Further, considering all modules, (5) is extended to the form of an equation set:

As a linear system, the coefficient matrix of unknowns can be replaced by , and (6) can be rewritten as

where is the total stiffness matrix. The expressions for the stiffness matrix will be discussed in the following section.

2.3. Analysis of Spatial Beam

The six-DOF motion of a module’s center is restricted by the deformation condition of the equivalent beam between two modules.

As Lu [20] has written, each side of the beam element i-j has six components of displacements and forces:

where are forces in the ith cross-section and are the corresponding bending moments. Similar expressions can be obtained for the jth cross-section. Generally, the forces and moments are assumed to be imposed on the center of the cross-section of the beam element.

The relationship between and is

where is the stiffness matrix that is used in (5), which is given in Lu et al. (2016) [16]. The specific instructions for the matrix are described by McGuire et al. (2000). The size of the matrix is . The matrix is affected by the size of each body, Young’s modulus, and the Poisson ratio. When there are floating bodies, the size of the total matrix in (7) will be .

3. Hydroelastic Response of the VLFS

3.1. Numerical Model

To illustrate this methodology, we use the model mentioned by Yokosuka [21]. This VLFS is a scaled model of the Mega-Float, constructed and developed for use in sheltered waters. The wave amplitude is 1 m, and the main parameters are listed in Table 1. The grid in the hydrodynamic calculation of the VLFS is .

Table 1: Details of the VLFS model.

Parameters
Units
Value

Length
L/m
300

Width
B/m
60

Depth
D/m
2

Draft
d/m
0.5

Young’s modulus
E/N/m²
1.19e10

Poisson ratio

0.13

Density
/kg/m³
256.25

Water depth
H/m
58.5

3.2. Displacement along the VLFS

Lu (2015) [20] proposed a method for the computation of the displacement along the VLFS based on spatial beam theory. In this paper, a more generalized method will be developed.

It can be seen from (9) that once the displacement has been solved, the section load can be computed. The displacement along a given spatial beam can be calculated based on the end interfaces’ displacements and section loads.

Equation (9) can be rewritten as

The end interface loads and can be obtained once the displacements and have been solved.

To get the displacement and section load defined at any position between the two end interfaces, (10) is modified as follows:

Then the displacement and section load defined at section can be computed by

Lu’s method (Lu, 2015) can only calculate the vertical displacement response. Equation (12) can provide the six-degree-of-freedom displacement along the VLFS, which is a big improvement over the original method proposed by Lu (2015) [20].

The present results are compared with experimental results [21] and numerically calculated data from three-dimensional hydroelastic theory [22], which are shown in Figure 2. In this part, the continuous VLFS is divided into five modules with 0-degree wave direction. Values calculated in this paper are denoted by “Present.” Experimental values [21] are denoted by “Experiment,” and numerically calculated data from three-dimensional hydroelasticity theory [22] are denoted by “Fu et al. 2007.” These results are in good agreement. The hydrodynamic responses of the model in 30- and 60-degree wave directions are calculated. The vertical displacement response of the model is compared with the results of three-dimensional hydroelasticity theory [2]. The results are listed in Figure 15.

Figure 2: Vertical response contrast of VLFS.

3.3. Bending Moment Distribution along the VLFS

To obtain the moment of the VLFS, Lu [20] proposed an equal force method based on the Darren Bell principle. It is equivalent to sets of forces being exerted on the center of each module for the VLFS. The whole problem can be described as follows: a free beam is subjected to a number of concentrated forces to reach static equilibrium and to obtain the moment distribution of the beam section. This is shown in Figure 3. The forces conform to the following formulae:

Figure 3: Schematic diagram of the equilibrium forces exerted on a free beam.

Lu adopted this method and computed the bending moment along the VLFS, where eight submodules were used, and the numerical results were compared to the values calculated by the hydroelastic theory [22]. It can be seen from Figure 4 that the results based on beam theory are in good agreement with the results calculated by Fu et al. Since the continuous distribution forces are simplified to be the concentrated force on the gravity center of the submodule, only the connection interfaces of the submodules are exactly right, which is explained in the next section; that is, eight submodules can obtain nine exact points; then the bending moment distribution along the VLFS can be computed by B-spline interpolation.

Figure 4: Moment distribution of the model.

To demonstrate the influence of the number of submodules, Figure 5 shows the comparison of bending moment distributions of five and eight submodules, respectively, in which the frequency is chosen as 0.6 rad/s. Here, FBP denotes the six exact values computed by the five submodules; FPS represents the bending moment distribution along the VLFS via B-spline interpolation; EBP means the nine exact values of the bending moment; and EPS represents the bending moment distribution along the VLFS via B-spline interpolation.

Figure 5: Moment distributions with different numbers of submodules.

It can be seen that there is a big difference between the bending moment distributions for five submodules and eight submodules. To overcome this problem, some improvements have been proposed in the following section, specifically to use fewer submodules to obtain high accuracy of the bending moment distribution.

3.4. High-Order Interpolation Method

In the multisubmodule beam theory of Lu (2015) [20], the continuous distribution pressure is treated as the concentrated force exerted on the gravity centers of the submodules. Therefore, the displacement as well as the bending moment distributions can be determined by four parameters, as seen in the following equations:

However, the continuous distribution pressure has to be taken into account when solving for the displacement and the bending moment, as seen in the following equations:

where is the vertical displacement at . , , , , , and are the displacement, rotation, equivalent bending moment and shear force, exact bending moment, and shear force, respectively, for the center of gravity of the i-th cross-section; is the continuous distribution pressure; and is the distance of the target point from the center of the i-th cross-section. Equations (14) through (16) explain why Lu’s method can only obtain several exact values.

In Lu’s original work, only the exact bending moment is used to perform the spline interpolation. However, it can be seen from (15) that, on the exact point, both the bending moment value and its partial derivative with respect to the variable x are known. Therefore, the higher-order interpolation scheme can be applied. Li (1986) showed that the interpolation curves can be cubic order if both the point value and its partial derivative at the two ending points are known; that is, and are known. Then, the interpolation curve can be expressed as

where

Since the current method takes into account both the bending moment and the shear force when solving the bending moment distribution along the VLFS, the interpolation accuracy must be greater than that of the original method of Lu. Figure 6 gives the values of the bending moment computed by five submodules and eight submodules, respectively. Here, FBHOI means the values are computed by the five submodules through the high-order interpolation method, and EBHOI means the values are computed by the eight submodules through the high-order interpolation method. It can be seen that the trend of the values computed by FBHOI is the same as that computed by EBHOI, which demonstrates that the HOI method is more accurate than the spline interpolation scheme. Besides, the comparison of EBS and EBHOI shows that if the number of submodules is enough, both the HOI method and the spline interpolation scheme can give reasonable results for the bending moment distribution along the VLFS.

Figure 6: Moment distributions computed by the high-order interpolation method.

Though the trend of the values computed by the HOI method is the same as for the spline interpolation scheme, there is still some difference in the absolute values of the bending moment distribution. Therefore, the displacement-to-displacement to moment (DDM) method is developed to improve the accuracy of the absolute values. Figure 7 gives the displacement distributions along the VLFS computed by five submodules and eight submodules, respectively; here, the wave frequencies are chosen as 0.3 rad/s and 0.6 rad/s. It can be seen that the displacement distributions along the VLFS computed by five submodules are nearly the same as those computed by eight submodules. The core idea of DDM is that the displacement distributions calculated by five submodules and eight submodules are considered the same, so the six DOF of eight modules can be converted by the six DOF of five modules, and then the bending moment can be calculated by eight submodules. In the HOI scheme, the bending moment and the shear force are taken into account, but the DDM scheme also considers the displacement information, which means this method is more accurate than the HOI scheme in application.

Figure 7: Displacement distributions with different submodules.

Figure 8 gives the bending moment distributions along the VLFS computed by the different interpolation schemes; the frequencies are chosen as 0.3 rad/s, 0.4 rad/s, 0.6 rad/s, and 0.7 rad/s. Here, FTEP means the eight bending moment points computed by the DDM scheme using five submodules, and FTEHOI represents the bending moment distributions along the VLFS computed by the DDM scheme using five submodules. It can be seen that the DDM scheme is better than the HOI method, and in realistic computations, the combination of five submodules and the DDM scheme can produce a reasonable result.

Figure 8: Bending moment distribution with different interpolation schemes.

In the actual situation, considering the computational efficiency, fewer modules are used. The DDM method is useful to obtain high-precision calculation results.

4. Hydroelastic Response with Moving Mass Points

4.1. Time Domain Hydroelastic Theory

To take into account the effects of moving mass points, the simulation must be carried out in the time domain. In terms of the convolution of body motions according to the Cummins equations [18] and Ogilive relations [19], for a multibody system (with bodies), the equations of motion can be rewritten in the following form:

where and are the mass and the added mass with infinite coefficient, respectively; is the delay function matrix; is the hydrostatic resilience matrix; and the over dot means the displacement and the partial derivative with respect to time; and is the wave exciting force.

The delay function and the added mass with infinite coefficient can be computed via the following equations:

where a and b are the frequency domain added mass and damping coefficient and is an arbitrarily chosen frequency.

Based on Lu’s theory, for the time domain hydroelastic simulations, (19) should include the stiffness matrix as follows:

where K is the stiffness matrix coefficient.

If there exists some other force, for example, moving point mass, then the above equations should be modified as follows:

where f is the additional external force.

Section 3 has shown that the combination of five submodules and the DDM scheme can provide reasonable results for both the displacement and bending moment distributions along the VLFS. Here, indirect time domain simulations are carried out, and the time domain results are transformed into the frequency domain. The comparisons of time domain (TD) values with the frequency domain (FD) results are shown in Figure 9. The results are in good agreement, but there are some numerical errors. This paper uses the Runge-Kutta method in time domain iteration. Normally, numerical methods have computational errors. In the process of using the IRF method, we need to truncate the delay function. This is another numerical error. For a number of floating bodies, the truncation term will increase, so the error will increase. The vertical response along the VLFS is calculated by (12), and the term causes calculation error. When multiplied by a coefficient matrix, the whole displacement response causes calculation error.

Figure 9: Comparisons of the time domain values with frequency domain results.

4.2. Consideration of the Moving Point Mass

To describe the position of the moving point mass, two coordinate systems are introduced in Figure 10. These are the space fixed coordinate system OXYZ and the body fixed coordinate system oxyz. With the origin located at in coordinate system OXYZ, the position of the moving point mass is in the body fixed coordinate system. Therefore, the position of the moving point mass in the space fixed coordinate system can be written as

Figure 10: Demonstration of the absolute acceleration of the moving point mass.

If the relative velocity of the moving point mass with respect to the VLFS is taken as , and keeping connected with the VLFS, then the acceleration of the moving point mass in the inertial coordinate system can be computed by the following equation (Shen, 2003; Shi min, 2001):

where is the rotational velocity of oxyz, is the absolute acceleration, is the transport acceleration, is the relative acceleration, and is the Coriolis acceleration.

For the linear case, (24) can be simplified as follows (Guo-ping, 1989):

where

Therefore, the external force exerted on the VLFS coming from the moving point mass can be computed by

where

4.3. Hydroelastic Response of the VLFS with a Moving Point Mass

Here, we use a nonrigorous method to study the hydroelastic response of the VLFS with a moving mass point.

Through the combination of time domain hydroelastic theory (22) and the external force computation equation of the moving point mass (29), the hydroelastic response of the VLFS with a moving mass point is investigated. If the moving point mass is located at the i-th submodule, then the motion equation of the VLFS can be computed as follows:

where is computed by (27).

The main purpose of this paper is to discuss the applications of the proposed method for the computation of the hydroelastic response of the VLFS, with or without the moving point mass. Therefore, only one wave condition (wave length ) is taken as an example. In the numerical simulations, the total VLFS is divided into five submodules, as shown in Figure 11; the moving point mass is taken as 10⁶ kg; the velocity of the moving point mass is 5 m/s; and the moving point mass begins to go on the VLFS at time 60 s.

Figure 11: Division of the total VLFS.

Figure 12 gives the time history of the displacement and bending moment of the VLFS, while Figures 13 and 14 give the contour of the displacement and bending moment of the VLFS with a moving point mass, respectively. It can be seen that before the moving point mass moves on the VLFS, the time variation of the VLFS is stable, and then the figure shows the influence of the moving point mass on the bending moment and the displacement.

Figure 12: Bending moment and vertical response amplitude in midpoint.

Figure 13: Contour of the displacement of the VLFS.

Figure 14: Contour of the bending moment of the VLFS.

Figure 15: The hydrodynamic response of VLFS in different wave direction.

5. Conclusion

This paper first describes a new method, based on both multibody hydrodynamics and the Euler-Bernoulli assumption, to calculate the response of the VLFS. Then two methods, the HOI method and DDM method, are proposed to improve the accuracy of the computation of the bending moment distributions along the VLFS. Via Fourier transformations, the method has been extended to time domain analysis and then used to investigate the hydroelastic response of the VLFS with a moving point mass on it.

Some test cases have been carried out. All the results are compared with the values obtained by three-dimensional hydroelastic theory and the model experiment [21]. The results show that the new method is accurate enough to deal with the hydroelastic response of the VLFS. Some initial attempts have been made to study the hydroelastic response of the VLFS with a moving load. Further investigation will be conducted to study the influence of the moving point mass on the displacement, as well as the bending moment distributions of the VLFS.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

References

J.-S. Wu and P.-Y. Shih, “Moving-load-induced vibrations of a moored floating bridge,” Computers and Structures, vol. 66, no. 4, pp. 435–461, 1998.

S. Fu, The Nonlinear Hydroelastic Analysis of Flexible Floating Body and Mooring Carrying Bridge (in Chinese) [Phd. thesis], Shanghai Jiaotong Univercity, Shanghai, China, 2005.

H. Endo, “The behavior of a VLFS and an airplane during takeoff/landing run in wave condition,” Marine Structures, vol. 13, no. 4-5, pp. 477–491, 2000.

H. Shin, H. Y. Lee, C. G. Lim et al., “Hydroelastic response simulation of a VLFS in irregular waves considering taking-off and landing impact,” in Proceedings of the 20th International Conference on Offshore Mechanics and Arctic Engineering, pp. 1–5, 2001.

C. V. Betts, R. E. D. Bishop, and W. G. Price, “The symmetric generalized fluid forces applied to a ship in a seaway,” RINA Supplementary Papers, vol. 119, 1977.

R. E. D. Bishop and W. G. Price, “Generalized anti-symmetry fluid forces applied to a ship in a seaway,” International Shipbuilding Progress, vol. 24, pp. 3–14, 1977.

Y. S. Wu, Hydroelasticity of Floating Bodies [Phd thesis], Brunel University, UK, 1984.

D. Wang, The analysis method of 3-D ship hydroelastic in time domain, China Ship Scientific Research Center, Wuxi,. China, 1996, 35–37.

X. Liu and S. Sakai, “Time domain analysis on the dynamic response of a flexible floating structure to waves,” Journal of Engineering Mechanics, vol. 128, no. 1, pp. 48–56, 2002.

S. Chakrabarti, “Time domain simulation of multiple-floating structures with nonlinear constraints,” in Proceedings of the Offshore Mechanics and Arctic Engineering (OMAE), ASME, 2001.

K. Jacobsen and G. F. Clauss, “Time-domain simulations of multi-body systems in deterministic wave trains,” in Proceedings of 25th International Conference on Offshore Mechanics and Arctic Engineering, OMAE 2006, pp. 375–384, June 2006.

Q. Shen, X. Chen, Y. Li et al., “An analysis method of the self-oscillation in time domain of a multi-body system with a floating base,” Ocean Engineering/Haiyang Gongcheng, vol. 20, no. 2, pp. 7–13, 2002.

J.-S. Wu and J.-J. Sheu, “An exact solution for a simplified model of the heave and pitch motions of a ship hull due to a moving load and a comparison with some experimental results,” Journal of Sound and Vibration, vol. 192, no. 2, pp. 495–520, 1996.

Q. Shen, Y. Li, and X. Chen, “Dynamic analysis of multibodies system with a floating base for rolling of ro-ro ship caused by wave and slip of heavy load,” Journal of Marine Science and Application, vol. 2, no. 2, pp. 17–24, 2003.

Q. Shen, Y. Li, and G. Wu, “Dynamic analysis of multi-body system with floating base in heavy equipment urgent unload process,” Journal of PLA University of Science and Technology (Natural Science), vol. 1, no. 9, 2003.

D. Lu, S. Fu, X. Zhang, F. Guo, and Y. Gao, “A method to estimate the hydroelastic behaviour of VLFS based on multi-rigid-body dynamics and beam bending,” Ships and Offshore Structures, pp. 1–9, 2016.

G. Miao, H. Ishida, and T. Saitoh, “Influence of gaps between multiple floating bodies on wave forces,” China Ocean Engineering, vol. 14, no. 4, pp. 407–422, 2000.

W. E. Cummins, “The impulse response function and ship motions,” Schiffstechnik, pp. 101–109, 1962.

T. Ogilvie, “Recent progress towards the understanding and prediction of ship motions,” in Proceedings of the Sixth Symposium on Naval Hydrodynamics, 1964.

D. Lu, Hydrodynamic Performance Study of Floating Multi-Body with Connector, Shanghai Jiao Tong University, Shanghai, China, 2016.

K. Yago and H. Endo, “Model experiment and numerical calculation of the hydroelastic behavior of matlike VLFS (VLFS '96),” in Proceedings of (VLFS '96), pp. 209–216, 1996.

S. Fu, T. Moan, X. Chen, and W. Cui, “Hydroelastic analysis of flexible floating interconnected structures,” Ocean Engineering, vol. 34, no. 11-12, pp. 1516–1531, 2007.