- Email: [email protected]
Stochastic hydroelastic analysis of a very large floating structure using pseudo-excitation method
Applied Ocean Research 48 (2014) 202–213
Contents lists available at ScienceDirect
Applied Ocean Research journal homepage: www.elsevier.com/locate/apor
Stochastic hydroelastic analysis of a very large floating structure using pseudo-excitation method Hai-tao Li, Qi-bin Wang, Zhi Zong ∗ , Lei Sun, Hui Liang State Key Laboratory of Structural Analysis for Industrial Equipment, School of Naval Architecture, Dalian University of Technology, Dalian 116024, PR China
a r t i c l e
i n f o
Article history: Received 11 April 2014 Received in revised form 1 August 2014 Accepted 1 August 2014 Keywords: VLFS Hydroelasticity Pseudo-excitation method Stationary stochastic response
a b s t r a c t This paper is concerned with the linear hydroelastic response of a pontoon-type very large floating structure (VLFS) in short-crested irregular waves. The linear potential theory is employed for the analysis of VLFS in frequency domain. To decouple the fluid–structure interaction, the higher-order boundary element method (HOBEM) combined with the finite element method (FEM) is adopted. VLFS is modeled as a Mindlin plate, and the mode superposition method is used to reduce the dimension of dynamic equations. The pseudo-excitation method (PEM) is adopted to analyze the stationary stochastic response of the floating structure. The efficiency of this new calculation scheme with the application of PEM is investigated in comparison with the conventional method for stochastic response by analyzing the computational complexity theoretically. Finally, the new calculation scheme is validated by comparing with the experimental data as well as the existing numerical results calculated in the conventional way. In addition, the efficiency of the present numerical approach is also testified which indicates that the proposed numerical scheme is time-saving. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction The very large floating structure (VLFS) is able to expand the living space for human beings. In contrast to the land reclamation from the sea, VLFS has less impact on the surrounding ecosystem and costs less. Because of the advantages above, VLFS has been applied to floating bridges, floating piers, floating fuel storage facilities, and the floating airport model Mega-Float, etc. [1]. Generally, the height of VLFS is only a few meters, which is much smaller than its horizontal dimensions. Therefore, VLFS is usually modeled as a floating plate for simplicity, and it is necessary to consider the hydroelasticity of the structure, which allows for the interaction between the floating elastic body and water waves. Most hydroelastic analyses of VLFS are carried out in regular waves in frequency domain. For linear hydroelastic analyses, the hydrodynamic loads and the dynamic response of the floating structure can be dealt with separately, i.e., firstly the hydrodynamic pressure on the floating structure is generally solved by the boundary element method (BEM), and then the response of the floating structure under the pressure is calculated. Several methods have been used to deal with the deformation of the structure, including the Galerkin method [2,3] combined with the “generalized mode” proposed by Wu [4] and Newman [5], and the mode
∗ Corresponding author. Tel.: +86 411 84707694. E-mail address: rcmh [email protected] (Z. Zong). 0141-1187/$ – see front matter © 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.apor.2014.08.003
superposition method consisting of the wet mode approach [6] and the dry mode approach [7,8], of which the dry mode is obtained by the finite element method (FEM). Different from the methods above, Meylen [9] proposed a variational equation for a floating thin plate subjected to wave forces, which combines the variational equation of the plate with the free-surface Green function method together. BEM can generate an unsymmetrical full matrix, which is difficult to deal with. In addition, irregular frequencies are present if the free-surface Green function is used. Therefore, some other methods are proposed for solving the fluid motion around VLFS, including FEM [10,11] and the hybrid finite/infinite element method [12]. However, it is difficult to analyze the hydroelasticity of VLFS using BEM or FEM for the huge horizontal dimensions of the floating structure. Therefore, analytical methods, and semianalytical methods using eigenfunctions have been widely used for the study of pontoon-type VLFSs of simple geometries and small drafts [13–16], as these methods are less time-consuming and more efficient in the parametric study of the hydroelastic response. In recent years, several computer codes for the hydroelastic analysis of VLFS have been developed, and a detailed comparative study of the simulation results from them was carried out by Riggs et el. [17]. In general, the floating bodies in actual engineering are of complex geometries; in this case, the conventional BEM is generally employed for the hydrodynamic computation of VLFS, but the required computations are extremely huge. To deal with the problem in the computation, some studies of efficient computation
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
have been carried out. Wu et el. [18] proposed the double composite singularity distribution method for the hydroelastic analysis of a floating structure which has two planes of symmetry and nearly 75% of the total equations can be reduced. Wang et el. [19] introduced a cut-off criterion for the Green function and its derivatives to convert full matrices formed by BEM into sparse matrices, and then applied an iterative sparse solver to the calculation of sparse matrices to reduce the CPU time. Kashiwagi [3] adopted bi-cubic B-spline functions to represent the unknown pressure, and employed “relative similarity relations” to evaluate the coefficient matrix formed by the pressure-distribution method, which can reduce the computational time drastically with a small number of unknowns required only. Recently, another two fast algorithms, namely the pre-corrected Fast Fourier Transform (pFFT) method and the fast multipole method (FMM), have been applied to the acceleration of BEM and make it possible to study the hydrodynamics of large floating structures of general geometries. The pFFT method has been applied to the analysis of the Mobile Offshore Base (MOB) [20,21] and can reduce the computational time and memory cost for BEM from the order O(N2 ) to O(N log N). Utsunomiya et el. [22,23] have developed another form of the higher-order boundary element method (HOBEM) accelerated by FMM, which can reduce the storage requirement and CPU time to O(N) and O(N log N), respectively. Some investigations have been carried out on the stochastic hydroelastic analysis of VLFS for a reliable design, which will be subjected to irregular waves during its whole service life. Hamamoto [24] analyzed the stochastic response of a circular floating island subjected to wind-induced waves under both long-term and shortterm descriptions of loadings, of which the short-term description of external loadings was represented by spectral density functions. Chen et al. [25,26] carried out the research on the second order response of VLFS induced by coupling of first-order wave potentials or the membrane force of the plate in multidirectional irregular waves. Ertekin et el. [27] investigated the hydroelastic response of a pontoon-type VLFS sheltered by a breakwater in regular and irregular waves based on the linear Green–Naghdi theory and the eigenfunction-expansion matching method. Using the same method as Ref. [12], Miyajima et el. [28] analyzed the hydroelastic response of the Mega-Float Phase-II model in short-crested irregular waves and compared their results with the measured values. Refs. [27,28] adopted the formula for linear systems of single degree of freedom (SDOF) to evaluate the response spectrum of VLFS, which means that the cross-correlation between responses of any two points of the floating structure is neglected, and only approximate results can be obtained in this way. Recently, Papaioannou et el. [29] studied the linear stochastic response of VLFS subjected to a directional wave spectrum by the stationary random vibration theory, which is an accurate method to the stochastic analysis of systems of multiple degrees of freedom (MDOF) without any cross-correlation terms neglected. However, the response spectra of VLFS are quite difficult to calculate due to the computations of O(N3 ) by the stationary random vibration theory, in which case the required number of unknowns N is huge. In recent years, Lin et el. [30,31] have developed the pseudo-excitation method (PEM), which is efficient in dealing with stochastic responses of linear systems of large degrees of freedom, and has been applied to stochastic analyses of long-span bridges, high-speed trains [32], and the fatigue damage of the deepwater riser [33], etc. In this study, we adopt PEM to study the linear hydroelasticity of a pontoon-type VLFS subjected to short-crested irregular waves in frequency domain. VLFS is modeled as a Mindlin plate with allowance for the effect of transverse shear deformation and rotary inertia. The fluid motion is solved by HOBEM, and the deformation of the floating structure is calculated by the mode superposition
203
Fig. 1. The schematic diagram of the fluid–structure interaction problem.
method, of which dry modes of the structure obtained by FEM are chosen. A directional wave spectrum is used for the short-term description of short-crested irregular waves. The efficiency of this calculation scheme with the application of PEM is discussed theoretically and verified by a numerical example at the end of this paper. Finally, the accuracy of the present method is also verified by comparing our numerical results with the experimental data as well as the existing numerical results obtained by the stationary random vibration theory. Some meaningful conclusions about the present calculation scheme are proposed. 2. Hydroelastic analysis of VLFS In this section, the linear hydroelasticity of VLFS in regular waves is analyzed in frequency domain. The schematic diagram of the fluid–structure interaction problem is shown in Fig. 1. The pontoontype VLFS is defined with length L, width B, height T, and modeled as a floating plate assumed to be flat with free edges. A zero-draft assumption is adopted for simplicity. The symbols ˝, Sf , Sb and Sd denote the fluid domain, the free surface, the bottom surface of the structure and the flat seabed, respectively. Three-dimensional Cartesian coordinate system o-xyz is established with o-xy plane coinciding with the undisturbed free surface and z-axis orienting positively upwards. 2.1. Equations of motion for the plate The height of VLFS is much smaller than its length and width, so the floating structure is usually modeled as a Kirchhoff plate which does not allow for transverse shear deformation. It is accurate to use the Kirchhoff plate theory when the plate is homogeneous and thin, in which case the effect of shear deformation is of no significance. However, if VLFS consists of different kinds of materials which behaves like sandwich plates, the effect of shear deformation cannot be neglected. In addition, when analyzing the wave-induced transverse shear force of VLFS, Gao et al. [8] found that the maximum shear force obtained by the Mindlin plate theory is about 10% larger than the one by the Kirchhoff plate theory. Therefore, in this study we model VLFS as an isotropic Mindlin plate for a more accurate evaluation of shear forces of the structure. The displacements of a Mindlin plate are given as [34] UP (x, y, z, t) = zy (x, y, t)
(1a)
VP (x, y, z, t) = −zx (x, y, t)
(1b)
WP (x, y, z, t) = W (x, y, t)
(1c)
where UP , VP , and WP denote the displacements along x-axis, y-axis, and z-axis, respectively. x , y represent the rotations about x-axis and y-axis, respectively, and W is the vertical displacement of the neutral surface of the plate.
204
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
A two-dimensional 8-node isoparametric element is used to build the FEM model. Then the stiffness matrix [K]Nf ×Nf and the consistent mass matrix [M]Nf ×Nf are obtained by standard FEM procedures, where Nf denotes the number of displacements of the plate. The free boundary condition should be satisfied, namely the bending moment, the transverse shear force and the twisting moment must vanish at the free edges of the plate. The steady state of the structural motion is considered, so the displacements of the plate can be expressed as x (x, y, t) = R[x (x, y)e−iωt ] −iωt
y (x, y, t) = R[y (x, y)e
(2a) (2b)
]
W (x, y, t) = R[w(x, y)e−iωt ]
(2c)
where √ x , y and w are the complex displacements of the plate. i = −1 is the imaginary unit, and ω is the circular frequency. R denotes the real part of a complex expression. Regardless of the structural damping, dynamic equations for the structure are given as
(−ω2 [M] + [K]) d where
d
Nf
=
=
F
d1 , d2 , . . .dN
(3)
T
and
F
Nf
are global vectors of
the complex displacement and the complex load, respectively, of which the displacement vector of the ith node is given as di =
T
xi , yi , wi . The subscript N denotes the number of nodes of the plate and satisfies Nf = 3N in the rest of this paper. Assuming that the free vibration of the plate has a timeharmonic solution with the frequency ωi , the natural modes of the plate can be obtained by the following eigenvalue simultaneous equations:
(−ωi2 [M] + [K])ci
=0
(4)
where ωi and ci are the ith natural frequency and dry mode vector of the Mindlin plate, respectively. 0 denotes the zero vector. 2.2. Equations of motion for the fluid domain The fluid is assumed to be ideal and incompressible, and its motion is irrotational so that the velocity potential ˚ exists. The linear potential theory is employed based on the assumption of small fluid and structural motions. The velocity potential and the dynamic pressure on the plate can be written as ˚(x, y, z, t) = R[(x, y, z)e−iωt ] −iωt
P(x, y, t) = R[p(x, y)e
]
(5b)
(6)
where I (x, y, z), D (x, y, z) and R (x, y, z) denote the incident wave potential, the diffraction potential and the radiation potential, respectively. The incident wave potential has the following form according to the linear wave theory: I (x, y, z) = −
igA cosh k(z + h) ik(x cos +y sin ) e ω cosh kh
∇ 2 = 0 in ˝
(7)
where g and A represent the gravitational acceleration and the wave amplitude, respectively. h and represent the water depth and the propagation direction of the incident wave, respectively, as illustrated in Fig. 1. k denotes the wave number for finite water depth which satisfies the linear wave dispersion relationship ω2 = gk tanh kh.
(8a)
ω2
∂ on Sf = g ∂z
(8b)
∂ = −iωw(x, y) on Sb ∂z
(8c)
∂ = 0 on Sd ∂z
(8d)
lim
√ r
r→∞
∂ − ik ∂r
( − I ) = 0 on S∞
(8e)
where r denotes the horizontal distance between the origin and the referred point, and S∞ represents the artificial fluid boundary at infinity. Eq. (8b) is the linearized free-surface boundary condition, which implies that the normal velocity of the fluid is equal to that of the free surface and the pressure on the free surface is atmospheric. Eq. (8c) is the plate surface condition and implies that there is no gap between the plate and the surface of the fluid domain. Eq. (8d) is the seabed boundary condition which implies the impermeability of the seabed. Eq. (8e) is the Sommerfeld radiation condition for D and R as r→ ∞ which indicates the propagation of diffracted and radiated waves is outgoing at infinity. By substituting Eq. (6) into Eq. (8c), we obtain another form of the plate surface condition, i.e.,
∂D ∂I =− ∂z ∂z
on Sb
(9a)
∂R = −iωw(x, y) on Sb ∂z
(9b)
The free-surface Green function for finite water depth [35] which satisfies the linearized free-surface boundary condition, the seabed boundary condition and the radiation condition at infinity is used, i.e.,
G(x, x0 ) = 2i
(5a)
where (x, y, z) and p(x, y) are the complex velocity potential and the complex dynamic pressure, respectively. The complex velocity potential consists of three parts based on the linear potential theory, i.e., (x, y, z) = I (x, y, z) + D (x, y, z) + R (x, y, z)
The velocity potential must satisfy Laplace’s equation and the following boundary conditions:
k02 − 2 h(k02 − 2 ) +
∞
+ 4
n=1
kn2 + 2 h(kn2
(1)
cosh k0 (z + h) cosh k0 (z0 +h)H0 (k0 R)
+ 2 ) −
cos kn (z + h) cos kn (z0 + h)K0 (kn R) (10)
where x = (x, y, z) and x0 = (x0 , y0 , z0 ) are the field point and the source point, respectively. = ω2 /g is the wave number for deep water. k0 and kn are the real roots of ω2 = gk0 tanh k0 h and ω2 = − gkn tan kn h, respectively. R represents the horizontal dis(1) tance between the source point and the field point. H0 (x) is the Hankel function of the first kind of order zero, and K0 (x) is the modified Bessel function of the second kind of order zero. Applying the free-surface Green function and Laplace’s equation to Green’s second identity, we obtain the following boundary integral equation:
C(x)(x)+
(x0 ) Sb
∂G(x, x0 ) ds = ∂z
G(x, x0 ) Sb
∂(x0 ) ds (11) ∂z
where C(x) is the solid angle at the field point, and equal to 2, 3, 3.5 if the field point is located at the bottom surface, the edges, the corners of the plate, respectively.
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
By discretizing the plate surface into a finite number of twodimensional 8-node isoparametric elements, we obtain the velocity potential of one point, i.e., =
8
Ni i = N
(12)
e
i=1
[A] D
[A] R
= [B]
= [B]
∂D ∂z ∂R ∂z
= −[B]
∂I ∂z
R
= [A]−1 [B]
(13b)
∂R ∂z
= −iω[A]−1 [B] {w}
p(x, y) = iω(x, y, 0) − gw(x, y)
(15)
where is the density of water. And the vector of the pressure induced by the radiation potential and the vertical motion of the plate is {pW } = ω2 [A]−1 [B] {ω} − g {ω}
(16)
Then the element added mass matrix [ma ], the element wave damping matrix [ca ] and the element hydrostatic stiffness matrix [ka ] are obtained by the integration of pW over an element, i.e., [ma ] = [S]R([a])
(17a)
[ca ] = ω[S]I([a])
(17b)
where [S]8×8 =
(17c)
T
Se
N N ds is an integration over an element.
[a]8×8 is a submatrix of [A]−1 [B] corresponding to the locations of 8 potentials of an element in {R }. I denotes the imaginary part of a complex expression. The element vector of the exciting force is obtained by
FS
e
NT N ds
= iω Se
I
e
+ D
{d} =
e
q
i ci = [c]
(19)
i=1
where { }q = { 1 , 2 , . . ., q }T consists of complex amplitudes of the first q dry modes of the plate, and [c]Nf ×q = [c1 , c2 , . . ., cq ] contains the first q dry mode vectors obtained in Eq. (4). After substituting [Ma ], [Ca ], [Ka ] and {FS } into Eq. (3), another form of simultaneous equations for fluid–structure interaction in frequency domain is obtained, i.e.,
(−ω2 [M] − ω2 [Ma ] − iω[Ca ] + [K] + [Ka ]) d
=
FS
(20)
By substituting Eq. (19) into Eq. (20), and multiplying both sides of Eq. (20) by [c]T , we obtain the ultimate form of simultaneous equations for fluid–structure interaction, of which the dimension is reduced by the mode superposition method, i.e.,
[c]T (−ω2 [M] − ω2 [Ma ] − iω[Ca ] + [K] + [Ka ])[c]
= [c]T FS
(21)
(14)
where {w}N = {w1 , w2 , . . ., wN }T contains the vertical displacement at each point of the plate. Then the pressure distribution on the bottom surface of the plate is obtained from linearized unsteady Bernoulli’s equation, i.e.,
[ka ] = g[S]
The mode superposition method is used to solve the displacement field of the plate, i.e.,
(13a)
The entries of [A]N×N and [B]N×N can be evaluated by numerical integration methods. When the source point and the field point are located in different elements, the boundary integral can be evaluated by the Gauss–Legendre quadrature rule directly. But if both the source and the field points are located in the same element, the Green function and its derivatives approach infinity at the rates of 1/R and 1/R2 , respectively; for this case, the triangular polar coordinate transformation proposed by Li et al. [36] is adopted to evaluate the boundary integral. By substituting Eq. (9b) into Eq. (13b), we obtain
is obtained by assembling the subvector in Eq. (18). It should be noted that the entries of the above matrices and vector, which are corresponding to the rotational displacements in {d}, namely x and y , are equal to zero as no hydrodynamic forces are assumed to be induced by rotational displacements of the plate. 2.3. Dynamic equations for fluid-structure interaction
where N = {N1 , N2 , . . ., N8 } is the shape function vector of the 8-node isoparametric element, and {}e = {1 , 2 , . . ., 8 }T is the potential vector that contains potentials at 8 nodes of an element. By substituting Eqs. (9a) and (12) into Eq. (11), we obtain the discretized form of the boundary integral equation [6], i.e.,
205
(18)
where the subscript e denotes the component of an element. Finally the global added mass matrix [Ma ]Nf ×Nf , the global wave damping matrix [Ca ]Nf ×Nf , and the global hydrostatic stiffness matrix [Ka ]Nf ×Nf are constructed by assembling the submatrices in
Eqs. (17a)–(17c), and the global vector of the exciting force FS
Nf
3. Stochastic analysis of VLFS 3.1. Directional wave spectrum Generally, irregular wind waves are assumed to be stationary, ergodic and can be described by a zero mean Gaussian process during a short period. The characteristics of short-crested irregular waves can be represented by a directional wave spectrum S(ω, ), which represents the distribution of wave energy on wave frequencies and propagation directions. The response of a linear system satisfies the Gaussian distribution if the input can be described by a Gaussian process. Therefore, the frequency characteristics of the response induced by short-crested irregular waves can be described as a directional spectrum. In this section, the linear stochastic response of VLFS is investigated in frequency domain. The relationship between a directional wave spectrum S(ω, ) and a one-dimensional frequency spectrum S(ω) can be expressed as S(ω, ) = S(ω)D(|ω)
(22)
where D(|ω) is the directional spreading function representthe directional distribution of wave energy, and it satisfies
ing
D(|ω) d = 1. The free-surface elevation of irregular wind waves (x, y, t) can be obtained by the linear superposition of infinite independent monochromatic waves of various frequencies and propagation directions, i.e., −
(x, y, t) = R
∞ ∞
i[kn (x cos m +y sin m )−ωn t+εmn ]
amn e
(23)
m=1 n=1
where kn , m , and ωn are the wave number, the wave direction and the circular frequency of a monochromatic wave, respectively. εmn
206
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
denotes the random phase distributed uniformly in [0, 2], and amn is the amplitude of a monochromatic wave and can be obtained by
time-harmonic excitation, which is referred to as the pseudoexcitation, is given as
amn =
x(t) =
(24)
2S(ω, )ıωı
where ıω and ı denote small intervals around the wave frequency ω and the wave direction , respectively. We assume the independence of the directional distribution on wave frequencies, and propagation directions of waves are assumed to range from −/2 to /2 around the mean wave direction . The directional spreading function proposed by Pierson et el. [37] is adopted, i.e.,
D( | ω) = D() =
⎧2 ⎨ cos2 ( − ) | − | ≤
2 | − | > 2
⎩0
(25)
S(ω, ) d
(26)
−/2
In this study, the two-parameter Bretschneider–Mitsuyasu spectrum that describes the frequency characteristic of fully developed wind waves [38] is adopted, i.e., SBM (ω) =
−4 2 0.257H1/3 T1/3
ω −5 2
ω exp −1.03 T1/3 2
−4
(27)
where H1/3 and T1/3 are the significant wave height and the significant wave period, respectively.
For a linear MDOF system, the relationship between a timeharmonic excitation with the frequency ω and the corresponding steady-state response can be described as
(28)
where [H(ω)] is the transfer function matrix for a linear MDOF system. Then, the response spectrum matrix and cross-spectrum matrices of a stationary random process can be obtained in the stationary random vibration theory, i.e., [Syy (ω)] = [H(ω)]∗ [Sxx (ω)][H(ω)]T
(29a)
[Sxy (ω)] = [Sxx (ω)][H(ω)]T
(29b)
∗
[Syx (ω)] = [H(ω)] [Sxx (ω)]
y(t) = HT (ω)x(t) = HT (ω) Sxx (ω) e−iωt
(31)
where HT (ω) is the transfer function for a linear SDOF system. We can obtain the auto-spectrum of the response and crossspectra between the excitation and the response by
2
∗
(29c)
where [Sxx (ω)] and [Syy (ω)] denote spectrum matrices of the excitation and the response, respectively. [Sxy (ω)] and [Syx (ω)] are cross-spectrum matrices between the excitation and the response. The superscript * denotes complex conjugate. Eqs. (29a)–(29c) have simple forms and are comparatively convenient for engineering applications, which is the conventional way to analyze the stochastic response. However, if the number of degrees of freedom M is large, the conventional method will be very time-consuming due to the O(M3 ) multiplications of the traditional algorithm for matrix multiplication. Therefore, PEM is applied to the stochastic analysis of VLFS in this study due to its high efficiency in dealing with stochastic responses of linear systems of large degrees of freedom. The stationary random vibration of a linear SDOF system is taken as an example to explain the basic principle of PEM. A
(32a)
x˜ (t) y(t) = Sxx (ω)HT (ω) = Sxy (ω)
(32b)
y(t)∗x(t) = HT (ω)∗ Sxx (ω) = Syx (ω)
(32c)
The process involving Eqs. (30)–(32c) is known as PEM for a linear SDOF system, and the response spectrum matrix and crossspectrum matrices between the excitation and the response of a linear MDOF system can be obtained using PEM in the similar way, i.e., [Syy (ω)] = [Sxy (ω)] = [Syx (ω)] =
∗
y˜ (t)
∗
x˜ (t)
3.2. Pseudo-excitation method [30,31]
{y(t)} = [H(ω)] x(t)
(30)
where Sxx (ω) is the auto-spectrum of the stationary random excitation x(t). Then the corresponding pseudo-response is obtained by substix(t) for x(t) in the formula y(t) = HT (ω) x(t), i.e., tuting
/2
S(ω) =
Sxx (ω) e−iωt
y˜ (t)∗ y˜ (t) = HT (ω) Sxx (ω) = Syy (ω)
Therefore, the relationship between the one-dimensional frequency spectrum and the directional wave spectrum can be written as
(33a)
T
y˜ (t)
∗
y˜ (t)
T
y˜ (t)
(33b)
T
x˜ (t)
(33c)
where x˜ (t) and y˜ (t) are vectors of the pseudo-excitation and the pseudo-response, respectively. It is obvious that the results obtained from Eqs. (33a)–(33c) are identical to those from Eqs. (29a)–(29c), which demonstrates that PEM is equivalent to the stationary random vibration theory when dealing with the stationary stochastic response of a linear MDOF system. Moreover, if the auto-spectrum of the response alone is desired, the auto-spectrum of the ith response can be obtained by [Syy (ω)]i,i =
∗
y˜ (t)
i
y˜ (t)
i
(34)
where y(t) is the ith entry of the pseudo-response vector. i
Eq. (26) assumes that components of the wave spectrum in different wave directions are independent of each other. As the transfer function from wave loads to the response of the floating structure is independent of wave directions as well, the total response spectrum can be obtained by the superposition of response spectra induced by waves of different directions, which satisfies Eq. (26) as well. Therefore, firstly we study the displacement spectrum of VLFS in single wave direction with PEM, and then the total displacement spectrum is obtained by the superposition of components in all directions. 3.3. Pseudo-excitation method for VLFS The component of the response spectrum with the frequency ω and the wave direction is considered. Based on the zero-draft assumption of VLFS, the relationship between the incident wave potential and the free-surface elevation is I (x, y, 0) = HI , (ω) (x, y) where HI , (ω) = −ig/ω is the transfer function from to I .
(35)
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
Incident wave potentials at any two points, namely the ith point (xi , yi ) and the jth point (xj , yj ) of the plate are obtained from Eq. (7), i.e., igA ik(x I (xi , yi , 0) = − e i ω I (xj , yj , 0) = −
igA ik(x e j ω
cos +yi sin )
(36a)
cos +yj sin )
(36b)
So the relationship between incident wave potentials at these two points can be written as ik(xij cos +yij sin )
I (xj , yj , 0) = I (xi , yi , 0) e
(37)
where xij = xj − xi , yij = yj − yi , and eik(xij cos +yij sin ) indicates that a phase-lag exists between incident wave potentials at two different points. The phase-lag is caused by the propagation of water waves, which is referred to as the wave passage effect. As VLFS is generally several kilometers long, a large phase difference exists between excitations of the two points with a large spatial distance. Therefore, the wave passage effect must be taken into account for the stochastic analysis of VLFS. Considering that a linear MDOF system is subjected to a waveinduced stationary random excitation in the form of F0 (t) at the origin of the coordinate system, the excitations of other points can be given as
FE (t)
⎧ ⎫ F1 (t) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ F2 (t) ⎪ ⎬
=
N
⎪ ⎪ ⎪ ⎪ ⎩
.. . FN (t)
⎪ ⎪ ⎪ ⎪ ⎭
=
⎧ ⎫ F0 (t − t1 ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ F0 (t − t2 ) ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎩
.. . F0 (t − tN )
⎪ ⎪ ⎪ ⎪ ⎭
FE (t)
=
N
. ⎪ ⎪ ⎪ ⎪ . ⎪ ⎪ ⎪ ⎩ . ⎪ ⎭
(38)
N
=
where H∂I /∂z,I = k tanh kh is the transfer function from I to Then applying the pseudo-excitation vector of the incident wave potential to Eq. (41), we obtain the following pseudo-excitation vector:
∂I ∂z
=
(ω, )
I (ω, )
N
H∂I /∂z,I
(42)
N
vector of the diffraction potential The pseudo-excitation D (ω, )
is obtained by substituting Eq. (42) into Eq. (13a).
Substituting the element pseudo-excitation vectors of the incident wave potential and the diffraction potential into Eq.(18), weobtain the element vector of the pseudo exciting force F˜S (ω, ) and
FS (ω, )
sequently the global one
by assembling
each element. The pseudo-amplitude vector of dry modes
FS (ω, )
e
F˜S (ω, )
e
of
(ω, ) is
into Eq. (21). Then we back-
the pseudo-displacement vector
d(ω, ) and sequently the dis-
placement spectrum matrix in the wave direction , i.e.,
∗ T d(ω, ) d(ω, ) Nf
(43)
Nf
Then, we substitute the results of Eq. (43) into Eq. (26) and obtain the displacement spectrum matrix [Sdd (ω)], and the jth order moment of the auto-spectrum of the ith displacement in the global displacement vector mi,j can be obtained by
∞
ωj [Sdd (ω)]i,i dω
(44)
0
SFF (ω) e−iωt
.. . eik(xN
(41)
∂I /∂z.
mi,j =
⎧ ik(x1 cos +y1 sin ) ⎫ e ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ eik(x2 cos +y2 sin ) ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎩
= H∂I /∂z,I I z=0
[Sdd (ω, )]Nf ×Nf =
(39)
As the time-harmonic factor e−iωt in a pseudo-excitation and its conjugate eiωt always come in pairs and are offset after multiplication, the time-harmonic factor is extracted from the pseudo-excitation and the pseudo-response afterwards. In this study, the directional wave spectrum is assumed to be uniform at each point of the plate, so the pseudo-excitation vector of the incident wave potential in the wave direction , which allows for the wave passage effect, can be written as
I (ω, )
∂I ∂z
substitute the pseudo-amplitude vector into Eq. (19) and obtain
eiωtN
The derivative of I in the z direction is
obtained by substituting
where t1 , t2 , . . ., tN denote time-lags between excitations of the origin and other points. SFF (ω) e−iωt for Substituting the pseudo-excitation F0 (t) = F0 (t) in Eq. (38), where SFF (ω) is the auto-spectrum of F0 (t), we obtain the pseudo-excitation of each point considering the wave passage effect, i.e.,
⎧ iωt1 ⎫ e ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ eiωt2 ⎪ ⎪ ⎨ ⎬
207
cos +yN sin )
⎪ ⎪ ⎪ ⎪ ⎭
S
(ω)D() HI , (ω) (40)
where k(xi cos + yi sin ) denotes the phase-lag between the ith point and the origin, and is equal to ωti as shown in Eq. (39). It is obvious that the spectrum matrix of the incident wave potential in the direction can be obtained by substituting Eq. (40) into Eq. (33a), which has the same form as that in Ref. [29].
where [Sdd (ω)]i,i is the ith diagonal entry of [Sdd (ω)]. The variance of the displacement can be obtained by setting j = 0. 3.4. The comparison of computational efficiency between the present calculation scheme and the conventional method When the traditional BEM is used to analyze the hydrodynamics of floating bodies, the computational time for factorization solvers and iterative solvers is O(N3 ) and O(N2 ), respectively, which is huge for the analysis of VLFS. Therefore, some fast algorithms, such as the pFFT method and FMM, have been applied to the acceleration of BEM to make it possible to solve the hydroelastic response of the floating structure. However, there still exists a problem when dealing with the stochastic response of VLFS in the conventional way, as studied in Ref. [29], i.e., one has to calculate Eq. (29a) twice to obtain the load spectrum and the displacement spectrum, which needs the computations of O(2Nf3 ) using the traditional algorithm for matrix multiplication. In this case, the computational time for creating response spectrum matrices is several times more than that for solving the discrete boundary integral equations by LU factorization (with the computations of O(N3 /3)), and even more than that for the same solution using iterative methods or accelerated BEM. Some favorite math libraries, including the Basic Linear Algebra Subprograms (BLAS) [39] and the Intel Math Kernel Library (MKL) implementation of BLAS [40], provide highly efficient algorithms
208
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
Table 1 Parameters of the VLFS model. Parameter
Symbol
Unit
Value
Total length Total width Total height Young’s modulus Poisson’s ratio Plate density Water density Water depth
L B T E p p h
m m m N/m2
300 60 2 1.19 × 1010 0.13 256.25 1025 58.5
kg/m3 kg/m3 m
for matrix multiplication. The computational time for multiplication of two square matrices of order M is O(M3 ) by routine GEMM involved in both libraries, yet it is still huge for problems with large number of unknowns. So the major concern about the fast analysis of VLFS in irregular waves is to reduce the operations for creating response spectrum matrices, which can be achieved with the application of PEM. The formation of response spectrum matrices using PEM just requires the matrix multiplication of two vectors once with the computations of O(Nf2 ), as shown in Eq. (33a), and even O(Nf ) if the auto-spectrum is required only, as described in Eq. (34). During the process of the stochastic analysis, Eqs. (29a), (33a) or (34) have to be executed for many discrete frequencies to ensure the accuracy of spectral curves and the corresponding variances. Therefore, the present calculation scheme using PEM has a great advantage over the conventional method in stochastic analyses of floating bodies with large degrees of freedom. As the computational time for creating response spectrum matrices in the conventional way is greater than that for solving the boundary integral equations, which are the most time-consuming parts in the stochastic analysis of VLFS, we conclude that more than half of the total computational time can be reduced by the present calculation scheme. 4. Numerical results and discussion 4.1. Verification of the accuracy of the present calculation scheme The accuracy of the present calculation scheme is verified by comparing our numerical results with a group of experimental data as well as other numerical results. The parameters of the VLFS model used here are identical to those in Ref. [29], as listed in Table 1. 4.1.1. Verification of the present calculation scheme by experimental data Due to the lack of experiment data on the stochastic response of a VLFS model, in this part the accuracy of the present calculation scheme is verified by comparing our numerical results with the experimental data of Yago et al. [41], which provides ratios between amplitudes of the vertical displacement and the wave in regular
waves of various wavelengths and wave directions. However, it should be noted that the present method is used for the hydroelastic analysis of VLFS in short-crested irregular waves, rather than that in regular waves. Hence, we compare the response of VLFS calculated by the present method with the experimental data in the case that the floating structure is subjected to a one-dimensional frequency spectrum S(ω), as long-crested irregular waves can be generated by a linear superposition of infinite independent monochromatic waves with a uniform wave direction. The process for the stochastic analysis in single wave direction, as described in Eqs. (40)–(43), is employed for this verification. For a given wave frequency ω and wave direction , the amplitudes of a monochromatic wave A(ω) and the vertical displacement w(ω, ) of one point can be expressed as A(ω) =
2S(ω)ıω
(45a)
w(ω, ) = 2Sww (ω, )ıω
(45b)
where S(ω) is the one-dimensional frequency spectrum as shown in Eq. (27), and Sww (ω, ) is the auto-spectrum of the vertical displacement obtained from Eq. (43). As a consequence, the ratio between amplitudes of the vertical displacement and the wave can be obtained by
w(ω, ) A(ω)
!
=
Sww (ω, ) S(ω)
(46)
As the test model in Ref. [41] is similar to the present one, the non-dimensional displacements measured in the experiment should coincide with those calculated by the present method theoretically. Therefore, we compare the calculated amplitude ratios along the longitudinal centerline of the present VLFS model with the experimental data. In this part, four wave directions, namely = 0◦ , 30◦ , 60◦ and 90◦ , and five wavelength to structure length ratios, namely /L = 0.2, 0.4, 0.6, 0.8 and 1.0 (the corresponding wave frequencies are 1.01 rad/s, 0.72 rad/s, 0.58 rad/s, 0.48 rad/s and 0.42 rad/s, respectively), are chosen for a thorough comparison. 500 8-node elements with 6 m × 6 m are adopted to discretize the plate, and the first 30 dry modes of the structure are used. The significant wave height and wave period are H1/3 = 2 m and T1/3 = 6.3 s, respectively. The comparative results are exhibited in Figs. 2–5, in which solid lines and scatters represent the numerical results calculated by the present method and the experimental data, respectively. It is clear to see that the agreement between our numerical results and the experimental data is good for all cases, except for a little difference in the case of /L = 0.4 and = 60◦ ; however, it is evident that the tendency of the numerical result in this case matches well with the experimental data. Based on the assumption of the independence of spectrum components in different wave directions, the total response spectrum can be obtained by the superposition of response spectra induced
Fig. 2. The comparison of the vertical displacement amplitude along the longitudinal centerline of VLFS for = 0◦ .
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
209
Fig. 3. The comparison of the vertical displacement amplitude along the longitudinal centerline of VLFS for = 30◦ .
Fig. 4. The comparison of the vertical displacement amplitude along the longitudinal centerline of VLFS for = 60◦ .
Fig. 5. The comparison of the vertical displacement amplitude along the longitudinal centerline of VLFS for = 90◦ .
by waves of different directions. Therefore, due to the excellent agreement between our numerical results and the experimental data for various wave frequencies and directions, it is plausible that the present calculation scheme is of enough accuracy for the analysis of VLFS subjected to short-crested irregular waves. 4.1.2. Verification of the present calculation scheme by other numerical results In this part, the accuracy of the present calculation scheme is verified by comparing our numerical results with those in Ref. [29] in the case that VLFS is subjected to a directional wave spectrum. For a more intuitive comparison, the present VLFS model, as shown in Table 1, and the stochastic wave conditions are identical to those in Ref. [29]: 2000 8-node elements with 3 m × 3 m are adopted to discretize the plate; the significant wave height and wave period are H1/3 = 2 m and T1/3 = 6.3 s, respectively; four mean wave directions, namely = 0◦ , 30◦ , 60◦ and 90◦ are chosen to investigate the effect of incident mean wave directions on the hydroelastic response of VLFS; four corner points (P1–P4) and the center point (P5), as shown in Fig. 6, are selected to study the stochastic response of the floating structure. The only difference is that the first 60 dry modes are used here for a more accurate analysis of the response
Fig. 6. Positions of 5 selected points.
of the floating structure, while the first 30 dry modes are adopted in the reference. Figs. 7–10 show the spectral curves of vertical displacements at 5 selected points for four mean wave directions. It is clear to see that the results in this study are in good agreement with those in Ref. [29], and the frequency of the peak point in the spectral curve obtained by the present method is consistent well with that in the reference in each case. However, the values of vertical displacement spectra in this study are a few smaller; in the case of = 30◦ , the frequency of the second peak point in the spectral curve of P3
210
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
Fig. 7. The comparison of response spectra of vertical displacements at 5 selected points for = 0◦ .
Fig. 8. The comparison of response spectra of vertical displacements at 5 selected points for = 30◦ .
Fig. 9. The comparison of response spectra of vertical displacements at 5 selected points for = 60◦ .
Fig. 10. The comparison of response spectra of vertical displacements at 5 selected points for = 90◦ .
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
211
Fig. 11. The comparison of standard deviations of vertical displacements for = 0◦ . The left result is obtained from Ref. [29], and the right is calculated by the present method.
Fig. 12. The comparison of standard deviations of vertical displacements for = 30◦ . The left result is obtained from Ref. [29], and the right is calculated by the present method.
Fig. 13. The comparison of standard deviations of vertical displacements for = 60◦ . The left result is obtained from Ref. [29], and the right is calculated by the present method.
is about 1.15 rad/s in this study compared with 1.4 rad/s in the reference; regarding the case of = 60◦ , a large discrepancy on the frequency of the peak point in the spectral curve of P1 exists, and the peak value of P1 is greater than that of P2 in this study, which is opposite to the result of the reference. The probable reason is that different BEM models cause the aforementioned discrepancies, rather than the inaccuracy of the present method, as PEM is equivalent to the stationary random vibration theory which is employed in the reference.
In Figs. 11–14, the distributions of standard deviations of vertical displacements for four mean wave directions are presented. The left four figures denote the results given in Ref. [29], and the right four are obtained by the present method. It is clear that the distribution of standard deviations calculated by the present method agrees well with that of the reference in each case, and the maximum discrepancy of standard deviations is less than 10%, which is caused by the results that the spectra of vertical displacements in this study are smaller, as demonstrated in Figs. 7–10.
212
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
Fig. 14. The comparison of standard deviations of vertical displacements for = 90◦ . The left result is obtained from Ref. [29], and the right is calculated by the present method. Table 2 The computational time for different element sizes. Case
Number of nodes
(a) (b) (c) (d)
1621 3556 6241 8929
Present method
Conventional method
LU factorization
Matrix multiplication
LU factorization
Matrix multiplication
6s 59 s 5.3 min 16.4 min
1s 5s 16 s 41 s
6s 1 min 5.3 min 16.9 min
4.9 min 49 min 4.6 h 12.8 h
4.2. Verification of the efficiency of the present calculation scheme To testify the high efficiency of the present calculation scheme with the application of PEM, a comparison is conducted between the computational time of the present method and the conventional method. The conventional method used in Ref. [29] is adopted, of which the routine GEMM from Intel MKL is used for matrix multiplication. The traditional HOBEM is employed in both methods, and the boundary integral equations are solved by LU factorization, for which the routine from Intel MKL is adopted as well. The VLFS model and the wave spectrum in this numerical example are the same as the last one, while the case of ω = 1.0 rad/s and = 0◦ is considered only, and the first 30 dry modes of the floating structure are used. Four cases of different element sizes, namely (a) 6 m × 6 m, (b) 4 m × 4 m, (c) 3 m × 3 m and (d) 2.5 m × 2.5 m, are executed. Here the computational time of the most time-consuming processes, i.e., the solution of the boundary integral equations by LU factorization and the formation of response spectra by matrix multiplication, are investigated. The CPU is AMD Opteron Processor 6272 with the speed of a single core 2.10 GHz, and programs for these two methods are executed serially, namely only one core is used. The computational time in different cases is listed in Table 2. In Table 2, the high efficiency of the present calculation scheme with the application of PEM is demonstrated. In each case, it is evident that the computational time for matrix multiplication in the conventional method is much more than that spent in the present method, especially in Case (d) that 12.8 h have to be spent on matrix multiplication in the conventional method, compared with 41 s in the present method. In addition, we can see that the CPU time for matrix multiplication in the conventional method is several times more than that spent on LU factorization in all cases. The reason for the above results is that the computations for LU factorization are O(N3 /3), and the computations for matrix multiplication are O(Nf2 ) in the present method but O(2Nf3 ) in the conventional method. The present VLFS model is far smaller and simpler than that used in actual engineering, and this is the reason why the computational time spent on LU factorization is less in these four cases. However,
when the number of unknowns of VLFS is huge, as 284,929 nodes used in Ref. [22], it is plausible that the present method will spend much less time than the conventional method in the stochastic analysis of VLFS.
5. Conclusions In this study, a new calculation scheme for the stochastic hydroelastic analysis of VLFS subjected to short-crested irregular waves is developed with the application of PEM. The VLFS is modeled as a Mindlin plate, and HOBEM combined with FEM are adopted to decouple the fluid–structure interaction in frequency domain. A directional wave spectrum is used for the short-term description of short-crested irregular waves. The present calculation scheme is verified by three examples, and the results demonstrate its accuracy and high efficiency in dealing with the linear stochastic hydroelastic analysis of VLFS in irregular waves. It is noted that the computations for creating response spectrum matrices using the conventional stationary random vibration theory can occupy more than half of the total computations for the stochastic analysis of VLFS, which can be cut by an order of magnitude by the present calculation scheme, and even two orders of magnitude if the auto-spectrum of the response is required only. Therefore, we conclude that more than half of the total computational time spent on the stochastic analysis of VLFS in the conventional way can be reduced by the present method. During the process of stochastic analyses, the response induced by independent monochromatic waves of many discrete frequencies has to be calculated separately to obtain accurate spectral curves and corresponding variances, so it is proposed to give a fast assessment on the hydroelastic response of VLFS subjected to irregular waves using the present calculation scheme. The wave spectrum at each point of VLFS is uniform in this study. As the length of VLFS may be thousands of meters long, different parts of the structure may be subjected to various wave spectra at the same time. Therefore, multiple wave spectra acting
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
on the floating body which takes inhomogeneous sea states into consideration can be considered for further study. Acknowledgements This research was supported by the National Natural Science Foundation of China (Grant No. 51221961 and 51309051), the National Key Basic Research Special Foundation of China (Grant No. 2013CB036101 and 2010CB832704), and High-Tech Ship Research Projects sponsored by Ministry of Industry and Information Technology of the People’s Republic of China. References [1] Watanabe E, Utsunomiya T, Wang CM. Hydroelastic analysis of pontoon-type VLFS: a literature survey. Eng Struct 2004;26(2):245–56. [2] Wu C, Utsunomiya T, Watanabe E. Application of Galerkin’s method in wave response analysis of flexible floating plates. In: Proceedings of the 6th international offshore and polar engineering conference. 1996. p. 307–14. [3] Kashiwagi M. A B-spline Galerkin scheme for calculating the hydroelastic response of a very large floating structure in waves. J Mar Sci Technol 1998;3(1):37–49. [4] Wu YS. Hydroelasticity of floating bodies. UK: Brunel University; 1984 [Ph.D. dissertation]. [5] Newman JN. Wave effects on deformable bodies. Appl Ocean Res 1994;16(1):47–59. [6] Hamamoto T, Fujita K. Wet-mode superposition for evaluating the hydroelastic response of floating structures with arbitrary shape. In: Proceedings of the 12th international offshore and polar engineering conference. 2002. p. 290–7. [7] Fu SX, Moan T, Chen XJ, Cui WC. Hydroelastic analysis of flexible floating interconnected structures. Ocean Eng 2007;34(11–12):1516–31. [8] Gao RP, Tay ZY, Wang CM, Koh CG. Hydroelastic response of very large floating structure with a flexible line connection. Ocean Eng 2011;38(17–18):1957–66. [9] Meylan MH. A variational equation for the wave forcing of floating thin plates. Appl Ocean Res 2001;23(4):195–206. [10] Bai KJ, Yoo BS, Kim JW. A localized finite-element analysis of a floating runway in a harbor. Mar Struct 2001;14(1–2):89–102. [11] Kyoung JH, Hong SY, Kim BW, Cho SK. Hydroelastic response of a very large floating structure over a variable bottom topography. Ocean Eng 2005;32(17–18):2040–52. [12] Seto H, Ohta M, Ochi M, Kawakado S. Integrated hydrodynamic-structural analysis of very large floating structures (VLFS). Mar Struct 2005;18(2):181–200. [13] Ohkusu M, Namba Y. Hydroelastic analysis of a large floating structure. J Fluids Struct 2004;19(4):543–55. [14] Wu C, Watanabe E, Utsunomiya T. An eigenfunction expansion-matching method for analyzing the wave-induced responses of an elastic floating plate. Appl Ocean Res 1995;17(5):301–10. [15] Kim JW, Ertekin RC. An eigenfunction-expansion method for predicting hydroelastic behavior of a shallow-draft VLFS. In: Proceedings of the 2nd international conference on hydroelasticity in marine technology. 1998. p. 47–59. [16] Ohmatsu S. Numerical calculation method for the hydroelastic response of a pontoon-type very large floating structure close to a breakwater. J Mar Sci Technol 2000;5(4):147–60. [17] Riggs HR, Suzuki H, Ertekin RC, Kim JW, Iijima K. Comparison of hydroelastic computer codes based on the ISSC VLFS benchmark. Ocean Eng 2008;35(7):589–97. [18] Wu YS, Wang DY, Riggs HR, Ertekin RC. Composite singularity distribution method with application to hydroelasticity. Mar Struct 1993;6(2–3):143–63.
213
[19] Wang SQ, Ertekin RC, Riggs HR. Computationally efficient techniques in the hydroelasticity analysis of very large floating structures. Comput Struct 1997;62(4):603–10. [20] Kring D, Korsmeyer T, Singer J, White J. Analyzing mobile offshore bases using accelerated boundary-element methods. Mar Struct 2000;13(4–5):301–13. [21] Newman JN, Lee CH. Boundary-element methods in offshore structure analysis. J Offshore Mech Arct Eng 2002;124(2):81–9. [22] Utsunomiya T, Watanabe E, Nishimura N. Fast multipole algorithm for wave diffraction/radiation problems and its application to VLFS in variable water depth and topography. In: Proceedings of OMAE’01, 20th international conference on offshore mechanics and arctic engineering. 2001. OMAE2001/OSU-5202. [23] Utsunomiya T, Watanabe E. Accelerated higher order boundary element method for wave diffraction/radiation problems and its applications. In: Proceedings of the 12th international offshore and polar engineering conference. 2002. p. 305–12. [24] Hamamoto T. Stochastic fluid–structure interaction of large circular floating islands during wind waves and seaquakes. Probab Eng Mech 1995;10(4):209–24. [25] Chen XJ, Jensen JJ, Cui WC, Tang XF. Hydroelastic analysis of a very large floating plate with large deflections in stochastic seaway. Mar Struct 2004;17(6):435–54. [26] Chen XJ, Moan T, Fu SX, Cui WC. Second-order hydroelastic analysis of a floating plate in multidirectional irregular waves. Int J Nonlinear Mech 2006;41(10):1206–18. [27] Ertekin RC, Kim JW, Xia DW. Hydroelastic response of a mat-type, floating runway near a breakwater in irregular seas. In: OCEANS ‘99 MTS/IEEE. Riding the Crest into the 21st Century. 1999. p. 839–47, http://dx.doi.org/10.1109/OCEANS.1999.804985. [28] Miyajima S, Seto H, Ohta M. Hydroelastic responses of the Mega-Float Phase-II model in waves. In: Proceedings of the 12th international offshore and polar engineering conference. 2002. p. 298–304. [29] Papaioannou I, Gao RP, Rank E, Wang CM. Stochastic hydroelastic analysis of pontoon-type very large floating structures considering directional wave spectrum. Probab Eng Mech 2013;33:26–37. [30] Lin JH. A fast CQC algorithm of psd matrices for random seismic responses. Comput Struct 1992;44(3):683–7, http://dx.doi.org/10.1016/ 0045-7949(92)90401-K. [31] Lin JH, Zhang YH. Seismic random vibration of long-span structures. In: de Silva C, editor. Vibration and shock handbook. Florida: CRC Press; 2005 [chapter 30]. [32] Lin JH, Zhang YH, Zhao Y. Pseudo excitation method and some recent developments. Procedia Eng 2011;14:2453–8. [33] Sun L, Liu CF, Zong Z, Dong XL. Fatigue damage analysis of the deepwater riser from VIV using pseudo-excitation method. Mar Struct 2014;37:86–110. [34] Petyt M. Introduction to finite element vibration analysis. Cambridge: Cambridge University Press; 1998. [35] Wehausen JV, Laitone EV. Surface waves. Encyclopedia of physics, vol. 9. Berlin: Springer-Verlag; 1960. p. 446–778. [36] Li HB, Han GM, Mang HA. A new method for evaluating singular integrals in stress analysis of solids by the direct boundary element method. Int J Numer Methods Eng 1985;21(11):2071–98. [37] Pierson WJ, Tuttell JJ, Woolley JA. The theory of the refraction of a short crested Gaussian sea surface with application to the Northern New Jersey Coast. In: Proceedings of 3rd conference on coastal engineering. 1952. p. 86–108. [38] Mitsuyasu H. On the growth of spectrum of wind generated waves. 2: Spectral shapes of wind waves at final fetch. In: Proceedings of 17th Japanese conference on coastal engineering. 1970. p. 1–7 [in Japanese]. [39] Blackford LS, Demmel J, Dongarra J, Duff I, Hammarling S, Henry G, et al. An updated set of basic linear algebra subprograms (BLAS). ACM Trans Math Soft 2002;28(2):135–51. [40] Intel Corporation. Intel math kernel library reference manual 10.1; 2008. [41] Yago K, Endo H. On the hydroelastic response of box-shaped floating structure with shallow draft (tank test with large scale model). J Soc Naval Arch Jpn 1996;180:341–52 [in Japanese].
Contents lists available at ScienceDirect
Applied Ocean Research journal homepage: www.elsevier.com/locate/apor
Stochastic hydroelastic analysis of a very large floating structure using pseudo-excitation method Hai-tao Li, Qi-bin Wang, Zhi Zong ∗ , Lei Sun, Hui Liang State Key Laboratory of Structural Analysis for Industrial Equipment, School of Naval Architecture, Dalian University of Technology, Dalian 116024, PR China
a r t i c l e
i n f o
Article history: Received 11 April 2014 Received in revised form 1 August 2014 Accepted 1 August 2014 Keywords: VLFS Hydroelasticity Pseudo-excitation method Stationary stochastic response
a b s t r a c t This paper is concerned with the linear hydroelastic response of a pontoon-type very large floating structure (VLFS) in short-crested irregular waves. The linear potential theory is employed for the analysis of VLFS in frequency domain. To decouple the fluid–structure interaction, the higher-order boundary element method (HOBEM) combined with the finite element method (FEM) is adopted. VLFS is modeled as a Mindlin plate, and the mode superposition method is used to reduce the dimension of dynamic equations. The pseudo-excitation method (PEM) is adopted to analyze the stationary stochastic response of the floating structure. The efficiency of this new calculation scheme with the application of PEM is investigated in comparison with the conventional method for stochastic response by analyzing the computational complexity theoretically. Finally, the new calculation scheme is validated by comparing with the experimental data as well as the existing numerical results calculated in the conventional way. In addition, the efficiency of the present numerical approach is also testified which indicates that the proposed numerical scheme is time-saving. © 2014 Elsevier Ltd. All rights reserved.
1. Introduction The very large floating structure (VLFS) is able to expand the living space for human beings. In contrast to the land reclamation from the sea, VLFS has less impact on the surrounding ecosystem and costs less. Because of the advantages above, VLFS has been applied to floating bridges, floating piers, floating fuel storage facilities, and the floating airport model Mega-Float, etc. [1]. Generally, the height of VLFS is only a few meters, which is much smaller than its horizontal dimensions. Therefore, VLFS is usually modeled as a floating plate for simplicity, and it is necessary to consider the hydroelasticity of the structure, which allows for the interaction between the floating elastic body and water waves. Most hydroelastic analyses of VLFS are carried out in regular waves in frequency domain. For linear hydroelastic analyses, the hydrodynamic loads and the dynamic response of the floating structure can be dealt with separately, i.e., firstly the hydrodynamic pressure on the floating structure is generally solved by the boundary element method (BEM), and then the response of the floating structure under the pressure is calculated. Several methods have been used to deal with the deformation of the structure, including the Galerkin method [2,3] combined with the “generalized mode” proposed by Wu [4] and Newman [5], and the mode
∗ Corresponding author. Tel.: +86 411 84707694. E-mail address: rcmh [email protected] (Z. Zong). 0141-1187/$ – see front matter © 2014 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.apor.2014.08.003
superposition method consisting of the wet mode approach [6] and the dry mode approach [7,8], of which the dry mode is obtained by the finite element method (FEM). Different from the methods above, Meylen [9] proposed a variational equation for a floating thin plate subjected to wave forces, which combines the variational equation of the plate with the free-surface Green function method together. BEM can generate an unsymmetrical full matrix, which is difficult to deal with. In addition, irregular frequencies are present if the free-surface Green function is used. Therefore, some other methods are proposed for solving the fluid motion around VLFS, including FEM [10,11] and the hybrid finite/infinite element method [12]. However, it is difficult to analyze the hydroelasticity of VLFS using BEM or FEM for the huge horizontal dimensions of the floating structure. Therefore, analytical methods, and semianalytical methods using eigenfunctions have been widely used for the study of pontoon-type VLFSs of simple geometries and small drafts [13–16], as these methods are less time-consuming and more efficient in the parametric study of the hydroelastic response. In recent years, several computer codes for the hydroelastic analysis of VLFS have been developed, and a detailed comparative study of the simulation results from them was carried out by Riggs et el. [17]. In general, the floating bodies in actual engineering are of complex geometries; in this case, the conventional BEM is generally employed for the hydrodynamic computation of VLFS, but the required computations are extremely huge. To deal with the problem in the computation, some studies of efficient computation
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
have been carried out. Wu et el. [18] proposed the double composite singularity distribution method for the hydroelastic analysis of a floating structure which has two planes of symmetry and nearly 75% of the total equations can be reduced. Wang et el. [19] introduced a cut-off criterion for the Green function and its derivatives to convert full matrices formed by BEM into sparse matrices, and then applied an iterative sparse solver to the calculation of sparse matrices to reduce the CPU time. Kashiwagi [3] adopted bi-cubic B-spline functions to represent the unknown pressure, and employed “relative similarity relations” to evaluate the coefficient matrix formed by the pressure-distribution method, which can reduce the computational time drastically with a small number of unknowns required only. Recently, another two fast algorithms, namely the pre-corrected Fast Fourier Transform (pFFT) method and the fast multipole method (FMM), have been applied to the acceleration of BEM and make it possible to study the hydrodynamics of large floating structures of general geometries. The pFFT method has been applied to the analysis of the Mobile Offshore Base (MOB) [20,21] and can reduce the computational time and memory cost for BEM from the order O(N2 ) to O(N log N). Utsunomiya et el. [22,23] have developed another form of the higher-order boundary element method (HOBEM) accelerated by FMM, which can reduce the storage requirement and CPU time to O(N) and O(N log N), respectively. Some investigations have been carried out on the stochastic hydroelastic analysis of VLFS for a reliable design, which will be subjected to irregular waves during its whole service life. Hamamoto [24] analyzed the stochastic response of a circular floating island subjected to wind-induced waves under both long-term and shortterm descriptions of loadings, of which the short-term description of external loadings was represented by spectral density functions. Chen et al. [25,26] carried out the research on the second order response of VLFS induced by coupling of first-order wave potentials or the membrane force of the plate in multidirectional irregular waves. Ertekin et el. [27] investigated the hydroelastic response of a pontoon-type VLFS sheltered by a breakwater in regular and irregular waves based on the linear Green–Naghdi theory and the eigenfunction-expansion matching method. Using the same method as Ref. [12], Miyajima et el. [28] analyzed the hydroelastic response of the Mega-Float Phase-II model in short-crested irregular waves and compared their results with the measured values. Refs. [27,28] adopted the formula for linear systems of single degree of freedom (SDOF) to evaluate the response spectrum of VLFS, which means that the cross-correlation between responses of any two points of the floating structure is neglected, and only approximate results can be obtained in this way. Recently, Papaioannou et el. [29] studied the linear stochastic response of VLFS subjected to a directional wave spectrum by the stationary random vibration theory, which is an accurate method to the stochastic analysis of systems of multiple degrees of freedom (MDOF) without any cross-correlation terms neglected. However, the response spectra of VLFS are quite difficult to calculate due to the computations of O(N3 ) by the stationary random vibration theory, in which case the required number of unknowns N is huge. In recent years, Lin et el. [30,31] have developed the pseudo-excitation method (PEM), which is efficient in dealing with stochastic responses of linear systems of large degrees of freedom, and has been applied to stochastic analyses of long-span bridges, high-speed trains [32], and the fatigue damage of the deepwater riser [33], etc. In this study, we adopt PEM to study the linear hydroelasticity of a pontoon-type VLFS subjected to short-crested irregular waves in frequency domain. VLFS is modeled as a Mindlin plate with allowance for the effect of transverse shear deformation and rotary inertia. The fluid motion is solved by HOBEM, and the deformation of the floating structure is calculated by the mode superposition
203
Fig. 1. The schematic diagram of the fluid–structure interaction problem.
method, of which dry modes of the structure obtained by FEM are chosen. A directional wave spectrum is used for the short-term description of short-crested irregular waves. The efficiency of this calculation scheme with the application of PEM is discussed theoretically and verified by a numerical example at the end of this paper. Finally, the accuracy of the present method is also verified by comparing our numerical results with the experimental data as well as the existing numerical results obtained by the stationary random vibration theory. Some meaningful conclusions about the present calculation scheme are proposed. 2. Hydroelastic analysis of VLFS In this section, the linear hydroelasticity of VLFS in regular waves is analyzed in frequency domain. The schematic diagram of the fluid–structure interaction problem is shown in Fig. 1. The pontoontype VLFS is defined with length L, width B, height T, and modeled as a floating plate assumed to be flat with free edges. A zero-draft assumption is adopted for simplicity. The symbols ˝, Sf , Sb and Sd denote the fluid domain, the free surface, the bottom surface of the structure and the flat seabed, respectively. Three-dimensional Cartesian coordinate system o-xyz is established with o-xy plane coinciding with the undisturbed free surface and z-axis orienting positively upwards. 2.1. Equations of motion for the plate The height of VLFS is much smaller than its length and width, so the floating structure is usually modeled as a Kirchhoff plate which does not allow for transverse shear deformation. It is accurate to use the Kirchhoff plate theory when the plate is homogeneous and thin, in which case the effect of shear deformation is of no significance. However, if VLFS consists of different kinds of materials which behaves like sandwich plates, the effect of shear deformation cannot be neglected. In addition, when analyzing the wave-induced transverse shear force of VLFS, Gao et al. [8] found that the maximum shear force obtained by the Mindlin plate theory is about 10% larger than the one by the Kirchhoff plate theory. Therefore, in this study we model VLFS as an isotropic Mindlin plate for a more accurate evaluation of shear forces of the structure. The displacements of a Mindlin plate are given as [34] UP (x, y, z, t) = zy (x, y, t)
(1a)
VP (x, y, z, t) = −zx (x, y, t)
(1b)
WP (x, y, z, t) = W (x, y, t)
(1c)
where UP , VP , and WP denote the displacements along x-axis, y-axis, and z-axis, respectively. x , y represent the rotations about x-axis and y-axis, respectively, and W is the vertical displacement of the neutral surface of the plate.
204
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
A two-dimensional 8-node isoparametric element is used to build the FEM model. Then the stiffness matrix [K]Nf ×Nf and the consistent mass matrix [M]Nf ×Nf are obtained by standard FEM procedures, where Nf denotes the number of displacements of the plate. The free boundary condition should be satisfied, namely the bending moment, the transverse shear force and the twisting moment must vanish at the free edges of the plate. The steady state of the structural motion is considered, so the displacements of the plate can be expressed as x (x, y, t) = R[x (x, y)e−iωt ] −iωt
y (x, y, t) = R[y (x, y)e
(2a) (2b)
]
W (x, y, t) = R[w(x, y)e−iωt ]
(2c)
where √ x , y and w are the complex displacements of the plate. i = −1 is the imaginary unit, and ω is the circular frequency. R denotes the real part of a complex expression. Regardless of the structural damping, dynamic equations for the structure are given as
(−ω2 [M] + [K]) d where
d
Nf
=
=
F
d1 , d2 , . . .dN
(3)
T
and
F
Nf
are global vectors of
the complex displacement and the complex load, respectively, of which the displacement vector of the ith node is given as di =
T
xi , yi , wi . The subscript N denotes the number of nodes of the plate and satisfies Nf = 3N in the rest of this paper. Assuming that the free vibration of the plate has a timeharmonic solution with the frequency ωi , the natural modes of the plate can be obtained by the following eigenvalue simultaneous equations:
(−ωi2 [M] + [K])ci
=0
(4)
where ωi and ci are the ith natural frequency and dry mode vector of the Mindlin plate, respectively. 0 denotes the zero vector. 2.2. Equations of motion for the fluid domain The fluid is assumed to be ideal and incompressible, and its motion is irrotational so that the velocity potential ˚ exists. The linear potential theory is employed based on the assumption of small fluid and structural motions. The velocity potential and the dynamic pressure on the plate can be written as ˚(x, y, z, t) = R[(x, y, z)e−iωt ] −iωt
P(x, y, t) = R[p(x, y)e
]
(5b)
(6)
where I (x, y, z), D (x, y, z) and R (x, y, z) denote the incident wave potential, the diffraction potential and the radiation potential, respectively. The incident wave potential has the following form according to the linear wave theory: I (x, y, z) = −
igA cosh k(z + h) ik(x cos +y sin ) e ω cosh kh
∇ 2 = 0 in ˝
(7)
where g and A represent the gravitational acceleration and the wave amplitude, respectively. h and represent the water depth and the propagation direction of the incident wave, respectively, as illustrated in Fig. 1. k denotes the wave number for finite water depth which satisfies the linear wave dispersion relationship ω2 = gk tanh kh.
(8a)
ω2
∂ on Sf = g ∂z
(8b)
∂ = −iωw(x, y) on Sb ∂z
(8c)
∂ = 0 on Sd ∂z
(8d)
lim
√ r
r→∞
∂ − ik ∂r
( − I ) = 0 on S∞
(8e)
where r denotes the horizontal distance between the origin and the referred point, and S∞ represents the artificial fluid boundary at infinity. Eq. (8b) is the linearized free-surface boundary condition, which implies that the normal velocity of the fluid is equal to that of the free surface and the pressure on the free surface is atmospheric. Eq. (8c) is the plate surface condition and implies that there is no gap between the plate and the surface of the fluid domain. Eq. (8d) is the seabed boundary condition which implies the impermeability of the seabed. Eq. (8e) is the Sommerfeld radiation condition for D and R as r→ ∞ which indicates the propagation of diffracted and radiated waves is outgoing at infinity. By substituting Eq. (6) into Eq. (8c), we obtain another form of the plate surface condition, i.e.,
∂D ∂I =− ∂z ∂z
on Sb
(9a)
∂R = −iωw(x, y) on Sb ∂z
(9b)
The free-surface Green function for finite water depth [35] which satisfies the linearized free-surface boundary condition, the seabed boundary condition and the radiation condition at infinity is used, i.e.,
G(x, x0 ) = 2i
(5a)
where (x, y, z) and p(x, y) are the complex velocity potential and the complex dynamic pressure, respectively. The complex velocity potential consists of three parts based on the linear potential theory, i.e., (x, y, z) = I (x, y, z) + D (x, y, z) + R (x, y, z)
The velocity potential must satisfy Laplace’s equation and the following boundary conditions:
k02 − 2 h(k02 − 2 ) +
∞
+ 4
n=1
kn2 + 2 h(kn2
(1)
cosh k0 (z + h) cosh k0 (z0 +h)H0 (k0 R)
+ 2 ) −
cos kn (z + h) cos kn (z0 + h)K0 (kn R) (10)
where x = (x, y, z) and x0 = (x0 , y0 , z0 ) are the field point and the source point, respectively. = ω2 /g is the wave number for deep water. k0 and kn are the real roots of ω2 = gk0 tanh k0 h and ω2 = − gkn tan kn h, respectively. R represents the horizontal dis(1) tance between the source point and the field point. H0 (x) is the Hankel function of the first kind of order zero, and K0 (x) is the modified Bessel function of the second kind of order zero. Applying the free-surface Green function and Laplace’s equation to Green’s second identity, we obtain the following boundary integral equation:
C(x)(x)+
(x0 ) Sb
∂G(x, x0 ) ds = ∂z
G(x, x0 ) Sb
∂(x0 ) ds (11) ∂z
where C(x) is the solid angle at the field point, and equal to 2, 3, 3.5 if the field point is located at the bottom surface, the edges, the corners of the plate, respectively.
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
By discretizing the plate surface into a finite number of twodimensional 8-node isoparametric elements, we obtain the velocity potential of one point, i.e., =
8
Ni i = N
(12)
e
i=1
[A] D
[A] R
= [B]
= [B]
∂D ∂z ∂R ∂z
= −[B]
∂I ∂z
R
= [A]−1 [B]
(13b)
∂R ∂z
= −iω[A]−1 [B] {w}
p(x, y) = iω(x, y, 0) − gw(x, y)
(15)
where is the density of water. And the vector of the pressure induced by the radiation potential and the vertical motion of the plate is {pW } = ω2 [A]−1 [B] {ω} − g {ω}
(16)
Then the element added mass matrix [ma ], the element wave damping matrix [ca ] and the element hydrostatic stiffness matrix [ka ] are obtained by the integration of pW over an element, i.e., [ma ] = [S]R([a])
(17a)
[ca ] = ω[S]I([a])
(17b)
where [S]8×8 =
(17c)
T
Se
N N ds is an integration over an element.
[a]8×8 is a submatrix of [A]−1 [B] corresponding to the locations of 8 potentials of an element in {R }. I denotes the imaginary part of a complex expression. The element vector of the exciting force is obtained by
FS
e
NT N ds
= iω Se
I
e
+ D
{d} =
e
q
i ci = [c]
(19)
i=1
where { }q = { 1 , 2 , . . ., q }T consists of complex amplitudes of the first q dry modes of the plate, and [c]Nf ×q = [c1 , c2 , . . ., cq ] contains the first q dry mode vectors obtained in Eq. (4). After substituting [Ma ], [Ca ], [Ka ] and {FS } into Eq. (3), another form of simultaneous equations for fluid–structure interaction in frequency domain is obtained, i.e.,
(−ω2 [M] − ω2 [Ma ] − iω[Ca ] + [K] + [Ka ]) d
=
FS
(20)
By substituting Eq. (19) into Eq. (20), and multiplying both sides of Eq. (20) by [c]T , we obtain the ultimate form of simultaneous equations for fluid–structure interaction, of which the dimension is reduced by the mode superposition method, i.e.,
[c]T (−ω2 [M] − ω2 [Ma ] − iω[Ca ] + [K] + [Ka ])[c]
= [c]T FS
(21)
(14)
where {w}N = {w1 , w2 , . . ., wN }T contains the vertical displacement at each point of the plate. Then the pressure distribution on the bottom surface of the plate is obtained from linearized unsteady Bernoulli’s equation, i.e.,
[ka ] = g[S]
The mode superposition method is used to solve the displacement field of the plate, i.e.,
(13a)
The entries of [A]N×N and [B]N×N can be evaluated by numerical integration methods. When the source point and the field point are located in different elements, the boundary integral can be evaluated by the Gauss–Legendre quadrature rule directly. But if both the source and the field points are located in the same element, the Green function and its derivatives approach infinity at the rates of 1/R and 1/R2 , respectively; for this case, the triangular polar coordinate transformation proposed by Li et al. [36] is adopted to evaluate the boundary integral. By substituting Eq. (9b) into Eq. (13b), we obtain
is obtained by assembling the subvector in Eq. (18). It should be noted that the entries of the above matrices and vector, which are corresponding to the rotational displacements in {d}, namely x and y , are equal to zero as no hydrodynamic forces are assumed to be induced by rotational displacements of the plate. 2.3. Dynamic equations for fluid-structure interaction
where N = {N1 , N2 , . . ., N8 } is the shape function vector of the 8-node isoparametric element, and {}e = {1 , 2 , . . ., 8 }T is the potential vector that contains potentials at 8 nodes of an element. By substituting Eqs. (9a) and (12) into Eq. (11), we obtain the discretized form of the boundary integral equation [6], i.e.,
205
(18)
where the subscript e denotes the component of an element. Finally the global added mass matrix [Ma ]Nf ×Nf , the global wave damping matrix [Ca ]Nf ×Nf , and the global hydrostatic stiffness matrix [Ka ]Nf ×Nf are constructed by assembling the submatrices in
Eqs. (17a)–(17c), and the global vector of the exciting force FS
Nf
3. Stochastic analysis of VLFS 3.1. Directional wave spectrum Generally, irregular wind waves are assumed to be stationary, ergodic and can be described by a zero mean Gaussian process during a short period. The characteristics of short-crested irregular waves can be represented by a directional wave spectrum S(ω, ), which represents the distribution of wave energy on wave frequencies and propagation directions. The response of a linear system satisfies the Gaussian distribution if the input can be described by a Gaussian process. Therefore, the frequency characteristics of the response induced by short-crested irregular waves can be described as a directional spectrum. In this section, the linear stochastic response of VLFS is investigated in frequency domain. The relationship between a directional wave spectrum S(ω, ) and a one-dimensional frequency spectrum S(ω) can be expressed as S(ω, ) = S(ω)D(|ω)
(22)
where D(|ω) is the directional spreading function representthe directional distribution of wave energy, and it satisfies
ing
D(|ω) d = 1. The free-surface elevation of irregular wind waves (x, y, t) can be obtained by the linear superposition of infinite independent monochromatic waves of various frequencies and propagation directions, i.e., −
(x, y, t) = R
∞ ∞
i[kn (x cos m +y sin m )−ωn t+εmn ]
amn e
(23)
m=1 n=1
where kn , m , and ωn are the wave number, the wave direction and the circular frequency of a monochromatic wave, respectively. εmn
206
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
denotes the random phase distributed uniformly in [0, 2], and amn is the amplitude of a monochromatic wave and can be obtained by
time-harmonic excitation, which is referred to as the pseudoexcitation, is given as
amn =
x(t) =
(24)
2S(ω, )ıωı
where ıω and ı denote small intervals around the wave frequency ω and the wave direction , respectively. We assume the independence of the directional distribution on wave frequencies, and propagation directions of waves are assumed to range from −/2 to /2 around the mean wave direction . The directional spreading function proposed by Pierson et el. [37] is adopted, i.e.,
D( | ω) = D() =
⎧2 ⎨ cos2 ( − ) | − | ≤
2 | − | > 2
⎩0
(25)
S(ω, ) d
(26)
−/2
In this study, the two-parameter Bretschneider–Mitsuyasu spectrum that describes the frequency characteristic of fully developed wind waves [38] is adopted, i.e., SBM (ω) =
−4 2 0.257H1/3 T1/3
ω −5 2
ω exp −1.03 T1/3 2
−4
(27)
where H1/3 and T1/3 are the significant wave height and the significant wave period, respectively.
For a linear MDOF system, the relationship between a timeharmonic excitation with the frequency ω and the corresponding steady-state response can be described as
(28)
where [H(ω)] is the transfer function matrix for a linear MDOF system. Then, the response spectrum matrix and cross-spectrum matrices of a stationary random process can be obtained in the stationary random vibration theory, i.e., [Syy (ω)] = [H(ω)]∗ [Sxx (ω)][H(ω)]T
(29a)
[Sxy (ω)] = [Sxx (ω)][H(ω)]T
(29b)
∗
[Syx (ω)] = [H(ω)] [Sxx (ω)]
y(t) = HT (ω)x(t) = HT (ω) Sxx (ω) e−iωt
(31)
where HT (ω) is the transfer function for a linear SDOF system. We can obtain the auto-spectrum of the response and crossspectra between the excitation and the response by
2
∗
(29c)
where [Sxx (ω)] and [Syy (ω)] denote spectrum matrices of the excitation and the response, respectively. [Sxy (ω)] and [Syx (ω)] are cross-spectrum matrices between the excitation and the response. The superscript * denotes complex conjugate. Eqs. (29a)–(29c) have simple forms and are comparatively convenient for engineering applications, which is the conventional way to analyze the stochastic response. However, if the number of degrees of freedom M is large, the conventional method will be very time-consuming due to the O(M3 ) multiplications of the traditional algorithm for matrix multiplication. Therefore, PEM is applied to the stochastic analysis of VLFS in this study due to its high efficiency in dealing with stochastic responses of linear systems of large degrees of freedom. The stationary random vibration of a linear SDOF system is taken as an example to explain the basic principle of PEM. A
(32a)
x˜ (t) y(t) = Sxx (ω)HT (ω) = Sxy (ω)
(32b)
y(t)∗x(t) = HT (ω)∗ Sxx (ω) = Syx (ω)
(32c)
The process involving Eqs. (30)–(32c) is known as PEM for a linear SDOF system, and the response spectrum matrix and crossspectrum matrices between the excitation and the response of a linear MDOF system can be obtained using PEM in the similar way, i.e., [Syy (ω)] = [Sxy (ω)] = [Syx (ω)] =
∗
y˜ (t)
∗
x˜ (t)
3.2. Pseudo-excitation method [30,31]
{y(t)} = [H(ω)] x(t)
(30)
where Sxx (ω) is the auto-spectrum of the stationary random excitation x(t). Then the corresponding pseudo-response is obtained by substix(t) for x(t) in the formula y(t) = HT (ω) x(t), i.e., tuting
/2
S(ω) =
Sxx (ω) e−iωt
y˜ (t)∗ y˜ (t) = HT (ω) Sxx (ω) = Syy (ω)
Therefore, the relationship between the one-dimensional frequency spectrum and the directional wave spectrum can be written as
(33a)
T
y˜ (t)
∗
y˜ (t)
T
y˜ (t)
(33b)
T
x˜ (t)
(33c)
where x˜ (t) and y˜ (t) are vectors of the pseudo-excitation and the pseudo-response, respectively. It is obvious that the results obtained from Eqs. (33a)–(33c) are identical to those from Eqs. (29a)–(29c), which demonstrates that PEM is equivalent to the stationary random vibration theory when dealing with the stationary stochastic response of a linear MDOF system. Moreover, if the auto-spectrum of the response alone is desired, the auto-spectrum of the ith response can be obtained by [Syy (ω)]i,i =
∗
y˜ (t)
i
y˜ (t)
i
(34)
where y(t) is the ith entry of the pseudo-response vector. i
Eq. (26) assumes that components of the wave spectrum in different wave directions are independent of each other. As the transfer function from wave loads to the response of the floating structure is independent of wave directions as well, the total response spectrum can be obtained by the superposition of response spectra induced by waves of different directions, which satisfies Eq. (26) as well. Therefore, firstly we study the displacement spectrum of VLFS in single wave direction with PEM, and then the total displacement spectrum is obtained by the superposition of components in all directions. 3.3. Pseudo-excitation method for VLFS The component of the response spectrum with the frequency ω and the wave direction is considered. Based on the zero-draft assumption of VLFS, the relationship between the incident wave potential and the free-surface elevation is I (x, y, 0) = HI , (ω) (x, y) where HI , (ω) = −ig/ω is the transfer function from to I .
(35)
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
Incident wave potentials at any two points, namely the ith point (xi , yi ) and the jth point (xj , yj ) of the plate are obtained from Eq. (7), i.e., igA ik(x I (xi , yi , 0) = − e i ω I (xj , yj , 0) = −
igA ik(x e j ω
cos +yi sin )
(36a)
cos +yj sin )
(36b)
So the relationship between incident wave potentials at these two points can be written as ik(xij cos +yij sin )
I (xj , yj , 0) = I (xi , yi , 0) e
(37)
where xij = xj − xi , yij = yj − yi , and eik(xij cos +yij sin ) indicates that a phase-lag exists between incident wave potentials at two different points. The phase-lag is caused by the propagation of water waves, which is referred to as the wave passage effect. As VLFS is generally several kilometers long, a large phase difference exists between excitations of the two points with a large spatial distance. Therefore, the wave passage effect must be taken into account for the stochastic analysis of VLFS. Considering that a linear MDOF system is subjected to a waveinduced stationary random excitation in the form of F0 (t) at the origin of the coordinate system, the excitations of other points can be given as
FE (t)
⎧ ⎫ F1 (t) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ F2 (t) ⎪ ⎬
=
N
⎪ ⎪ ⎪ ⎪ ⎩
.. . FN (t)
⎪ ⎪ ⎪ ⎪ ⎭
=
⎧ ⎫ F0 (t − t1 ) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ F0 (t − t2 ) ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎩
.. . F0 (t − tN )
⎪ ⎪ ⎪ ⎪ ⎭
FE (t)
=
N
. ⎪ ⎪ ⎪ ⎪ . ⎪ ⎪ ⎪ ⎩ . ⎪ ⎭
(38)
N
=
where H∂I /∂z,I = k tanh kh is the transfer function from I to Then applying the pseudo-excitation vector of the incident wave potential to Eq. (41), we obtain the following pseudo-excitation vector:
∂I ∂z
=
(ω, )
I (ω, )
N
H∂I /∂z,I
(42)
N
vector of the diffraction potential The pseudo-excitation D (ω, )
is obtained by substituting Eq. (42) into Eq. (13a).
Substituting the element pseudo-excitation vectors of the incident wave potential and the diffraction potential into Eq.(18), weobtain the element vector of the pseudo exciting force F˜S (ω, ) and
FS (ω, )
sequently the global one
by assembling
each element. The pseudo-amplitude vector of dry modes
FS (ω, )
e
F˜S (ω, )
e
of
(ω, ) is
into Eq. (21). Then we back-
the pseudo-displacement vector
d(ω, ) and sequently the dis-
placement spectrum matrix in the wave direction , i.e.,
∗ T d(ω, ) d(ω, ) Nf
(43)
Nf
Then, we substitute the results of Eq. (43) into Eq. (26) and obtain the displacement spectrum matrix [Sdd (ω)], and the jth order moment of the auto-spectrum of the ith displacement in the global displacement vector mi,j can be obtained by
∞
ωj [Sdd (ω)]i,i dω
(44)
0
SFF (ω) e−iωt
.. . eik(xN
(41)
∂I /∂z.
mi,j =
⎧ ik(x1 cos +y1 sin ) ⎫ e ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ eik(x2 cos +y2 sin ) ⎪ ⎬ ⎪ ⎪ ⎪ ⎪ ⎩
= H∂I /∂z,I I z=0
[Sdd (ω, )]Nf ×Nf =
(39)
As the time-harmonic factor e−iωt in a pseudo-excitation and its conjugate eiωt always come in pairs and are offset after multiplication, the time-harmonic factor is extracted from the pseudo-excitation and the pseudo-response afterwards. In this study, the directional wave spectrum is assumed to be uniform at each point of the plate, so the pseudo-excitation vector of the incident wave potential in the wave direction , which allows for the wave passage effect, can be written as
I (ω, )
∂I ∂z
substitute the pseudo-amplitude vector into Eq. (19) and obtain
eiωtN
The derivative of I in the z direction is
obtained by substituting
where t1 , t2 , . . ., tN denote time-lags between excitations of the origin and other points. SFF (ω) e−iωt for Substituting the pseudo-excitation F0 (t) = F0 (t) in Eq. (38), where SFF (ω) is the auto-spectrum of F0 (t), we obtain the pseudo-excitation of each point considering the wave passage effect, i.e.,
⎧ iωt1 ⎫ e ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ eiωt2 ⎪ ⎪ ⎨ ⎬
207
cos +yN sin )
⎪ ⎪ ⎪ ⎪ ⎭
S
(ω)D() HI , (ω) (40)
where k(xi cos + yi sin ) denotes the phase-lag between the ith point and the origin, and is equal to ωti as shown in Eq. (39). It is obvious that the spectrum matrix of the incident wave potential in the direction can be obtained by substituting Eq. (40) into Eq. (33a), which has the same form as that in Ref. [29].
where [Sdd (ω)]i,i is the ith diagonal entry of [Sdd (ω)]. The variance of the displacement can be obtained by setting j = 0. 3.4. The comparison of computational efficiency between the present calculation scheme and the conventional method When the traditional BEM is used to analyze the hydrodynamics of floating bodies, the computational time for factorization solvers and iterative solvers is O(N3 ) and O(N2 ), respectively, which is huge for the analysis of VLFS. Therefore, some fast algorithms, such as the pFFT method and FMM, have been applied to the acceleration of BEM to make it possible to solve the hydroelastic response of the floating structure. However, there still exists a problem when dealing with the stochastic response of VLFS in the conventional way, as studied in Ref. [29], i.e., one has to calculate Eq. (29a) twice to obtain the load spectrum and the displacement spectrum, which needs the computations of O(2Nf3 ) using the traditional algorithm for matrix multiplication. In this case, the computational time for creating response spectrum matrices is several times more than that for solving the discrete boundary integral equations by LU factorization (with the computations of O(N3 /3)), and even more than that for the same solution using iterative methods or accelerated BEM. Some favorite math libraries, including the Basic Linear Algebra Subprograms (BLAS) [39] and the Intel Math Kernel Library (MKL) implementation of BLAS [40], provide highly efficient algorithms
208
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
Table 1 Parameters of the VLFS model. Parameter
Symbol
Unit
Value
Total length Total width Total height Young’s modulus Poisson’s ratio Plate density Water density Water depth
L B T E p p h
m m m N/m2
300 60 2 1.19 × 1010 0.13 256.25 1025 58.5
kg/m3 kg/m3 m
for matrix multiplication. The computational time for multiplication of two square matrices of order M is O(M3 ) by routine GEMM involved in both libraries, yet it is still huge for problems with large number of unknowns. So the major concern about the fast analysis of VLFS in irregular waves is to reduce the operations for creating response spectrum matrices, which can be achieved with the application of PEM. The formation of response spectrum matrices using PEM just requires the matrix multiplication of two vectors once with the computations of O(Nf2 ), as shown in Eq. (33a), and even O(Nf ) if the auto-spectrum is required only, as described in Eq. (34). During the process of the stochastic analysis, Eqs. (29a), (33a) or (34) have to be executed for many discrete frequencies to ensure the accuracy of spectral curves and the corresponding variances. Therefore, the present calculation scheme using PEM has a great advantage over the conventional method in stochastic analyses of floating bodies with large degrees of freedom. As the computational time for creating response spectrum matrices in the conventional way is greater than that for solving the boundary integral equations, which are the most time-consuming parts in the stochastic analysis of VLFS, we conclude that more than half of the total computational time can be reduced by the present calculation scheme. 4. Numerical results and discussion 4.1. Verification of the accuracy of the present calculation scheme The accuracy of the present calculation scheme is verified by comparing our numerical results with a group of experimental data as well as other numerical results. The parameters of the VLFS model used here are identical to those in Ref. [29], as listed in Table 1. 4.1.1. Verification of the present calculation scheme by experimental data Due to the lack of experiment data on the stochastic response of a VLFS model, in this part the accuracy of the present calculation scheme is verified by comparing our numerical results with the experimental data of Yago et al. [41], which provides ratios between amplitudes of the vertical displacement and the wave in regular
waves of various wavelengths and wave directions. However, it should be noted that the present method is used for the hydroelastic analysis of VLFS in short-crested irregular waves, rather than that in regular waves. Hence, we compare the response of VLFS calculated by the present method with the experimental data in the case that the floating structure is subjected to a one-dimensional frequency spectrum S(ω), as long-crested irregular waves can be generated by a linear superposition of infinite independent monochromatic waves with a uniform wave direction. The process for the stochastic analysis in single wave direction, as described in Eqs. (40)–(43), is employed for this verification. For a given wave frequency ω and wave direction , the amplitudes of a monochromatic wave A(ω) and the vertical displacement w(ω, ) of one point can be expressed as A(ω) =
2S(ω)ıω
(45a)
w(ω, ) = 2Sww (ω, )ıω
(45b)
where S(ω) is the one-dimensional frequency spectrum as shown in Eq. (27), and Sww (ω, ) is the auto-spectrum of the vertical displacement obtained from Eq. (43). As a consequence, the ratio between amplitudes of the vertical displacement and the wave can be obtained by
w(ω, ) A(ω)
!
=
Sww (ω, ) S(ω)
(46)
As the test model in Ref. [41] is similar to the present one, the non-dimensional displacements measured in the experiment should coincide with those calculated by the present method theoretically. Therefore, we compare the calculated amplitude ratios along the longitudinal centerline of the present VLFS model with the experimental data. In this part, four wave directions, namely = 0◦ , 30◦ , 60◦ and 90◦ , and five wavelength to structure length ratios, namely /L = 0.2, 0.4, 0.6, 0.8 and 1.0 (the corresponding wave frequencies are 1.01 rad/s, 0.72 rad/s, 0.58 rad/s, 0.48 rad/s and 0.42 rad/s, respectively), are chosen for a thorough comparison. 500 8-node elements with 6 m × 6 m are adopted to discretize the plate, and the first 30 dry modes of the structure are used. The significant wave height and wave period are H1/3 = 2 m and T1/3 = 6.3 s, respectively. The comparative results are exhibited in Figs. 2–5, in which solid lines and scatters represent the numerical results calculated by the present method and the experimental data, respectively. It is clear to see that the agreement between our numerical results and the experimental data is good for all cases, except for a little difference in the case of /L = 0.4 and = 60◦ ; however, it is evident that the tendency of the numerical result in this case matches well with the experimental data. Based on the assumption of the independence of spectrum components in different wave directions, the total response spectrum can be obtained by the superposition of response spectra induced
Fig. 2. The comparison of the vertical displacement amplitude along the longitudinal centerline of VLFS for = 0◦ .
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
209
Fig. 3. The comparison of the vertical displacement amplitude along the longitudinal centerline of VLFS for = 30◦ .
Fig. 4. The comparison of the vertical displacement amplitude along the longitudinal centerline of VLFS for = 60◦ .
Fig. 5. The comparison of the vertical displacement amplitude along the longitudinal centerline of VLFS for = 90◦ .
by waves of different directions. Therefore, due to the excellent agreement between our numerical results and the experimental data for various wave frequencies and directions, it is plausible that the present calculation scheme is of enough accuracy for the analysis of VLFS subjected to short-crested irregular waves. 4.1.2. Verification of the present calculation scheme by other numerical results In this part, the accuracy of the present calculation scheme is verified by comparing our numerical results with those in Ref. [29] in the case that VLFS is subjected to a directional wave spectrum. For a more intuitive comparison, the present VLFS model, as shown in Table 1, and the stochastic wave conditions are identical to those in Ref. [29]: 2000 8-node elements with 3 m × 3 m are adopted to discretize the plate; the significant wave height and wave period are H1/3 = 2 m and T1/3 = 6.3 s, respectively; four mean wave directions, namely = 0◦ , 30◦ , 60◦ and 90◦ are chosen to investigate the effect of incident mean wave directions on the hydroelastic response of VLFS; four corner points (P1–P4) and the center point (P5), as shown in Fig. 6, are selected to study the stochastic response of the floating structure. The only difference is that the first 60 dry modes are used here for a more accurate analysis of the response
Fig. 6. Positions of 5 selected points.
of the floating structure, while the first 30 dry modes are adopted in the reference. Figs. 7–10 show the spectral curves of vertical displacements at 5 selected points for four mean wave directions. It is clear to see that the results in this study are in good agreement with those in Ref. [29], and the frequency of the peak point in the spectral curve obtained by the present method is consistent well with that in the reference in each case. However, the values of vertical displacement spectra in this study are a few smaller; in the case of = 30◦ , the frequency of the second peak point in the spectral curve of P3
210
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
Fig. 7. The comparison of response spectra of vertical displacements at 5 selected points for = 0◦ .
Fig. 8. The comparison of response spectra of vertical displacements at 5 selected points for = 30◦ .
Fig. 9. The comparison of response spectra of vertical displacements at 5 selected points for = 60◦ .
Fig. 10. The comparison of response spectra of vertical displacements at 5 selected points for = 90◦ .
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
211
Fig. 11. The comparison of standard deviations of vertical displacements for = 0◦ . The left result is obtained from Ref. [29], and the right is calculated by the present method.
Fig. 12. The comparison of standard deviations of vertical displacements for = 30◦ . The left result is obtained from Ref. [29], and the right is calculated by the present method.
Fig. 13. The comparison of standard deviations of vertical displacements for = 60◦ . The left result is obtained from Ref. [29], and the right is calculated by the present method.
is about 1.15 rad/s in this study compared with 1.4 rad/s in the reference; regarding the case of = 60◦ , a large discrepancy on the frequency of the peak point in the spectral curve of P1 exists, and the peak value of P1 is greater than that of P2 in this study, which is opposite to the result of the reference. The probable reason is that different BEM models cause the aforementioned discrepancies, rather than the inaccuracy of the present method, as PEM is equivalent to the stationary random vibration theory which is employed in the reference.
In Figs. 11–14, the distributions of standard deviations of vertical displacements for four mean wave directions are presented. The left four figures denote the results given in Ref. [29], and the right four are obtained by the present method. It is clear that the distribution of standard deviations calculated by the present method agrees well with that of the reference in each case, and the maximum discrepancy of standard deviations is less than 10%, which is caused by the results that the spectra of vertical displacements in this study are smaller, as demonstrated in Figs. 7–10.
212
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
Fig. 14. The comparison of standard deviations of vertical displacements for = 90◦ . The left result is obtained from Ref. [29], and the right is calculated by the present method. Table 2 The computational time for different element sizes. Case
Number of nodes
(a) (b) (c) (d)
1621 3556 6241 8929
Present method
Conventional method
LU factorization
Matrix multiplication
LU factorization
Matrix multiplication
6s 59 s 5.3 min 16.4 min
1s 5s 16 s 41 s
6s 1 min 5.3 min 16.9 min
4.9 min 49 min 4.6 h 12.8 h
4.2. Verification of the efficiency of the present calculation scheme To testify the high efficiency of the present calculation scheme with the application of PEM, a comparison is conducted between the computational time of the present method and the conventional method. The conventional method used in Ref. [29] is adopted, of which the routine GEMM from Intel MKL is used for matrix multiplication. The traditional HOBEM is employed in both methods, and the boundary integral equations are solved by LU factorization, for which the routine from Intel MKL is adopted as well. The VLFS model and the wave spectrum in this numerical example are the same as the last one, while the case of ω = 1.0 rad/s and = 0◦ is considered only, and the first 30 dry modes of the floating structure are used. Four cases of different element sizes, namely (a) 6 m × 6 m, (b) 4 m × 4 m, (c) 3 m × 3 m and (d) 2.5 m × 2.5 m, are executed. Here the computational time of the most time-consuming processes, i.e., the solution of the boundary integral equations by LU factorization and the formation of response spectra by matrix multiplication, are investigated. The CPU is AMD Opteron Processor 6272 with the speed of a single core 2.10 GHz, and programs for these two methods are executed serially, namely only one core is used. The computational time in different cases is listed in Table 2. In Table 2, the high efficiency of the present calculation scheme with the application of PEM is demonstrated. In each case, it is evident that the computational time for matrix multiplication in the conventional method is much more than that spent in the present method, especially in Case (d) that 12.8 h have to be spent on matrix multiplication in the conventional method, compared with 41 s in the present method. In addition, we can see that the CPU time for matrix multiplication in the conventional method is several times more than that spent on LU factorization in all cases. The reason for the above results is that the computations for LU factorization are O(N3 /3), and the computations for matrix multiplication are O(Nf2 ) in the present method but O(2Nf3 ) in the conventional method. The present VLFS model is far smaller and simpler than that used in actual engineering, and this is the reason why the computational time spent on LU factorization is less in these four cases. However,
when the number of unknowns of VLFS is huge, as 284,929 nodes used in Ref. [22], it is plausible that the present method will spend much less time than the conventional method in the stochastic analysis of VLFS.
5. Conclusions In this study, a new calculation scheme for the stochastic hydroelastic analysis of VLFS subjected to short-crested irregular waves is developed with the application of PEM. The VLFS is modeled as a Mindlin plate, and HOBEM combined with FEM are adopted to decouple the fluid–structure interaction in frequency domain. A directional wave spectrum is used for the short-term description of short-crested irregular waves. The present calculation scheme is verified by three examples, and the results demonstrate its accuracy and high efficiency in dealing with the linear stochastic hydroelastic analysis of VLFS in irregular waves. It is noted that the computations for creating response spectrum matrices using the conventional stationary random vibration theory can occupy more than half of the total computations for the stochastic analysis of VLFS, which can be cut by an order of magnitude by the present calculation scheme, and even two orders of magnitude if the auto-spectrum of the response is required only. Therefore, we conclude that more than half of the total computational time spent on the stochastic analysis of VLFS in the conventional way can be reduced by the present method. During the process of stochastic analyses, the response induced by independent monochromatic waves of many discrete frequencies has to be calculated separately to obtain accurate spectral curves and corresponding variances, so it is proposed to give a fast assessment on the hydroelastic response of VLFS subjected to irregular waves using the present calculation scheme. The wave spectrum at each point of VLFS is uniform in this study. As the length of VLFS may be thousands of meters long, different parts of the structure may be subjected to various wave spectra at the same time. Therefore, multiple wave spectra acting
H.-t. Li et al. / Applied Ocean Research 48 (2014) 202–213
on the floating body which takes inhomogeneous sea states into consideration can be considered for further study. Acknowledgements This research was supported by the National Natural Science Foundation of China (Grant No. 51221961 and 51309051), the National Key Basic Research Special Foundation of China (Grant No. 2013CB036101 and 2010CB832704), and High-Tech Ship Research Projects sponsored by Ministry of Industry and Information Technology of the People’s Republic of China. References [1] Watanabe E, Utsunomiya T, Wang CM. Hydroelastic analysis of pontoon-type VLFS: a literature survey. Eng Struct 2004;26(2):245–56. [2] Wu C, Utsunomiya T, Watanabe E. Application of Galerkin’s method in wave response analysis of flexible floating plates. In: Proceedings of the 6th international offshore and polar engineering conference. 1996. p. 307–14. [3] Kashiwagi M. A B-spline Galerkin scheme for calculating the hydroelastic response of a very large floating structure in waves. J Mar Sci Technol 1998;3(1):37–49. [4] Wu YS. Hydroelasticity of floating bodies. UK: Brunel University; 1984 [Ph.D. dissertation]. [5] Newman JN. Wave effects on deformable bodies. Appl Ocean Res 1994;16(1):47–59. [6] Hamamoto T, Fujita K. Wet-mode superposition for evaluating the hydroelastic response of floating structures with arbitrary shape. In: Proceedings of the 12th international offshore and polar engineering conference. 2002. p. 290–7. [7] Fu SX, Moan T, Chen XJ, Cui WC. Hydroelastic analysis of flexible floating interconnected structures. Ocean Eng 2007;34(11–12):1516–31. [8] Gao RP, Tay ZY, Wang CM, Koh CG. Hydroelastic response of very large floating structure with a flexible line connection. Ocean Eng 2011;38(17–18):1957–66. [9] Meylan MH. A variational equation for the wave forcing of floating thin plates. Appl Ocean Res 2001;23(4):195–206. [10] Bai KJ, Yoo BS, Kim JW. A localized finite-element analysis of a floating runway in a harbor. Mar Struct 2001;14(1–2):89–102. [11] Kyoung JH, Hong SY, Kim BW, Cho SK. Hydroelastic response of a very large floating structure over a variable bottom topography. Ocean Eng 2005;32(17–18):2040–52. [12] Seto H, Ohta M, Ochi M, Kawakado S. Integrated hydrodynamic-structural analysis of very large floating structures (VLFS). Mar Struct 2005;18(2):181–200. [13] Ohkusu M, Namba Y. Hydroelastic analysis of a large floating structure. J Fluids Struct 2004;19(4):543–55. [14] Wu C, Watanabe E, Utsunomiya T. An eigenfunction expansion-matching method for analyzing the wave-induced responses of an elastic floating plate. Appl Ocean Res 1995;17(5):301–10. [15] Kim JW, Ertekin RC. An eigenfunction-expansion method for predicting hydroelastic behavior of a shallow-draft VLFS. In: Proceedings of the 2nd international conference on hydroelasticity in marine technology. 1998. p. 47–59. [16] Ohmatsu S. Numerical calculation method for the hydroelastic response of a pontoon-type very large floating structure close to a breakwater. J Mar Sci Technol 2000;5(4):147–60. [17] Riggs HR, Suzuki H, Ertekin RC, Kim JW, Iijima K. Comparison of hydroelastic computer codes based on the ISSC VLFS benchmark. Ocean Eng 2008;35(7):589–97. [18] Wu YS, Wang DY, Riggs HR, Ertekin RC. Composite singularity distribution method with application to hydroelasticity. Mar Struct 1993;6(2–3):143–63.
213
[19] Wang SQ, Ertekin RC, Riggs HR. Computationally efficient techniques in the hydroelasticity analysis of very large floating structures. Comput Struct 1997;62(4):603–10. [20] Kring D, Korsmeyer T, Singer J, White J. Analyzing mobile offshore bases using accelerated boundary-element methods. Mar Struct 2000;13(4–5):301–13. [21] Newman JN, Lee CH. Boundary-element methods in offshore structure analysis. J Offshore Mech Arct Eng 2002;124(2):81–9. [22] Utsunomiya T, Watanabe E, Nishimura N. Fast multipole algorithm for wave diffraction/radiation problems and its application to VLFS in variable water depth and topography. In: Proceedings of OMAE’01, 20th international conference on offshore mechanics and arctic engineering. 2001. OMAE2001/OSU-5202. [23] Utsunomiya T, Watanabe E. Accelerated higher order boundary element method for wave diffraction/radiation problems and its applications. In: Proceedings of the 12th international offshore and polar engineering conference. 2002. p. 305–12. [24] Hamamoto T. Stochastic fluid–structure interaction of large circular floating islands during wind waves and seaquakes. Probab Eng Mech 1995;10(4):209–24. [25] Chen XJ, Jensen JJ, Cui WC, Tang XF. Hydroelastic analysis of a very large floating plate with large deflections in stochastic seaway. Mar Struct 2004;17(6):435–54. [26] Chen XJ, Moan T, Fu SX, Cui WC. Second-order hydroelastic analysis of a floating plate in multidirectional irregular waves. Int J Nonlinear Mech 2006;41(10):1206–18. [27] Ertekin RC, Kim JW, Xia DW. Hydroelastic response of a mat-type, floating runway near a breakwater in irregular seas. In: OCEANS ‘99 MTS/IEEE. Riding the Crest into the 21st Century. 1999. p. 839–47, http://dx.doi.org/10.1109/OCEANS.1999.804985. [28] Miyajima S, Seto H, Ohta M. Hydroelastic responses of the Mega-Float Phase-II model in waves. In: Proceedings of the 12th international offshore and polar engineering conference. 2002. p. 298–304. [29] Papaioannou I, Gao RP, Rank E, Wang CM. Stochastic hydroelastic analysis of pontoon-type very large floating structures considering directional wave spectrum. Probab Eng Mech 2013;33:26–37. [30] Lin JH. A fast CQC algorithm of psd matrices for random seismic responses. Comput Struct 1992;44(3):683–7, http://dx.doi.org/10.1016/ 0045-7949(92)90401-K. [31] Lin JH, Zhang YH. Seismic random vibration of long-span structures. In: de Silva C, editor. Vibration and shock handbook. Florida: CRC Press; 2005 [chapter 30]. [32] Lin JH, Zhang YH, Zhao Y. Pseudo excitation method and some recent developments. Procedia Eng 2011;14:2453–8. [33] Sun L, Liu CF, Zong Z, Dong XL. Fatigue damage analysis of the deepwater riser from VIV using pseudo-excitation method. Mar Struct 2014;37:86–110. [34] Petyt M. Introduction to finite element vibration analysis. Cambridge: Cambridge University Press; 1998. [35] Wehausen JV, Laitone EV. Surface waves. Encyclopedia of physics, vol. 9. Berlin: Springer-Verlag; 1960. p. 446–778. [36] Li HB, Han GM, Mang HA. A new method for evaluating singular integrals in stress analysis of solids by the direct boundary element method. Int J Numer Methods Eng 1985;21(11):2071–98. [37] Pierson WJ, Tuttell JJ, Woolley JA. The theory of the refraction of a short crested Gaussian sea surface with application to the Northern New Jersey Coast. In: Proceedings of 3rd conference on coastal engineering. 1952. p. 86–108. [38] Mitsuyasu H. On the growth of spectrum of wind generated waves. 2: Spectral shapes of wind waves at final fetch. In: Proceedings of 17th Japanese conference on coastal engineering. 1970. p. 1–7 [in Japanese]. [39] Blackford LS, Demmel J, Dongarra J, Duff I, Hammarling S, Henry G, et al. An updated set of basic linear algebra subprograms (BLAS). ACM Trans Math Soft 2002;28(2):135–51. [40] Intel Corporation. Intel math kernel library reference manual 10.1; 2008. [41] Yago K, Endo H. On the hydroelastic response of box-shaped floating structure with shallow draft (tank test with large scale model). J Soc Naval Arch Jpn 1996;180:341–52 [in Japanese].