## Abstract

In this paper, we use the lattice Boltzmann method with the Bhatnagar–Gross–Krook linear collision operator to study the flow physics induced by a rigid lamina undergoing moderately large harmonic oscillations in a viscous fluid. We propose a refill procedure for the hydrodynamic quantities in the lattice sites that are in the vicinity of the oscillating lamina. The numerically estimated flow field is used to compute the complex hydrodynamic function that describes the added mass and hydrodynamic damping experienced by the lamina. Results of the numerical simulations are validated against theoretical predictions for small amplitude vibrations and experimental and numerical findings for moderately large oscillations.

## 1. Introduction

In recent years, considerable efforts have been devoted to the problem of flexible cantilevers undergoing vibrations within fluid environments in view of their multiple applications, including cooling devices [1], biomimetic underwater propulsion systems [2], flapping wings in micro air vehicles [3], atomic force microscopy [4] and energy harvesting via smart materials [5]. In these analyses, cantilevers are generally modelled as slender beams with thin rectangular cross sections. Under these assumptions, the beam is considered as infinitely thin in the direction of vibration and, for low flexural modes, the flow of the encompassing fluid is assumed to be two-dimensional in the beam cross-section plane (see [4,6,7]).

Most of the available analytical efforts are based on the classical works in Tuck [8] and Sader [4]. Therein, closed-form solutions based on linearized Navier–Stokes equations are derived for oscillating cantilevers of circular and thin rectangular cross section immersed in a viscous fluid. Moderate to large oscillations are traditionally described by using the so-called Keulegan–Carpenter (KC) number in conjunction with Morison’s formula for circular cross sections [9]. Numerical and experimental studies on cantilevers of rectangular cross sections are presented in Bidkar *et al.* [6] and Aureli & Porfiri [7].

In recent investigations [10,11], the lattice Boltzmann method (LBM) has been employed in the analysis of oscillating cantilevers. The LBM is a numerical method for computational fluid dynamics (CFD) that solves Boltzmann’s kinetic equation for the particle distribution functions (see [12,13]). The method is intrinsically mesoscopic, as opposed to conventional methods based on macroscopic conservation laws. The LBM is particularly effective in simulating fluid flows involving complicated boundaries and/or complex phenomena characterizing dissipative systems far from equilibrium (see [14,15]). Its distinctive features as a computational solver are: (i) the approach is intrinsically unsteady with particle distribution functions evolving towards a local equilibrium, which is a function of the macroscopic hydrodynamic field, and (ii) information travels along straight lines defined by the (constant) particle velocities associated with the lattice, rather than along space–time dependent material lines defined by the flow speed.

In Colosqui *et al.* [10], high-frequency oscillations of nanocantilevers are investigated. In these cases, Reynolds numbers are very low; namely, . As a consequence, the effects of fluid inertia are negligible and the force acting on the cantilever is largely due to the pressure. In Shi & Sader [11], a modified linear LBM is proposed to analyse harmonic oscillatory flows in the frequency domain. This work focuses on small displacements and low-frequency ranges, for which convective effects are irrelevant. For rectangular cantilevers, this corresponds to oscillatory Reynolds numbers of the order of –10.

In this paper, we use the LBM to study a rigid lamina undergoing transverse oscillations in a viscous fluid. We implement the LBM with the linear collision operator, that is, the lattice Bhatnagar–Gross–Krook operator (see [13]), in a range of oscillatory Reynolds numbers of approximately –200. Therefore, the focus of this study is the intermediate range of frequencies, when compared with the existing LBM literature. In addition, we consider finite oscillation amplitudes and consistently retain the full convective term in the fluid-governing equations. In this problem, vorticity convection and diffusion have a dramatic influence on the forces experienced by the lamina, resulting in nonlinear inertial and damping effects as demonstrated in Bidkar *et al.* [6] and Aureli & Porfiri [7]. We propose a novel refill procedure within the LBM simulation framework to accurately capture the flow physics in the vicinity of the lamina edges and properly describe the complex vorticity patterns in the flow. The main practical goal of this work is to estimate the force acting on the lamina and explore its dependence on the oscillation frequency and amplitude. We use the complex hydrodynamic function introduced in Sader [4] to describe the added mass and the hydrodynamic damping experienced by the oscillating body as a function of the Reynolds and KC numbers (see [4, 6–9]).

## 2. Problem statement and numerical method

### (a) Governing equations

