Propulsive performance of a flapping plate near a free surface

Journal of Fluids and Structures 65 (2016) 411–432

Contents lists available at ScienceDirect

Journal of Fluids and Structures journal homepage: www.elsevier.com/locate/jfs

Propulsive performance of a flapping plate near a free surface Meng-Hsuan Chung Department of Naval Architecture and Ocean Engineering, National Kaohsiung Marine University, No. 142, Haijhuan Rd., Nanzih District, Kaohsiung City 811, Taiwan, ROC

a r t i c l e i n f o

abstract

Article history: Received 3 November 2015 Received in revised form 29 June 2016 Accepted 1 July 2016

The propulsive performance, i.e., the time-averaged thrust coefficient or the propulsive efficiency, of a flapping flat plate advancing near an otherwise quiescent free surface (liquid–gas interface) with Re of 1000, Fr of 0.2 and 0.8, and various submergence depths is numerically investigated by employing an adaptive Cartesian cut-cell/level-set method. The flapping kinematics parameters excluding the pitch-leading-heave phase angle were fixed as those commonly seen in literature. Results show that for submergence depth larger than the heave amplitude, the propulsive performance peaks at a smaller pitchleading-heave phase angle with a shallower submergence for Fr of 0.2 but at the same phase angle for Fr of 0.8. Proximity to the free surface enhances the peak propulsive performance for Fr of 0.2 but the influence is minor for Fr of 0.8. The propulsive performance with Fr of 0.2 increases with decreasing chord-normalized submergence depth for the pitch-leading-heave phase angle smaller than 100°. The trend is reversed for the pitch-leading-heave phase angle larger than 100°. However, the propulsive performance with Fr of 0.8 hardly depends on the chord-normalized submergence depth. For submergence depth equal to the heave amplitude, the temporal variation in the thrust coefficient exhibits characteristics inherently different from those with other submergence depths for Fr of 0.2. Also, the time-averaged thrust coefficient exhibits a unique variation with the pitch-leading-heave phase angle. However, the various characteristics of the propulsive performance are similar to those with other submergence depths for Fr of 0.8. For submergence depth smaller than the heave amplitude and Fr of 0.2, the propulsive performance gains much from exposure of the upper surface of the plate to the gas phase. The efficiency enhancement is linked to the weakening of the leading edge vortices. A second harmonic with significant amplitude is found in the upstream wave for Fr of 0.2 with a typical pitch-leading-heave phase angle. & 2016 Elsevier Ltd. All rights reserved.

Keywords: Flapping plate Free surface Submergence depth Froude number Pitch-leading-heave phase angle Adaptive Cartesian cut-cell/level-set method

1. Introduction In designing and studying man-made vehicles, locomotion (flight or swimming) of natural species is a good source of inspiration and rule of guidance for the engineers and scientists in the field of aerodynamics or hydrodynamics. For these vehicles based on the science of biomimetics, the propulsion system plays a significant role. For locomotion in water, propulsion by flapping hydrofoil is very common. It is similar to the thunniform mode of fish swimming (Breder, 1926; Lindsey, 1978) which is by far the most efficient locomotion mode in nature evolving in the aquatic environment, capable of maintaining high cruising speeds for long periods. For this mode of locomotion, the swimmer’s body keeps nearly straight E-mail addresses: [email protected], [email protected] http://dx.doi.org/10.1016/j.jfluidstructs.2016.07.003 0889-9746/& 2016 Elsevier Ltd. All rights reserved.

412

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

and significant lateral movements occur only at the caudal fin (producing more than 90% of the thrust) and at the area near the narrow peduncle. On the other hand, the flapping motion, combining with the rowing motion (Blake, 1983), can well model the oscillatory movements of the pectoral fins in the labriform mode of fish swimming (Breder, 1926; Lindsey, 1978) which is suitable for maneuvering and stabilization (Videler, 1993). For the flapping motion, the following simple pitch-and-plunge (or pitch-and-heave) motion is most considered in the literature.

yh ( t ) = h1 sin ( 2πft ),

(1)

α ( t ) = α0 + α1 sin ( 2πft + δ ),

(2)

where yh(t) and α(t) are the instantaneous heave position of the pitch axis and angle of attack respectively, h1 the heave (or plunge) amplitude, α0 the mean pitch angle, α1 the amplitude of the sinusoidal pitch angle variation, f the flapping frequency and δ the phase angle of pitch leading heave. The Reynolds number, Re, is defined as Uc/ν where U represents the free stream velocity, c the foil chord, and ν the kinematic viscosity of the fluid. There has been a lot of experimental works on underwater oscillating (mostly flapping) foils (Triantafyllou et al., 2004). The most comprehensive numerical investigations were given by Ramamurti and Sandberg (2001), who solved the Reynolds Averaged Navier-Stokes (RANS) equations using a finite element method on an unstructured grid for the NACA 0012 foil. The predicted thrust coefficient versus δ is generally consistent with the experiments of Anderson (1996) at Re of 1.1
103, but there is significant disagreement of higher-frequency thrust coefficients with the experimental results of Koochesfahani (1987) at Re of 1.2
104. For inviscid flows the thrust and power coefficients reach their minima and the propulsive efficiency its maxima approximately at δ of 90° for a low-frequency and small-amplitude flapping foil , as indicated in the review paper by Platzer et al. (2008). For viscous flows in the above regimes of Reynolds number, the maximum efficiency also occurs at δ of about 90° and the leading edge vortex (LEV) is beneficial to the thrust enhancement but harmful to the propulsive efficiency (Platzer et al., 2008). In contrast to the inviscid flow, the maximal thrust coefficient and the peak efficiency occur at approximately the same phase angle (Ramamurti and Sandberg, 2001). It is interesting to investigate the effects of free surface proximity on the pitch-leading-heave phase angle with the peak performance and the relation between the propulsive efficiency and the LEV. Hydrodynamics of foil oscillating beneath a free surface has attracted less attention as compared to that in unbounded fluid. Two additional parameters, i.e., Froude number Fr (≡U/(gc)1/2 with g denoting the gravitational acceleration) and depth of submergence, have to be considered in this problem setup. Isshiki (1982a) improved the linearized potential flow theory for oscillating foil propulsors in waves (Wu, 1972; Wu and Chwang, 1975) by including the free surface effect and applied it to the investigation of a passive-type hydrofoil propulsor (Isshiki, 1982b; Isshiki and Murakami, 1984). Isshiki (1982b) pointed out a possibility of an optimum depth of submergence for the wave devouring efficiency when the foil moves without oscillation. Similarly, Grue et al. (1988) studied a flapping flat plate moving with a prescribed forward speed in calm water or waves. They found that, for moderate values of Uf/g, the free surface strongly influences the vortex wake and the forces upon the foil. However, at some lower wave numbers there were significant discrepancies between their theory and experimental results of Isshiki and Murakami (1984). They attributed these inconsistencies to not considering non-linear effects and modeling only partial free-surface effects in the theory. Zhu et al. (2006) and Cleaver et al. (2013) studied a plunging foil beneath a free surface, using the potential flow theory and experimental test respectively. Cleaver et al. (2013) found that the drag reduction as a function of the dimensionless plunging frequency departures more from the predictions by the inviscid linear theory (Garrick, 1937) with decreasing submergence depth. Vortex lock-in was reported for the first time in experiments. Xu and Wu (2013), via the linearized potential flow theory, analyzed wave radiation and diffraction of a surge-heave-pitch hydrofoil advancing in waves. The effects of submergence depth were not discussed, only knowing that it was set to some values larger than 0.6c. Further, propulsive performance was not of their concern. Filippas and Belibassakis (2014), based on the boundary-element-method numerical results for two depths of submergence (1.5c and 2.5c), concluded that the thrust coefficient is more reduced when a flapping foil advances closer to a calm free surface, essentially due to development of wave resistance. Although some effects of the submergence depth and Froude number have been investigated in the above literature, the former was always limited to large values for not touching/piercing the free surface. Further, all the theoretical analyses were based on the potential-flow assumption, causing difficulty or even impossibility to resolve complicated viscous flows, e.g., those due to large-amplitude motion of foil. In recent years, De Silva and Yamaguchi (2012) used the commercial software FLUENT to numerically study the potential of wave energy exploitation by an active-type flapping-foil propulsor with Re of 5
107. There were only few choices of submergence depth (41c) and no systematic analysis was conducted. To gain deep insight into the fluid dynamics mechanism associated with a flapping foil advancing in waves, the study on a flapping foil advancing near an otherwise quiescent free surface should be first conducted. In a similar respect, Kajtar and Monaghan (2012) applied the Smoothed Particle Hydrodynamics (SPH) method to study the two-dimensional fish-swimming (mimicked by three linked relative-angle-changing rigid bodies) performance near a free surface. The energy consumption per unit distance is lowest when the fish swims within one third of the fish length from the free surface if wave breaking is not severe, as compared with deeper swimming. This trend is inconsistent with that of the thrust coefficient in inviscid flows mentioned above (Filippas and Belibassakis, 2014). Finally, the influences of Froude number were not discussed.

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

