- Email: [email protected]
Fully nonlinear simulation of resonant motion of liquid confined between floating structures
Computers & Fluids 44 (2011) 89–101
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Fully nonlinear simulation of resonant motion of liquid confined between floating structures C.Z. Wang a, G.X. Wu b,⇑, B.C. Khoo a a b
Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore Department of Mechanical Engineering, University College London, Torrington Place, London, WC1E 7JE, UK
a r t i c l e
i n f o
Article history: Received 10 May 2010 Received in revised form 9 November 2010 Accepted 15 December 2010 Available online 30 December 2010 Keywords: Resonance Floating bodies Potential theory with fully nonlinear boundary conditions High order finite element method
a b s t r a c t A finite element based numerical method is employed to analyze the resonant oscillations of the liquid confined within multiple floating bodies based on fully nonlinear wave theory. The velocity potentials at each time step are obtained through the finite element method (FEM) with quadratic shape functions. The matrix equation of the FEM is solved through an iteration. The waves at the open boundary are absorbed by the method of combination of the damping zone method and the Sommerfeld–Orlanski equation. Numerical examples are given by floating two rectangular cylinders, two wedge-shaped cylinders and two semi elliptic cylinders undergoing forced oscillations at the resonant frequencies. The numerical results are compared with the first and second order solutions by previous study, which showed that the first order resonance happens at the odd order natural frequencies x2i1 ði ¼ 1; 2; Þ and second order resonance at the half of even order x2i =2ði ¼ 1; 2; Þ for antisymmetric motions and happen at x2i ði ¼ 1; 2; Þ and x2i =2ði ¼ 1; 2; Þ for symmetric motions, and it is found that they are in good agreement in smaller amplitude motions within a sufficiently long period of time. However, difference begins to appear as time increases further. This happens even in smaller amplitude motions. In all the calculated cases, the results always become periodic in time eventually. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Resonance is one of the important features in nature. For a mechanical system in periodic motion, it occurs at a frequency when the inertial term is cancelled by the stiffness term. As a result its response to external excitation at such a frequency can become exceedingly large. One of the well known examples in fluid mechanics is sloshing in a tank [1–3]. Another typical example is the motion of the liquid confined between two floating bodies, as considered by Wang and Wu [4] for two dimensional rectangular cylinders, by Sun et al. [5] for two barges, by Malenica et al. [6], and Chen and Fang [7] for water confined by two ships with varying gap between their near sides. When resonance does occur, it has some serious implications to the safety and operation of a mechanic system. Thus to accurately predict resonant behaviour is hugely important. The analysis of the resonant behaviour is usually made based on the perturbation theory, starting with the linear theory. The linear theory usually captures the values of natural frequencies and the overall behaviour of the system at resonance quite well. It shows that when the excitation frequency x is equal to one of the natural ⇑ Corresponding author. E-mail address: [email protected] (G.X. Wu). 0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2010.12.014
frequencies xi ; i ¼ 1; 2; 3; :::, resonance will occur. Though subharmonics, one would then expect that resonance in the second order result would occur when x ¼ xi =2; i ¼ 1; 2; 3; :::. In fact, when the excitation has frequency components xej ; j ¼ 1; 2; 3; :::, resonance in the second order result could occur when xej1 þ xej2 ¼ xi . These have been confirmed in the analyses by Wu [3] for a sloshing tank and by Wang and Wu [4] for the liquid motion between two bodies floating on the water surface. While the perturbation theory can capture some important features of resonance, it does not always provide accurate results quantitatively. The main reason for this is that the perturbation theory is based on small parameter expansion. At resonance, the motion becomes very large. In fact, the motion can become infinite based on the linear theory when there is no damping in the system. This clearly contradicts with the foundation of the perturbation theory. Thus an appropriate approach at resonance is to use the fully nonlinear theory. The nonlinear theory has been used previously for the resonant behaviour of a sloshing tank [8,2]. Very little has, however, been done using the fully nonlinear theory for the motion of the liquid confined between two floating bodies. This has in fact important implications to the marine structures, such as a catamaran or ships very close to each during operation [5]. The present paper uses a high order finite element method to analyze the motion of the liquid confined between two floating
90
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
cylinders. While the method can be used for general shapes, cases considered include twin rectangles, twin wedges and twin semi ellipses. Extensive simulations are made at the first and the second order resonant frequencies. The agreement and disagreement between the perturbation theory and the fully nonlinear theory are discussed.
@/ 1 p ¼ q þ jr/j2 þ gy @t 2
ð8Þ
where q is the fluid density. The hydrodynamic force acting on the body can be expressed as *
F¼
Z
*
p n ds
ð9Þ
S
2. Mathematical formulation
~¼ M
Zb
*
pð~ r n Þds
ð10Þ
Sb
2.1. Governing equation and boundary conditions
*
We consider the problem of motion of the liquid confined between two floating bodies, as shown in Fig. 1. A Cartesian coordinate system oxy is defined in which the x-axis coincides with the undisturbed free surface, and the y-axis is positive upward. The free surface and the floating body surface are denoted as Sf and Sb respectively. The bottom of the liquid is assumed to be a horizontal plane at y = h and is denoted by Sbot . The fluid is assumed incompressible and inviscid, and the flow is irrotational. The liquid motion can then be described by a velocity potential /. It satisfies the Laplace equation in the fluid domain 8, or
r2 / ¼ 0
in 8
ð1Þ
The Eulerian form of the kinematic and dynamic boundary conditions on the free surface y ¼ g can be written as
@/ @ g @/ @ g ¼ 0 on Sf @y @t @x @x @/ 1 þ g g þ jr/j2 ¼ 0 on Sf @t 2
ð2Þ ð3Þ
where g is the acceleration due to gravity. Eqs. (2) and (3) can also be written as in the Lagrangian form
d/ 1 ¼ g g þ r/r/ dt 2 dx @/ dy @/ ¼ ; ¼ dt @x dt @y
ð4Þ ð5Þ
on Sb ; Sc and Sbot
ð6Þ
where V n is the velocity of the surface along the inward normal vector ~ n. In addition, the initial condition can be expressed as
/ðx; y ¼ n; t ¼ 0Þ ¼ /0 ðxÞ;
gðx; t ¼ 0Þ ¼ nðxÞ
r2 /t ¼ 0
ð7Þ
2.2. Evaluation of hydrodynamic forces Once the velocity potential has been found, the pressure in the fluid can be determined by the Bernoulli equation
ð11Þ
On the fixed boundary it satisfies
@/t ¼0 @n
ð12Þ
On the free surface /t is given by the Bernoulli equation as
1 /t ¼ gy r/r/ 2
On the rigid surface, the boundary condition is
@/ ¼ Vn @n
where r ¼ x i þy~ j is the position vector from the origin of the system to the point on the body surface. In Eqs. (9) and (10), the computation of integration of 1 jr/j2 þ gy over the body surface does not present real diffi2 culty. There are however some difficulties associated with ou/ot. This is the partial temporal derivative which is taken when the spatial point is fixed. This can be problematic when the body is in motion. It may be calculated by following a point on the body surface and using the finite difference method. However, it is often found that the force obtained this method has sawteeth behaviour, especially in the case of large amplitude motion. In order to avoid this, we treat /t as another harmonic function. The procedure can be summarized below. Assume the translational velocity of the body at the origin is ~ V and the angular velocity about the axis passing through the origin ~. In the fluid domain, the time derivative / satisfies the Laplace is X t equation
ð13Þ
On the moving boundary, the condition can be written as (Wu [9])
h i h i @ _ ~_ ~ ~ ~ @ r/ ~ @ ~ ~ ð/t Þ ¼ ~ þX r ðV r/Þ V þX r n V: @n @n @n
ð14Þ
~ indicates the derivative with respect to where the dot over ~ V and X time. If the acceleration is known, as in the forced motion for example, the above problem can then be solved to obtain /t . For a freely floating body, the acceleration will depend on the force which in turn depends on the acceleration. To overcome this mutual dependence, an auxiliary function can be introduced [9]. The acceleration can then be found before the pressure is known.
3. Finite element discretization and mesh generation We shall adopt the finite element method to solve the velocity potential problem at each time step. Quadrilateral parametric element (see Fig. 2) is used. The shape functions defined in a local coordinate system ~ n ¼ ðn; gÞ corresponding to element e with eight nodes may be expressed as
9 ðeÞ Ni ðn; gÞ ¼ 14 ð1 þ ni nÞð1 gi gÞð1 þ ni n þ gi gÞ ði ¼ 1; 2; 3; 4Þ > > = ðeÞ Ni ðn; gÞ ¼ 12 ð1 n2 Þð1 þ gi gÞ ði ¼ 5; 7Þ > > ; ðeÞ ði ¼ 6; 8Þ Ni ðn; gÞ ¼ 12 ð1 g2 Þð1 þ ni nÞ ð15Þ Fig. 1. Sketch of the computational domain.
This gives
91
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
The Galerkin method is, on the other had, suitable for general cases but needs to solve another matrix equation at each time step. When higher element is used, the derivatives can be obtained by differentiating the shape function directly. This has in fact been used by Wang and Khoo [16] and Wang and Wu [4]. We write n @/ X @Nj ¼ /j ; @x @x j¼1
n @/ X @Nj ¼ /j @y @y j¼1
ð19Þ
where
2
3
@N j 4 @x @N j @y
Fig. 2. Meshes for (a) two semielliptic cylinders; (b) two wedge-shaped cylinders. 8 X
x¼
ðeÞ
xk N k ;
y¼
k¼1
8 X
ðeÞ
yk N k
k¼1
within the element. The velocity potential / in the fluid domain can be expressed as
/ðx; yÞ ¼
n X
@n @x @g
@y @n @y @g
#1 2 @Nj 3 4 @n 5 @N j @g
ð20Þ
within each element. The second order derivatives @ 2 /=@x2 ,@ 2 /=@x@y and @ 2 /=@y2 can be obtained in a similar way through the Jacobian. The fourth order Runge–Kutta method is adopted for the integration with respect to time to update the wave elevation and the potential on the free surface. Let P denote one of these three variables ðx; y; /Þ. If Pn is the value at time step nDt, the new value at t ¼ ðn þ 1ÞDt can be obtained through the following equation
Pnþ1 ¼ Pn þ
Dt ðk1 þ 3k2 þ 3k3 þ k4 Þ 8
ð21Þ
where
/J N J ðn; gÞ
ð16Þ
J¼1
ki ¼
where /J is the potential at node J and n is the number of nodes. N j ðeÞ in the equation is equal to Nk if node J in the global system is node k of element e, and ðx; yÞ 2 e. Otherwise N j ¼ 0. Through the Galerkin method, we have
Z Z
5¼
" @x
i1 d i1 ði ¼ 1; 2; 3; 4Þ P nDt þ Dt; Q nþ 3 dt 3
and
Q n ¼ Pn ; 1
Q nþ3 ¼ Pn þ
2
r /Ni d8 ¼ 0
ð17Þ
8
using Green’s identity and the boundary conditions, we can obtain the following matrix equations
½Kf/g ¼ fFg
Dt k1 ; 3
2 1 Q nþ3 ¼ Pn þ Dt k1 þ k2 ; 3 Q nþ1 ¼ Pn þ Dtðk1 k2 þ k3 Þ
ð18Þ
where
4.2. Remeshing and smoothing
f/g ¼ ½/1 ; /2 ; . . . ; /I ; . . . ; /n T ðI R Sp Þ; Z Z rN I rNJ d8 ðI R Sp &J R Sp Þ; K IJ ¼
When the simulation is over a substantial period of time, the nodes on the free surface may cluster and cause elements to be distorted. In order to avoid this, nodes on the free surface should be rearranged every several time steps. A B-spline based remeshing technique [14] is used in the present paper. A smoothing technique called the energy method is also adopted to avoid the sawteeth problem on the free surface [14].
FI ¼
Z Sn
8
N I fn dS
Z Z
rNI
8
n X
J¼1 j2Sp
ðfp ÞJ rNJ d8 ðI R Sp Þ
Sp in above equations represents the Dirichlet boundary on which the potential, denoted by fp, are known, and Sn represents the Neumann boundary on which the normal derivative of the potential, denoted by fn, are known. The integration in the equation is calculated element by element with respect to ðn; gÞ through the Jacobian. Once the coefficients are found, the matrix equations are then solved through an iteration based on the conjugate gradient method with a symmetric successive over relaxation (SSOR) preconditioner. 4. Numerical procedures 4.1. Evaluation of derivatives of potential with respect to coordinates There are first and second order derivatives of the potential in Eqs. (4), (5), and (14), respectively. Approaches such as difference method [10–12] and Galerkin method [13–15] have been used associated with linear elements. The difference method is efficient but is only suitable for wall-sided bodies and structured meshes.
4.3. Numerical radiation conditions For long time simulations, an appropriate radiation condition should be imposed on the boundary Sc to minimize the wave reflection. Here, we use a combination of an artificial damping zone and Sommerfeld–Orlanski condition. The damping zone method is also called a sponge layer or an artificial beach. The artificial damping zone is an absorbing boundary. Outgoing waves will be first absorbed through artificial damping terms added to free surface boundary conditions upon their arrival at the damping zone and the most remanding waves will propagate out of the computational domain through the Sommerfeld–Orlanski condition [17]. Cointe et al. [18] adopted the damping zone technique for their 2D numerical wave tank simulations and this method was also used by Tanizawa [19] in fully nonlinear 2D simulations. In the damping zone, the following modified absorbing free surface boundary condition is used
92
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
Table 1 Number of elements and nodes for meshes one and two.
Mesh one Mesh two
NF1
NF2
ND
NB
NH
NE
NN
50 70
32 36
8 12
12 16
18 24
2580 4608
8117 14,337
d/ 1 ¼ g g þ r/r/ mðxÞ/ dt 2 dx @/ dy @/ ¼ ; ¼ mðxÞy dt @x dt @y
ð22Þ ð23Þ
Fig. 5. Histories of waves at the right side of cylinder one at x ¼ x02 and A = 0.05d.
with
(
mðxÞ ¼
xx0 2
ax 0
k
x0 6 x 6 x1 ¼ x0 þ bk 0 < x < x0 or x > x1
where x is the wave frequency and k is the linear wave length. The damping zone starts from point x0 and extends for a length of bk. Parameters a and b control the strength and the length of the damping zone respectively, and they may be obtained through numerical experiments. 5. Numerical results 5.1. Two rectangular cylinders in forced motions
Fig. 3. Wave runups along cylinder one, (a) left hand side, (b) right hand side at x ¼ x02 and A = 0.0125d; —— mesh one; – – – mesh two.
Wang and Wu [4] have undertaken some detailed analysis for the liquid motion confined between twin rectangular cylinders floating on the free surface based on the perturbation analysis. They confirmed that when the motion frequency of the cylinders is at one the natural frequencies very large transient wave motion between cylinders can occur based on the linear theory. They further showed that when the motion frequency of the cylinders is at half of one of the natural frequencies the second order transient motion of the liquid can also be very large. This analysis has raised the question of whether these behaviours are a result of the perturbation theory. This is because the perturbation theory is based on small motion assumption and the results obtained from this theory become inaccurate quantitatively when the motion is big. This question will be discussed in detail by the fully nonlinear theory below.
Fig. 4. Wave runups along cylinder one, (a) left hand side, (b) right hand side at x ¼ x02 and A = 0.0125d; —— DT = t/100; – – –DT = t/200.
Table 2 Number of elements and nodes for meshes three and four.
Mesh three Mesh four
NF1
NF2
ND
NB
NH
NE
NN
60 80
36 42
8 12
12 16
18 24
3048 5232
9573 16,261
Fig. 6. Comparisons of waves at the right side of cylinder one; (a) A = 0.0125d; (b) A = 0.025d; (c) A = 0.05d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
93
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
Fig. 7. Comparisons of forces on cylinder one between the fully nonlinear results and the second order solutions; (a), (b) A = 0.0125d; (c), (d) A = 0.025d; (e), (f) A = 0.05d; —— fully nonlinear; linear; linear + second order. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
We assume that the two cylinders are identical. The cylinder on the left hand side is called cylinder one and the one on the right is called cylinder two. The breadth at mean water level and draught of each cylinder are 2b and d respectively. The centre lines of the two cylinders are located at x = Lcy and Lcy, respectively. In the following simulations xn is defined as
Fig. 8. Comparison of waves at the right side of cylinder one at A = 0.0125d and x ¼ x02 within two hundred periods.
xn ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi np g; 2ðLcy bÞ
n ¼ 1; 2; . . .
ð24Þ
This should be a good estimate for resonant frequencies when the cylinders are rectangular and the confined liquid resembles the sloshing in a tank (see, for example, Wu [3]). In particular, as discussed by Wu [3], when n is an odd number, the liquid motion is anti symmetrical and when n is an even number, the motion is symmetrical. For other cylinders such as twin wedges and twin semi
Fig. 9. Comparison of waves at the right side of cylinder one at x ¼ x02 within two hundred periods.
94
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
ellipses considered below, the resonant frequencies are expected to further differ from Eq. (24). The two cylinder are assumed to undergo the same vertical motion defined below
X ¼ A sin xt
h = 10d, and d = b and Lcy = 8b. We choose A = 0.0125d. Because the problem is symmetric about x = 0, the resonant behaviour may occur only when x is around x2n in Eq. (24). The real resonant
ð25Þ
where x and A are the motion frequency and amplitude, respectively. To avoid starting the motion suddenly with a non zero velocity, the displacement function in this equation could be modified initially and then returns to the normal one after a short period of time. For the results given below, the water depth is taken as
Fig. 13. Comparisons of forces on cylinder one between the fully nonlinear results and the second order solutions at x ¼ 0:5x02 and A = 0.05d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
Fig. 10. Histories of waves and forces for cylinder one at x ¼ x02 and A = 0.2d.
Fig. 14. Comparison of waves between the fully nonlinear result and the second order solution at x ¼ 0:5x02 and A = 0.2d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
Fig. 11. Wave profiles at t/T = 11, 11.05, 11.1, . . . , 12 at x ¼ 0:5x02 and A = 0.2d.
Fig. 12. Comparison of waves between the fully nonlinear result and the second order solution at x ¼ 0:5x02 and A = 0.05d.
Fig. 15. Comparisons of forces between the fully nonlinear results and the second order solutions at x ¼ 0:5x02 and A = 0.2d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
frequency x0n may be different from xn. We denote x0n ¼ C xn , where C is a constant. In fact, through detailed numerical tests, Wang and Wu [4] have found that first order resonance occurs when C = 1.01 or x ¼ x02 ¼ 1:01x2 . The control surfaces Sc at both ends are located at distance about four times wavelength away the nearest cylinder. Two meshes are used to test the numerical convergence. The details are given in Table 1, in which NF1 is the segment number along the free surface on the left of cylinder one or the right of cylinder two, NF2 the segments on the free surface between the two cylinders, ND and NB the segments along the vertical and horizontal faces of one cylinder respectively, NH the segment on the both control surfaces Sc at the far ends, NE and NN are the total number of elements and the total number of nodes in the whole fluid domain respectively. The results from the meshes are given in Fig. 3. It can be seen that they are graphically indistinguishable over the entire simulations of 200 cycles. Further convergence study with respect to time step is made using mesh one. Fig. 4 gives results with Dt ¼ T=100 and Dt = T/200 and the wave histories from the two time steps are in good agreement. These tests show that mesh one and Dt = T/100 can provide convergent results for this case. Further convergence study at a larger amplitude of A = 0.05d is also made. The details of the meshes are given in Table 2 and the result
95
Fig. 16. History of the second order wave at the right side of cylinder one at x ¼ 0:5x02 .
is given in Fig. 5. It can be seen that mesh three and Dt = T/200 give the convergent results. We now undertake detailed analysis of the nonlinear behaviour of the liquid motion at the natural frequency x = 1.01x2 by comparing the results with those obtained from the perturbation theory [4]. Three amplitudes A = 0.0125d, 0.025d and 0.05d have been chosen and meshes one with Dt = T/100, three with Dt = T/ 150 and three with Dt = T/200 have bee used for these amplitudes, respectively. The wave runups along the right side of cylinder one are given in Fig. 6. It can be seen that the linear and the combined linear and second order results are in very good agreement with those from the fully nonlinear theory at A = 0.0125d. When the amplitude is increased to A = 0.025d, there is some visible discrepancy between the results from the linear theory and the fully nonlinear theory, and the combined linear the second order theory give a better agreement. When the amplitude is increased further to A = 0.05d, the differences between the results from the perturbation theory and the fully nonlinear theory become obvious. The differences appear not only in the magnitudes of the peaks and troughs, but also in the phase. Comparisons are also made for the hydrodynamic forces without the contribution of initial buoyancy. Results from the perturbation theory and the fully nonlinear theory are given in Fig. 7. As in the case of wave elevation, the results are in good agreement at smaller motion amplitude and discrepancy appears when it increases. It is also observed that there are double peaks in the vertical force within one cycle at larger motion amplitude, which implies the high harmonic nonlinear effects. Figs. 6 and 7 show that when A increases the nondimensionalized nonlinear results are much smaller than those from the perturbation theory. In fact, the latter is expected to increase further due to the near resonance behaviour. At smaller A, however, the results are in good agreement over the entire period of simulation. It is important to realize, however, that this should be put in the context that the smaller the A is, the longer period will be over which the results from the perturbation theory and the fully nonlinear theory will be in good agreement. It is expected that as time progresses the discrepancy will develop even at small A. Thus we run the simulations over two hundred periods at A = 0.0125d. Fig. 8 shows that after about 40 periods, significant difference appears between the results from the perturbation theory and the fully nonlinear theory. The envelope of the wave amplitude from the nonlinear theory begins to oscillate at a lower frequency. As the confined liquid is linked to the open domain, some radiation
Fig. 17. Histories of waves and forces at x ¼ 0:5x02 and A = 0.2d.
Fig. 18. Histories of waves at the right side of cylinder one at (a) x ¼ 0:5x02 and (b) x ¼ 0:5x02 in sway in same direction (A = 0.05d).
96
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
damping is expected. Thus it is expected that the results from the perturbation theory will not increase indefinitely and an envelope of the amplitude oscillation will develop eventually, although the frequency will be much lower. Results from the nonlinear theory at A = 0.0125d, 0.025d and 0.05d are given in Fig. 9 for comparison. It can be seen that the nondimensionalized amplitude of the envelope decreases when A increases, while its frequency increases. Further simulation is made at even larger A = 0.2d. Fig. 10 gives the wave elevation and force. The result shows that after a very relatively short transition period, they become more or less periodic. All these results show that at a small A, the fully nonlinear results exhibit the resonant behaviour of the linear theory over a long period of time, but the results will eventually depart form those the perturbation theory when time progresses further. As large A, the resonant behaviour is far less obvious due to the highly nonlinear effects. Fig. 11 gives the development of the wave profile between t = 11T and t = 12T with an interval of 0.05T. It is seen that the wave between the two cylinders is much larger than those outside the two cylinders.
Following the work of Wu [3], Wang and Wu [4] also found that at x ¼ 0:5x02 , the linear theory does not show any resonance but the second order result does. We choose A = 0.05d and run the fully nonlinear simulation at x ¼ 0:5x02 . The results are shown in Figs. 12 and 13. The comparison shows that the results are in excellent agreement within those obtained from the combined linear and the second order theory for whole 40 cycles in the figure. This gives more evidence that the second order resonance does reflect the physical reality to large extent in some cases. Simulation is then made at A = 0.2d and the results are given in Figs. 14 and 15. The difference between results from the fully nonlinear theory and the perturbation theory becomes evident. This is obviously due to the same reason as that discussed right above. We run the simulation over a longer time of two hundred periods. Fig. 16 clearly shows the resonant behaviour of the second order wave elevation at x ¼ 0:5x02 (k is the linear wave number in the figure), which has been discussed in detail by Wang and Wu [4]. Fig. 17 gives the wave elevation at the right side of cylinder one and force from the nonlinear theory within 200 periods. The results become periodic quickly and there is no obvious resonant behaviour due to the nonlinear effects. We now consider some cases of cylinders in swaying motions. Two cases are simulated. Case one is the two cylinders moving in the same direction along x-axis with X ¼ A sin xt. Case two is the
Fig. 19. Histories of waves and forces for cylinder one at x ¼ 0:5x02 in sway in same direction; – – – A = 0.05d; —— A = 0.1d; . . . . . . A = 0.2d. Fig. 21. Numerical tests to determine the constant C.
Fig. 20. Histories of waves and forces for cylinder one at x ¼ x02 in sway in opposite directions; —— A = 0.0125d; A = 0.025d; A = 0.05d. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 22. Histories of waves at (a) x = Lcy 3b and (b) x = Lcy + 3b at x ¼ x02 for twin wedges in heave;—— A = 0.0125d; A = 0.025d; A = 0.05d.
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
two cylinders moving in opposite directions along the x-axis, which means that the motion for cylinder one is given as X ¼ A sin xt and X ¼ A sin xt for cylinder two. Simulations are made at the first and half of second resonant frequencies, or x01 and 0:5x02 (x01 1:068x1 and x02 1:01x2 ) in Case one. Fig. 18 shows the waves at the right side of cylinder one at A = 0.05d. Waves corresponding to both frequencies reach periodic state relatively early within two hundred periods. The nondimensionalized maximum value of wave at the right side is about 20.6 at around t/ T = 22.15. The peak is about 15.3 when the wave becomes periodic. The wave at the left side is, however, about 2.8 only. The wave at half of the second resonant frequency 0:5x02 has a more evident nonlinear effect, which may also be seen by comparisons of waves at the right side of cylinder one and forces on cylinder one at A = 0.05d, 0.1d and 0.2d in Fig. 19. Case two is similar to the vertical motion. The flow is symmetric about x = 0, and the first- and second order resonances happen at x ¼ x02 and x ¼ 0:5x02 respectively. We give results at x ¼ x02 only. Fig. 20 shows waves and forces at A = 0.0125d, 0.025d and 0.05d. The overall behaviour of the results is very similar to that in the heave motion. 5.2. Two wedge-shaped cylinders in forced motions A significant difference between this case and the rectangular cylinders is that the distance between the near sides of the two cylinders varies along the water depth. During the motion the distance between the two cylinders at the wave surface level also
97
varies. We notice that this distance is an important parameter for the natural frequency. Thus when its variation is included in the nonlinear model, the natural frequency changes with time, which makes the resonant effect less profound. On the other hand, when the body is non-wall-sided, it will have more significant nonlinear effect on the force. We consider a case in which the deadrise angle of each of the two identical wedges is 600 and Lcy = 16b (2b is the breath of the wedge at the calm water level when the body is stationary). The best approximation of the natural frequency is chosen by numerical tests based on the linear theory. Fig. 21 shows that the amplitudes of wave and forces on cylinder one are the largest when C = 1.06, which implies the resonance may occur at x ¼ x02 1:06x2 . We first consider three cases with amplitudes A = 0.0125d, 0.025d and 0.05d at x = 1.06x2 in the same vertical motion as the previous rectangular cylinders. On the free surface, there are 60 segments on the left side of cylinder one and the right side of cylinders two, 30 segments between the two cylinders and 16 segments on each cylinder surface and 18 segments on the control surface Sc, respectively, which corresponds to a total of 2700 elements and 8469 nodes. The time intervals for A = 0.0125d, 0.025d and 0.05d are T/100, T/150 and T/200, respectively. The histories of waves at x ¼ Łcy 3b are given in Fig. 22. It is noticed that all wave amplitudes increase at the earlier stage of the simulations and then reaches periodic state more quickly than the rectangular cases shown in Fig. 9. As mentioned above, one should note that
Fig. 23. Comparisons of waves between the fully nonlinear results and the second order solutions for twin wedges in heave; (a), (b) A = 0.0125d; (c), (d) A = 0.025d; (e), (f) A = 0.05d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
98
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
Fig. 24. Comparisons of forces on cylinder one between the fully nonlinear results and the second order solutions for twin wedges in heave; (a), (b) A = 0.0125d; (c), (d) A = 0.025d; (e), (f) A = 0.05d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
the wedges are non-wall-sided bodies. During the vertical motion, the distance between the near sides of the two wedges changes. Because of this variation of the configuration of the fluid domain confined between the two bodies, the natural frequency obtained
approximately based on the mean domain through the linear theory would have less significant effect. However, this effect is still not negligible. Fig. 22 shows that the wave peaks at the right side x ¼ Łcy þ 3b (confined between the wedges) is about four times
Fig. 25. Histories of waves at x = Lcy + 3b and x ¼ 0:5x02 for twin wedges in heave; (a) A = 0.05d; (b) A = 0.1d; (b) A = 0.2d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
Fig. 26. Histories of waves at x = Lcy + 3b and forces on cylinder one at x ¼ 0:5x02 for twin wedges in heave.
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
99
that at the left side x ¼ Łcy 3b (open water) at A = 0.0125d and nearly six times at A = 0.05d. Comparisons are made between the fully nonlinear results and the linear and linear plus second order solutions at A = 0.0125d, 0.025d and 0.05d are shown in Figs. 23 and 24. The waves correspond to two locations: one at x ¼ Łcy þ 3b and the other at x = 0, which is the middle point between the two cylinders. For the case of A = 0.0125d, all the results are generally in good agreement even though the fully nonlinear result is slightly smaller than the linear and linear plus second order. For the largest A = 0.05d, the differences between them increase. Unlike the rectangular cases, the linear and linear plus second order solutions in these two cases at this first resonant frequency do not always increase with time. They become periodic, or their peaks or troughs become constants within the period of simulation. This suggests that the resonance may have less effect when a body has flare. It is interesting to notice that the above water flare is not included in the fully linear theory. Thus it is quite likely that the linear theory could overestimate the resonant effect for the flare bodies. Fig. 27. Histories of waves at x = Lcy + 3b and x ¼ x02 for twin wedges in sway in same direction; (a) A = 0.05d; (b) A = 0.1d; (c) A = 0.2d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
Fig. 28. Histories of waves at x = Lcy + 3b and forces on cylinder one at x ¼ x02 for twin wedges in heave; – – – A = 0.05d; —— A = 0.1d; . . . . . . A = 0.2d.
Fig. 29. Histories of waves at x = Lcy + 3b and forces on cylinder one at x ¼ 0:5x02 for twin wedges in heave; – – – A = 0.1d; —— A = 0.2d; . . . . . . A = 0.4d.
Fig. 30. Histories of waves at (a) x = Lcy 3b (b) x = Lcy + 3b and (c) x = 0 at x ¼ x02 for twin wedges in sway in opposite directions; – – – A = 0.0125d; —— A = 0.025d; . . . . . . A = 0.05d.
Fig. 31. Histories of forces on cylinder one at x ¼ x02 for twin wedges in sway in opposite directions; – – – A = 0.0125d; —— A = 0.025d; . . . . . . A = 0.05d.
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
sults from the fully nonlinear theory at A = 0.05d, 0.1d and 0.2d are given in Fig. 28. These figures show that the nonlinear effect does not seem to be as strong as the heave motion. Results at 0:5x02 are given in Fig. 29. The nonlinear effect becomes more significant. In
4 2
h/A
Further simulations at x ¼ 0:5x02 are made. The results of waves at the right side of cylinder one at A = 0.05d, A = 0.1d and A = 0.2d together with the linear and linear plus second order solutions are given in Fig. 25, long after they have reached the periodic state. The comparisons with the linear plus second order show a good agreement at A = 0.05d. Discrepancy is evident at A = 0.1d and is even more at A = 0.2d due to stronger nonlinearity. Further comparisons between waves and forces of the fully nonlinear results at A = 0.1d, A = 0.2d and A = 0.4d are given in Fig. 26 and the nonlinear effect is strong. Simulations in horizontal motions are also made at the first and half of second resonant frequency x01 ðx01 1:16x1 , which is obtained in a similar way as in Fig. 21) and 0:5x02 (x02 1:06x2 ) when the wedges move in the same direction. Fig. 27 shows the wave histories at x = Lcy + 3b and x = 0 with x ¼ x01 and at three amplitudes A = 0.05d, 0.1d and 0.2d. It is seen that the differences between the fully nonlinear results and the linear and linear plus solutions are small even for a larger amplitude A = 0.2d. Further re-
0 -2 -4 0
5
10
15
20
25
30
35
190
195
200
t/T 1
Fx/rgdA
100
0 -1 0
5
10
15
20
25
30
35
190
192
194
196
198
200
192
194
196
198
200
Fy/rgdA
t/T 1 0 -1 0
5
10
15
20
25
30
35
190
t/T Fig. 34. Histories of waves at x = Lcy + 1.5b and forces cylinder one at x ¼ 0:5x02 for twin semi elliptic cylinders in heave; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
Fig. 32. Histories of waves and forces on cylinder one at x ¼ 0:5x02 for twin wedges in sway in opposite directions; – – – A = 0.1d; —— A = 0.2d; . . . . . . A = 0.4d.
Fig. 35. Histories of second order wave at x = Lcy + 1.5b and x ¼ 0:5x02 for twin semielliptic cylinders in heave.
Fig. 33. Histories of waves at x = Lcy + 1.5b and forces on cylinder one at x ¼ x02 for twin semielliptic cylinders in heave; – – – A = 0.005d; —— A = 0.01d; . . . . . . A = 0.02d.
Fig. 36. Comparison of linear and second order solutions at x ¼ 0:5x02 and A = 0.05d for twin semielliptic cylinders in heave.
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
particular, the time history curves have departed further away from the sinusoidal shape. Simulations of horizontal motions when two wedges move in the opposite directions at the second and the half of second frequencies x02 and 0:5x02 are also made. The result for x ¼ x02 and x ¼ 0:5x02 are given in Figs. 30–32, respectively. This is a symmetrical case and behaviour of the results is therefore similar to that in the heave motion. 5.3. Two semi elliptic cylinders in forced motions The previous examples are for bodies with no curvature. Here, we consider the case of two semi elliptic cylinders. b = 0.75 is the semi minor axis of each cylinder in the x direction and d = 1/b is the semi major axis in the y direction. The horizontal distance between the centre lines of the two cylinders is Lcy = 8b. Simulations are made for vertical motion at x ¼ x02 and x ¼ 0:5x02 respectively 0 x2 1:01x2 . The wave elevation at xc ¼ Lcy þ 1:5b, and forces on cylinder one with x ¼ x02 at three amplitudes A = 0.005d, 0.01d & 0.02d are given in Fig. 33. It should be noted although the body has curvature, it is wall sided at the mean water line, i.e., the tangential of the body surface at its intersection with the free surface is in the vertical direction. The variation of the distance between the near sides of the cylinders at the water line is smaller than that of the wedges. Thus the effect of the natural frequency is less significant than that in the case of two wedges. This is reflected in Fig. 33 by the longer period over which the amplitude continues to increase and the large amplitude when the motion becomes periodic. The figure also shows that at the periodic state, the nonlinear effects are clearly very important. The results at the half of second resonant frequency are given in Fig. 34. In the figure, the linear and linear plus second order solutions are also given for comparison. It is seen that results from these methods are in good agreement for both the wave elevation and the force during the first thirty periods and difference between them begins to develop. It can be observed that the nonlinear effects on both the wave elevation and the force are strong. Based on the second order theory, the second order solution in this case may increase to become very large, as discussed in Wu [3]. This however will not be infinite, as discussed in Wang and Wu [4] due to the radiation damping. This is reflected in the history of the second order wave at x ¼ Lcy þ 1:5b in Fig. 35. Compared with the case of two rectangles in Fig. 16, the increase of the wave elevation stops much earlier. It is also interesting to see in Fig. 36 that the second order wave becomes larger than the linear one after t/ T = 14.8. This implies that full nonlinear theory may be needed. 6. Conclusions A higher order finite element method has been used to analyze the fully nonlinear motion of the liquid combined between floating twin-cylinders in forced oscillations. The velocity potential at each time step is obtained through solving matrix equations based on an iteration method. The fourth order Runge–Kutta method is used to update the wave elevation and potential on the free surface at each
101
time step. The radiation condition is imposed through a combination of the damping zone method and the Sommerfeld–Orlanski equation. Extensive simulations have been made for twin rectangular, wedge-shaped and semi elliptic cylinders in vertical and horizontal motions. It confirms the first order and the second resonant behaviour at small body motion amplitude, predicted by the linear and the second order theory. The wave elevation amplitude increases over a long period of time in the perturbation theory but the increase will stop much earlier in the fully nonlinear theory. This becomes more evident when the body motion amplitude increases. It is also found that the resonant effect is very significant for the rectangular cylinders. When the body has flare or curvature at the water line, the resonant effect is less significant even though it may be still important. This may suggest that the linear theory could over predict the resonant effect for bodies with flare or curvature. References [1] Faltinsen OM. A numerical non-linear method of sloshing in tanks with twodimensional flow. J Ship Res 1978;22:193–202. [2] Faltinsen OM, Timokha AN. Adaptive multimodal approach to nonlinear sloshing in a rectangular tank. J Fluid Mech 2001;432:167–200. [3] Wu GX. Second order resonance of sloshing in a tank. Ocean Eng 2007;34:2345–9. [4] Wang CZ, Wu GX. Analysis of second order resonance in wave interactions with floating bodies through a finite element method. Ocean Eng 2008;35:717–26. [5] Sun L, Taylor PH, Eatock Taylor R. First and second order wave effects in narrow gaps between moored vessels marine operations specialty symposium, Singapore; 2008. p. 113–24. [6] Malenica S, Orozco JM, Chen XB. Some aspects of multibody interactions in seakeeping. In: Proceedings of 15th international offshore and polar engineering conference, Seoul, Korea; 2005. [7] Chen GR, Fang MC. Hydrodynamic interactions between two ships advancing in waves. Ocean Eng 2001;28:1053–78. [8] Wu GX, Ma QW, Eatock Taylor R. Numerical simulation of sloshing waves in a 3D tank based on a finite element method. Appl Ocean Res 1998;20:337–55. [9] Wu GX, Eatock Taylor R. The coupled finite element and boundary element analysis of nonlinear interactions between waves and bodies. Ocean Eng 2003;30:387–400. [10] Ma QW, Wu GX, Eatock Taylor R. Finite element simulation of fully nonlinear interaction between vertical cylinders and steep waves. Part 1: methodology and numerical procedure. Int J Num Methods Fluids 2001;36:265–85. [11] Ma QW, Wu GX, Eatock Taylor R. Finite element simulation of fully nonlinear interaction between vertical cylinders and steep waves. Part 2: numerical results and validation. Int J Num Methods Fluids 2001;36:287–308. [12] Wang CZ, Wu GX. Time domain analysis of second order wave diffraction by an array of vertical cylinders. J Fluids Struct 2007;23:605–31. [13] Wu GX, Eatock Taylor R. Finite element analysis of two dimensional nonlinear transient water waves. Appl Ocean Res 1994;16:363–72. [14] Wang CZ, Wu GX. An unstructured mesh based finite element simulation of wave interactions with non-wall-sided bodies. J Fluids Struct 2006;22:441–61. [15] Wang CZ, Wu GX, Drake KR. Interactions between fully nonlinear water wave and non-wall-sided 3D structures. Ocean Eng 2007;34:1182–96. [16] Wang CZ, Khoo BC. Finite element analysis of two-dimensional nonlinear sloshing problems in random excitations. Ocean Eng 2005;32:107–33. [17] Orlanski I. A simple boundary condition for unbounded hyperbolic flows. J Comput Phys 1976;21:251–69. [18] Cointe R, Geyer P, King B, Molin B, Tramoni M. Nonlinear and linear motions of a rectangular barge in a perfect fluid. In: Proc. of the 18th symp on naval hydro, Ann Arbor, Michigan; 1990. p. 85–98. [19] Tanizawa K. Long time fully nonlinear simulation of floating body motions with artificial damping zone. J Soc Naval Archit Jpn 1996;180:311–9.
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Fully nonlinear simulation of resonant motion of liquid confined between floating structures C.Z. Wang a, G.X. Wu b,⇑, B.C. Khoo a a b
Department of Mechanical Engineering, National University of Singapore, 10 Kent Ridge Crescent, Singapore 119260, Singapore Department of Mechanical Engineering, University College London, Torrington Place, London, WC1E 7JE, UK
a r t i c l e
i n f o
Article history: Received 10 May 2010 Received in revised form 9 November 2010 Accepted 15 December 2010 Available online 30 December 2010 Keywords: Resonance Floating bodies Potential theory with fully nonlinear boundary conditions High order finite element method
a b s t r a c t A finite element based numerical method is employed to analyze the resonant oscillations of the liquid confined within multiple floating bodies based on fully nonlinear wave theory. The velocity potentials at each time step are obtained through the finite element method (FEM) with quadratic shape functions. The matrix equation of the FEM is solved through an iteration. The waves at the open boundary are absorbed by the method of combination of the damping zone method and the Sommerfeld–Orlanski equation. Numerical examples are given by floating two rectangular cylinders, two wedge-shaped cylinders and two semi elliptic cylinders undergoing forced oscillations at the resonant frequencies. The numerical results are compared with the first and second order solutions by previous study, which showed that the first order resonance happens at the odd order natural frequencies x2i1 ði ¼ 1; 2; Þ and second order resonance at the half of even order x2i =2ði ¼ 1; 2; Þ for antisymmetric motions and happen at x2i ði ¼ 1; 2; Þ and x2i =2ði ¼ 1; 2; Þ for symmetric motions, and it is found that they are in good agreement in smaller amplitude motions within a sufficiently long period of time. However, difference begins to appear as time increases further. This happens even in smaller amplitude motions. In all the calculated cases, the results always become periodic in time eventually. Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction Resonance is one of the important features in nature. For a mechanical system in periodic motion, it occurs at a frequency when the inertial term is cancelled by the stiffness term. As a result its response to external excitation at such a frequency can become exceedingly large. One of the well known examples in fluid mechanics is sloshing in a tank [1–3]. Another typical example is the motion of the liquid confined between two floating bodies, as considered by Wang and Wu [4] for two dimensional rectangular cylinders, by Sun et al. [5] for two barges, by Malenica et al. [6], and Chen and Fang [7] for water confined by two ships with varying gap between their near sides. When resonance does occur, it has some serious implications to the safety and operation of a mechanic system. Thus to accurately predict resonant behaviour is hugely important. The analysis of the resonant behaviour is usually made based on the perturbation theory, starting with the linear theory. The linear theory usually captures the values of natural frequencies and the overall behaviour of the system at resonance quite well. It shows that when the excitation frequency x is equal to one of the natural ⇑ Corresponding author. E-mail address: [email protected] (G.X. Wu). 0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2010.12.014
frequencies xi ; i ¼ 1; 2; 3; :::, resonance will occur. Though subharmonics, one would then expect that resonance in the second order result would occur when x ¼ xi =2; i ¼ 1; 2; 3; :::. In fact, when the excitation has frequency components xej ; j ¼ 1; 2; 3; :::, resonance in the second order result could occur when xej1 þ xej2 ¼ xi . These have been confirmed in the analyses by Wu [3] for a sloshing tank and by Wang and Wu [4] for the liquid motion between two bodies floating on the water surface. While the perturbation theory can capture some important features of resonance, it does not always provide accurate results quantitatively. The main reason for this is that the perturbation theory is based on small parameter expansion. At resonance, the motion becomes very large. In fact, the motion can become infinite based on the linear theory when there is no damping in the system. This clearly contradicts with the foundation of the perturbation theory. Thus an appropriate approach at resonance is to use the fully nonlinear theory. The nonlinear theory has been used previously for the resonant behaviour of a sloshing tank [8,2]. Very little has, however, been done using the fully nonlinear theory for the motion of the liquid confined between two floating bodies. This has in fact important implications to the marine structures, such as a catamaran or ships very close to each during operation [5]. The present paper uses a high order finite element method to analyze the motion of the liquid confined between two floating
90
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
cylinders. While the method can be used for general shapes, cases considered include twin rectangles, twin wedges and twin semi ellipses. Extensive simulations are made at the first and the second order resonant frequencies. The agreement and disagreement between the perturbation theory and the fully nonlinear theory are discussed.
@/ 1 p ¼ q þ jr/j2 þ gy @t 2
ð8Þ
where q is the fluid density. The hydrodynamic force acting on the body can be expressed as *
F¼
Z
*
p n ds
ð9Þ
S
2. Mathematical formulation
~¼ M
Zb
*
pð~ r n Þds
ð10Þ
Sb
2.1. Governing equation and boundary conditions
*
We consider the problem of motion of the liquid confined between two floating bodies, as shown in Fig. 1. A Cartesian coordinate system oxy is defined in which the x-axis coincides with the undisturbed free surface, and the y-axis is positive upward. The free surface and the floating body surface are denoted as Sf and Sb respectively. The bottom of the liquid is assumed to be a horizontal plane at y = h and is denoted by Sbot . The fluid is assumed incompressible and inviscid, and the flow is irrotational. The liquid motion can then be described by a velocity potential /. It satisfies the Laplace equation in the fluid domain 8, or
r2 / ¼ 0
in 8
ð1Þ
The Eulerian form of the kinematic and dynamic boundary conditions on the free surface y ¼ g can be written as
@/ @ g @/ @ g ¼ 0 on Sf @y @t @x @x @/ 1 þ g g þ jr/j2 ¼ 0 on Sf @t 2
ð2Þ ð3Þ
where g is the acceleration due to gravity. Eqs. (2) and (3) can also be written as in the Lagrangian form
d/ 1 ¼ g g þ r/r/ dt 2 dx @/ dy @/ ¼ ; ¼ dt @x dt @y
ð4Þ ð5Þ
on Sb ; Sc and Sbot
ð6Þ
where V n is the velocity of the surface along the inward normal vector ~ n. In addition, the initial condition can be expressed as
/ðx; y ¼ n; t ¼ 0Þ ¼ /0 ðxÞ;
gðx; t ¼ 0Þ ¼ nðxÞ
r2 /t ¼ 0
ð7Þ
2.2. Evaluation of hydrodynamic forces Once the velocity potential has been found, the pressure in the fluid can be determined by the Bernoulli equation
ð11Þ
On the fixed boundary it satisfies
@/t ¼0 @n
ð12Þ
On the free surface /t is given by the Bernoulli equation as
1 /t ¼ gy r/r/ 2
On the rigid surface, the boundary condition is
@/ ¼ Vn @n
where r ¼ x i þy~ j is the position vector from the origin of the system to the point on the body surface. In Eqs. (9) and (10), the computation of integration of 1 jr/j2 þ gy over the body surface does not present real diffi2 culty. There are however some difficulties associated with ou/ot. This is the partial temporal derivative which is taken when the spatial point is fixed. This can be problematic when the body is in motion. It may be calculated by following a point on the body surface and using the finite difference method. However, it is often found that the force obtained this method has sawteeth behaviour, especially in the case of large amplitude motion. In order to avoid this, we treat /t as another harmonic function. The procedure can be summarized below. Assume the translational velocity of the body at the origin is ~ V and the angular velocity about the axis passing through the origin ~. In the fluid domain, the time derivative / satisfies the Laplace is X t equation
ð13Þ
On the moving boundary, the condition can be written as (Wu [9])
h i h i @ _ ~_ ~ ~ ~ @ r/ ~ @ ~ ~ ð/t Þ ¼ ~ þX r ðV r/Þ V þX r n V: @n @n @n
ð14Þ
~ indicates the derivative with respect to where the dot over ~ V and X time. If the acceleration is known, as in the forced motion for example, the above problem can then be solved to obtain /t . For a freely floating body, the acceleration will depend on the force which in turn depends on the acceleration. To overcome this mutual dependence, an auxiliary function can be introduced [9]. The acceleration can then be found before the pressure is known.
3. Finite element discretization and mesh generation We shall adopt the finite element method to solve the velocity potential problem at each time step. Quadrilateral parametric element (see Fig. 2) is used. The shape functions defined in a local coordinate system ~ n ¼ ðn; gÞ corresponding to element e with eight nodes may be expressed as
9 ðeÞ Ni ðn; gÞ ¼ 14 ð1 þ ni nÞð1 gi gÞð1 þ ni n þ gi gÞ ði ¼ 1; 2; 3; 4Þ > > = ðeÞ Ni ðn; gÞ ¼ 12 ð1 n2 Þð1 þ gi gÞ ði ¼ 5; 7Þ > > ; ðeÞ ði ¼ 6; 8Þ Ni ðn; gÞ ¼ 12 ð1 g2 Þð1 þ ni nÞ ð15Þ Fig. 1. Sketch of the computational domain.
This gives
91
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
The Galerkin method is, on the other had, suitable for general cases but needs to solve another matrix equation at each time step. When higher element is used, the derivatives can be obtained by differentiating the shape function directly. This has in fact been used by Wang and Khoo [16] and Wang and Wu [4]. We write n @/ X @Nj ¼ /j ; @x @x j¼1
n @/ X @Nj ¼ /j @y @y j¼1
ð19Þ
where
2
3
@N j 4 @x @N j @y
Fig. 2. Meshes for (a) two semielliptic cylinders; (b) two wedge-shaped cylinders. 8 X
x¼
ðeÞ
xk N k ;
y¼
k¼1
8 X
ðeÞ
yk N k
k¼1
within the element. The velocity potential / in the fluid domain can be expressed as
/ðx; yÞ ¼
n X
@n @x @g
@y @n @y @g
#1 2 @Nj 3 4 @n 5 @N j @g
ð20Þ
within each element. The second order derivatives @ 2 /=@x2 ,@ 2 /=@x@y and @ 2 /=@y2 can be obtained in a similar way through the Jacobian. The fourth order Runge–Kutta method is adopted for the integration with respect to time to update the wave elevation and the potential on the free surface. Let P denote one of these three variables ðx; y; /Þ. If Pn is the value at time step nDt, the new value at t ¼ ðn þ 1ÞDt can be obtained through the following equation
Pnþ1 ¼ Pn þ
Dt ðk1 þ 3k2 þ 3k3 þ k4 Þ 8
ð21Þ
where
/J N J ðn; gÞ
ð16Þ
J¼1
ki ¼
where /J is the potential at node J and n is the number of nodes. N j ðeÞ in the equation is equal to Nk if node J in the global system is node k of element e, and ðx; yÞ 2 e. Otherwise N j ¼ 0. Through the Galerkin method, we have
Z Z
5¼
" @x
i1 d i1 ði ¼ 1; 2; 3; 4Þ P nDt þ Dt; Q nþ 3 dt 3
and
Q n ¼ Pn ; 1
Q nþ3 ¼ Pn þ
2
r /Ni d8 ¼ 0
ð17Þ
8
using Green’s identity and the boundary conditions, we can obtain the following matrix equations
½Kf/g ¼ fFg
Dt k1 ; 3
2 1 Q nþ3 ¼ Pn þ Dt k1 þ k2 ; 3 Q nþ1 ¼ Pn þ Dtðk1 k2 þ k3 Þ
ð18Þ
where
4.2. Remeshing and smoothing
f/g ¼ ½/1 ; /2 ; . . . ; /I ; . . . ; /n T ðI R Sp Þ; Z Z rN I rNJ d8 ðI R Sp &J R Sp Þ; K IJ ¼
When the simulation is over a substantial period of time, the nodes on the free surface may cluster and cause elements to be distorted. In order to avoid this, nodes on the free surface should be rearranged every several time steps. A B-spline based remeshing technique [14] is used in the present paper. A smoothing technique called the energy method is also adopted to avoid the sawteeth problem on the free surface [14].
FI ¼
Z Sn
8
N I fn dS
Z Z
rNI
8
n X
J¼1 j2Sp
ðfp ÞJ rNJ d8 ðI R Sp Þ
Sp in above equations represents the Dirichlet boundary on which the potential, denoted by fp, are known, and Sn represents the Neumann boundary on which the normal derivative of the potential, denoted by fn, are known. The integration in the equation is calculated element by element with respect to ðn; gÞ through the Jacobian. Once the coefficients are found, the matrix equations are then solved through an iteration based on the conjugate gradient method with a symmetric successive over relaxation (SSOR) preconditioner. 4. Numerical procedures 4.1. Evaluation of derivatives of potential with respect to coordinates There are first and second order derivatives of the potential in Eqs. (4), (5), and (14), respectively. Approaches such as difference method [10–12] and Galerkin method [13–15] have been used associated with linear elements. The difference method is efficient but is only suitable for wall-sided bodies and structured meshes.
4.3. Numerical radiation conditions For long time simulations, an appropriate radiation condition should be imposed on the boundary Sc to minimize the wave reflection. Here, we use a combination of an artificial damping zone and Sommerfeld–Orlanski condition. The damping zone method is also called a sponge layer or an artificial beach. The artificial damping zone is an absorbing boundary. Outgoing waves will be first absorbed through artificial damping terms added to free surface boundary conditions upon their arrival at the damping zone and the most remanding waves will propagate out of the computational domain through the Sommerfeld–Orlanski condition [17]. Cointe et al. [18] adopted the damping zone technique for their 2D numerical wave tank simulations and this method was also used by Tanizawa [19] in fully nonlinear 2D simulations. In the damping zone, the following modified absorbing free surface boundary condition is used
92
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
Table 1 Number of elements and nodes for meshes one and two.
Mesh one Mesh two
NF1
NF2
ND
NB
NH
NE
NN
50 70
32 36
8 12
12 16
18 24
2580 4608
8117 14,337
d/ 1 ¼ g g þ r/r/ mðxÞ/ dt 2 dx @/ dy @/ ¼ ; ¼ mðxÞy dt @x dt @y
ð22Þ ð23Þ
Fig. 5. Histories of waves at the right side of cylinder one at x ¼ x02 and A = 0.05d.
with
(
mðxÞ ¼
xx0 2
ax 0
k
x0 6 x 6 x1 ¼ x0 þ bk 0 < x < x0 or x > x1
where x is the wave frequency and k is the linear wave length. The damping zone starts from point x0 and extends for a length of bk. Parameters a and b control the strength and the length of the damping zone respectively, and they may be obtained through numerical experiments. 5. Numerical results 5.1. Two rectangular cylinders in forced motions
Fig. 3. Wave runups along cylinder one, (a) left hand side, (b) right hand side at x ¼ x02 and A = 0.0125d; —— mesh one; – – – mesh two.
Wang and Wu [4] have undertaken some detailed analysis for the liquid motion confined between twin rectangular cylinders floating on the free surface based on the perturbation analysis. They confirmed that when the motion frequency of the cylinders is at one the natural frequencies very large transient wave motion between cylinders can occur based on the linear theory. They further showed that when the motion frequency of the cylinders is at half of one of the natural frequencies the second order transient motion of the liquid can also be very large. This analysis has raised the question of whether these behaviours are a result of the perturbation theory. This is because the perturbation theory is based on small motion assumption and the results obtained from this theory become inaccurate quantitatively when the motion is big. This question will be discussed in detail by the fully nonlinear theory below.
Fig. 4. Wave runups along cylinder one, (a) left hand side, (b) right hand side at x ¼ x02 and A = 0.0125d; —— DT = t/100; – – –DT = t/200.
Table 2 Number of elements and nodes for meshes three and four.
Mesh three Mesh four
NF1
NF2
ND
NB
NH
NE
NN
60 80
36 42
8 12
12 16
18 24
3048 5232
9573 16,261
Fig. 6. Comparisons of waves at the right side of cylinder one; (a) A = 0.0125d; (b) A = 0.025d; (c) A = 0.05d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
93
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
Fig. 7. Comparisons of forces on cylinder one between the fully nonlinear results and the second order solutions; (a), (b) A = 0.0125d; (c), (d) A = 0.025d; (e), (f) A = 0.05d; —— fully nonlinear; linear; linear + second order. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
We assume that the two cylinders are identical. The cylinder on the left hand side is called cylinder one and the one on the right is called cylinder two. The breadth at mean water level and draught of each cylinder are 2b and d respectively. The centre lines of the two cylinders are located at x = Lcy and Lcy, respectively. In the following simulations xn is defined as
Fig. 8. Comparison of waves at the right side of cylinder one at A = 0.0125d and x ¼ x02 within two hundred periods.
xn ¼
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi np g; 2ðLcy bÞ
n ¼ 1; 2; . . .
ð24Þ
This should be a good estimate for resonant frequencies when the cylinders are rectangular and the confined liquid resembles the sloshing in a tank (see, for example, Wu [3]). In particular, as discussed by Wu [3], when n is an odd number, the liquid motion is anti symmetrical and when n is an even number, the motion is symmetrical. For other cylinders such as twin wedges and twin semi
Fig. 9. Comparison of waves at the right side of cylinder one at x ¼ x02 within two hundred periods.
94
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
ellipses considered below, the resonant frequencies are expected to further differ from Eq. (24). The two cylinder are assumed to undergo the same vertical motion defined below
X ¼ A sin xt
h = 10d, and d = b and Lcy = 8b. We choose A = 0.0125d. Because the problem is symmetric about x = 0, the resonant behaviour may occur only when x is around x2n in Eq. (24). The real resonant
ð25Þ
where x and A are the motion frequency and amplitude, respectively. To avoid starting the motion suddenly with a non zero velocity, the displacement function in this equation could be modified initially and then returns to the normal one after a short period of time. For the results given below, the water depth is taken as
Fig. 13. Comparisons of forces on cylinder one between the fully nonlinear results and the second order solutions at x ¼ 0:5x02 and A = 0.05d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
Fig. 10. Histories of waves and forces for cylinder one at x ¼ x02 and A = 0.2d.
Fig. 14. Comparison of waves between the fully nonlinear result and the second order solution at x ¼ 0:5x02 and A = 0.2d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
Fig. 11. Wave profiles at t/T = 11, 11.05, 11.1, . . . , 12 at x ¼ 0:5x02 and A = 0.2d.
Fig. 12. Comparison of waves between the fully nonlinear result and the second order solution at x ¼ 0:5x02 and A = 0.05d.
Fig. 15. Comparisons of forces between the fully nonlinear results and the second order solutions at x ¼ 0:5x02 and A = 0.2d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
frequency x0n may be different from xn. We denote x0n ¼ C xn , where C is a constant. In fact, through detailed numerical tests, Wang and Wu [4] have found that first order resonance occurs when C = 1.01 or x ¼ x02 ¼ 1:01x2 . The control surfaces Sc at both ends are located at distance about four times wavelength away the nearest cylinder. Two meshes are used to test the numerical convergence. The details are given in Table 1, in which NF1 is the segment number along the free surface on the left of cylinder one or the right of cylinder two, NF2 the segments on the free surface between the two cylinders, ND and NB the segments along the vertical and horizontal faces of one cylinder respectively, NH the segment on the both control surfaces Sc at the far ends, NE and NN are the total number of elements and the total number of nodes in the whole fluid domain respectively. The results from the meshes are given in Fig. 3. It can be seen that they are graphically indistinguishable over the entire simulations of 200 cycles. Further convergence study with respect to time step is made using mesh one. Fig. 4 gives results with Dt ¼ T=100 and Dt = T/200 and the wave histories from the two time steps are in good agreement. These tests show that mesh one and Dt = T/100 can provide convergent results for this case. Further convergence study at a larger amplitude of A = 0.05d is also made. The details of the meshes are given in Table 2 and the result
95
Fig. 16. History of the second order wave at the right side of cylinder one at x ¼ 0:5x02 .
is given in Fig. 5. It can be seen that mesh three and Dt = T/200 give the convergent results. We now undertake detailed analysis of the nonlinear behaviour of the liquid motion at the natural frequency x = 1.01x2 by comparing the results with those obtained from the perturbation theory [4]. Three amplitudes A = 0.0125d, 0.025d and 0.05d have been chosen and meshes one with Dt = T/100, three with Dt = T/ 150 and three with Dt = T/200 have bee used for these amplitudes, respectively. The wave runups along the right side of cylinder one are given in Fig. 6. It can be seen that the linear and the combined linear and second order results are in very good agreement with those from the fully nonlinear theory at A = 0.0125d. When the amplitude is increased to A = 0.025d, there is some visible discrepancy between the results from the linear theory and the fully nonlinear theory, and the combined linear the second order theory give a better agreement. When the amplitude is increased further to A = 0.05d, the differences between the results from the perturbation theory and the fully nonlinear theory become obvious. The differences appear not only in the magnitudes of the peaks and troughs, but also in the phase. Comparisons are also made for the hydrodynamic forces without the contribution of initial buoyancy. Results from the perturbation theory and the fully nonlinear theory are given in Fig. 7. As in the case of wave elevation, the results are in good agreement at smaller motion amplitude and discrepancy appears when it increases. It is also observed that there are double peaks in the vertical force within one cycle at larger motion amplitude, which implies the high harmonic nonlinear effects. Figs. 6 and 7 show that when A increases the nondimensionalized nonlinear results are much smaller than those from the perturbation theory. In fact, the latter is expected to increase further due to the near resonance behaviour. At smaller A, however, the results are in good agreement over the entire period of simulation. It is important to realize, however, that this should be put in the context that the smaller the A is, the longer period will be over which the results from the perturbation theory and the fully nonlinear theory will be in good agreement. It is expected that as time progresses the discrepancy will develop even at small A. Thus we run the simulations over two hundred periods at A = 0.0125d. Fig. 8 shows that after about 40 periods, significant difference appears between the results from the perturbation theory and the fully nonlinear theory. The envelope of the wave amplitude from the nonlinear theory begins to oscillate at a lower frequency. As the confined liquid is linked to the open domain, some radiation
Fig. 17. Histories of waves and forces at x ¼ 0:5x02 and A = 0.2d.
Fig. 18. Histories of waves at the right side of cylinder one at (a) x ¼ 0:5x02 and (b) x ¼ 0:5x02 in sway in same direction (A = 0.05d).
96
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
damping is expected. Thus it is expected that the results from the perturbation theory will not increase indefinitely and an envelope of the amplitude oscillation will develop eventually, although the frequency will be much lower. Results from the nonlinear theory at A = 0.0125d, 0.025d and 0.05d are given in Fig. 9 for comparison. It can be seen that the nondimensionalized amplitude of the envelope decreases when A increases, while its frequency increases. Further simulation is made at even larger A = 0.2d. Fig. 10 gives the wave elevation and force. The result shows that after a very relatively short transition period, they become more or less periodic. All these results show that at a small A, the fully nonlinear results exhibit the resonant behaviour of the linear theory over a long period of time, but the results will eventually depart form those the perturbation theory when time progresses further. As large A, the resonant behaviour is far less obvious due to the highly nonlinear effects. Fig. 11 gives the development of the wave profile between t = 11T and t = 12T with an interval of 0.05T. It is seen that the wave between the two cylinders is much larger than those outside the two cylinders.
Following the work of Wu [3], Wang and Wu [4] also found that at x ¼ 0:5x02 , the linear theory does not show any resonance but the second order result does. We choose A = 0.05d and run the fully nonlinear simulation at x ¼ 0:5x02 . The results are shown in Figs. 12 and 13. The comparison shows that the results are in excellent agreement within those obtained from the combined linear and the second order theory for whole 40 cycles in the figure. This gives more evidence that the second order resonance does reflect the physical reality to large extent in some cases. Simulation is then made at A = 0.2d and the results are given in Figs. 14 and 15. The difference between results from the fully nonlinear theory and the perturbation theory becomes evident. This is obviously due to the same reason as that discussed right above. We run the simulation over a longer time of two hundred periods. Fig. 16 clearly shows the resonant behaviour of the second order wave elevation at x ¼ 0:5x02 (k is the linear wave number in the figure), which has been discussed in detail by Wang and Wu [4]. Fig. 17 gives the wave elevation at the right side of cylinder one and force from the nonlinear theory within 200 periods. The results become periodic quickly and there is no obvious resonant behaviour due to the nonlinear effects. We now consider some cases of cylinders in swaying motions. Two cases are simulated. Case one is the two cylinders moving in the same direction along x-axis with X ¼ A sin xt. Case two is the
Fig. 19. Histories of waves and forces for cylinder one at x ¼ 0:5x02 in sway in same direction; – – – A = 0.05d; —— A = 0.1d; . . . . . . A = 0.2d. Fig. 21. Numerical tests to determine the constant C.
Fig. 20. Histories of waves and forces for cylinder one at x ¼ x02 in sway in opposite directions; —— A = 0.0125d; A = 0.025d; A = 0.05d. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Fig. 22. Histories of waves at (a) x = Lcy 3b and (b) x = Lcy + 3b at x ¼ x02 for twin wedges in heave;—— A = 0.0125d; A = 0.025d; A = 0.05d.
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
two cylinders moving in opposite directions along the x-axis, which means that the motion for cylinder one is given as X ¼ A sin xt and X ¼ A sin xt for cylinder two. Simulations are made at the first and half of second resonant frequencies, or x01 and 0:5x02 (x01 1:068x1 and x02 1:01x2 ) in Case one. Fig. 18 shows the waves at the right side of cylinder one at A = 0.05d. Waves corresponding to both frequencies reach periodic state relatively early within two hundred periods. The nondimensionalized maximum value of wave at the right side is about 20.6 at around t/ T = 22.15. The peak is about 15.3 when the wave becomes periodic. The wave at the left side is, however, about 2.8 only. The wave at half of the second resonant frequency 0:5x02 has a more evident nonlinear effect, which may also be seen by comparisons of waves at the right side of cylinder one and forces on cylinder one at A = 0.05d, 0.1d and 0.2d in Fig. 19. Case two is similar to the vertical motion. The flow is symmetric about x = 0, and the first- and second order resonances happen at x ¼ x02 and x ¼ 0:5x02 respectively. We give results at x ¼ x02 only. Fig. 20 shows waves and forces at A = 0.0125d, 0.025d and 0.05d. The overall behaviour of the results is very similar to that in the heave motion. 5.2. Two wedge-shaped cylinders in forced motions A significant difference between this case and the rectangular cylinders is that the distance between the near sides of the two cylinders varies along the water depth. During the motion the distance between the two cylinders at the wave surface level also
97
varies. We notice that this distance is an important parameter for the natural frequency. Thus when its variation is included in the nonlinear model, the natural frequency changes with time, which makes the resonant effect less profound. On the other hand, when the body is non-wall-sided, it will have more significant nonlinear effect on the force. We consider a case in which the deadrise angle of each of the two identical wedges is 600 and Lcy = 16b (2b is the breath of the wedge at the calm water level when the body is stationary). The best approximation of the natural frequency is chosen by numerical tests based on the linear theory. Fig. 21 shows that the amplitudes of wave and forces on cylinder one are the largest when C = 1.06, which implies the resonance may occur at x ¼ x02 1:06x2 . We first consider three cases with amplitudes A = 0.0125d, 0.025d and 0.05d at x = 1.06x2 in the same vertical motion as the previous rectangular cylinders. On the free surface, there are 60 segments on the left side of cylinder one and the right side of cylinders two, 30 segments between the two cylinders and 16 segments on each cylinder surface and 18 segments on the control surface Sc, respectively, which corresponds to a total of 2700 elements and 8469 nodes. The time intervals for A = 0.0125d, 0.025d and 0.05d are T/100, T/150 and T/200, respectively. The histories of waves at x ¼ Łcy 3b are given in Fig. 22. It is noticed that all wave amplitudes increase at the earlier stage of the simulations and then reaches periodic state more quickly than the rectangular cases shown in Fig. 9. As mentioned above, one should note that
Fig. 23. Comparisons of waves between the fully nonlinear results and the second order solutions for twin wedges in heave; (a), (b) A = 0.0125d; (c), (d) A = 0.025d; (e), (f) A = 0.05d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
98
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
Fig. 24. Comparisons of forces on cylinder one between the fully nonlinear results and the second order solutions for twin wedges in heave; (a), (b) A = 0.0125d; (c), (d) A = 0.025d; (e), (f) A = 0.05d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
the wedges are non-wall-sided bodies. During the vertical motion, the distance between the near sides of the two wedges changes. Because of this variation of the configuration of the fluid domain confined between the two bodies, the natural frequency obtained
approximately based on the mean domain through the linear theory would have less significant effect. However, this effect is still not negligible. Fig. 22 shows that the wave peaks at the right side x ¼ Łcy þ 3b (confined between the wedges) is about four times
Fig. 25. Histories of waves at x = Lcy + 3b and x ¼ 0:5x02 for twin wedges in heave; (a) A = 0.05d; (b) A = 0.1d; (b) A = 0.2d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
Fig. 26. Histories of waves at x = Lcy + 3b and forces on cylinder one at x ¼ 0:5x02 for twin wedges in heave.
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
99
that at the left side x ¼ Łcy 3b (open water) at A = 0.0125d and nearly six times at A = 0.05d. Comparisons are made between the fully nonlinear results and the linear and linear plus second order solutions at A = 0.0125d, 0.025d and 0.05d are shown in Figs. 23 and 24. The waves correspond to two locations: one at x ¼ Łcy þ 3b and the other at x = 0, which is the middle point between the two cylinders. For the case of A = 0.0125d, all the results are generally in good agreement even though the fully nonlinear result is slightly smaller than the linear and linear plus second order. For the largest A = 0.05d, the differences between them increase. Unlike the rectangular cases, the linear and linear plus second order solutions in these two cases at this first resonant frequency do not always increase with time. They become periodic, or their peaks or troughs become constants within the period of simulation. This suggests that the resonance may have less effect when a body has flare. It is interesting to notice that the above water flare is not included in the fully linear theory. Thus it is quite likely that the linear theory could overestimate the resonant effect for the flare bodies. Fig. 27. Histories of waves at x = Lcy + 3b and x ¼ x02 for twin wedges in sway in same direction; (a) A = 0.05d; (b) A = 0.1d; (c) A = 0.2d; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
Fig. 28. Histories of waves at x = Lcy + 3b and forces on cylinder one at x ¼ x02 for twin wedges in heave; – – – A = 0.05d; —— A = 0.1d; . . . . . . A = 0.2d.
Fig. 29. Histories of waves at x = Lcy + 3b and forces on cylinder one at x ¼ 0:5x02 for twin wedges in heave; – – – A = 0.1d; —— A = 0.2d; . . . . . . A = 0.4d.
Fig. 30. Histories of waves at (a) x = Lcy 3b (b) x = Lcy + 3b and (c) x = 0 at x ¼ x02 for twin wedges in sway in opposite directions; – – – A = 0.0125d; —— A = 0.025d; . . . . . . A = 0.05d.
Fig. 31. Histories of forces on cylinder one at x ¼ x02 for twin wedges in sway in opposite directions; – – – A = 0.0125d; —— A = 0.025d; . . . . . . A = 0.05d.
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
sults from the fully nonlinear theory at A = 0.05d, 0.1d and 0.2d are given in Fig. 28. These figures show that the nonlinear effect does not seem to be as strong as the heave motion. Results at 0:5x02 are given in Fig. 29. The nonlinear effect becomes more significant. In
4 2
h/A
Further simulations at x ¼ 0:5x02 are made. The results of waves at the right side of cylinder one at A = 0.05d, A = 0.1d and A = 0.2d together with the linear and linear plus second order solutions are given in Fig. 25, long after they have reached the periodic state. The comparisons with the linear plus second order show a good agreement at A = 0.05d. Discrepancy is evident at A = 0.1d and is even more at A = 0.2d due to stronger nonlinearity. Further comparisons between waves and forces of the fully nonlinear results at A = 0.1d, A = 0.2d and A = 0.4d are given in Fig. 26 and the nonlinear effect is strong. Simulations in horizontal motions are also made at the first and half of second resonant frequency x01 ðx01 1:16x1 , which is obtained in a similar way as in Fig. 21) and 0:5x02 (x02 1:06x2 ) when the wedges move in the same direction. Fig. 27 shows the wave histories at x = Lcy + 3b and x = 0 with x ¼ x01 and at three amplitudes A = 0.05d, 0.1d and 0.2d. It is seen that the differences between the fully nonlinear results and the linear and linear plus solutions are small even for a larger amplitude A = 0.2d. Further re-
0 -2 -4 0
5
10
15
20
25
30
35
190
195
200
t/T 1
Fx/rgdA
100
0 -1 0
5
10
15
20
25
30
35
190
192
194
196
198
200
192
194
196
198
200
Fy/rgdA
t/T 1 0 -1 0
5
10
15
20
25
30
35
190
t/T Fig. 34. Histories of waves at x = Lcy + 1.5b and forces cylinder one at x ¼ 0:5x02 for twin semi elliptic cylinders in heave; —— fully nonlinear; – – – linear; . . . . . . linear + second order.
Fig. 32. Histories of waves and forces on cylinder one at x ¼ 0:5x02 for twin wedges in sway in opposite directions; – – – A = 0.1d; —— A = 0.2d; . . . . . . A = 0.4d.
Fig. 35. Histories of second order wave at x = Lcy + 1.5b and x ¼ 0:5x02 for twin semielliptic cylinders in heave.
Fig. 33. Histories of waves at x = Lcy + 1.5b and forces on cylinder one at x ¼ x02 for twin semielliptic cylinders in heave; – – – A = 0.005d; —— A = 0.01d; . . . . . . A = 0.02d.
Fig. 36. Comparison of linear and second order solutions at x ¼ 0:5x02 and A = 0.05d for twin semielliptic cylinders in heave.
C.Z. Wang et al. / Computers & Fluids 44 (2011) 89–101
particular, the time history curves have departed further away from the sinusoidal shape. Simulations of horizontal motions when two wedges move in the opposite directions at the second and the half of second frequencies x02 and 0:5x02 are also made. The result for x ¼ x02 and x ¼ 0:5x02 are given in Figs. 30–32, respectively. This is a symmetrical case and behaviour of the results is therefore similar to that in the heave motion. 5.3. Two semi elliptic cylinders in forced motions The previous examples are for bodies with no curvature. Here, we consider the case of two semi elliptic cylinders. b = 0.75 is the semi minor axis of each cylinder in the x direction and d = 1/b is the semi major axis in the y direction. The horizontal distance between the centre lines of the two cylinders is Lcy = 8b. Simulations are made for vertical motion at x ¼ x02 and x ¼ 0:5x02 respectively 0 x2 1:01x2 . The wave elevation at xc ¼ Lcy þ 1:5b, and forces on cylinder one with x ¼ x02 at three amplitudes A = 0.005d, 0.01d & 0.02d are given in Fig. 33. It should be noted although the body has curvature, it is wall sided at the mean water line, i.e., the tangential of the body surface at its intersection with the free surface is in the vertical direction. The variation of the distance between the near sides of the cylinders at the water line is smaller than that of the wedges. Thus the effect of the natural frequency is less significant than that in the case of two wedges. This is reflected in Fig. 33 by the longer period over which the amplitude continues to increase and the large amplitude when the motion becomes periodic. The figure also shows that at the periodic state, the nonlinear effects are clearly very important. The results at the half of second resonant frequency are given in Fig. 34. In the figure, the linear and linear plus second order solutions are also given for comparison. It is seen that results from these methods are in good agreement for both the wave elevation and the force during the first thirty periods and difference between them begins to develop. It can be observed that the nonlinear effects on both the wave elevation and the force are strong. Based on the second order theory, the second order solution in this case may increase to become very large, as discussed in Wu [3]. This however will not be infinite, as discussed in Wang and Wu [4] due to the radiation damping. This is reflected in the history of the second order wave at x ¼ Lcy þ 1:5b in Fig. 35. Compared with the case of two rectangles in Fig. 16, the increase of the wave elevation stops much earlier. It is also interesting to see in Fig. 36 that the second order wave becomes larger than the linear one after t/ T = 14.8. This implies that full nonlinear theory may be needed. 6. Conclusions A higher order finite element method has been used to analyze the fully nonlinear motion of the liquid combined between floating twin-cylinders in forced oscillations. The velocity potential at each time step is obtained through solving matrix equations based on an iteration method. The fourth order Runge–Kutta method is used to update the wave elevation and potential on the free surface at each
101
time step. The radiation condition is imposed through a combination of the damping zone method and the Sommerfeld–Orlanski equation. Extensive simulations have been made for twin rectangular, wedge-shaped and semi elliptic cylinders in vertical and horizontal motions. It confirms the first order and the second resonant behaviour at small body motion amplitude, predicted by the linear and the second order theory. The wave elevation amplitude increases over a long period of time in the perturbation theory but the increase will stop much earlier in the fully nonlinear theory. This becomes more evident when the body motion amplitude increases. It is also found that the resonant effect is very significant for the rectangular cylinders. When the body has flare or curvature at the water line, the resonant effect is less significant even though it may be still important. This may suggest that the linear theory could over predict the resonant effect for bodies with flare or curvature. References [1] Faltinsen OM. A numerical non-linear method of sloshing in tanks with twodimensional flow. J Ship Res 1978;22:193–202. [2] Faltinsen OM, Timokha AN. Adaptive multimodal approach to nonlinear sloshing in a rectangular tank. J Fluid Mech 2001;432:167–200. [3] Wu GX. Second order resonance of sloshing in a tank. Ocean Eng 2007;34:2345–9. [4] Wang CZ, Wu GX. Analysis of second order resonance in wave interactions with floating bodies through a finite element method. Ocean Eng 2008;35:717–26. [5] Sun L, Taylor PH, Eatock Taylor R. First and second order wave effects in narrow gaps between moored vessels marine operations specialty symposium, Singapore; 2008. p. 113–24. [6] Malenica S, Orozco JM, Chen XB. Some aspects of multibody interactions in seakeeping. In: Proceedings of 15th international offshore and polar engineering conference, Seoul, Korea; 2005. [7] Chen GR, Fang MC. Hydrodynamic interactions between two ships advancing in waves. Ocean Eng 2001;28:1053–78. [8] Wu GX, Ma QW, Eatock Taylor R. Numerical simulation of sloshing waves in a 3D tank based on a finite element method. Appl Ocean Res 1998;20:337–55. [9] Wu GX, Eatock Taylor R. The coupled finite element and boundary element analysis of nonlinear interactions between waves and bodies. Ocean Eng 2003;30:387–400. [10] Ma QW, Wu GX, Eatock Taylor R. Finite element simulation of fully nonlinear interaction between vertical cylinders and steep waves. Part 1: methodology and numerical procedure. Int J Num Methods Fluids 2001;36:265–85. [11] Ma QW, Wu GX, Eatock Taylor R. Finite element simulation of fully nonlinear interaction between vertical cylinders and steep waves. Part 2: numerical results and validation. Int J Num Methods Fluids 2001;36:287–308. [12] Wang CZ, Wu GX. Time domain analysis of second order wave diffraction by an array of vertical cylinders. J Fluids Struct 2007;23:605–31. [13] Wu GX, Eatock Taylor R. Finite element analysis of two dimensional nonlinear transient water waves. Appl Ocean Res 1994;16:363–72. [14] Wang CZ, Wu GX. An unstructured mesh based finite element simulation of wave interactions with non-wall-sided bodies. J Fluids Struct 2006;22:441–61. [15] Wang CZ, Wu GX, Drake KR. Interactions between fully nonlinear water wave and non-wall-sided 3D structures. Ocean Eng 2007;34:1182–96. [16] Wang CZ, Khoo BC. Finite element analysis of two-dimensional nonlinear sloshing problems in random excitations. Ocean Eng 2005;32:107–33. [17] Orlanski I. A simple boundary condition for unbounded hyperbolic flows. J Comput Phys 1976;21:251–69. [18] Cointe R, Geyer P, King B, Molin B, Tramoni M. Nonlinear and linear motions of a rectangular barge in a perfect fluid. In: Proc. of the 18th symp on naval hydro, Ann Arbor, Michigan; 1990. p. 85–98. [19] Tanizawa K. Long time fully nonlinear simulation of floating body motions with artificial damping zone. J Soc Naval Archit Jpn 1996;180:311–9.