We consider a rigid lamina of length *L* undergoing harmonic oscillations within a quiescent fluid environment. The system is described in a Cartesian reference frame, with basis unit vectors ** i** and

**along the**

*j**x*- and

*y*-axis, respectively. In the two-dimensional framework, the lamina is modelled as a line segment with zero thickness. The oscillation amplitude and radian frequency are denoted with

*A*and

*ω*, respectively. At time

*t*, the region occupied by the lamina is the domain . The lamina represents a moving boundary for the encompassing fluid that is assumed to be incompressible and Newtonian. Body forces are neglected and the Navier–Stokes equations governing the velocity field

**are 2.1where D**

*u***/D**

*u**t*is the material derivative of the velocity,

*ρ*is the fluid constant mass density and

**T**=−

*p*

**I**+2

*μ*

**D**is the Cauchy stress tensor. Here,

**I**is the identity tensor,

*p*and

*μ*denote the pressure and the dynamic viscosity, respectively, and

**D**=Sym∇

**is the symmetric part of the velocity gradient. Boundary conditions are imposed on the velocity field**

*u***in the form on and**

*u***→**

*u***0**as .

Non-dimensionalization of equation (2.1) and its boundary conditions shows that two parameters govern the flow characteristics (see also [16]). In this work, consistently with Aureli & Porfiri [7], we introduce the frequency parameter *β*=*ρωL*^{2}/(2*πμ*) and the non-dimensional amplitude of oscillation *ε*=*A*/*L* that is directly related to the KC number *κ*=2*πε*. Note that a pertinent oscillatory Reynolds number is defined accordingly as (see also [6]).

### (b) Lattice Boltzmann implementation

In this work, we consider a two-dimensional uniformly spaced lattice characterized by nine discrete directions (including the zero vector identifying particles at rest), namely D2 (see [13] and figure 1*a*). Physical lengths, time and mass density are non-dimensionalized with the lattice spacing *δ*_{L}, the simulation time step *δ*_{t} and the fluid mass density scale *δ*_{ρ} selected to be equal to the fluid mass density *ρ*, respectively. For a generic physical variable *a*, we let , where the underline bar denotes non-dimensional lattice quantities and *δ*_{a} is the corresponding dimensional scaling constant. The lattice Boltzmann equation (LBE) is as follows (see [12,13]):
2.2In equation (2.2), *f*_{i}, with *i*=0,…,8, is a probability density function, or LB population, and represents the probability of finding a fluid particle at site ** x** and time that is moving along the

*i*-th lattice direction with a discrete speed

*c*_{i}. The left-hand side of equation (2.2) represents the advection, or streaming, of fluid particles moving along the discrete directions of the lattice during the considered time step, whereas the right-hand side stands for particle collisions in the given lattice site

**, represented through a relaxation towards local (Maxwellian) equilibrium on a non-dimensional time scale**

*x**τ*=1/

*ϖ*. The discrete local equilibria in equation (2.2) are given by a second-order expansion in the Mach number of a local Maxwellian, 2.3where is the lattice speed of sound and

*w*

_{i}are the weights for the different lattice speeds (see [13]). Each computational time step, or cycle, includes one streaming and one collision process, providing the values of in equation (2.2). Once the LB populations are determined, the non-dimensional fluid density and flow velocity are given by and . A Chapman–Enskog analysis [18] of the LBE shows that, in the limit of long wavelengths and low frequency, the LBE reproduces exactly the Navier–Stokes equations for quasi-incompressible flows, with an ideal equation of state and a non-dimensional kinematic viscosity . Therefore, the incompressibility condition in equation (2.1) is only approximately satisfied.

The rows immediately above and below the lamina do not participate in this standard streaming process because the lamina acts as a wall for the particles’ motion. The motion of the oscillating lamina within the background lattice is reconstructed by means of the Filippova and Hänel enhanced procedure [17,19]. We consider the nodes and located on the rows immediately above and below the lamina, respectively (figure 1*b*). The location on the boundary denotes the intersections of the lamina with the lattice link between and . The length fraction of the intersected link between and is . The probability density functions of the nodes surrounding the lamina are recovered by using the linear interpolation proposed in Filippova & Hänel [19], that is,
2.4where is the post-collision population, subscript denotes the direction from a wall node to a fluid node , and subscripts *i* and identify opposite directions. Once the populations are computed according to equation (2.4), the standard collision process takes place for all the lattice sites in the computational domain. In equation (2.4), is a fictitious equilibrium density function given by
2.5In equations (2.4) and (2.5), the quantities *χ* and are determined as follows:
2.6
Note that the updating process is repeated for the nodes below the lamina by exchanging `b` with `f` and *i* with in equation (2.4). The lamina position and, consequently, the distances between the lamina and the two adjacent rows of nodes change with time; therefore, the value of Δ is updated at each time step.

Figure 1*a* illustrates the implemented computational domain and boundary conditions. To reduce the computational complexity, a symmetry boundary condition is imposed with respect to the oscillation axis. The symmetry boundary condition enforces a free-slip wall, that is, the normal component of distributions entering the boundary is reflected while the component parallel to the boundary is left unchanged. In Bidkar *et al.* [6] and Aureli & Porfiri [7], it is shown that vortex shedding at the sharp edges of the lamina has a dramatic influence on the global hydrodynamic field. Adequately capturing such structures requires that the lamina’s motion involves several lattice sites. To this aim, we propose an ad hoc devised refill procedure for the lattice sites overcome by the oscillating lamina. The procedure is based on the evidence that, in the parameter ranges explored in the present work, the local velocity of both the lamina and the fluid is very small. In particular, the maximum value of the lamina speed is in lattice units, which is much below the compressibility limit, that is, , where is the local Mach number. Therefore, we impose that the nodes `b` overcome by the oscillating lamina have the same velocity of the lamina itself and that they share the density of the nearest neighbours that reside on the same side of the lamina (figure 1*b*). Once the density and velocity fields are determined, the LB populations of such nodes are set to their equilibrium values stemming from equation (2.3). As a result, the lamina is not ‘transparent’ to the fluid, that is, no-penetration conditions are enforced (see [20]), and no actual ‘communication’ is performed between the fluid above and below the lamina, except for the interactions due to vortices at the lamina edges.

### (c) Force resultant evaluation

We evaluate the force resultant per unit length ** F** exerted by the fluid on the oscillating lamina through a control area approach as discussed in Prince

*et al.*[21]. On account of symmetry, the force resultant in the

**direction vanishes. Application of the divergence theorem on the integrated form of equation (2.1), along with the incompressibility condition, yields the following expression for the force resultant per unit length in the**

*i***direction 2.7where represents the control area surrounding the oscillating lamina, is its boundary, with unit normal vector**

*j***, d**

*n**l*represents the line element on and ⊗ denotes tensor product. Equation (2.7) is implemented in the LBM scheme by replacing integrals with summations involving lattice densities and velocities over pertinent lattice nodes and by replacing the stress tensor

**T**in equation (2.7) with its lattice counterpart [22] 2.8where is the non-equilibrium part of the density function, that is, . We select the lattice control area as the closest box surrounding the lamina (figure 1

*b*); therefore, prior to the use of symmetry conditions, its extent is given by nodes

^{2}, where is the lamina’s length in lattice units, that is,

*L*/