413

ρ1 , μ1 ρ2 , μ2 h U

g

yh(t)

α(t)

c Pitch axis

w

Fig. 1. Definition of the problem. U: inflow velocity, g: gravitational acceleration, c: plate length, w: plate thickness, h: distance from the pitch axis to the quiescent free surface when the pitch axis is located at zero heave position, yh(t): instantaneous heave position of the pitch axis, α(t): instantaneous angle of attack, ρi: mass density of the i-th fluid phase, μi: dynamic viscosity of the i-th fluid phase.

1.1. Research objectives Kajtar and Monaghan (2012) treated the self-propulsion scenario to understand undulation-like swimming performance. However, a large amount of natural swimmers and most man-made vehicles have oscillating parts as major propulsors, which can be often regarded hydrodynamically separated from the main body. The study of a flapping foil, as in the present work, thus would be an important issue for these applications. A computational fluid dynamics research code based on the Cartesian cut cell/level set approach with adaptive mesh refinement (AMR) was used to study the free surface flow past a flapping flat plate with Re of 1000, h1/c of 0.5, α0 of 0°, α1 of 20°, and the Strouhal number St (defined as fc/U) of 0.2. Fig. 1 shows the schematics of the present problem. With the present choice of Re, we avoided the ambiguity in higher-Re regimes mentioned in Mittal (2004) and acquired as high as possible resolution considering the limit of the laminar code. The present choice of St matches the numerical results for the optimal propulsive efficiency with Re of 1100 reported by Ramamurti and Sandberg (2001) and Guglielmini and Blondeaux (2004). The choice of α1 of 20° ensures generation of the LEV, which is crucial to the hydrodynamic performance of a flapping foil. The pitch axis is located one third chord aft the foil leading edge (L.E.). Because the sectional shape of the foil is not of our concern at present, a flat-plate foil with a thickness-to-chord ratio (w) of 1/20 and rounded ends was used. For an advancing flapping foil, the present work is the first study featuring viscous flow and very close proximity to an otherwise quiescent free surface. We aim to examine the effects of changing the Froude number, normalized submergence depth, h/c, and phase lead of pitch over heave, δ, on the time-averaged thrust coefficient and the propulsive efficiency.

2. Governing equations 2.1. Single-phase fluid flow In the present work two dimensional flow field of two immiscible incompressible fluids is calculated. The solutions to both phases, denoted as phase 1 (gas) and phase 2 (liquid), are obtained simultaneously. If ρ2, U, and c are used as the characteristic mass density, velocity, and length respectively, the nondimensional equations of motion are given by the incompressible Navier–Stokes equations

(3)

∇⋅u = 0, ⎞ ∂u 1 1⎛ 1 ⎡ + ( u⋅∇) u = − 2 j + ⎜ −∇p + ∇⋅⎣ μ ∇u + ∇uT ⎤⎦ ⎟, ⎠ ∂t ρ⎝ Fr Re

(

)

(4)

where u¼ u i þ v j is the fluid velocity vector, p ¼ p(x,t) the static pressure, ρ ¼ ρ(x,t) the fluid density, μ ¼ μ(x,t) the viscosity ratio of the local fluid to the liquid phase, x ¼x iþy j the position vector, and t the time. The superscript “T” represents transpose of matrix. Vector i (j) is the unit vector directed toward the positive x (y) coordinate axis. The Reynolds number is defined as Re≡ρ2Uc/μ2. All physical quantities below are made dimensionless by ρ2, U, and c unless stated otherwise. 2.2. Free surface capturing For immiscible fluids the density and viscosity are constant following the fluid particle, hence rapidly jump at the interface. Severe numerical diffusion is inevitable when capturing this jump by an Eulerian approach. In the present work, the interface is instead captured through the use of level set function (Sussman et al., 1994). Let ϕ denote the level set function with positive value in phase 1 and negative value in phase 2, then the interface can be represented by the zero level set of ϕ . The initial values of ϕ are prescribed as the signed normal distance to the interface. Following the evolutionary equation,

∂ϕ + ( u⋅∇) ϕ = 0, ∂t

(5)

the zero level set will move in the same way as the interface does. Further, due to smoothness of ϕ, the above equation can

414

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

uB = 1.5 u1 - 0.5 u2, vB = 0, φ B = φ 1

(-Lx, Ly)

lb

ls

(xhigh, ystill) y

x

Sponge Buffer layer layer ls

uB prescribed according to the flapping kinematics and φ B = φ 1

lb

uB = u1+Δu, vB = v1, φB = φstill

uB = 1, vB = 0, φ B = φstill

Buffer Sponge layer layer (xlow, ystill)

(Lx, Ly)

2 1

(-Lx, -Ly)

uB = 1.5 u1 - 0.5 u2, vB = 0, φB = φ1

B

(Lx, -Ly)

Fig. 2. Schematic of computational domain, boundary conditions of the fluid velocity and level set function, and other numerical parameters. The subscript “B” denotes the boundary node. The subscript “1” and “2” denote respectively the nearest and second nearest interior solution nodes in the direction normal to the boundary. Example positions of these nodes are schematically illustrated on the bottom boundary. Δu¼ uniform correction of u to fulfill the global conservation of mass. ϕstill ¼initial value of the level set function in the still fluid.

be solved without problem in numerical stability. Then the fluid density and dynamic viscosity can be obtained algebraically as smooth functions of ϕ (Chung, 2015). In order to damp outgoing disturbances generated by solid bodies, numerical sponge layers are introduced close to the upstream and downstream boundaries of the computational domain. That is, the evolutionary equation is modified by adding a damping term the x-direction strength of which grows with the distance from local to some prescribed position, xlow or xhigh in Fig. 2 (Chung, 2015). To avoid boundary conditions conflicting with the residual disturbances therein, we further introduce a buffer layer where the value of ϕ is fixed to the local initial value (Fig. 2). The success of interface capturing is guaranteed only if ϕ remains a distance function. However, this requirement usually can not be met in the evolution of flow field. Hence ϕ must be frequently reinitialized to recover to a distance function. The reinitialization method of Sussman et al. (1994) was chosen as in Chung (2015).

