Transverse harmonic oscillations of laminae in viscous fluids: a lattice Boltzmann study

Giacomo Falcucci, Matteo Aureli, Stefano Ubertini, Maurizio Porfiri


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 Embedded Image are very low; namely, Embedded Image. 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 Embedded Image–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 Embedded Image–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, 69]).

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 j along the 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 Embedded Image. 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 u are Embedded Image2.1where Du/Dt is the material derivative of the velocity, ρ is the fluid constant mass density and T=−pI+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∇u is the symmetric part of the velocity gradient. Boundary conditions are imposed on the velocity field u in the form Embedded Image on Embedded Image and u0 as Embedded Image.

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 β=ρωL2/(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 Embedded Image (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 1a). 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 Embedded Image, 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]): Embedded Image2.2In equation (2.2), fi, 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 Embedded Image that is moving along the i-th lattice direction with a discrete speed ci. 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 x, represented through a relaxation towards local (Maxwellian) equilibrium on a non-dimensional time scale τ=1/ϖ. The discrete local equilibria Embedded Image in equation (2.2) are given by a second-order expansion in the Mach number of a local Maxwellian, Embedded Image2.3where Embedded Image is the lattice speed of sound and wi 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 Embedded Image in equation (2.2). Once the LB populations are determined, the non-dimensional fluid density and flow velocity are given by Embedded Image and Embedded Image. 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 Embedded Image and a non-dimensional kinematic viscosity Embedded Image. Therefore, the incompressibility condition in equation (2.1) is only approximately satisfied.

Figure 1.

(a) Sketch of the computational domain with the symmetry boundary condition. (b) Boundary treatment method, adapted from Mei et al. [17]; for a given node f, the lamina intersects the link connecting the node f with b. The intersection is represented by w. The node b, with neighbouring node Graphic, is refilled in terms of velocity and LB populations when it is overcome by the oscillating lamina.

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 Embedded Image and Embedded Image located on the rows immediately above and below the lamina, respectively (figure 1b). The location Embedded Image on the boundary denotes the intersections of the lamina with the lattice link between Embedded Image and Embedded Image. The length fraction of the intersected link between Embedded Image and Embedded Image is Embedded Image. 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, Embedded Image2.4where Embedded Image is the post-collision population, subscript Embedded Image denotes the direction from a wall node Embedded Image to a fluid node Embedded Image, and subscripts i and Embedded Image 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), Embedded Image is a fictitious equilibrium density function given by Embedded Image2.5In equations (2.4) and (2.5), the quantities χ and Embedded Image are determined as follows: Embedded Image2.6 Note that the updating process is repeated for the nodes below the lamina by exchanging b with f and i with Embedded Image 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 1a 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 Embedded Image in lattice units, which is much below the compressibility limit, that is, Embedded Image, where Embedded Image 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 Embedded Image that reside on the same side of the lamina (figure 1b). 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 i 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 j direction Embedded Image2.7where Embedded Image represents the control area surrounding the oscillating lamina, Embedded Image is its boundary, with unit normal vector n, dl represents the line element on Embedded Image and ⊗ denotes tensor product. Equation (2.7) is implemented in the LBM scheme by replacing integrals with summations involving lattice densities Embedded Image and velocities Embedded Image over pertinent lattice nodes and by replacing the stress tensor T in equation (2.7) with its lattice counterpart Embedded Image [22] Embedded Image2.8where Embedded Image is the non-equilibrium part of the density function, that is, Embedded Image. We select the lattice control area Embedded Image as the closest box surrounding the lamina (figure 1b); therefore, prior to the use of symmetry conditions, its extent is given by Embedded Image nodes2, where Embedded Image 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 Embedded Image 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 Embedded Image. As a consequence, the time scale is determined as Embedded Image, 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 Fy, as computed via equation (2.7), to the domain size. For domains larger than Embedded Image nodes, force time histories converge satisfactorily to a common value. Therefore, we select a domain size of approximately Embedded Image nodes.

Figure 2.

Convergence study on the influence of the domain size. Force per unit length and displacement are given in dimensional units, for the case L=10.2 mm, ρ=1000 kg m−3 and ν= 10−6 m2 s−1, as in Aureli & Porfiri [7]. The thick line represents the position of the lamina. Size of the domain: (a) Graphic, black solid line; (b) Graphic, black dashed line; (c) Graphic, black dashed-dotted line; (d) Graphic, red solid line; (e) Graphic, red dashed line; (f) Graphic, red dashed-dotted line; (g) Graphic, blue solid line; and (h) Graphic, blue dashed line. Transient effects are damped out after the first oscillation cycle. (Online version in colour.)

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 Embedded Image 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.

Figure 3.

Representative flow fields for ε=0.065 in the case (ac) Graphic and (df) β=200 Graphic. Quantities are given in dimensional units, as in figure 2. Velocity magnitude is scaled in the range [−15,15] mm s−1 in (a) and (d); relative pressure is scaled in the range [0.98,1.02]×105 Pa in (b) and (e); and out-of-plane vorticity is scaled in the range [−50,50] s−1 in (c) and (f). Vorticity patterns are highly comparable with those observed in the flow visualization experiments reported in Aureli & Porfiri [7] for the cases Graphic and Graphic. Red and blue areas correspond to maximum and minimum values, respectively. (Online version in colour.)

(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 Fy in equation (2.7) as Embedded Image3.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 4ag 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 β.

Figure 4.

Comparison of Re[Θ] (black markers) and −Im[Θ] (red markers) for the hydrodynamic function as determined with the present LBM simulations (filled circles), representative results (open markers) in Aureli & Porfiri [7], theoretical predictions (asterisks) from Sader [4] at ε=0, and numerical findings (plus symbols) from Bidkar et al. [6]. Dashed lines in panels (ag) are linear fits of LBM simulation results, whose slopes (filled inverted triangles) and intercepts (right-facing triangles) are shown in (h). Data in (i) and (j) are presented in aggregated form. Representative data from Aureli & Porfiri [7] in (i) and (j) refer to: β=27 (open squares, κ=2πε=0.5,0.9,1.5), β=60 (diamonds, κ=2πε=0.4,0.7,1.1) and β=165 (open triangles, κ=2πε=0.3,0.5,0.7). Data shown are computed by adopting the affine dependence of Θ on the KC number proposed in Aureli & Porfiri [7] and using the experimentally determined slopes Δ(β) reported therein and the observed values of κ at the tip of the oscillating beam. (Online version in colour.)

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 Embedded Image–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.


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.



View Abstract