*δ*

_{L}. Note that, in the force estimation procedure, the control area occupies a fixed region of space and its position is only instantaneously updated when the lamina overcomes a row of nodes.

## 3. Results and discussion

### (a) Representative flow fields

We conduct LBM simulations in the parameter ranges *ε*={0.02,0.05,0.065,0.075,0.1} and *β*={20,50,100,150,200,250,300}. The dimensions of the computational domain are selected to guarantee an adequate resolution for the fluid field. Fine grid spacing, which is necessary to resolve vortex shedding from the lamina edges, and the large computational domain, which is crucial to minimize perturbations and reflections from the boundaries, are balanced with the computational effort. For a fixed pair (*ε*,*β*), we conduct a convergence study to determine a good compromise between resolution and computational time, by varying the oscillation amplitude in terms of spanned lattice sites and the domain size. This procedure implies that the lattice length of the lamina is varied to account for the varying value of *ε*.

We fix the oscillation amplitude to 15 lattice sites to accurately capture the vortices at the lamina edges. This selection implies the value for the length scale *δ*_{L}=*εL*/15. Additionally, to enhance numerical stability (see [13,17]), we let *τ*=3/4, thus yielding *ϖ*≈1.33 and . As a consequence, the time scale is determined as , to ensure Reynolds number similarity. Figure 2 shows a representative time history of the lamina position and the fluid force resultant, displaying the sensitivity of the force *F*_{y}, as computed via equation (2.7), to the domain size. For domains larger than nodes, force time histories converge satisfactorily to a common value. Therefore, we select a domain size of approximately nodes.