3. Numerical method 3.1. Fluid flow solver To solve the above governing equations of fluid flow, we adopt the adaptive Cartesian cut-cell/level-set method (Chung, 2015) which is an improved version of Chung (2013) in view of the nominally second-order accuracy in space and time and the upgraded refinement strategy. The method is briefly introduced as follows. The governing equations, Eqs. (3)–(5), constitute a nonlinear and coupled system of partial differential equations. A cell-centered collocated finite volume Cartesian grid method with the cut cell approach is applied to discretize this set of coupled equations. The coupling between Eqs. (3) and (4) is solved by a pressure-free projection method. Eqs. (4) and (5) are numerically advanced in time using the Crank-Nicholson scheme to generate intermediate velocities, which are in turn corrected by solving the pressure Poisson equation derived from Eq. (3). The u, v, and ϕ equations are solved in a fully coupled way. Implementation details of boundary conditions on embedded solid-body surface, cell merging, and interpolation strategies for irregular cells can be found in Chung (2006); the treatment of newly-exposed cells is the same as in Chung (2015). To solve the reinitialization equation, the second-order Runge–Kutta method is adopted for time integration and the second-order ENO scheme (Osher and Sethian, 1988) to approximate the first-order derivatives in space. For each solution variable (u, v, p, or ϕ), the discretized equations for all cells form a system of algebraic equations with a sparse coefficient matrix. The restarted Generalized Minimal RESidual method (Saad and Schultz, 1986), called GMRES(m), in the free open-source package of libraries, SPARSKIT Version 2, was selected to solve the system. 3.2. Propulsive performance parameters The two common parameters are the thrust coefficient, CT, and propulsive efficiency, η, which are defined respectively as

CT = 2FT /(p2 U 2c )

(6)

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

415

and

η=

FT U P

FT =

∮S ⎜⎝ pnx + Re1 ∂∂uns ny ⎟⎠ ds

(7)

where

⎛

⎞

(8)

and

P=

⎡ ⎛

⎞

⎛

⎞⎤

∮S ⎢⎣ u ⎜⎝ pnx + Re1 ∂∂uns ny ⎟⎠ + v ⎜⎝ pny + Re1 ∂∂uns nx ⎟⎠ ⎥⎦ ds

(9)

are respectively the instantaneous thrust force and mechanical power input to the ambient fluid by the hydrofoil motion. In Eqs. (8) and (9), nx and ny are respectively the x and y components of the unit vector normal to a differential arc length, ds, along the surface of hydrofoil, S, pointing toward fluid; us is the velocity component parallel to ds. The operator ⋅ denotes time averaging over one flapping period, T≡1/f, defined as

⋅ ≡

1 T

∫t

t0+ T

⋅ dt .

(10)

0

Note that η ¼〈CT〉 /〈CP〉 in the present nondimensionalization system where the power coefficient CP≡2P/(ρ2U c). For the cut-cell method, the solid boundary frequently moves across the mesh cell faces such that, from the current to next time step, some mesh cell centers become solution nodes and some others turn to be not involved in the solution. Mesh-scale numerical oscillations in the pressure signal thus appear. The mesh-scale oscillation however will not influence the physical carrier signal, which varies with much lower frequencies compared with the former. For visual clearness, we smooth out these numerical noises before any time record of force is plotted. 3

4. Results and discussion In all the following simulations, ρ2/ρ1 ¼815.7 and μ2/μ1 ¼63.88, unless stated otherwise. The settings to determine the time-step size and the convergence criterion for solution of the algebraic equation system is exactly the same as in Chung (2015). 4.1. Computational domain and other numerical parameters The computational domain, boundary conditions of the fluid velocity and level set function, and other numerical parameters are shown in Fig. 2. The computational domain is a rectangle centered at (0, 0) with sizes of 2Lx and 2Ly in the x and y directions respectively, excluding the flat plate. The equilibrium position of the pitching center or pitch axis is at the origin. A horizontal uniform flow enters with unit speed at the inlet plane. The fluid is not allowed to flow across the top/bottom boundary and the slip velocity is linearly extrapolated out using the first two nearest solution nodes. On the outlet plane, the values of u and v simply take those of the nearest solution node. However, u is further corrected uniformly to fulfill the global conservation of mass (Chung, 2015). The boundary condition of pressure comes from the equation to correct the cellface-center velocity. The projection of this equation on the direction normal to the boundary shows that the normal derivative of pressure is constituted by the velocity correction and gravitational acceleration terms. On the inlet/outlet plane, the level set function keeps constant at the initial value in the still fluid (Chung, 2015). On the top/bottom slip boundary, the value of ϕ takes that of the nearest solution node (ϕB ¼ ϕ1 as illustrated in Fig. 2), i.e., we apply ∂ϕ/∂n¼ 0 where n is the coordinate along the direction normal to the boundary surface. On the plate surface, the noslip boundary condition causes any contact line (intersection of the liquid-gas interface with the solid surface) fixed relative to the solid surface. However, the fluid flow surrounding a given contact line will later guide some other liquid particle to touch the solid surface at another location, forming a new contact line. Continuous emergence of a series of contact lines brings about the dynamic contact line. Because the thermodynamic equilibrium among the three phases (solid, liquid, and gas) is not considered for contact lines in the present work, the position of contact line on the solid surface is determined by the reasonable condition that the liquid-gas interface orthogonally intersects with the solid surface, i.e., ∂ϕ/∂n ¼0, which is the same as applied on the slip boundary. In all the following simulations, Lx ¼ 48 and Ly ¼32, which are the same as that determined in studying a transversely vibrating circular cylinder beneath a free surface (Chung, 2015). Due to the higher Re in the present study than in Chung (2015), the former needs a numerical sponge layer longer than that used in the latter, i.e., 7. Through a series of tests, ls and Cdamp, were determined to be 31 and 0.01 respectively such that the outgoing disturbances can be damped out without incurring significant inward reflection near the inner end of the sponge layer. Refer to Chung (2015) for the meaning of the

416

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

0.35

0.12

0.30

0.10

η



0.25 0.08 0.20 0.06 0.15 0.04

0.02

0.10

0

1

2

N

3

4

0.05

Fig. 3. Test of grid independency for the predicted o CT 4 (square) and η (circle) with Fr of 0.8, h/c of 0.5, and δ of 90°. “N” denotes the grid resolution index, with 1, 2, and 3 representing coarse, medium, and fine respectively. Two thickness-to-chord ratios, w of 0.125 (hollow) and 0.05 (solid), were considered.

damping coefficient, Cdamp. The length of buffer layer, lb, is 1. 4.2. Test of grid independency We tested grid independency by simulating the case with Fr of 0.8, h/c of 0.5, and δ of 90° using three grids whose smallest mesh sizes, Δmin, are 1/32, 1/64, and 1/128 in order. Besides the flow structure dependent strategies for AMR (Chung, 2015), three region dependent criteria are imposed to further reduce the grid size to a manageable scale: (1) the smallest allowable mesh size is Δmin within distance of 0.2 of the plate surface; (2) the mesh size is fixed at Δmin near the free surface; (3) the smallest allowable mesh sizes in the rectangular regions corresponding to the x-direction range of [  Lx, 2], [2, 6], [6, 10], and [10, Lx] in order are 2, 4, 8, and 16 Δmin. Two thickness-to-chord ratios, w of 0.125 and 0.05, were considered. Results show that both the time-averaged thrust coefficient and the propulsive efficiency differ little between the second and third grids, irrespective of the plate thickness (Fig. 3). For the thicker plate, the first-grid predictions deviate significantly from those of the other two grids. All the following simulation results were thus obtained with the second grid resolution unless stated otherwise. 4.3. Validation Three validation cases, which are relevant to the present study in terms of geometry and kinematics, were tested. For free surface involved flows, the present numerical method had been validated against previous results for the water entry problems of a wedge and a circular cylinder (Chung, 2013). The results demonstrated the capability of the method to accurately predict the hydrodynamic pressure and force exerted on the solid-body surface. For further validation, a water exit problem was simulated. 4.3.1. Impulsive start of a normal flat plate at Re of 1000 The time histories of the drag coefficient, CD ¼ CT, in the whole time span except the early stage agree well among various contributions as shown in Fig. 4. The present prediction with w of 0.025, which was obtained using the third grid resolution (Δmin of 1/128), matches those with w of 0.023 in Wang and Eldredge (2013) and w of 0.02 with square end in Hejlesen et al. (2015). Combining the present result with w of 0.05 and that with w of 0 in Koumoutsakos and Shiels (1996), we can observe strong effects of plate thickness in the early stage of impulsive start. 4.3.2. Constant angular velocity pitch-up motion There is no heave motion for this problem configuration. The pitch axis is located at the L.E. and the pitch angle, α(t), increases with t from 0° to 90°, linearly with a pitch rate ω in the whole span of the pitch-up course except at the early and final stages, where the smoothed variation (Eldredge et al., 2009) replaces the linear one. Granlund et al. (2010) have experimentally studied the case with w of 0.02, the reduced frequency K (≡cω/2U) of 0.2, and Re of 2
104 in a facility with the blockage ratio of 1/8. Wang and Eldredge (2013) numerically studied a similar problem with w of 0.023 and Re of 1000 by the viscous vortex particle method (VVPM). We performed a simulation with w of 0.02 and Re of 1000, using Δmin of 1/ 256. The drag and lift coefficients as a function of pitch angle for the three studies are similar in the global trend as shown in Fig. 5. The lift coefficient, CL, is defined as 2FL/(ρ2U2c) where FL is the instantaneous lift force. Numerical results of Liang (2009) showed strong influence of the Reynolds number on the flow structures, i.e., the vortex structure became stronger and more coherent with increasing Re, in the regime of 50 r Re r 10000 for K of 0.1. Because the present reduced

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

