- Email: [email protected]
Wave–structure interaction using coupled structured–unstructured finite element meshes
Applied Ocean Research 25 (2003) 63–77 www.elsevier.com/locate/apor
Wave –structure interaction using coupled structured –unstructured finite element meshes M.S. Turnbulla, A.G.L. Borthwickb,*, R. Eatock Taylorb a HR Wallingford Ltd, Howbery Park, Wallingford, Oxon OX10 8BA, UK Department of Engineering Science, University of Oxford, Parks Road, Oxford OX1 3PJ, UK
b
Received 31 August 2002; accepted 3 April 2003
Abstract The interaction of inviscid gravity waves with submerged fixed horizontal structures is modelled in two-dimensions using a finite element numerical wave tank. An adaptive hybrid coupled mesh is utilised that tracks the free surface through vertical movement of the free surface nodes in a Lagrange – Eulerian scheme. The hybrid mesh consists of a combination of structured and Voronoi unstructured meshes, with the submerged structure located in the Voronoi mesh region. Validation tests include free and forced sloshing in a rectangular tank, regular progressive wave propagation in a flume, and regular wave loading on a horizontal cylinder. The results are found to be in close agreement with analytical potential theory solutions for waves of small amplitude. Non-linear effects are noted for steeper waves. The wave-induced force components on the horizontal cylinder match the expected results from Ogilvie’s [J Fluid Mech 16 (1963) 451] linear theory and Vada’s [J Fluid Mech 174 (1987) 23] second order model. Wave interactions with a pair of submerged horizontal cylinders spaced at half the wavelength of undisturbed regular waves are examined. q 2003 Elsevier Ltd. All rights reserved. Keywords: Finite elements; Wave tank; Cylinder; Forces
1. Introduction Steep free surface waves give rise to highly non-linear forces on fixed offshore structures, and can cause the phenomenon of ringing, which has been the subject of considerable recent research (see e.g. Chaplin et al. [1]). Theoretical modelling of the interaction of steep gravity waves with complicated offshore structures is a challenging task because both the kinematic and dynamic free surface boundary conditions are non-linear and the flow domain has an irregular outer boundary surrounding an interior obstacle or obstacles. To deal with the awkward geometries encountered by offshore engineers in practice, it is necessary to develop efficient and accurate numerical models that can cope with steep wave non-linearity and the presence of submerged obstacles. Longuet-Higgins and Cokelet [2] assumed periodicity and hence simulated steep Stokes waves in infinite depth * Corresponding author. Tel.: þ 44-1865-273047; fax: þ 44-1865273010. E-mail address: [email protected] (A.G.L. Borthwick). 0141-1187/03/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0141-1187(03)00032-4
using a boundary element method. Vinje and Brevig [3] extended this approach to finite depth with floating bodies. Cao et al. [4] used a time stepping technique with the Desingularised Boundary Integral Equation Method (DBIEM) to study various non-linear free surface problems including the non-linear waves generated by free surface pressure disturbances. The same method was used by Celebi et al. [5] to investigate non-linear wave interactions with a stationary vertical cylinder. Wang et al. [6] developed a multi-subdomain approach that significantly improved the performance of their 2D boundary element numerical wave tank. Wu and Eatock Taylor [7,8] developed a finite element numerical wave tank using structured meshes that stretched between the bed and liquid free surface. Wu and Eatock Taylor demonstrated that the finite element method could be as efficient computationally as the boundary element method for inviscid gravity wave calculations. The present paper implements the finite element formulation derived by Wu and Eatock Taylor [7,8] on a combination of unstructured Voronoi and structured meshes.
64
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
George [9] has reviewed the generation of finite element meshes, and broadly classifies the different approaches as enumerative, functional mapping, numerical mapping, hierarchical, block decomposition, advancingfront and Voronoi –Delaunay’s construction. Of the foregoing, the hierarchical and latter two unstructured approaches are becoming popular due to their inherent robustness, efficiency and ease of adaptation. For example, Yerry and Shephard [10,11] and Greaves and Borthwick [12] have developed 2D quadtree and threedimensional octree hierarchical mesh generators. Greaves et al. [13] solved free surface sloshing in a tank on a finite element quadtree mesh. Zhu et al. [14] have developed a viscous flow numerical wave tank using triangularised quadtree finite elements where the free surface has been tracked in a Lagrangian manner. Lo [15] proposed the advancing front technique that was extended by Peraire et al. [16]. Lee and Hobbs [17] describe an advancing front method for the creation of adaptive finite element meshes of triangles and quadrilaterals. Chan [18] solved a variety of free surface flow problems including water sloshing in a tank, solitary wave propagation, and flow over a submerged hydrofoil using a finite volume solution of the Navier – Stokes equations on unstructured 2D meshes. In an alternative approach to unstructured mesh generation, Lawson [19] proposed an algorithm based on edge swapping to create Delaunay triangulations. A more efficient Voronoi –Delaunay algorithm was proposed by Bowyer [20] and Watson [21]. Mavriplis [22] used Voronoi meshes adaptively for viscous flow in two dimensions. Further work on the generation of Voronoi meshes was carried out by Borouchaki and George [23]. Rivara and Inostroza [24] compared two longest side bisection techniques for the automatic refinement of Delaunay triangulations. Xu et al. [25] presented an automatic coarsening procedure. As far as the present authors are aware, unstructured Voronoi meshes have not yet been applied to an inviscid numerical wave tank. This paper aims to use a flexible mesh generation strategy with a Voronoi unstructured mesh applied in the vicinity of a submerged cylinder and structured meshes elsewhere in regions away from cylinder. The validation tests have been selected to check that the model properly simulates free surface waves of increasing steepness, and gives accurate results for a submerged horizontal cylinder under regular wave loading. The technique, though limited to two-dimensions in this paper, is readily extendible to three-dimensions and to more complicated structural configurations typical of offshore structures.
[20]. An initial mesh is created, consisting of a square divided diagonally into two triangular elements and having a total of four nodes within which the computational domain fits. Each of these elements has an associated circumcircle that intersects its nodes. Nodes are inserted that delineate the boundary of the flow domain. Each new node is checked, and if it lies within a circumcircle, then the element is removed, and new elements created to occupy the space left behind by connecting the surrounding nodes to the node under consideration. Following this process, any elements created outside the flow domain boundary are deleted. A recursive procedure is then invoked. Each element is checked in turn. If its area exceeds a prescribed criterion a new node is inserted such that it makes an equilateral triangle when joined to the shortest side of the element. Overlapping nodes are coalesced. This procedure is repeated until no more changes occur to the mesh. The element limiting area criterion is as follows. For a uniform mesh, the maximum area of any element is set to a constant value throughout the domain. For a variable density mesh, an analytical formula gives the desired length of a line segment, l; anywhere in the domain. Consider the domain shown in Fig. 1, which is to have uniform node spacing on its upper and lower outer boundaries, and variable spacing on its lateral outer boundaries. This domain could represent a submerged horizontal cylinder under waves, where a fine mesh would be required near the free surface and close to the cylinder. Point P within the domain is located a radial distance r from the inner circular boundary and a distance y measured upwards from the lower outer boundary. The distance between the upper and lower outer boundaries at P is given by h: R is a prescribed radial distance measured outward from the boundary of the circle. Let s; b1 and b2 be lengths of line segments at the upper outer, bottom outer, and inner cylinder boundaries, respectively. Also, d is the number of line segments on the lateral outer boundaries. If b2 is smaller than either s or b1 ; then the line segment l is set equal to the smaller of la and lb ; otherwise, l is set equal to the larger of la and lb ; where f2 A b1 b la ¼ s and lb ¼ b1 2 ð1Þ s b1
2. Voronoi mesh generator The Voronoi method used herein generates a mesh using the node insertion technique proposed by Bowyer
Fig. 1. Definition sketch.
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
65
Fig. 2. Voronoi meshes generated around a circle in a square domain.
where A¼
log10 ðf1B þ 1Þ ; log10 ð2Þ
f1 ¼ 1 2 ðy=hÞ; B ¼ ðs=b1 Þðd=ðd21ÞÞ þ 1 and f2 ¼ 1 2 ðr=RÞ: The upper limit on the size of an element is given by q1 l2 where q1 is a prescribed parameter. During initial mesh generation, the following rules are invoked to improve mesh quality. † Before a new node is incorporated into the mesh, it is checked whether it lies closer than a prescribed distance q2 l to an existing node. If it is closer than q2 l; the hypothetical node is not inserted. † If the smallest side of an element is larger than q3 l then no node is inserted into that element, regardless of its size. This rule prevents nodes being created too far from the boundary initially. The values of q1 ; q2 and q3 can all be adjusted in order to optimize the quality of the mesh produced. After the initial mesh has been created, smoothing in four iterative sweeps is then performed to improve mesh quality.
Each node is repositioned so that it lies approximately at the average location of the nodes to which it is connected. In the test cases considered herein, the mesh parameters, q1 ; q2 and q3 ; have been set to 0.6, 0.5 and 10, respectively, after extensive numerical tests (Turnbull, [26]). Fig. 2 shows the results of the Voronoi method applied to a square domain containing a circle. The analytical formulas (1) were used to specify the mesh density, with R equal to twice the circle radius. For the uniform square, s ¼ b1 : In Fig. 2(a), the mesh is uniform, with 32 segments on the circumference of the circle for b2 ¼ s: The mesh in Fig. 2(b) is finer in the region around the circle which has 64 segments for b2 ¼ s=2: The converse occurs in Fig. 2(c) where there are only 16 segments around the circle for b2 ¼ 2s and the mesh is coarser in the vicinity of the circle. In Fig. 2(d), the mesh is uniform, but the circle is not at the centre of the domain. In all these cases, the Voronoi method has produced elements of acceptable quality, with the correct mesh density where specified. Fig. 3 shows several examples of Voronoi meshes for a domain of width equal to twice its height, and containing a circular inner boundary. Fig. 3(a) presents the mesh
66
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
the kind of mesh that would be used for finite element simulations, with finer mesh density close to the free surface and near the cylinder, and coarser mesh density towards the bottom of the domain, where velocities are expected to be smaller.
3. Governing equations Consider an incompressible, inviscid liquid undergoing irrotational flow in two-dimensions. Mass conservation may be expressed by 72 f ¼ 0;
ð2Þ
where f is the velocity potential, and 7 is the gradient vector. The non-linear free surface kinematic boundary condition is
›f ›h ›f ›h ¼ þ ; ›z ›t ›x ›x
ð3Þ
where h is the elevation of the free surface vertically above still water level, x is distance measured horizontally along the free surface, z is distance vertically above the free surface and t is time. The non-linear dynamic free surface boundary condition is
›f 1 þ 7f7f ¼ 0 ð4Þ 2 ›t in which g is gravitational acceleration. At solid boundaries, such as the fixed bed, lateral walls and submerged obstacle surfaces, gh þ
›f ¼ q; ð5Þ ›n where q is the normal velocity at the boundary and n is normal to the boundary and is directed out of the fluid domain. At open boundaries, the Sommerfeld radiation condition is ›f ›f ¼ 2c ›t ›x where c is the wave celerity.
ð6Þ
4. Finite element formulation Fig. 3. Meshes for domains with circular inner boundary and varying mesh density.
generated when the upper boundary is straight, whereas Fig. 3(b) –(d) depict meshes where the upper boundary is sinusoidal. These meshes are representative of those that would be produced during the finite element simulation of waves passing over a submerged horizontal cylinder. In each case the cylinder has 32 segments and b2 ¼ 0:5: The number of segments vertically, d is 16. In Fig. 3(a) and (b), s ¼ b1 ¼ 1:0; in Fig. 3(c) s ¼ 0:5 and b1 ¼ 1:0; in Fig. 3(d), s ¼ 0:5 and b1 ¼ 2:0: Fig. 3(d) shows
4.1. Discretisation The continuous velocity potential, f; is expressed in terms of nodal potentials fi and corresponding linear shape functions Ni ðx; zÞ as
f¼
N X
fi Ni ðx; zÞ
ð7Þ
i¼1
where N is the total number of nodes in the mesh. In the finite element mesh, each element is triangular with three corner nodes.
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
At each time step, a solution of the Laplace equation is obtained using the Galerkin technique. Wu and Eatock Taylor [8] integrated the product of Eq. (2) and the nodal shape functions over the domain, V; to obtain ð ð ð ›f dS 2 Ni 72 f dV ¼ Ni 7Ni 7f dV ¼ 0; ð8Þ ›n V S V
system: 0 P X B lkx lkx B B k¼1 B B P BX @ lkx lkz k¼1
where S is the boundary of V and i ¼ 1; …; N: Wu and Eatock Taylor then substituted f with its nodal approximation and replaced ›f=›n with q to give
67
8 1 P X 8 9 > ›f > ›f > > C> lkx k > > > C> < = < C ›x k¼1 k¼1 ›l C ¼ C > > P P C> X X > ›f : ›f > ; > > > lkz lkz A ›z lkz k > : ›l P X
lkx lkz
k¼1
k¼1
9 > > > > > = > > > > > ;
ð13Þ
where S1 and S2 are the surfaces where f and ›f=›n are specified. The terms on the right hand side of the equation are evaluated using the boundary conditions. In matrix form, the system of equations is,
in which lkx and lkz are the x and z components of lk divided by its magnitude of lk : The method is suitable for any kind of unstructured mesh. A higher order calculation is used to estimate the vertical velocity component at nodes on the free surface. The method requires that each free surface node and a prescribed number of nodes below it are situated along a vertical line. This restriction was incorporated within the Voronoi mesh generator. Implementation is as follows. Let z be the vertical distance upwards from the mean free surface. For a free surface node with elevation z1 ; the three nodes vertically below it are located at elevations, whose positions are given by
½A{f} ¼ {B}
z2 ¼ z1 2 Dz
ð R
7Ni
¼2
N X
fj 7Nj dRljS1
j¼1
ð R
7Ni
N X
fj 7Nj dRlj[S1 þ
j¼1
ð S2
Ni q dS
ð9Þ
ð10Þ
where {f} is a vector containing nodal velocity potentials and ½A and {B} have the following coefficients ð 7Ni 7Nj dR or Aij ¼ R
ð11Þ
Aij ¼ 1
if i ¼ j and i [ S1
Aij ¼ 0
if ði [ S1 or j [ S1 Þ and i – j
and Bi ¼ 2
N X j¼1
fj
ð
R
7Ni ·7Nj dR lj[s1 þ
ð S2
Ni q dS ð12Þ
or Bi ¼ fi if i [ S1 : In the above, fi is specified on the free surface ði [ S1 Þ as so does not need to be calculated there. The relations following the word ‘or’ in Eqs. (11) and (12) allow the known values to be transferred from the vector to the vector {f}: Gaussian elimination is used to solve matrix Eq. (10).
z3 ¼ z1 2 aDz
z4 ¼ z1 2 bDz
ð14Þ
where Dz is the spacing between the uppermost pair of nodes, and a and b are real numbers. The velocity potential is expressed as a power series, z 2 z1 z 2 z1 2 f ¼ A0 þ A1 þ· · · þ A2 Dz Dz z 2 z1 nv 21 ; ð15Þ þ An21 Dz in which nv is the number of nodes in the vertical that are utilised in estimating the vertical velocity component. Eq. (15) can be expressed as nv simultaneous equations by inserting values for f and z at the free surface node and the nv 2 1 nodes vertically beneath it. The vertical velocity component at the free surface is calculated by differentiating Eq. (15) and then taking z to be equal to z1 : When nv ¼ 3 and the nodes are unevenly spaced in the vertical, the vertical velocity component at the free surface becomes: w¼
›f ða2 2 1Þf1 2 a2 f2 þ f3 ¼ : ›z aða 2 1ÞDz
4.2. Surface velocity calculation
For a ¼ 2; the nodes are evenly spaced, and
The least-squares method used to estimate the horizontal free surface velocity component is implemented as follows. Consider an arbitrary free surface node that is connected to p underlying nodes in the finite element mesh. Let lk denote the position vector of the k th node relative to the arbitrary free surface node and connected to it. The potential gradient vector is obtained from the following least-squares equation
w¼
3f1 2 4f2 þ f3 : 2Dz
ð16Þ
ð17Þ
For evenly spaced nodes when nv ¼ 5; the following higher order formula results: w¼
25f1 2 48f2 þ 36f3 2 16f4 þ 3f5 12Dz
ð18Þ
68
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
4.3. Updating the free surface
5. Results
A mixed Euler –Lagrange approach is used where the free surface nodes are restricted to move solely in the vertical direction. The free surface elevation is updated using a fourth-order Runge – Kutta integration of the free surface kinematic boundary condition, expressed as
5.1. Free and forced sloshing in a rectangular tank
›h ›h ¼w2u : ð19Þ ›t ›x Similarly, the free surface potential is updated using the dynamic free surface boundary condition: ›f 1 ›h ¼ 2gh 2 ðu2 þ w2 Þ þ w : 2 ›t ›t
ð20Þ
Fig. 4 shows initial entirely Voronoi and combined Voronoi-structured meshes used to check mesh convergence for free sloshing waves of amplitude to depth ratio, a=h ¼ 0:1: Three levels of coarseness were considered: coarse, with 32 elements along the upper free surface boundary, and 16 vertically along the lateral walls; medium, with 64 free surface and 32 vertical wall elements; and fine, with 128 free surface and 64 lateral wall elements. The structured mesh is created using a linear stretching between the free surface and bed (Wu and Eatock Taylor, [7]). The Voronoi meshes are symmetric about the vertical centre of the tank, and contain greatest cell variations at about mid-depth. Fig. 5
Fig. 4. Initial meshes: ap ¼ 0:1:
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
69
Fig. 5. Time history of the free surface elevation at the centre of the tank. Free sloshing: ap ¼ 0:1: Voronoi mesh. Dotted line: 32 by 16; dashed line: 64 by 32; solid line: 128 by 64.
Fig. 6. Time history of the free surface elevation at the centre of the tank. Free sloshing: ap ¼ 0:001: Solid line: Voronoi; thin line: Voronoi-structured; dashed line: second order analytical.
Fig. 7. Time history of the free surface elevation at the centre of the tank. Free sloshing: ap ¼ 0:1: Solid line: Voronoi; thin line: Voronoi-structured; dashed line: second order analytical.
70
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
shows the motions of the free surface elevation at the tank centre simulated using the three combined Voronoistructured meshes. The 32 by 16 mesh gives slightly diffusive results with the wave crests tending to be underpredicted, in comparison with the corresponding results obtained using the finer meshes. Simulations using the 64 by 32 and 128 by 64 meshes are in almost exact agreement, indicating that mesh convergence has been achieved. It should pffiffiffiffibe noted that the non-dimensional time step, Dtp ¼ Dt g=h; was set to be Dtp ¼ 0:1; 0.05 and 0.025 for each mesh, respectively, in order of its fineness. The entirely Voronoi mesh scheme gave results that were indistinguishable from those using the combined mesh approach. Mesh and time step convergence occur for a 64 by 32 fully Voronoi mesh with non-dimensional time step, Dtp ¼ 0:05: Figs. 6 and 7 demonstrate the increasing non-linear behaviour of the free surface as the relative wave amplitude ap ¼ a=h is increased from 0.001 to 0.1. For small amplitude free sloshing where ap ¼ 0:001; the free surface motions are identical to the analytical solution from second order potential theory (Wu and Eatock Taylor, [7]) which itself is very close in this case to linear theory. For larger amplitude sloshing where ap ¼ 0:1; the free surface waves are of longer period and have higher extreme maxima than those from the analytical solution. The non-linear free surface motions are in excellent agreement with those obtained using alternative (fully non-linear) pseudospectral and finite element numerical approaches by Chern et al. [27] and Wu and Eatock Taylor [7], respectively. An asymmetric sloshing wave is next considered. It was first utilised in a DnV study (Nestega˚rd, [28]) that compared results from several industrial and academic numerical codes for fully non-linear free surface sloshing. The tank has length 160 m and still water depth 70 m. At time t ¼ 0; the free surface elevation is given by " 2 # " 2 # x x exp 2 ð21Þ hðx; 0Þ ¼ as 1 2 bs gs where as ¼ 12 m, bs ¼ 53 m and gs ¼ 76 m. Fig. 8 shows the initial combined Voronoi-structured mesh. Table 1 presents the free surface elevation, and surface velocity components at x ¼ 60 m and t ¼ 9.2 s obtained by the participants in the DnV study and using the present method with unstructured and hybrid meshes. It is clear that the present scheme gives results in close agreement with those of the majority of the participants. Fig. 9 shows the timedependent free surface motions at x ¼ 60 m (corresponding to xp ¼ 0:8571 where xp ¼ x=h). After non-dimensional time tp ¼ 10; the most significant component appears to be third-order. Horizontally forced sloshing in a base-excited tank has also been considered. Fig. 10 shows the bursting characteristics of the free surface motions at the left wall when the excitation frequency v is close to the lowest natural sloshing
Fig. 8. Initial meshes for asymmetric sloshing wave.
frequency v0 : The sloshing motions exhibit a beating effect, caused by interactions between the natural sloshing frequency of the liquid and the forcing frequency. There is close agreement evident between the numerical simulation and the linear analytical solution. 5.2. Regular progressive wave propagation in a flume For the remaining cases, a combined mesh was utilised, with a Voronoi mesh occupying from xp ¼ 9.5 to 10.5, with a regular mesh either side. The whole domain
Table 1 Results for 2D asymmetric sloshing at x ¼ 60 m, t ¼ 9:2 s [28] Surface elevation, h (m)
u-velocity (m/s)
v-velocity (m/s)
23.803 23.860 23.815 23.759 23.820 23.803 23.720 23.811 23.814 23.810 23.854 23.867
22.456 22.480 22.423 22.411 22.417 22.417 22.480 22.411 22.453 22.240 22.449 22.408
20.363 20.560 20.577 20.602 20.580 20.572 20.690 20.550 20.579 20.560 20.583 20.583
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
71
Fig. 9. Time history of the free surface elevation at xp ¼ 0:8571: Solid line: Voronoi; dashed line: Voronoi-structured.
Fig. 10. Time history of free surface elevation at left wall of base-excited tank. Excitation frequency, v ¼ 1:1 v0 : Solid line: Voronoi; thin line: Voronoistructured; dashed line: analytical.
Fig. 11. Time history of free surface elevation at xp ¼ 10:0: Solid line: Voronoi-regular; dashed: sigma transform model.
72
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
Fig. 12. Space profile of free surface elevation at tp ¼ 135:9: Solid line: Voronoi-regular; dashed: sigma transform model.
occupied from xp ¼ 0 to xp ¼ 30: At the left hand end, the effect of a paddle wave generator was simulated, whereas at the right hand end a damping zone was used in conjunction with a radiation condition to prevent any wave reflection. Further details are given by Turnbull [26] and Turnbull et al., [29]. Fig. 11 shows a portion of the time history of the free surface elevation at xp ¼ 10; at the location where a submerged horizontal cylinder is to be placed later, for a wave of non-dimensional amplitude ap ¼ 0:001 and non-dimensional frequency vp ¼ 1:8495 (Section 5.3). The amplitude builds up slowly while tp , 20; and then monotonically increases to the steady state value which is approximately reached by tp ¼ 100: The wave is smallamplitude, and its temporal profile agrees almost exactly with linear wave theory. Fig. 12 shows the spatial profile of the free surface along the tank at time tp ¼ 135:9: Results from a sigma-transformed finite element model (Turnbull et al. [29]) are also plotted, and are in close agreement with the profiles predicted using the hybrid Voronoi-structured mesh. (It should be noted that the sigma-transformed finite element model utilises a linearly stretched mesh between the bed and free surface, and is limited to non-overturning non-breaking waves in the absence of submerged bodies.)
amplitude ap ¼ 0:001 and non-dimensional frequency vp ¼ 1:8495: The corresponding results for the case without the cylinder are also shown in Figs. 14 and 15. The time-dependent wave motions in the presence of the
5.3. Regular wave loading on a single submerged horizontal cylinder In this case, a submerged horizontal cylinder of nondimensional radius apc ¼ 0:06 is inserted at xp ¼ 10 with its centre 0.12 below the still water free surface. This case corresponds to Case E of Chaplin [30], for which experimental results were obtained. The overall mesh structure is the same as above, other than the region close to the cylinder, and is shown in Fig. 13. The mesh contains a total of 20,728 elements, of which 14,848 are in the regular regions and 5880 are in the Voronoi space. Fig. 14 shows a portion of the time-dependent free surface elevation history at xp ¼ 10; and Fig. 15 shows part of the spatial free surface distribution along the tank at tp ¼ 135:9 for waves of non-dimensional
Fig. 13. Voronoi section of mesh for progressive wave simulations.
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
73
Fig. 14. Time history of free surface elevation at cylinder ðxp ¼ 10:0Þ: Solid line: with cylinder; dashed line: without cylinder.
Fig. 15. Space profile of free surface elevation in region of cylinder ðtp ¼ 135:9Þ: Solid line: with cylinder; dashed line: without cylinder.
cylinder are nearly identical to those in its absence for xp , 9; but show significant difference in the vicinity of xp ¼ 10: This is in agreement with the theoretical result of Dean [31] who showed to first-order that the
reflection coefficient of a submerged cylinder is zero. The presence of the cylinder disturbs the local wave crests and troughs, causing the free surface to undulate slightly, as is evident in Fig. 15. There is a marked effect
Fig. 16. (a) Time history of force component Fa : Solid line: x-component; dashed line: z-component. (b) Final segment of time history of force component Fa : Solid line: x-component; dashed line: z-component.
74
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
Fig. 17. (a) Time history of force component Fb : Solid line: x-component; dashed line: z-component. (b) Final segment of time history of force component Fb : Solid line: x-component; dashed line: z-component.
on the wave phase, but the amplitude remains the same; waves that pass the cylinder experience a lag of about 158 in comparison with waves in the absence of the cylinder.
Table 2 Predicted, analytical (Ogilvie, 1963) and experimental (Chaplin, 1984) mean vertical force component Predicted
Ogilvie (1963)
Chaplin (1984)
F z =a2 0.193
F z =a2 0.190
F z =a2 0.187
The force on the submerged horizontal cylinder is calculated by integrating the hydrodynamic pressure 2 2 ! ›f r ›f ›f p ¼ 2r 2 ð22Þ þ ›t 2 ›x ›z Table 3 Predicted, analytical (Ogilvie, 1963) and experimental (Chaplin, 1984) first order Fourier oscillatory force components Predicted
Predicted
Ogilvie (1963)
Ogilvie (1963)
Chaplin (1984)
Chaplin (1984)
lFxð1Þ l=a 0.0581
lFzð1Þ l=a 0.0581
lFxð1Þ l=a 0.0576
lFzð1Þ l=a 0.0576
lFxð1Þ l=a 0.0577
lFzð1Þ l=a 0.0577
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77 Table 4 Predicted, analytical (Vada, 1987) and experimental (Chaplin, 1984) second order Fourier oscillatory force components Predicted
Predicted
Vada (1987)
Vada (1987)
Chaplin (1984)
Chaplin (1984)
lFxð2Þ l=a2 0.293
lFzð2Þ l=a2 0.307
lFxð2Þ l=a2 0.31
lFzð2Þ l=a2 0.31
lFxð2Þ l=a2 0.39
lFzð2Þ l=a2 0.39
Fig. 18. Voronoi section of mesh for progressive wave simulations with two cylinders.
to give F ¼ F a þ Fb
ð23Þ
where Fa ¼ 2 r
ð Sc
›f ›t
nx nz
oscillation by a non-dimensional time equal to about 100. The oscillations are predominantly at the wave frequency. Fb exhibits a significant mean vertical force component, unlike Fa : Ogilvie’s [32] np parameter for this case is np ¼ p2 p v =g ¼ 3:42065; and so nac ¼ 0:205 and 2ne ¼ 0:821 (where e is the distance between the axis of the cylinder and the near free surface). The mean vertical component of F given in Table 2 matches Ogilvie’s [32] analytical and Chaplin’s [30] values. Fourier analysis of segments of the time histories of the x and z components of F (after steady state had been reached) was used to determine the first and second order Fourier oscillatory force components. The predicted non-dimensional first order Fourier oscillatory components on the cylinder are given in Table 3, and are in reasonable agreement with the results from Ogilvie [32] and Chaplin [30]. Almost the same results were achieved using a mesh with approximately half the number of nodes. Table 4 presents the predicted non-dimensional second order Fourier oscillatory force components, compared with results from the analysis by Vada [33] and the experiments of Chaplin [30]. The overall agreement suggests that the finite element hybrid mesh numerical model should be suitable for solving a range of non-linear wave structure interaction problems. 5.4. Regular wave interaction with a pair of submerged horizontal cylinders
! dSc
ð24Þ
and
! ! rð ›f 2 ›f 2 nx F ¼2 dSc þ 2 Sc ›x ›z nz b
75
ð25Þ
where the integrals are with respect to the distance around the cylinder surface, Sc : Figs. 16 and 17 show the time-dependent behaviour of the x- and z- components of the force vectors, Fa ; and Fb : The amplitudes of the force components are initially zero, then increase after a non-dimensional time equal to 20 approximately, and eventually build up to a steady state
Fig. 18 shows the physical mesh obtained for regular waves interacting with a pair of submerged horizontal cylinders, each of non-dimensional radius apc ¼ 0:06; at non-dimensional time tp ¼ 135:9 (corresponding to the previous case). The upwave cylinder is centred at xp ¼ 10 with its centre 0.12 below the still water free surface. The centres of the two cylinders are located a non-dimensional distance 0.91842 apart (equal to half the wavelength of undisturbed progressive waves). Fig. 19 shows the spatial variation in free surface elevation in the region of the pair of submerged cylinders at tp ¼ 135:9: The free surface upwave of the first cylinder is identical to the undisturbed regular wave profile and that ahead of the single submerged
Fig. 19. Space profile of free surface elevation in region of cylinders ðtp ¼ 135:9Þ: Solid line: two cylinders; dashed line: one cylinder; dotted line: no cylinder.
76
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
horizontal cylinder in Fig. 15. The effect of the pair of cylinders is to cause a further phase shift in the wave profile after the regular waves have passed over the cylinders. This effect is significantly larger than for the single submerged cylinder case.
6. Conclusions This paper has described a 2D finite element scheme for predicting wave-induced forces on a submerged horizontal body. The finite element mesh comprises a combination of Voronoi unstructured and regular structured mesh regions, the former being used in the vicinity of the cylinder. The hybrid mesh is adaptive and uses a Lagrange – Eulerian scheme to track the vertical locations of free surface nodes. The free sloshing validation tests indicate that almost identical agreement is achieved using entirely Voronoi and combined Voronoi-structured grids. The results confirm the influence of non-linear effects on steep sloshing reported by Wu and Eatock Taylor [7] and Chern et al. [27]. The present hybrid mesh scheme predicts velocity components under an asymmetric sloshing wave that closely match corresponding results from other codes, reported by Nestega˚rd [28]. Progressive regular waves are simulated in the absence of a submerged body, and close agreement obtained with fully non-linear simulations using an alternative sigmatransformed finite element model developed by Turnbull et al. [29]. Mean, first order oscillatory and second order oscillatory forces on a submerged horizontal cylinder are calculated, and found to be in close agreement with analytical results from Ogilvie [32] and Vada [33] as well as experimental measurements by Chaplin [30]. The effect of a second submerged horizontal cylinder is to accentuate the phase shift experienced by regular waves as they pass over a single submerged cylinder.
Acknowledgements The authors gratefully acknowledge the financial support provided by the UK Engineering and Physical Sciences Research Council (EPSRC) Grant GR/M56401. The authors would also like to thank Professor G.X. Wu of University College London who provided the original finite element code.
Appendix A
Notation list A Aij
intermediate quantity for calculation of l coefficients of matrix ½A
½A B Bi {B} F Fa Fb F z lFxð1Þ l lFxð2Þ l N Ni P R S Sc S1 S2 a ap ac apc b1 b2 c d e f1 f2 g gp h l la lb lk n nv nx nz p q q1 q2 q3 r s t
matrix for calculation of {f} intermediate quantity for calculation of l coefficients of vector {B} vector for calculation of {f} total force vector on cylinder component of force vector on cylinder component of force vector on cylinder mean vertical force on cylinder first order horizontal oscillatory force component on cylinder second order horizontal oscillatory force component on cylinder total number of nodes in the finite element mesh shape function a point in the domain a prescribed radial distance measured outward from the boundary of the circle fluid boundary distance around cylinder surface section of boundary where f is specified section of boundary where ›f=›n is specified initial amplitude for symmetrical sloshing case; amplitude of progressive wave non-dimensional amplitude ðap ¼ a=hÞ radius of cylinder non-dimensional radius of cylinder ðapc ¼ ac =hÞ length of line segment at bottom of domain length of line segment on inner circle wave celerity number of line segments on the lateral outer boundaries distance between cylinder centre and undisturbed free surface intermediate quantity for calculation of l intermediate quantity for calculation of l acceleration due to gravity non-dimensionalised gravity ðgp ¼ 1Þ: water depth; distance between lower and upper boundaries at P length of line segment length of line segment length of line segment vector used in least squared velocity calculation normal to the boundary number of nodes in the vertical used for velocity calculation x-component of normal vector z-component of normal vector pressure velocity normal to the domain mesh quality parameter mesh quality parameter mesh quality parameter radial distance from boundary of circle length of line segment at the top of the domain time
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
tp u w x y z V a as
b bs gs h n np r f fi {f} v vp
pffiffiffiffi non-dimensional time, tp ¼ t g=h horizontal velocity component vertical velocity component horizontal distance distance measured upwards from the lower outer boundary vertical distance fluid domain real number used for vertical velocity calculation parameter for definition of initial free surface for asymmetric sloshing case real number used for vertical velocity calculation parameter for definition of initial free surface for asymmetric sloshing case parameter for definition of initial free surface for asymmetric sloshing case free surface elevation Ogilvie parameter non-dimensional Ogilvie parameter, np ¼ nh fluid density velocity potential nodal velocity potential vector of nodal velocity potentials wave frequency for travelling wave pffiffiffiffi case non-dimensional time, vp ¼ v h=g
References [1] Chaplin JR, Rainey RCT, Yemm RW. Ringing of a vertical cylinder in waves. J Fluid Mech 1997;350:119 –47. [2] Longuet-Higgins MS, Cokelet ED. The deformation of steep surface waves on water, I: a numerical method of computation. Proc R Soc Lond, A 1976;350:1–25. [3] Vinje T, Brevig P. Non linear ship motions. Proceedings of the Third International Symposium Numerical Ship Hydrodynamics, Paris (Chairmen: Dern JC, Haussling HJ). Bassin d’Essais des Carenes, France, 1981. p. 257 –68. [4] Cao Y, Schultz WW, Beck RF. Three dimensional unsteady computations of non-linear waves caused by underwater disturbance. Proceedings of the 18th Symposium on Naval Hydrodynamics, Ann Arbor, MI, 1990. p. 417–27. [5] Celebi MS, Kim MH. Non linear wave body interactions in a numerical wave tank. Proceedings of the 12th International Workshop on Water Waves and Floating Bodies, France, 1997. [6] Wang P, Yao Y, Tulin MP. An efficient numerical tank for non-linear water waves, based on the multi-subdomain approach with BEM. Int J Numer Meth Fluids 1995;20:1315 –36. [7] Wu GX, Eatock Taylor R. Finite element analysis of twodimensional non-linear transient water waves. Appl Ocean Res 1994;16:363 –72. [8] Wu GX, Eatock Taylor R. Time stepping solutions of the two dimensional non-linear wave radiation problem. Ocean Engng 1995; 22(8):785–98. [9] George PL. Automatic mesh generation. Paris: Wiley; 1991.
77
[10] Yerry MA, Shephard MS. A modified quadtree approach to finite element mesh generation. IEEE Comput Graph Appl 1983;3(1): 39–46. [11] Yerry MA, Shephard MS. Automatic 3D mesh generation by the modified octree technique. Int J Numer Meth Engng 1984;20(11): 1965– 90. [12] Greaves DM, Borthwick AGL. Hierarchical tree-based finite element mesh generation. Int J Numer Meth Engng 1999;45:447–71. [13] Greaves DM, Borthwick AGL, Wu GX, Eatock Taylor R. A moving boundary finite element method for fully non-linear wave simulations. J Ship Res 1997;3:181– 94. [14] Zhu G, Borthwick AGL, Eatock Taylor R. A finite element model of interaction between viscous free surface waves and submerged cylinders. Ocean Engng 2001;28(8):989–1008. [15] Lo SH. A new mesh generation for arbitrary planar domains. Int J Numer Meth Engng 1985;21:1403. [16] Peraire J, Vahdati M, Morgan K, Zienkiewicz OC. Adaptive remeshing for compressible flow computations. J Comput Phys 1987;72(2):449–66. [17] Lee CK, Hobbs RE. Automatic adaptive finite element mesh generation over arbitrary two-dimensional domain using advancing front technique. Comput Struct 1999;71(1):9–34. [18] Chan CT. Computation of flows by the finite volume method as applied to unstructured meshes. PhD Thesis, University of London (Imperial College), 1997. [19] Lawson CL. Generation of a triangular grid with application to contour plotting. California Institute of Technology, Jet Propulsion Lab, Technical Memo 299, 1972. [20] Bowyer A. Computing Dirichlet tessellations. Comput J 1981;24:162. [21] Watson DF. Computing the n-dimensional Delaunay tessellation with application to Voronoi polytropes. Comput J 1981;24(2):167–72. [22] Mavriplis DJ. Adaptive mesh generation for viscous flows using Delaunay triangulation. J Comput Phys 1990;90(2):271–91. [23] Borouchaki H, George PL. Aspects of 2D Delaunay mesh generation. Int J Numer Meth Engng 1997;40(11):1957– 75. [24] Rivara MC, Inostroza P. Using longest side bisection techniques for the automatic refinement of Delauney triangulations. Int J Numer Meth Engng 1997;40(4):581 –97. [25] Xu X, Pain CC, Goddard AJH, de Oliveira CRE. An automatic adaptive meshing technique for Delaunay triangulations. Comput Meth Appl Mech Engng 1998;161(3–4):297– 303. [26] Turnbull MS. The numerical modelling of steep waves interacting with structures. DPhil Thesis, University of Oxford, UK, 1999. [27] Chern MJ, Borthwick AGL, Eatock Taylor R. A pseudospectral stransformation model of 2-D non-linear waves. J Fluids Struct 1999; 13:607–30. [28] Nestega˚rd A. Comparative study of fully non-linear wave simulation programs. DNV Report No 94-2041, 1994. [29] Turnbull MS, Borthwick AGL, Eatock Taylor R. Numerical wave tank based on a s-transformed finite element viscous flow solver. Int J Numer Meth Fluids, in press. [30] Chaplin JR. Nonlinear forces on a horizontal cylinder beneath waves. J Fluid Mech 1984;147:449 –64. [31] Dean WR. On the reflection of surface waves by a submerged cylinder. Proc Camb Phil Soc 1948;44:483– 91. [32] Ogilvie TF. First- and second-order forces on a cylinder submerged under a free surface. J Fluid Mech 1963;16:451–72. [33] Vada T. A numerical solution of the second-order wave-diffraction problem for a submerged cylinder of arbitrary shape. J Fluid Mech 1987;174:23 –37.
Wave –structure interaction using coupled structured –unstructured finite element meshes M.S. Turnbulla, A.G.L. Borthwickb,*, R. Eatock Taylorb a HR Wallingford Ltd, Howbery Park, Wallingford, Oxon OX10 8BA, UK Department of Engineering Science, University of Oxford, Parks Road, Oxford OX1 3PJ, UK
b
Received 31 August 2002; accepted 3 April 2003
Abstract The interaction of inviscid gravity waves with submerged fixed horizontal structures is modelled in two-dimensions using a finite element numerical wave tank. An adaptive hybrid coupled mesh is utilised that tracks the free surface through vertical movement of the free surface nodes in a Lagrange – Eulerian scheme. The hybrid mesh consists of a combination of structured and Voronoi unstructured meshes, with the submerged structure located in the Voronoi mesh region. Validation tests include free and forced sloshing in a rectangular tank, regular progressive wave propagation in a flume, and regular wave loading on a horizontal cylinder. The results are found to be in close agreement with analytical potential theory solutions for waves of small amplitude. Non-linear effects are noted for steeper waves. The wave-induced force components on the horizontal cylinder match the expected results from Ogilvie’s [J Fluid Mech 16 (1963) 451] linear theory and Vada’s [J Fluid Mech 174 (1987) 23] second order model. Wave interactions with a pair of submerged horizontal cylinders spaced at half the wavelength of undisturbed regular waves are examined. q 2003 Elsevier Ltd. All rights reserved. Keywords: Finite elements; Wave tank; Cylinder; Forces
1. Introduction Steep free surface waves give rise to highly non-linear forces on fixed offshore structures, and can cause the phenomenon of ringing, which has been the subject of considerable recent research (see e.g. Chaplin et al. [1]). Theoretical modelling of the interaction of steep gravity waves with complicated offshore structures is a challenging task because both the kinematic and dynamic free surface boundary conditions are non-linear and the flow domain has an irregular outer boundary surrounding an interior obstacle or obstacles. To deal with the awkward geometries encountered by offshore engineers in practice, it is necessary to develop efficient and accurate numerical models that can cope with steep wave non-linearity and the presence of submerged obstacles. Longuet-Higgins and Cokelet [2] assumed periodicity and hence simulated steep Stokes waves in infinite depth * Corresponding author. Tel.: þ 44-1865-273047; fax: þ 44-1865273010. E-mail address: [email protected] (A.G.L. Borthwick). 0141-1187/03/$ - see front matter q 2003 Elsevier Ltd. All rights reserved. doi:10.1016/S0141-1187(03)00032-4
using a boundary element method. Vinje and Brevig [3] extended this approach to finite depth with floating bodies. Cao et al. [4] used a time stepping technique with the Desingularised Boundary Integral Equation Method (DBIEM) to study various non-linear free surface problems including the non-linear waves generated by free surface pressure disturbances. The same method was used by Celebi et al. [5] to investigate non-linear wave interactions with a stationary vertical cylinder. Wang et al. [6] developed a multi-subdomain approach that significantly improved the performance of their 2D boundary element numerical wave tank. Wu and Eatock Taylor [7,8] developed a finite element numerical wave tank using structured meshes that stretched between the bed and liquid free surface. Wu and Eatock Taylor demonstrated that the finite element method could be as efficient computationally as the boundary element method for inviscid gravity wave calculations. The present paper implements the finite element formulation derived by Wu and Eatock Taylor [7,8] on a combination of unstructured Voronoi and structured meshes.
64
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
George [9] has reviewed the generation of finite element meshes, and broadly classifies the different approaches as enumerative, functional mapping, numerical mapping, hierarchical, block decomposition, advancingfront and Voronoi –Delaunay’s construction. Of the foregoing, the hierarchical and latter two unstructured approaches are becoming popular due to their inherent robustness, efficiency and ease of adaptation. For example, Yerry and Shephard [10,11] and Greaves and Borthwick [12] have developed 2D quadtree and threedimensional octree hierarchical mesh generators. Greaves et al. [13] solved free surface sloshing in a tank on a finite element quadtree mesh. Zhu et al. [14] have developed a viscous flow numerical wave tank using triangularised quadtree finite elements where the free surface has been tracked in a Lagrangian manner. Lo [15] proposed the advancing front technique that was extended by Peraire et al. [16]. Lee and Hobbs [17] describe an advancing front method for the creation of adaptive finite element meshes of triangles and quadrilaterals. Chan [18] solved a variety of free surface flow problems including water sloshing in a tank, solitary wave propagation, and flow over a submerged hydrofoil using a finite volume solution of the Navier – Stokes equations on unstructured 2D meshes. In an alternative approach to unstructured mesh generation, Lawson [19] proposed an algorithm based on edge swapping to create Delaunay triangulations. A more efficient Voronoi –Delaunay algorithm was proposed by Bowyer [20] and Watson [21]. Mavriplis [22] used Voronoi meshes adaptively for viscous flow in two dimensions. Further work on the generation of Voronoi meshes was carried out by Borouchaki and George [23]. Rivara and Inostroza [24] compared two longest side bisection techniques for the automatic refinement of Delaunay triangulations. Xu et al. [25] presented an automatic coarsening procedure. As far as the present authors are aware, unstructured Voronoi meshes have not yet been applied to an inviscid numerical wave tank. This paper aims to use a flexible mesh generation strategy with a Voronoi unstructured mesh applied in the vicinity of a submerged cylinder and structured meshes elsewhere in regions away from cylinder. The validation tests have been selected to check that the model properly simulates free surface waves of increasing steepness, and gives accurate results for a submerged horizontal cylinder under regular wave loading. The technique, though limited to two-dimensions in this paper, is readily extendible to three-dimensions and to more complicated structural configurations typical of offshore structures.
[20]. An initial mesh is created, consisting of a square divided diagonally into two triangular elements and having a total of four nodes within which the computational domain fits. Each of these elements has an associated circumcircle that intersects its nodes. Nodes are inserted that delineate the boundary of the flow domain. Each new node is checked, and if it lies within a circumcircle, then the element is removed, and new elements created to occupy the space left behind by connecting the surrounding nodes to the node under consideration. Following this process, any elements created outside the flow domain boundary are deleted. A recursive procedure is then invoked. Each element is checked in turn. If its area exceeds a prescribed criterion a new node is inserted such that it makes an equilateral triangle when joined to the shortest side of the element. Overlapping nodes are coalesced. This procedure is repeated until no more changes occur to the mesh. The element limiting area criterion is as follows. For a uniform mesh, the maximum area of any element is set to a constant value throughout the domain. For a variable density mesh, an analytical formula gives the desired length of a line segment, l; anywhere in the domain. Consider the domain shown in Fig. 1, which is to have uniform node spacing on its upper and lower outer boundaries, and variable spacing on its lateral outer boundaries. This domain could represent a submerged horizontal cylinder under waves, where a fine mesh would be required near the free surface and close to the cylinder. Point P within the domain is located a radial distance r from the inner circular boundary and a distance y measured upwards from the lower outer boundary. The distance between the upper and lower outer boundaries at P is given by h: R is a prescribed radial distance measured outward from the boundary of the circle. Let s; b1 and b2 be lengths of line segments at the upper outer, bottom outer, and inner cylinder boundaries, respectively. Also, d is the number of line segments on the lateral outer boundaries. If b2 is smaller than either s or b1 ; then the line segment l is set equal to the smaller of la and lb ; otherwise, l is set equal to the larger of la and lb ; where f2 A b1 b la ¼ s and lb ¼ b1 2 ð1Þ s b1
2. Voronoi mesh generator The Voronoi method used herein generates a mesh using the node insertion technique proposed by Bowyer
Fig. 1. Definition sketch.
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
65
Fig. 2. Voronoi meshes generated around a circle in a square domain.
where A¼
log10 ðf1B þ 1Þ ; log10 ð2Þ
f1 ¼ 1 2 ðy=hÞ; B ¼ ðs=b1 Þðd=ðd21ÞÞ þ 1 and f2 ¼ 1 2 ðr=RÞ: The upper limit on the size of an element is given by q1 l2 where q1 is a prescribed parameter. During initial mesh generation, the following rules are invoked to improve mesh quality. † Before a new node is incorporated into the mesh, it is checked whether it lies closer than a prescribed distance q2 l to an existing node. If it is closer than q2 l; the hypothetical node is not inserted. † If the smallest side of an element is larger than q3 l then no node is inserted into that element, regardless of its size. This rule prevents nodes being created too far from the boundary initially. The values of q1 ; q2 and q3 can all be adjusted in order to optimize the quality of the mesh produced. After the initial mesh has been created, smoothing in four iterative sweeps is then performed to improve mesh quality.
Each node is repositioned so that it lies approximately at the average location of the nodes to which it is connected. In the test cases considered herein, the mesh parameters, q1 ; q2 and q3 ; have been set to 0.6, 0.5 and 10, respectively, after extensive numerical tests (Turnbull, [26]). Fig. 2 shows the results of the Voronoi method applied to a square domain containing a circle. The analytical formulas (1) were used to specify the mesh density, with R equal to twice the circle radius. For the uniform square, s ¼ b1 : In Fig. 2(a), the mesh is uniform, with 32 segments on the circumference of the circle for b2 ¼ s: The mesh in Fig. 2(b) is finer in the region around the circle which has 64 segments for b2 ¼ s=2: The converse occurs in Fig. 2(c) where there are only 16 segments around the circle for b2 ¼ 2s and the mesh is coarser in the vicinity of the circle. In Fig. 2(d), the mesh is uniform, but the circle is not at the centre of the domain. In all these cases, the Voronoi method has produced elements of acceptable quality, with the correct mesh density where specified. Fig. 3 shows several examples of Voronoi meshes for a domain of width equal to twice its height, and containing a circular inner boundary. Fig. 3(a) presents the mesh
66
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
the kind of mesh that would be used for finite element simulations, with finer mesh density close to the free surface and near the cylinder, and coarser mesh density towards the bottom of the domain, where velocities are expected to be smaller.
3. Governing equations Consider an incompressible, inviscid liquid undergoing irrotational flow in two-dimensions. Mass conservation may be expressed by 72 f ¼ 0;
ð2Þ
where f is the velocity potential, and 7 is the gradient vector. The non-linear free surface kinematic boundary condition is
›f ›h ›f ›h ¼ þ ; ›z ›t ›x ›x
ð3Þ
where h is the elevation of the free surface vertically above still water level, x is distance measured horizontally along the free surface, z is distance vertically above the free surface and t is time. The non-linear dynamic free surface boundary condition is
›f 1 þ 7f7f ¼ 0 ð4Þ 2 ›t in which g is gravitational acceleration. At solid boundaries, such as the fixed bed, lateral walls and submerged obstacle surfaces, gh þ
›f ¼ q; ð5Þ ›n where q is the normal velocity at the boundary and n is normal to the boundary and is directed out of the fluid domain. At open boundaries, the Sommerfeld radiation condition is ›f ›f ¼ 2c ›t ›x where c is the wave celerity.
ð6Þ
4. Finite element formulation Fig. 3. Meshes for domains with circular inner boundary and varying mesh density.
generated when the upper boundary is straight, whereas Fig. 3(b) –(d) depict meshes where the upper boundary is sinusoidal. These meshes are representative of those that would be produced during the finite element simulation of waves passing over a submerged horizontal cylinder. In each case the cylinder has 32 segments and b2 ¼ 0:5: The number of segments vertically, d is 16. In Fig. 3(a) and (b), s ¼ b1 ¼ 1:0; in Fig. 3(c) s ¼ 0:5 and b1 ¼ 1:0; in Fig. 3(d), s ¼ 0:5 and b1 ¼ 2:0: Fig. 3(d) shows
4.1. Discretisation The continuous velocity potential, f; is expressed in terms of nodal potentials fi and corresponding linear shape functions Ni ðx; zÞ as
f¼
N X
fi Ni ðx; zÞ
ð7Þ
i¼1
where N is the total number of nodes in the mesh. In the finite element mesh, each element is triangular with three corner nodes.
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
At each time step, a solution of the Laplace equation is obtained using the Galerkin technique. Wu and Eatock Taylor [8] integrated the product of Eq. (2) and the nodal shape functions over the domain, V; to obtain ð ð ð ›f dS 2 Ni 72 f dV ¼ Ni 7Ni 7f dV ¼ 0; ð8Þ ›n V S V
system: 0 P X B lkx lkx B B k¼1 B B P BX @ lkx lkz k¼1
where S is the boundary of V and i ¼ 1; …; N: Wu and Eatock Taylor then substituted f with its nodal approximation and replaced ›f=›n with q to give
67
8 1 P X 8 9 > ›f > ›f > > C> lkx k > > > C> < = < C ›x k¼1 k¼1 ›l C ¼ C > > P P C> X X > ›f : ›f > ; > > > lkz lkz A ›z lkz k > : ›l P X
lkx lkz
k¼1
k¼1
9 > > > > > = > > > > > ;
ð13Þ
where S1 and S2 are the surfaces where f and ›f=›n are specified. The terms on the right hand side of the equation are evaluated using the boundary conditions. In matrix form, the system of equations is,
in which lkx and lkz are the x and z components of lk divided by its magnitude of lk : The method is suitable for any kind of unstructured mesh. A higher order calculation is used to estimate the vertical velocity component at nodes on the free surface. The method requires that each free surface node and a prescribed number of nodes below it are situated along a vertical line. This restriction was incorporated within the Voronoi mesh generator. Implementation is as follows. Let z be the vertical distance upwards from the mean free surface. For a free surface node with elevation z1 ; the three nodes vertically below it are located at elevations, whose positions are given by
½A{f} ¼ {B}
z2 ¼ z1 2 Dz
ð R
7Ni
¼2
N X
fj 7Nj dRljS1
j¼1
ð R
7Ni
N X
fj 7Nj dRlj[S1 þ
j¼1
ð S2
Ni q dS
ð9Þ
ð10Þ
where {f} is a vector containing nodal velocity potentials and ½A and {B} have the following coefficients ð 7Ni 7Nj dR or Aij ¼ R
ð11Þ
Aij ¼ 1
if i ¼ j and i [ S1
Aij ¼ 0
if ði [ S1 or j [ S1 Þ and i – j
and Bi ¼ 2
N X j¼1
fj
ð
R
7Ni ·7Nj dR lj[s1 þ
ð S2
Ni q dS ð12Þ
or Bi ¼ fi if i [ S1 : In the above, fi is specified on the free surface ði [ S1 Þ as so does not need to be calculated there. The relations following the word ‘or’ in Eqs. (11) and (12) allow the known values to be transferred from the vector to the vector {f}: Gaussian elimination is used to solve matrix Eq. (10).
z3 ¼ z1 2 aDz
z4 ¼ z1 2 bDz
ð14Þ
where Dz is the spacing between the uppermost pair of nodes, and a and b are real numbers. The velocity potential is expressed as a power series, z 2 z1 z 2 z1 2 f ¼ A0 þ A1 þ· · · þ A2 Dz Dz z 2 z1 nv 21 ; ð15Þ þ An21 Dz in which nv is the number of nodes in the vertical that are utilised in estimating the vertical velocity component. Eq. (15) can be expressed as nv simultaneous equations by inserting values for f and z at the free surface node and the nv 2 1 nodes vertically beneath it. The vertical velocity component at the free surface is calculated by differentiating Eq. (15) and then taking z to be equal to z1 : When nv ¼ 3 and the nodes are unevenly spaced in the vertical, the vertical velocity component at the free surface becomes: w¼
›f ða2 2 1Þf1 2 a2 f2 þ f3 ¼ : ›z aða 2 1ÞDz
4.2. Surface velocity calculation
For a ¼ 2; the nodes are evenly spaced, and
The least-squares method used to estimate the horizontal free surface velocity component is implemented as follows. Consider an arbitrary free surface node that is connected to p underlying nodes in the finite element mesh. Let lk denote the position vector of the k th node relative to the arbitrary free surface node and connected to it. The potential gradient vector is obtained from the following least-squares equation
w¼
3f1 2 4f2 þ f3 : 2Dz
ð16Þ
ð17Þ
For evenly spaced nodes when nv ¼ 5; the following higher order formula results: w¼
25f1 2 48f2 þ 36f3 2 16f4 þ 3f5 12Dz
ð18Þ
68
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
4.3. Updating the free surface
5. Results
A mixed Euler –Lagrange approach is used where the free surface nodes are restricted to move solely in the vertical direction. The free surface elevation is updated using a fourth-order Runge – Kutta integration of the free surface kinematic boundary condition, expressed as
5.1. Free and forced sloshing in a rectangular tank
›h ›h ¼w2u : ð19Þ ›t ›x Similarly, the free surface potential is updated using the dynamic free surface boundary condition: ›f 1 ›h ¼ 2gh 2 ðu2 þ w2 Þ þ w : 2 ›t ›t
ð20Þ
Fig. 4 shows initial entirely Voronoi and combined Voronoi-structured meshes used to check mesh convergence for free sloshing waves of amplitude to depth ratio, a=h ¼ 0:1: Three levels of coarseness were considered: coarse, with 32 elements along the upper free surface boundary, and 16 vertically along the lateral walls; medium, with 64 free surface and 32 vertical wall elements; and fine, with 128 free surface and 64 lateral wall elements. The structured mesh is created using a linear stretching between the free surface and bed (Wu and Eatock Taylor, [7]). The Voronoi meshes are symmetric about the vertical centre of the tank, and contain greatest cell variations at about mid-depth. Fig. 5
Fig. 4. Initial meshes: ap ¼ 0:1:
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
69
Fig. 5. Time history of the free surface elevation at the centre of the tank. Free sloshing: ap ¼ 0:1: Voronoi mesh. Dotted line: 32 by 16; dashed line: 64 by 32; solid line: 128 by 64.
Fig. 6. Time history of the free surface elevation at the centre of the tank. Free sloshing: ap ¼ 0:001: Solid line: Voronoi; thin line: Voronoi-structured; dashed line: second order analytical.
Fig. 7. Time history of the free surface elevation at the centre of the tank. Free sloshing: ap ¼ 0:1: Solid line: Voronoi; thin line: Voronoi-structured; dashed line: second order analytical.
70
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
shows the motions of the free surface elevation at the tank centre simulated using the three combined Voronoistructured meshes. The 32 by 16 mesh gives slightly diffusive results with the wave crests tending to be underpredicted, in comparison with the corresponding results obtained using the finer meshes. Simulations using the 64 by 32 and 128 by 64 meshes are in almost exact agreement, indicating that mesh convergence has been achieved. It should pffiffiffiffibe noted that the non-dimensional time step, Dtp ¼ Dt g=h; was set to be Dtp ¼ 0:1; 0.05 and 0.025 for each mesh, respectively, in order of its fineness. The entirely Voronoi mesh scheme gave results that were indistinguishable from those using the combined mesh approach. Mesh and time step convergence occur for a 64 by 32 fully Voronoi mesh with non-dimensional time step, Dtp ¼ 0:05: Figs. 6 and 7 demonstrate the increasing non-linear behaviour of the free surface as the relative wave amplitude ap ¼ a=h is increased from 0.001 to 0.1. For small amplitude free sloshing where ap ¼ 0:001; the free surface motions are identical to the analytical solution from second order potential theory (Wu and Eatock Taylor, [7]) which itself is very close in this case to linear theory. For larger amplitude sloshing where ap ¼ 0:1; the free surface waves are of longer period and have higher extreme maxima than those from the analytical solution. The non-linear free surface motions are in excellent agreement with those obtained using alternative (fully non-linear) pseudospectral and finite element numerical approaches by Chern et al. [27] and Wu and Eatock Taylor [7], respectively. An asymmetric sloshing wave is next considered. It was first utilised in a DnV study (Nestega˚rd, [28]) that compared results from several industrial and academic numerical codes for fully non-linear free surface sloshing. The tank has length 160 m and still water depth 70 m. At time t ¼ 0; the free surface elevation is given by " 2 # " 2 # x x exp 2 ð21Þ hðx; 0Þ ¼ as 1 2 bs gs where as ¼ 12 m, bs ¼ 53 m and gs ¼ 76 m. Fig. 8 shows the initial combined Voronoi-structured mesh. Table 1 presents the free surface elevation, and surface velocity components at x ¼ 60 m and t ¼ 9.2 s obtained by the participants in the DnV study and using the present method with unstructured and hybrid meshes. It is clear that the present scheme gives results in close agreement with those of the majority of the participants. Fig. 9 shows the timedependent free surface motions at x ¼ 60 m (corresponding to xp ¼ 0:8571 where xp ¼ x=h). After non-dimensional time tp ¼ 10; the most significant component appears to be third-order. Horizontally forced sloshing in a base-excited tank has also been considered. Fig. 10 shows the bursting characteristics of the free surface motions at the left wall when the excitation frequency v is close to the lowest natural sloshing
Fig. 8. Initial meshes for asymmetric sloshing wave.
frequency v0 : The sloshing motions exhibit a beating effect, caused by interactions between the natural sloshing frequency of the liquid and the forcing frequency. There is close agreement evident between the numerical simulation and the linear analytical solution. 5.2. Regular progressive wave propagation in a flume For the remaining cases, a combined mesh was utilised, with a Voronoi mesh occupying from xp ¼ 9.5 to 10.5, with a regular mesh either side. The whole domain
Table 1 Results for 2D asymmetric sloshing at x ¼ 60 m, t ¼ 9:2 s [28] Surface elevation, h (m)
u-velocity (m/s)
v-velocity (m/s)
23.803 23.860 23.815 23.759 23.820 23.803 23.720 23.811 23.814 23.810 23.854 23.867
22.456 22.480 22.423 22.411 22.417 22.417 22.480 22.411 22.453 22.240 22.449 22.408
20.363 20.560 20.577 20.602 20.580 20.572 20.690 20.550 20.579 20.560 20.583 20.583
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
71
Fig. 9. Time history of the free surface elevation at xp ¼ 0:8571: Solid line: Voronoi; dashed line: Voronoi-structured.
Fig. 10. Time history of free surface elevation at left wall of base-excited tank. Excitation frequency, v ¼ 1:1 v0 : Solid line: Voronoi; thin line: Voronoistructured; dashed line: analytical.
Fig. 11. Time history of free surface elevation at xp ¼ 10:0: Solid line: Voronoi-regular; dashed: sigma transform model.
72
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
Fig. 12. Space profile of free surface elevation at tp ¼ 135:9: Solid line: Voronoi-regular; dashed: sigma transform model.
occupied from xp ¼ 0 to xp ¼ 30: At the left hand end, the effect of a paddle wave generator was simulated, whereas at the right hand end a damping zone was used in conjunction with a radiation condition to prevent any wave reflection. Further details are given by Turnbull [26] and Turnbull et al., [29]. Fig. 11 shows a portion of the time history of the free surface elevation at xp ¼ 10; at the location where a submerged horizontal cylinder is to be placed later, for a wave of non-dimensional amplitude ap ¼ 0:001 and non-dimensional frequency vp ¼ 1:8495 (Section 5.3). The amplitude builds up slowly while tp , 20; and then monotonically increases to the steady state value which is approximately reached by tp ¼ 100: The wave is smallamplitude, and its temporal profile agrees almost exactly with linear wave theory. Fig. 12 shows the spatial profile of the free surface along the tank at time tp ¼ 135:9: Results from a sigma-transformed finite element model (Turnbull et al. [29]) are also plotted, and are in close agreement with the profiles predicted using the hybrid Voronoi-structured mesh. (It should be noted that the sigma-transformed finite element model utilises a linearly stretched mesh between the bed and free surface, and is limited to non-overturning non-breaking waves in the absence of submerged bodies.)
amplitude ap ¼ 0:001 and non-dimensional frequency vp ¼ 1:8495: The corresponding results for the case without the cylinder are also shown in Figs. 14 and 15. The time-dependent wave motions in the presence of the
5.3. Regular wave loading on a single submerged horizontal cylinder In this case, a submerged horizontal cylinder of nondimensional radius apc ¼ 0:06 is inserted at xp ¼ 10 with its centre 0.12 below the still water free surface. This case corresponds to Case E of Chaplin [30], for which experimental results were obtained. The overall mesh structure is the same as above, other than the region close to the cylinder, and is shown in Fig. 13. The mesh contains a total of 20,728 elements, of which 14,848 are in the regular regions and 5880 are in the Voronoi space. Fig. 14 shows a portion of the time-dependent free surface elevation history at xp ¼ 10; and Fig. 15 shows part of the spatial free surface distribution along the tank at tp ¼ 135:9 for waves of non-dimensional
Fig. 13. Voronoi section of mesh for progressive wave simulations.
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
73
Fig. 14. Time history of free surface elevation at cylinder ðxp ¼ 10:0Þ: Solid line: with cylinder; dashed line: without cylinder.
Fig. 15. Space profile of free surface elevation in region of cylinder ðtp ¼ 135:9Þ: Solid line: with cylinder; dashed line: without cylinder.
cylinder are nearly identical to those in its absence for xp , 9; but show significant difference in the vicinity of xp ¼ 10: This is in agreement with the theoretical result of Dean [31] who showed to first-order that the
reflection coefficient of a submerged cylinder is zero. The presence of the cylinder disturbs the local wave crests and troughs, causing the free surface to undulate slightly, as is evident in Fig. 15. There is a marked effect
Fig. 16. (a) Time history of force component Fa : Solid line: x-component; dashed line: z-component. (b) Final segment of time history of force component Fa : Solid line: x-component; dashed line: z-component.
74
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
Fig. 17. (a) Time history of force component Fb : Solid line: x-component; dashed line: z-component. (b) Final segment of time history of force component Fb : Solid line: x-component; dashed line: z-component.
on the wave phase, but the amplitude remains the same; waves that pass the cylinder experience a lag of about 158 in comparison with waves in the absence of the cylinder.
Table 2 Predicted, analytical (Ogilvie, 1963) and experimental (Chaplin, 1984) mean vertical force component Predicted
Ogilvie (1963)
Chaplin (1984)
F z =a2 0.193
F z =a2 0.190
F z =a2 0.187
The force on the submerged horizontal cylinder is calculated by integrating the hydrodynamic pressure 2 2 ! ›f r ›f ›f p ¼ 2r 2 ð22Þ þ ›t 2 ›x ›z Table 3 Predicted, analytical (Ogilvie, 1963) and experimental (Chaplin, 1984) first order Fourier oscillatory force components Predicted
Predicted
Ogilvie (1963)
Ogilvie (1963)
Chaplin (1984)
Chaplin (1984)
lFxð1Þ l=a 0.0581
lFzð1Þ l=a 0.0581
lFxð1Þ l=a 0.0576
lFzð1Þ l=a 0.0576
lFxð1Þ l=a 0.0577
lFzð1Þ l=a 0.0577
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77 Table 4 Predicted, analytical (Vada, 1987) and experimental (Chaplin, 1984) second order Fourier oscillatory force components Predicted
Predicted
Vada (1987)
Vada (1987)
Chaplin (1984)
Chaplin (1984)
lFxð2Þ l=a2 0.293
lFzð2Þ l=a2 0.307
lFxð2Þ l=a2 0.31
lFzð2Þ l=a2 0.31
lFxð2Þ l=a2 0.39
lFzð2Þ l=a2 0.39
Fig. 18. Voronoi section of mesh for progressive wave simulations with two cylinders.
to give F ¼ F a þ Fb
ð23Þ
where Fa ¼ 2 r
ð Sc
›f ›t
nx nz
oscillation by a non-dimensional time equal to about 100. The oscillations are predominantly at the wave frequency. Fb exhibits a significant mean vertical force component, unlike Fa : Ogilvie’s [32] np parameter for this case is np ¼ p2 p v =g ¼ 3:42065; and so nac ¼ 0:205 and 2ne ¼ 0:821 (where e is the distance between the axis of the cylinder and the near free surface). The mean vertical component of F given in Table 2 matches Ogilvie’s [32] analytical and Chaplin’s [30] values. Fourier analysis of segments of the time histories of the x and z components of F (after steady state had been reached) was used to determine the first and second order Fourier oscillatory force components. The predicted non-dimensional first order Fourier oscillatory components on the cylinder are given in Table 3, and are in reasonable agreement with the results from Ogilvie [32] and Chaplin [30]. Almost the same results were achieved using a mesh with approximately half the number of nodes. Table 4 presents the predicted non-dimensional second order Fourier oscillatory force components, compared with results from the analysis by Vada [33] and the experiments of Chaplin [30]. The overall agreement suggests that the finite element hybrid mesh numerical model should be suitable for solving a range of non-linear wave structure interaction problems. 5.4. Regular wave interaction with a pair of submerged horizontal cylinders
! dSc
ð24Þ
and
! ! rð ›f 2 ›f 2 nx F ¼2 dSc þ 2 Sc ›x ›z nz b
75
ð25Þ
where the integrals are with respect to the distance around the cylinder surface, Sc : Figs. 16 and 17 show the time-dependent behaviour of the x- and z- components of the force vectors, Fa ; and Fb : The amplitudes of the force components are initially zero, then increase after a non-dimensional time equal to 20 approximately, and eventually build up to a steady state
Fig. 18 shows the physical mesh obtained for regular waves interacting with a pair of submerged horizontal cylinders, each of non-dimensional radius apc ¼ 0:06; at non-dimensional time tp ¼ 135:9 (corresponding to the previous case). The upwave cylinder is centred at xp ¼ 10 with its centre 0.12 below the still water free surface. The centres of the two cylinders are located a non-dimensional distance 0.91842 apart (equal to half the wavelength of undisturbed progressive waves). Fig. 19 shows the spatial variation in free surface elevation in the region of the pair of submerged cylinders at tp ¼ 135:9: The free surface upwave of the first cylinder is identical to the undisturbed regular wave profile and that ahead of the single submerged
Fig. 19. Space profile of free surface elevation in region of cylinders ðtp ¼ 135:9Þ: Solid line: two cylinders; dashed line: one cylinder; dotted line: no cylinder.
76
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
horizontal cylinder in Fig. 15. The effect of the pair of cylinders is to cause a further phase shift in the wave profile after the regular waves have passed over the cylinders. This effect is significantly larger than for the single submerged cylinder case.
6. Conclusions This paper has described a 2D finite element scheme for predicting wave-induced forces on a submerged horizontal body. The finite element mesh comprises a combination of Voronoi unstructured and regular structured mesh regions, the former being used in the vicinity of the cylinder. The hybrid mesh is adaptive and uses a Lagrange – Eulerian scheme to track the vertical locations of free surface nodes. The free sloshing validation tests indicate that almost identical agreement is achieved using entirely Voronoi and combined Voronoi-structured grids. The results confirm the influence of non-linear effects on steep sloshing reported by Wu and Eatock Taylor [7] and Chern et al. [27]. The present hybrid mesh scheme predicts velocity components under an asymmetric sloshing wave that closely match corresponding results from other codes, reported by Nestega˚rd [28]. Progressive regular waves are simulated in the absence of a submerged body, and close agreement obtained with fully non-linear simulations using an alternative sigmatransformed finite element model developed by Turnbull et al. [29]. Mean, first order oscillatory and second order oscillatory forces on a submerged horizontal cylinder are calculated, and found to be in close agreement with analytical results from Ogilvie [32] and Vada [33] as well as experimental measurements by Chaplin [30]. The effect of a second submerged horizontal cylinder is to accentuate the phase shift experienced by regular waves as they pass over a single submerged cylinder.
Acknowledgements The authors gratefully acknowledge the financial support provided by the UK Engineering and Physical Sciences Research Council (EPSRC) Grant GR/M56401. The authors would also like to thank Professor G.X. Wu of University College London who provided the original finite element code.
Appendix A
Notation list A Aij
intermediate quantity for calculation of l coefficients of matrix ½A
½A B Bi {B} F Fa Fb F z lFxð1Þ l lFxð2Þ l N Ni P R S Sc S1 S2 a ap ac apc b1 b2 c d e f1 f2 g gp h l la lb lk n nv nx nz p q q1 q2 q3 r s t
matrix for calculation of {f} intermediate quantity for calculation of l coefficients of vector {B} vector for calculation of {f} total force vector on cylinder component of force vector on cylinder component of force vector on cylinder mean vertical force on cylinder first order horizontal oscillatory force component on cylinder second order horizontal oscillatory force component on cylinder total number of nodes in the finite element mesh shape function a point in the domain a prescribed radial distance measured outward from the boundary of the circle fluid boundary distance around cylinder surface section of boundary where f is specified section of boundary where ›f=›n is specified initial amplitude for symmetrical sloshing case; amplitude of progressive wave non-dimensional amplitude ðap ¼ a=hÞ radius of cylinder non-dimensional radius of cylinder ðapc ¼ ac =hÞ length of line segment at bottom of domain length of line segment on inner circle wave celerity number of line segments on the lateral outer boundaries distance between cylinder centre and undisturbed free surface intermediate quantity for calculation of l intermediate quantity for calculation of l acceleration due to gravity non-dimensionalised gravity ðgp ¼ 1Þ: water depth; distance between lower and upper boundaries at P length of line segment length of line segment length of line segment vector used in least squared velocity calculation normal to the boundary number of nodes in the vertical used for velocity calculation x-component of normal vector z-component of normal vector pressure velocity normal to the domain mesh quality parameter mesh quality parameter mesh quality parameter radial distance from boundary of circle length of line segment at the top of the domain time
M.S. Turnbull et al. / Applied Ocean Research 25 (2003) 63–77
tp u w x y z V a as
b bs gs h n np r f fi {f} v vp
pffiffiffiffi non-dimensional time, tp ¼ t g=h horizontal velocity component vertical velocity component horizontal distance distance measured upwards from the lower outer boundary vertical distance fluid domain real number used for vertical velocity calculation parameter for definition of initial free surface for asymmetric sloshing case real number used for vertical velocity calculation parameter for definition of initial free surface for asymmetric sloshing case parameter for definition of initial free surface for asymmetric sloshing case free surface elevation Ogilvie parameter non-dimensional Ogilvie parameter, np ¼ nh fluid density velocity potential nodal velocity potential vector of nodal velocity potentials wave frequency for travelling wave pffiffiffiffi case non-dimensional time, vp ¼ v h=g
References [1] Chaplin JR, Rainey RCT, Yemm RW. Ringing of a vertical cylinder in waves. J Fluid Mech 1997;350:119 –47. [2] Longuet-Higgins MS, Cokelet ED. The deformation of steep surface waves on water, I: a numerical method of computation. Proc R Soc Lond, A 1976;350:1–25. [3] Vinje T, Brevig P. Non linear ship motions. Proceedings of the Third International Symposium Numerical Ship Hydrodynamics, Paris (Chairmen: Dern JC, Haussling HJ). Bassin d’Essais des Carenes, France, 1981. p. 257 –68. [4] Cao Y, Schultz WW, Beck RF. Three dimensional unsteady computations of non-linear waves caused by underwater disturbance. Proceedings of the 18th Symposium on Naval Hydrodynamics, Ann Arbor, MI, 1990. p. 417–27. [5] Celebi MS, Kim MH. Non linear wave body interactions in a numerical wave tank. Proceedings of the 12th International Workshop on Water Waves and Floating Bodies, France, 1997. [6] Wang P, Yao Y, Tulin MP. An efficient numerical tank for non-linear water waves, based on the multi-subdomain approach with BEM. Int J Numer Meth Fluids 1995;20:1315 –36. [7] Wu GX, Eatock Taylor R. Finite element analysis of twodimensional non-linear transient water waves. Appl Ocean Res 1994;16:363 –72. [8] Wu GX, Eatock Taylor R. Time stepping solutions of the two dimensional non-linear wave radiation problem. Ocean Engng 1995; 22(8):785–98. [9] George PL. Automatic mesh generation. Paris: Wiley; 1991.
77
[10] Yerry MA, Shephard MS. A modified quadtree approach to finite element mesh generation. IEEE Comput Graph Appl 1983;3(1): 39–46. [11] Yerry MA, Shephard MS. Automatic 3D mesh generation by the modified octree technique. Int J Numer Meth Engng 1984;20(11): 1965– 90. [12] Greaves DM, Borthwick AGL. Hierarchical tree-based finite element mesh generation. Int J Numer Meth Engng 1999;45:447–71. [13] Greaves DM, Borthwick AGL, Wu GX, Eatock Taylor R. A moving boundary finite element method for fully non-linear wave simulations. J Ship Res 1997;3:181– 94. [14] Zhu G, Borthwick AGL, Eatock Taylor R. A finite element model of interaction between viscous free surface waves and submerged cylinders. Ocean Engng 2001;28(8):989–1008. [15] Lo SH. A new mesh generation for arbitrary planar domains. Int J Numer Meth Engng 1985;21:1403. [16] Peraire J, Vahdati M, Morgan K, Zienkiewicz OC. Adaptive remeshing for compressible flow computations. J Comput Phys 1987;72(2):449–66. [17] Lee CK, Hobbs RE. Automatic adaptive finite element mesh generation over arbitrary two-dimensional domain using advancing front technique. Comput Struct 1999;71(1):9–34. [18] Chan CT. Computation of flows by the finite volume method as applied to unstructured meshes. PhD Thesis, University of London (Imperial College), 1997. [19] Lawson CL. Generation of a triangular grid with application to contour plotting. California Institute of Technology, Jet Propulsion Lab, Technical Memo 299, 1972. [20] Bowyer A. Computing Dirichlet tessellations. Comput J 1981;24:162. [21] Watson DF. Computing the n-dimensional Delaunay tessellation with application to Voronoi polytropes. Comput J 1981;24(2):167–72. [22] Mavriplis DJ. Adaptive mesh generation for viscous flows using Delaunay triangulation. J Comput Phys 1990;90(2):271–91. [23] Borouchaki H, George PL. Aspects of 2D Delaunay mesh generation. Int J Numer Meth Engng 1997;40(11):1957– 75. [24] Rivara MC, Inostroza P. Using longest side bisection techniques for the automatic refinement of Delauney triangulations. Int J Numer Meth Engng 1997;40(4):581 –97. [25] Xu X, Pain CC, Goddard AJH, de Oliveira CRE. An automatic adaptive meshing technique for Delaunay triangulations. Comput Meth Appl Mech Engng 1998;161(3–4):297– 303. [26] Turnbull MS. The numerical modelling of steep waves interacting with structures. DPhil Thesis, University of Oxford, UK, 1999. [27] Chern MJ, Borthwick AGL, Eatock Taylor R. A pseudospectral stransformation model of 2-D non-linear waves. J Fluids Struct 1999; 13:607–30. [28] Nestega˚rd A. Comparative study of fully non-linear wave simulation programs. DNV Report No 94-2041, 1994. [29] Turnbull MS, Borthwick AGL, Eatock Taylor R. Numerical wave tank based on a s-transformed finite element viscous flow solver. Int J Numer Meth Fluids, in press. [30] Chaplin JR. Nonlinear forces on a horizontal cylinder beneath waves. J Fluid Mech 1984;147:449 –64. [31] Dean WR. On the reflection of surface waves by a submerged cylinder. Proc Camb Phil Soc 1948;44:483– 91. [32] Ogilvie TF. First- and second-order forces on a cylinder submerged under a free surface. J Fluid Mech 1963;16:451–72. [33] Vada T. A numerical solution of the second-order wave-diffraction problem for a submerged cylinder of arbitrary shape. J Fluid Mech 1987;174:23 –37.