Figure 3 shows velocity magnitude, relative pressure and vorticity fields for two representative cases of comparatively small and large non-dimensional frequency of oscillation; namely, for *β*=100 and *β*=200. The value of the non-dimensional amplitude of oscillation is kept fixed to the intermediate value of *ε*=0.065. The displayed snapshots correspond to a time instant when the lamina velocity attains a maximum. The position of the lamina is visually identified by the pressure jump across and through the high-velocity and -vorticity regions in correspondence with its edges. The case *β*=200, associated with a larger value of *ω*, shows a significant region of the domain characterized by rather large velocities. Analogously, extended regions characterized by strong vorticity can be observed in proximity to the lamina edges for the case *β*=200. This is consistent with the observation that, for larger *β*, viscous effects play a minor role in the vorticity transport mechanisms. As a consequence, the vorticity shed at the lamina edges is primarily convected into the field rather than diffused (see also [7]). Additionally, the expected singular behaviour of the pressure field in the vicinity of the lamina edges is clearly visible for *β*=200.

### (b) Extraction of the hydrodynamic function

By following Sader [4] and Aureli & Porfiri [7], in the frequency domain, we express the (dimensional) force per unit length *F*_{y} in equation (2.7) as
3.1Here, the hydrodynamic function *Θ* is a complex function of *ε* and *β*, through the dependence on *ω*. In equation (3.1), a superimposed hat denotes a phasor quantity. For each simulation, that is, for each value of the pair (*ε*,*β*), the force time history is fitted with a sine model at radian frequency *ω* by using a nonlinear least-squares method. The radian frequency *ω* is directly assigned in the fitting model from the value of *β*. After discarding approximately one cycle of oscillation related to the initial transient, two cycles of oscillation are analysed to obtain amplitude and phase of the force. Higher harmonics in the force signals are not observed in the simulations; therefore, purely harmonic fits provide an *a posteriori* justification of equation (3.1) and yield satisfactory results in estimating the amplitude and phase shift of the force resultant with respect to the lamina displacement.

Figure 4*a*–*g* displays the numerical results for Re[*Θ*] and −Im[*Θ*], describing the added mass and hydrodynamic damping effects exerted by the fluid on the oscillating lamina. For a fixed value of *β*, the present findings partially support the affine dependence for the real and imaginary part of *Θ* on the value of *ε*, generally showing an increase with the oscillation amplitude, as demonstrated through the superimposed linear fits. Additionally, linear fits suggest that the proposed method can potentially capture the theoretical limit in Sader [4] as *ε*→0. On the other hand, slopes from the linear fits seem to overestimate the experimentally determined slopes reported in Aureli & Porfiri [7] (therein referred to as the function Δ(*β*)), resulting in a significantly larger prediction for the nonlinear damping, especially at large *β*. Conversely, due to the approximately constant value of Re[*Θ*] over the range of *ε* studied, numerical predictions suggest that nonlinear inertial effects are not as prominent as nonlinear damping. Consistent with the findings in Aureli & Porfiri [7], added mass and hydrodynamic damping are generally larger as *β* decreases. Overall, a satisfactory agreement is observed among the trends identified by theoretical predictions [4] and studies on moderately large amplitudes [6,7] and those found with the present LB simulations over the full range of *ε* and *β*.

## 4. Conclusions

In this work, we proposed a numerical study based on the LBM for the flow physics analysis of a quiescent viscous fluid that is perturbed by a rigid oscillating lamina. To the best of our knowledge, this work represents the first LBM analysis of oscillating laminae with characteristic Reynolds numbers in the range of –200 undergoing moderately large oscillations characterized by considerable vortex convection. The results of the direct numerical simulations are compared with experimental and numerical findings in the literature, showing a satisfactory correlation in terms of both fluid force magnitude and phase shift with the lamina acceleration.

## Acknowledgements

This research was supported by the National Science Foundation under grant nos CMMI-0745753 and CMMI-0926791, the Office of Naval Research under grant no. N00014-10-1-0988 with Dr Y.D.S. Rajapakse as the program manager, and the Honors Center of Italian Universities. Views expressed herein are those of the authors and not of the funding agencies. The authors would like to gratefully acknowledge Prof. Sauro Succi and Dr Alessandro De Maio for useful discussions. G.F. would like to thank the Department of Mechanical and Aerospace Engineering of the Polytechnic Institute of New York University for hosting him as a visiting professor while this work was developed.

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society