417

8 Present (w = 0.05) Present (w = 0.025) Wang & Eldredge (2013) Hejlesen et al. (2015) Koumoutsakos and Shiels (1996)

CD

6

4

2

0

0

0.5

1

1.5

2

2.5

3

t Fig. 4. The time history of the drag coefficient for impulsive start of a normal flat plate at Re of 1000 as the first validation case. The present predictions with w of 0.05 and 0.025 use the second (medium) and third (fine) grid resolutions respectively.

5

4

4

3

3

CL

CD

5

2

2

1

Present Granlund et al. (2010) Wang and Eldredge (2013)

0

Present Granlund et al. (2010) Wang and Eldredge (2013)

1

0 0

30°

α

60°

90°

0

30°

α

60°

90°

Fig. 5. Drag and lift coefficients as a function of the pitch angle for constant angular velocity pitch-up motion as the second validation case.

frequency is comparable to 0.1, we can therefore argue that the discrepancy between the present result and that of Granlund et al. (2010) is mainly caused by the large difference in Reynolds number. The inconsistency between the present result and that of Wang and Eldredge (2013) would most possibly stem from adopting inherently different numerical methods to solve such a highly transient and separated flows. In our simulation, the plate keeps stationary for 0 r t r 2 and the pitch motion starts at t of 2. The small non-zero value of CD at α of 0° thus reflects the nearly steady-state drag coefficient for an uniform flow over an in-line flat plat.

4.3.3. Power-extraction flapping plate The pitch axis is located at c/3 aft the L.E. and the kinematic parameters are h1/c of 1, α0 of 0°, α1 of 76.3°, f of 0.14, and δ of 90°. For Re of 1100, the calculated power coefficient as a function of time exhibits a periodic oscillation with the frequency of 2f as shown in Fig. 6. The present prediction globally agrees with that of Tian et al. (2014) with the former exhibiting smaller amplitude and smoother variation than the latter. In the study of Tian et al. (2014), the mesh size is 0.0133 near the plate and ΔtE 0.0035. The present study adopted much higher resolutions in space and time, i.e., Δmin ¼0.00244 and Δt o0.001. Also, the differences in the results between the two studies could originate partly from the difference of the plate thickness: w ¼ 0.01 for the present and 0 for the study of Tian et al. (2014).

418

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

3 Present Tian et al. (2014)

CP

2

1

0

-1

1

1.5

2

t/T Fig. 6. Power coefficient as a function of time for a power-extraction flapping plate as the third validation case. The pitch axis is located at c/3 aft the L.E., h1 /c¼ 1, α0 ¼0°, α1 ¼76.3°, f¼0.14, and δ¼ 90°.

4.3.4. Water exit of a circular cylinder With quantities made dimensionless based on the exit velocity and the cylinder diameter, the cylinder center is initially (t¼0) stationary at depth of 3.5, accelerates vertically in a harmonic manner up to the unit velocity in the interval of 0 o t o 1, then rises with the unit velocity in the interval of 1 r t r 6. The cylinder top touches the calm water surface (located at y of 0) at t of 3.5. The equation of cylinder motion is

⎧ ⎡ ⎤ 1 ⎪ − 3 + 0.5 ⎢ t − 1 − sin ( πt ) ⎥ , 0 ≤ t < 1 ⎣ ⎦ yc (t ) = ⎨ π ⎪ ⎩ t−4 , 1≤t<6

(11)

where yc is the vertical coordinate of the cylinder center. Simulations were performed with Fr of 0.4627 and two Re’s of 500 and 1000. The computational domain is [  48, 48]
[ 28, 4] with Δmin of 1/64. As shown in Fig. 7, the calculated lift coefficients as a function of time for the two Re’s are nearly the same, close to that predicted by the Constrained Interpolation Profile (CIP) method (Zhu, 2006) and that measured by Miao (1989). The Incompressible Smoothed Particle Hydrodynamics (ISPH) method, employed by Bašić et al. (2014), predicted an unexplained weird sharp depression of the lift 11 Miao (1989) Zhu (2006) Basic et al (2014) Present (Re = 500) Present (Re = 1000)

10 9 8 7

CL

6 5 4 3 2 1 0 -1 -5

-4

-3

-2

-1

0

1

t

2

3

4

5

+

Fig. 7. Lift coefficient as a function of time for water exit of a circular cylinder as the fourth validation case. The abscissa is the dimensionless time based on the exit velocity and the cylinder radius, with the origin corresponding to the cylinder top touching the calm water surface. That is, t þ ¼ 2t  7.

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

0.4

0.5

+ +

+ #

# +

# +

# +

0.4 # +

#

# +

# +

+

0.3

#

0.2 +

+ #

+ #

+ #

+ #

# + # + # +

#

#

η



0.3

0.2

419

0.1

0.1

# #

0

0 -0.1 -0.1 30°

60°

δ

90°

120°

-0.2 30°

60°

δ

90°

120°

Fig. 8. o CT 4 and η as a function of δ for h/c of 0.3 ( þ), 0.4 (#), 0.5 (□), 0.7 (♢), 1.0 (∇), 1.5 (○), and 1 (solid line with symbol ) with Fr of 0.2.

coefficient when the cylinder top touches the calm water surface. As analyzed by Zhu (2006), the large oscillatory amplitude of the experimental lift coefficient (Miao, 1989) was caused by the transient effects of the hydraulic system, which was used to force the test cylinder out of the water. 4.4. Propulsive flapping plate near a free surface We performed simulations of a series of cases with the combinations of parameters covering two Froude numbers (Fr ¼0.2 and 0.8), seven normalized submergence depths (h/c¼0.3, 0.4, 0.5, 0.7, 1.0, 1.5, and 1), and various pitch-leadingheave phase angles (30° r δ r120°). 4.4.1. Fr ¼0.2 As a reference, the peak time-averaged thrust coefficient and propulsive efficiency for the deep-submergence case occur at δ of 100° and 90° respectively as shown in Fig. 8. For h Z h1, the peak performance, in terms of either oCT 4 or η, with shallower submergence occurs at smaller δ. The peak performance for any finite h/c is better than that for deep submergence. The propulsive performances for h 4h1 are better than those for h of 1 when δ o90°. The viscous computational results of Kajtar and Monaghan (2012) showed that the energy consumption per unit distance is lowest when a model fish swims within one third of the fish length from the free surface if wave breaking is not severe, as compared with deeper swimming. The consistency in the effect of free surface proximity on the efficiency, between the present and their results, confirms that the free surface proximity would be a positive effect if it can induce an additional propulsive power to compensate for the power loss due to the surface wave propagation. For h oh1, the values of oCT 4 and η are much larger than the corresponding ones with h Zh1 for any δ. To investigate the mechanisms behind the effects of h/c on the propulsive performance, we take as examples h/c of 0.3, 0.5, and 1 with δ of 70°. Fig. 9a illustrates that the much larger value of oCT 4 for h/c of 0.3 as compared with that for h/c of 0.5 and 1 is largely due to the value of CT for the former greatly exceeding those for the latter approximately during the subinterval of 7/16 r t’/T r 11/16 (spanning one quarter of the flapping period) in the downstroke of the heaving motion (1/4 r t’/T r 3/4). The great amount of excess of CT during this subinterval is in turn caused by some mechanisms explained below via pressure coefficient and vorticity distributions at typical timings, t’/T of 8/16, 9/16, and 10/16. Among the three values of h/c, the difference in the pattern of the pressure coefficient contours near the lower surface of the plate is much smaller than that near the upper surface at the same timing (Fig. 10). The pressure coefficient distribution on both the plate surfaces clearly reflects this point as well (Fig. 11). For the upper surface, a low-pressure zone moves along toward the trailing edge (T.E.) with increasing t’/T for h/c of 0.5 and 1 (Fig. 10). The low-pressure structure is more obvious for h/c of 1 as compared with 0.5 (Fig. 10). The pressure coefficient distribution on the upper surface exhibits a corresponding depression moving toward T.E. (Fig. 11). Inspection of vorticity contours at t’/T of 8/16 and 10/16 (Figs. 12 and 13) confirms that the low-pressure structure is exactly the tiny positive vortex very close to the upper surface and moving along toward the T. E. Actually, this vortex accompanies the negative LEV. For h/c of 0.3, a considerable portion of the upper (windward) surface is exposed to the gas at t’/T of 8/16 (Fig. 10). Because of the extremely large density ratio, this portion of the upper surface experiences much lower pressure than the wetted portion of the upper surface and the entire lower surface do (Fig. 11). This pressure difference, together with the angle of attack reaching the minimum at this instant, produces the substantial increase of the thrust coefficient as compared with larger h/c. Further, the reentry of the plate into the liquid makes some low-pressure liquid flow along the upper

420

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

2

5 4

++ + + + + + ++ + + ++ ++ + ++ + + + + ++ + + ++ + + ++ + + + + + + + + + + + + ++ + + + + ++ + + ++ + ++ 0 ++++

3

CP

CT

1

2 1 0 + + + -1 +

-1 0

1/4

2/4

t'/T

3/4

4/4

0

++ + + + + + + + ++ + ++ + + + + + + + + + + + + + + + + + + + + + + + + + ++ + + + +++ + +++ + ++ + + + + + + + + + + +

1/4

2/4

3/4

+

+

+

+

++

4/4

t'/T

Fig. 9. Time histories of (a) CT and (b) CP for δ of 70° and h/c of 0.3 ( þ), 0.5 (□), 0.7 (♢), 1.5 (○), and 1 (solid line) with Fr of 0.2. The time t’ is defined as (t–tref) where tref denotes some instant when the pitch axis reaches the zero-heave position during the upstroke of the heaving motion.

Fig. 10. Contours of pressure coefficient, 2(p–pstatic–p1), around the plate at selected times for Fr of 0.2, δ of 70°, and h/c of (a) 0.3, (b) 0.5, and (c) 1. The increment of contour level is 0.125. The far upstream pressure at y of 0 is assigned to p1. The hydrostatic component, pstatic, calculated as (ystill – y)/Fr2, has been subtracted for straightforward comparison of the non-hydrostatic component between different submergence depths. Thick lines indicate the free surfaces.

surface, causes the pressure coefficient on the wetted portion much lower than the upper-surface values for larger h/c (Fig. 11), and hence additionally contributing to the enhancement of the thrust coefficient. At the two total-submergence timings (t’/T¼9/16 and 10/16), the stronger free surface distortion causes the smaller effective angle of attack as compared with larger h/c. The boundary layer on the upper surface thus does not separate, attached to the plate surface all the way to the T.E. (as illustrated at t’/T of 10/16 in Fig. 14), hence inexistence of the low-pressure zone mentioned above. The uppersurface pressure coefficients at these two timings, similar to that on the wetted portion at t’/T of 8/16, are much smaller than

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

421

2

Pressure Coefficient

t'/T = 8/16 0 -2 -4 -6 -8

h/c = 0.3 h/c = 0.5

-10 -12

h/c = ∞ 0

0.2

0.4

s/c

0.6

0.8

1

2

Pressure Coefficient

t'/T = 9/16 0 -2 -4 -6 -8

h/c = 0.3 h/c = 0.5

-10 -12

h/c = ∞ 0

0.2

0.4

s/c

0.6

0.8

1

2

Pressure Coefficient

t'/T = 10/16 0 -2 -4 -6 -8

h/c = 0.3 h/c = 0.5

-10 -12

h/c = ∞ 0

0.2

0.4

s/c

0.6

0.8

1

Fig. 11. Distribution of pressure coefficient along the plate surface at selected times for Fr of 0.2, δ of 70°, and various h/c. Symbol “s” denotes the chordwise distance from the L.E. of the plate. Thin and thick lines indicate the lower and upper surfaces respectively. Refer to Fig. 10 for the definition of the pressure coefficient.

422

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

Fig. 12. Equi-vorticity lines (solid for positive and dotted for negative) at eight equally spaced consecutive instants, corresponding to t’/T of 1/8  1, from top to bottom, firstly the left and secondly the right columns, with Fr of 0.2, h/c of 0.5, and δ of 70°. The increment of vorticity between two adjacent lines is 2. The time t’ is defined as (t  tref) where tref denotes some instant when the pitch axis reaches the zero-heave position during the upstroke of the heaving motion. The thick solid line represents the free surface.

those for larger h/c on the most part of the upper surface (Fig. 11), hence much larger thrust coefficient (Fig. 9a). The numerical spike around t’/T of 3/8 in Fig. 9 is attributed to the high-pressure gas being trapped in between the plate and the free surface (Fig. 14). This unphysical pressure stems from not satisfying the law of mass conservation for the trapped gas, which is in turn due to numerical errors in predicting the free surface location. The spike affects little the time-averaged level because of the very short life of the trapped gas. The abrupt drop in the interval 1/2 o t’/T o 9/16 in Fig. 9a can be explained as below. The wetted area of the upper surface of the plate increases with t’/T increasing from 1/2 up to some instant when the split free surfaces meet together. The portion of exposure to the gas, hence the thrust force, is thus reduced with time increasing in this subinterval. After the free surfaces meet, the pressure field established around the totally immersed plate causes the recovery of the thrust force. Actually, the plate does not totally exit the liquid even though hoh1.

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

423

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

Fig. 13. Equi-vorticity lines (solid for positive and dotted for negative) at eight equally spaced consecutive instants, corresponding to t’/T of 1/8  1, from top to bottom, firstly the left and secondly the right columns, with h/c of 1 and δ of 70°. The increment of vorticity between two adjacent lines is 2. The time t’ is defined as (t  tref) where tref denotes some instant when the pitch axis reaches the zero-heave position during the upstroke of the heaving motion.

For h/c of 0.5, the algebraic pressure coefficient difference between the lower and upper surfaces, the lower part minus the upper one, is smaller than that for h/c of 1 through the chord at all the three timings. The thrust coefficient is thus smaller as compared with h/c of 1 at these timings (Fig. 9a). The time histories of CT and CP are very different from those for other h/c and unique to this h/c. The major reasons could be manifold. First, the liquid reentry of the L.E. part can generate extraordinarily small upper-surface pressure as illustrated for δ of 70° and h/c of 0.3 (Fig. 11). A similar phenomenon exists for any currently investigated δ with h/c of 0.3 because there is sufficient time to adjust the angle of attack of the plate such that the L.E. part pierces but not slams into the liquid upon reentry. For h/c of 0.5, i.e., h ¼ h1, the occurrence of similar situations however strongly depends on δ because there is no time to adjust the angle of attack. Second, the formation, movement, and intensity of the low-pressure positive vortex very close to the upper surface are highly correlated with δ.

424

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

Fig. 14. Equi-vorticity lines (solid for positive and dotted for negative) at eight equally spaced consecutive instants, corresponding to t’/T of 1/8  1, from top to bottom, firstly the left and secondly the right columns, with Fr of 0.2, h/c of 0.3, and δ of 70°. The increment of vorticity between two adjacent lines is 2. The time t’ is defined as (t  tref) where tref denotes some instant when the pitch axis reaches the zero-heave position during the upstroke of the heaving motion. The thick solid line represents the free surface. The arrow points out the negative vortex due to free surface distortion.

The combined effects produce distinctly different time histories of hydrodynamic coefficients from those for larger h/c even with the same δ. The unique behavior of oCT 4 as a function of δ (Fig. 8), as compared with other h/c, can be partly attributed to the combined effects. Another effect is related to the gas exposure of the plate surface. Specifically, oCT 4 for δ of 120° is higher than the corresponding expected level based on the performance curves for other h/c. The time history of CT for this δ (Fig. 15) reveals that oCT 4 is higher largely because the thrust coefficient is much higher than those of other h/c approximately within the interval 6/16 r t’/T r 7/16, during which a part of windward surface of the plate is exposed to the gas (Fig. 16). Because the T.E. is higher for δ of 120° (Fig. 16) than for δ of 70° (Fig. 12) during this time interval, the exposure to the gas occurs only with the former δ, causing extraordinarily high thrust coefficient. In other words, the mechanism involved is similar to that proposed above for h/c of 0.3 with δ of 70°.

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

425

2

CT

1

0

-1

0

1/4

2/4

3/4

4/4

t'/T Fig. 15. Time histories of CT for δ of 120° and h/c of 0.5 (□), 0.7 (♢), and 1 (solid line) with Fr of 0.2. The time t’ is defined as (t  tref) where tref denotes some instant when the pitch axis reaches the zero-heave position during the upstroke of the heaving motion.

As is the case with CT, the value of CP for h/c of 0.3 is higher than that for h/c of 0.5 during the subinterval of 7/16 r t’/T r 11/16 as shown in Fig. 9b, causing a higher oCP 4 for the former than for the latter. However, the resultant propulsive efficiency for h/c of 0.3 is significantly higher than that for h/c of 0.5 because the increase of oCT 4 can afford that of oCP 4 more than sufficiently. The efficiency enhancement with h/c of 0.3 is linked to the particular vortex structures as explained below. Accompanying the free surface distortion, a negative vortex forms under the lower surface of the plate and drifts downstream while the plate pierces the free surface up to the maximal heave position. This vortex eventually catches up and pairs with the positive LEV. The typical vortex structures in the whole process, clearly seen at t’/T of 1/8 and 2/8 in Fig. 14, would diminish whatever effects the positive LEV may have. Further, as previously mentioned, due to the stronger free surface distortion, the high-vorticity fluid is concentrated closer to the upper surface of the plate than that for h/c of 0.5, e.g. at t’/T of 5/8 in Figs. 14 and 12. The shedding of the negative LEV is thus weaker or even does not occur for h/c of 0.3 as compared with h/c of 0.5. Tuncer and Kaya (2005) and Young and Lei (2007) have demonstrated similar relationships between the efficiency and the vortex structures, though they studied the deep-submergence scenario and used higher Reynolds numbers than the present one. That is, the LEV shedding has an adverse effect on the propulsive efficiency. For h 4 h1, the effects of h/c on the propulsive performance are more significant in the small-δ regime as compared with the large-δ regime and very minor for 90° o δ o100°, where the peak performance occurs with h/c of 1. Both oCT 4 and η increase with decreasing h/c for δ o100° and increasing h/c for δ 4100°. 4.4.2. Fr ¼0.8 For h Z h1, the peak propulsive performance occurs at the same δ as in the deep-submergence case as shown in Fig. 17, in contrast to Fr of 0.2. The peak performance for any finite h 4 h1 is slightly poorer than that for deep submergence. Actually for any finite h 4 h1, the performance as a function of δ with Fr of 0.8 is much closer to that for deep submergence than with Fr of 0.2. The time history of CT, with h/c of 0.7 and δ of 70° for example, can give some insights into this observation. Large amount of difference in CT between the two Froude numbers is found during almost the whole downstroke of the heaving motion, approximately 3/8 o t’/T o 3/4, as shown in Fig. 18. During this time interval, the vortex structures for Fr of 0.2 (Fig. 19) are distinctly different from those for Fr of 0.8 (Fig. 20). The latter is more similar to those for h/c of 1 (Fig. 13) than the former, particularly at t’/T of 1/2 and 5/8. In other words, the influence of proximity to a free surface on the vortex structures, hence the hydrodynamic forces, is less significant when the free surface can distort such that the flow around the plate more resembles the deep-submergence one. The effects of h/c on the propulsive performance with a given δ are very minor for all cases with h 4 h1, dissimilar to the conclusion drawn for Fr of 0.2. For this regime of submergence depth, the propulsive performance is slightly degraded as compared with the deep-submergence performance, for δ values smaller than the one where the peak performance occurs (Fig. 17). For h o h1, some simulation cases (δ of 60°, 80°  120° for h/c of 0.3 and δ of 60°  90° for h/c of 0.4) failed to numerically converge because of strongly complicated free surface flow structures. Among cases with converged solution, the results with δ of 70° for h/c of 0.3 and δ of 100–120° for h/c of 0.4 could be questionable and hence not presented in

426

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

Fig. 16. Equi-vorticity lines (solid for positive and dotted for negative) at eight equally spaced consecutive instants, corresponding to t’/T of 1/8  1, from top to bottom, firstly the left and secondly the right columns, with Fr of 0.2, h/c of 0.5, and δ of 120°. The increment of vorticity between two adjacent lines is 2. The time t’ is defined as (t  tref) where tref denotes some instant when the pitch axis reaches the zero-heave position during the upstroke of the heaving motion. The thick solid line represents the free surface.

Fig. 17 because the above-mentioned unphysical high-pressure gas, trapped in between the plate and the free surface, survives for a long time. For reliable cases with positive time-averaged thrust coefficient, the propulsive performances for h/ c of 0.3 and 0.4 are respectively superior and comparable to that for h/c of 1.

4.4.3. Free surface waves To investigate effects of the plate motion on the free surface deformations or waves, we present in Fig. 21 instantaneous free surface profiles of four cases with δ of 70° at the same timing (highest heave position). The four cases correspond to Fr of 0.2 and 0.8 and two h/c of 0.3 and 0.7. Firstly, the amount of free surface deformation or wave amplitude is much larger

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

427

0.5

0.4

0.4 0.3

0.2

0.2

η



0.3

0.1

0.1

+ #

0 + #

0

+

+ #

-0.1

+ #

-0.1 30°

60°

90°

δ

120°

#

-0.2 30°

60°

δ

90°

Fig. 17. oCT 4 and η as a function of δ for h/c of 0.3 ( þ), 0.4 (#), 0.5 (□), 0.7 (♢), 1.0 (∇), 1.5 (○), and 1 (solid line with symbol

120°

) with Fr of 0.8.

0.6

0.4

CT

0.2

0

-0.2

-0.4

0

1/4

2/4

3/4

4/4

t'/T Fig. 18. Time histories of CT for h/c of 0.7 and δ of 70° with Fr of 0.2 (♢) and 0.8 (♦). Solid line: h/c of 1. The time t’ is defined as (t  tref) where tref denotes some instant when the pitch axis reaches the zero-heave position during the upstroke of the heaving motion.

for Fr of 0.8 than 0.2. Secondly, the free surface profiles upstream of the plate remain almost unchanged with time for the two values of h/c with Fr of 0.8. A hardly perturbed free surface is also seen downstream of the plate for Fr of 0.2 and h/c of 0.7. In the linear theory of small-amplitude free surface waves on a uniform current (Thompson, 1949), the dimensionless form of the dispersion relationship is

ω = ±k+

1 k tanh ( khY ) Fr

(12)

where ω is the radian frequency of the wave, k the wave number, and hY the depth of the liquid. The plus and minus signs apply to the following current (wave and current travel in the same direction, ω 4 0) and adverse current (wave and current travel in the opposite direction, ω o 0) respectively. We observed that the frequency of the free surface wave, fwv (¼ ω/2π), is an integer multiple of f for each case of the present study. With corresponding fwv and noting that hY ¼ 32þh/c, we can find the wave length, 2π/k, for given Fr and h/c by solving Eq. (12). The results for the four cases are listed in Table 1 along with the far-field wave lengths measured from Fig. 21. It is seen that the effect of h/c on the wave length is very minor for either theory or simulation, as expected from Eq. (12) if we note that hY changes little with h/c. For Fr of 0.2, the upstream wave propagates with two times of the flapping frequency (the second harmonic) in the simulation. The wave length is

428

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

Fig. 19. Equi-vorticity lines (solid for positive and dotted for negative) at eight equally spaced consecutive instants, corresponding to t’/T of 1/8  1, from top to bottom, firstly the left and secondly the right columns, with Fr of 0.2, h/c of 0.7, and δ of 70°. The increment of vorticity between two adjacent lines is 2. The time t’ is defined as (t  tref) where tref denotes some instant when the pitch axis reaches the zero-heave position during the upstroke of the heaving motion. The thick solid line represents the free surface.

shorter than predicted by the linear theory. The first harmonic is not necessarily absent in unbounded domain; it is possibly destroyed due to too small of a computational domain as compared with the linear-theory wave length ( E 87). Irregular ripples propagate downstream for h/c of 0.3 and no downstream propagating wave is observed for h/c of 0.7. From the results of the linear theory for Fr of 0.2, we can argue that the downstream wave lengths for f and 2f in the simulation could be too large to be felt with the present limited computational domain and the harmonics higher than 2f either nonlinearly interact to form ripples (for h/c of 0.3) or have neglectable amplitudes (for h/c of 0.7). For Fr of 0.8, both the theory and simulation predict no upstream propagating waves. The downstream waves propagate with the flapping frequency and, as is

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

429

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

0

1

2

3

Fig. 20. Equi-vorticity lines (solid for positive and dotted for negative) at eight equally spaced consecutive instants, corresponding to t’/T of 1/8  1, from top to bottom, firstly the left and secondly the right columns, with Fr of 0.8, h/c of 0.7, and δ of 70°. The increment of vorticity between two adjacent lines is 2. The time t’ is defined as (t  tref) where tref denotes some instant when the pitch axis reaches the zero-heave position during the upstroke of the heaving motion. The thick solid line represents the free surface.

the case with Fr of 0.2, a wave length shorter than predicted by the linear theory. The present distinctive finding, as compared with that for a circular cylinder (Chung, 2015), is the emergence of the second harmonic in the upstream wave for Fr of 0.2.

5. Conclusions We have numerically studied a flapping flat plate advancing near an otherwise quiescent free surface at Re of 1000 with various combinations of physical parameters. The flapping kinematics parameters, excluding the pitch-leading-heave phase

430

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

0.6 0.5

0.6

Fr = 0.2 h/c = 0.3

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1 -40 1

-30

0.9

Fr = 0.2 h/c = 0.7

-20

-10

0

10

20

30

40

0.1 -40 1

-30

0.9

Fr = 0.8 h/c = 0.7

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3 -40

-30

-20

-10

0

10

20

30

40

Fr = 0.8 h/c = 0.3

0.3 -40

-30

-20

-20

-10

0

10

20

30

40

-10

0

10

20

30

40

Fig. 21. Instantaneous free surface profiles of four selected cases with δ of 70°. The snapshots are taken when the pivot axis moves to the highest position. The plate is colored in black.

Table 1 Wave length of free surface wave for a flapping plate beneath a free surface with δ of 70°. The three columns under the “Linear theory” field correspond to the flapping frequency and the first two superharmonic ones respectively. The number in the bracket represents the index of harmonics observed in the simulation (1, 2, 3 denote f, 2f, 3f respectively). Fr

h/c

Upstream wave

Downstream wave

Linear theory f 2f 3f

Present

Linear theory f 2f 3f

Present

Irregular ripples – 9.3(1) 9.3(1)

0.2

0.3

87.3 19.6 7.4

8.2(2)

105.1 29.7 14.2

0.8

0.7 0.3 0.7

87.4 19.6 7.4 – – – – – –

8.2(2) – –

105.3 29.7 14.2 14.5 5.4 3.2 14.5 5.4 3.2

–: No wave.

angle, were fixed as those commonly seen in literature. Most conclusions are categorized according to whether the submergence depth exceeds the heave amplitude or not: (i) h 4 h1 For Fr of 0.2, either the time-averaged thrust coefficient or the propulsive efficiency peaks at a smaller pitch-leadingheave phase angle with a shallower submergence. For Fr of 0.8, the peak propulsive performance however occurs at the same phase angle for all the normalized submergence depths. Proximity to the free surface benefits the peak propulsive performance for Fr of 0.2 but the influence is minor for Fr of 0.8. For Fr of 0.2, the propulsive performance increases with decreasing normalized submergence depth for the pitchleading-heave phase angle smaller than 100°, similar to the conclusions previously drawn for a self-propelled model fish beneath a free surface; the trend is reversed for the pitch-leading-heave phase angle larger than 100°. For Fr of 0.8, the propulsive performance as a function of the pitch-leading-heave phase angle with any normalized submergence depth closely resembles that for deep submergence, though the performance is slightly degraded as compared with the deepsubmergence one for small pitch-leading-heave phase angles. The influence of proximity to a free surface on the vortex structures, hence the hydrodynamic forces, is less significant because of stronger free surface distortion, as compared with Fr of 0.2. (ii) h ¼ h1 For Fr of 0.2, the characteristics of temporal variation in the thrust coefficient or power coefficient are inherently

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

431

different from those for other h. This is because the pitch-leading-heave phase angle, or equivalently the angle of attack upon liquid reentry, plays a remarkable role (as compared with other h) in changing flow structures around the plate, especially near the upper surface. Because of the above mechanism and that the pitch-leading-heave phase angle strongly affects gas exposure area of the plate surface, the time-averaged thrust coefficient exhibits unique variations with the pitch-leading-heave phase angle. For Fr of 0.8, the various characteristics of the propulsive performance are similar to those for other h. (iii) h o h1 For Fr of 0.2, the exposure of the upper surface of the plate to the gas is the major mechanism for extraordinarily large timeaveraged thrust coefficients and propulsive efficiencies, as compared with deeper submergence. The second mechanism is the piercing reentry of the plate L.E. into the liquid which induces much lower top-surface pressure coefficient than those without reentry. The corresponding efficiency enhancement is linked to the weakening of the LEVs. As for the free surface waves induced by the flapping plate, the linear theory roughly supports some results of numerical simulation, with the simulated wave length shorter than the theory-predicted one. With a typical pitch-leading-heave phase angle, the upstream wave for Fr of 0.2 possesses a second harmonic with significant amplitude whether the submergence depth is larger than the heave amplitude or not. The present results have disclosed the effects of pitch-leading-heave phase angle on the propulsive performance with various environmental parameters, i.e., Froude number and normalized submergence depth. The findings would help to design the flapping kinematics and to get some insights into the appendage (e.g. tail) motions of cetacean-like swimmers or propulsors near the water surface.

References Anderson, J.M., 1996. Vortex Control for Efficient Propulsion. Massachusetts Institute of Technology, Cambridge, MA. Bašić, J., Degiuli, N., Werner, A., 2014. Simulation of water entry and exit of a circular cylinder using the ISPH method. Trans. FAMENA 38 (1), 45–62. Blake, R.W., 1983. Median and paired fin propulsion. In: Webb, P.W., Weihs, D. (Eds.), Fish Biomechanics, Praeger Publishers, New York. Breder, C.M., 1926. The locomotion of fishes. Zoologica 4, 159–256. Chung, M.-H., 2006. Cartesian cut cell approach for simulating incompressible flows with rigid bodies of arbitrary shape. Comput. Fluids 35 (6), 607–623. Chung, M.-H., 2013. An adaptive Cartesian cut-cell/level-set method to simulate incompressible two-phase flows with embedded moving solid boundaries. Comput. Fluids 71, 469–486. Chung, M.-H., 2015. Hydrodynamics of flow over a transversely oscillating circular cylinder beneath a free surface. J. Fluids Struct. 54, 27–73. Cleaver, D.J., Calderon, D.E., Wang, Z., Gursul, I., 2013. Periodically plunging foil near a free surface. Exp. Fluids 54, 1491. De Silva, L.W.A., Yamaguchi, H., 2012. Numerical study on active wave devouring propulsion. J. Mar. Sci. Technol. 17, 261–275. Eldredge, J.D., Wang, C.J., Ol, M.,2009. A computational study of a canonical pitch-up, pitch-down wing maneuver. AIAA 2009-3687. Filippas, E.S., Belibassakis, K.A., 2014. Hydrodynamic analysis of flapping-foil thrusters operating beneath the free surface and in waves. Eng. Anal. Bound. Elem. 41, 47–59. Garrick, I., 1937. Propulsion of a flapping and oscillating airfoil. NACA Rept. 567. Granlund, K., Ol, M., Garmann, D., Visbal, M.,Bernal, L.,2010. Experiments and computations on abstractions of perching. AIAA paper 2010-4943. Grue, J., Mo, A., Palm, E., 1988. Propulsion of a foil moving in water waves. J. Fluid Mech. 186, 393–417. Guglielmini, L., Blondeaux, P., 2004. Propulsive efficiency of oscillating foils. Eur. J. Mech. B Fluids 23 (2), 255–278, http://dx.doi.org/10.1016/j. euromechflu.2003.10.002. Hejlesen, M.M., Koumoutsakos, P., Leonard, A., Walther, J.H., 2015. Iterative Brinkman penalization for remeshed vortex methods. J Comput. Phys. 280, 547–562. Isshiki, H., 1982a. A theory of wave devouring propulsion (1st Report) – thrust generation by a linear wells turbine. J. Soc. Nav. Archit. Jpn. 151, 54–64. Isshiki, H., 1982b. A theory of wave devouring propulsion (2nd Report) – optimized foil motions for a passive-type wave devouring propulsor. J. Soc. Nav. Archit. Jpn. 152, 89–100. Isshiki, H., Murakami, M., 1984. A theory of wave devouring propulsion (4th Report) – a comparison between theory and experiment in case of passive-type hydrofoil propulsor. J. Soc. Nav. Archit. Jpn. 156, 102–114. Kajtar, J.B., Monaghan, J.J., 2012. On the swimming of fish like bodies near free and fixed boundaries. Eur. J. Mech. B Fluids 328, 1–13. Koochesfahani, M.M., 1987. Vortical patterns in the wake of an oscillating airfoil. In: Proceedings of the AIAA 25th Aerospace Sciences Meeting. Reno, NV. Koumoutsakos, P., Shiels, D., 1996. Simulations of the viscous flow normal to an impulsively started and uniformly accelerated flat plate. J. Fluid Mech. 328, 177–227. Liang, Y., 2009. Parametric study of a pitching flat plate at low Reynolds numbers. AIAA-2009-3688 Presented in 39th AIAA Fluid Dynamics Conference 22 25 June 2009, San Antonio, Texas. Lindsey, C.C., 1978. Form, function and locomotory habits in fish. Hoar, W.S., Randall, D.J. (Eds.), Fish Physiology, Vol. VII. , Academic Press, New York. Miao, G., 1989. Hydrodynamic Forces and Dynamic Responses of Circular Cylinders in Wave Zones. University of Trondheim, Trondheim, Norway. Mittal, R., 2004. Computational modeling in biohydrodynamics: trends, challenges, and recent advances. IEEE J. Ocean. Eng. 29 (3), 595–604. Osher, S., Sethian, J.A., 1988. Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulation. J. Comput. Phys. 79 (1), 12–49. Platzer, M.F., Jones, K.D., Young, J., Lai, J.C.S., 2008. Flapping-wing aerodynamics: progress and challenges. AIAA J. 46, 2136–2149. Ramamurti, R., Sandberg, W.C., 2001. Simulation of flow about flapping airfoils using a finite element incompressible flow solver. AIAA J. 39 (2), 253–260. Saad, Y., Schultz, M.H., 1986. GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 7, 856–869. Sussman, M., Smereka, P., Osher, S., 1994. A level set approach for computing solutions to incompressible two-phase flow. J. Comput. Phys. 114, 146–159. Thompson, P.D., 1949. The propagation of small surface disturbances through rotational flow. Ann. NY Acad. Sci. 51, 464–474. Tian, F.-B., Young, J., Lai, J.C.S., 2014. Improving power-extraction efficiency of a flapping plate: From passive deformation to active control. J. Fluids Struct. 51, 384–392. Triantafyllou, M.S., Techet, A.H., Hover, F.S., 2004. Review of experimental work in biomimetic foils. IEEE J. Ocean. Eng. 29 (3), 585–594. Tuncer, I.H., Kaya, M., 2005. Optimization of flapping airfoils for maximum thrust and propulsion efficiency. AIAA J. 43 (11), 2329–2341. Videler, J.J., 1993. Fish Swimming. Chapman & Hall, London. Wang, C., Eldredge, J.D., 2013. Low-order phenomenological modeling of leading-edge vortex formation. Theor. Comput. Fluid Dyn. 27, 577–598.

432

M.-H. Chung / Journal of Fluids and Structures 65 (2016) 411–432

Wu, T.Y., 1972. Extraction of flow energy by a wing oscillating in waves. J. Ship. Res. 14 (1), 66–78. Wu, T.Y., Chwang, A.T., 1975. Extraction of flow energy by fish and birds in a wavy stream. In: Wu, T.Y., Brokaw, C.J., Brennan, C. (Eds.), Swimming and Flying in Nature, Vol. II. , Plenum Press, New York, pp. 687–702. Xu, G.D., Wu, G.X., 2013. Hydrodynamics of a submerged hydrofoil advancing in waves. Appl. Ocean Res. 42, 70–78. Young, J., Lai, J.C.S., 2007. Mechanisms influencing the efficiency of oscillating airfoil propulsion. AIAA J. 45 (7), 1695–1702. Zhu, Q., Liu, Y., Yue, D.K.P., 2006. Dynamics of a three-dimensional oscillating foil near the free surface. AIAA J. 44 (12), 2997–3009. Zhu, X.Y., 2006. Application of the CIP method to strongly nonlinear wave-body interaction problems. Norwegian University of Science and Technology, Trondheim, Norway.