## Abstract

In this paper, we study flexural vibrations of two thin beams that are coupled through an otherwise quiescent viscous fluid. While most of the research has focused on isolated beams immersed in placid fluids, inertial and viscous hydrodynamic coupling is ubiquitous across a multitude of engineering and natural systems comprising arrays of flexible structures. In these cases, the distributed hydrodynamic loading experienced by each oscillating structure is not only related to its absolute motion but is also influenced by its relative motion with respect to the neighbouring structures. Here, we focus on linear vibrations of two identical beams for low Knudsen, Keulegan–Carpenter and squeeze numbers. Thus, we describe the fluid flow using unsteady Stokes hydrodynamics and we propose a boundary integral formulation to compute pertinent hydrodynamic functions to study the fluid effect. We validate the proposed theoretical approach through experiments on centimetre-size compliant cantilevers that are subjected to underwater base-excitation. We consider different geometric arrangements, beam interdistances and excitation frequencies to ascertain the model accuracy in terms of the relevant non-dimensional parameters.

## 1. Introduction

The accurate estimation of the hydrodynamic loading on flexible structures vibrating in liquids is of fundamental importance in the study of microelectromechanical systems, such as atomic force microscope probes [1–4] and radio frequency switches [5–9]; sensors and energy harvesting devices [10–16]; miniature propulsion systems across biology [17,18] and engineering [19–23]; and naval structures [24–26].

A particularly relevant problem entails the analysis of flexural vibrations of slender thin beams immersed in otherwise placid viscous incompressible fluids. For low Keulegan–Carpenter numbers, convective phenomena can be discarded and the distributed hydrodynamic loading on the structure is described through a hydrodynamic function that depends on the sole Reynolds number. Such function is estimated from the two-dimensional unsteady Stokes hydrodynamics of a thin lamina oscillating in the fluid. While the time-dependent Stokes equation is linear, it accounts for both the fluid-added mass and viscous damping. The real and imaginary parts of the hydrodynamic function are related thus to the added mass and the viscous damping produced by the encompassing fluid [27,28]. The role of the Keulegan–Carpenter number on the hydrodynamic function is discussed in [29–33], while compressibility effects are studied in [34].

The effect of the proximity to a solid surface can be described within the same theoretical framework, by accounting for the dependence of the hydrodynamic function on the ratio between the gap and the beam width. A boundary integral formulation for the two-dimensional unsteady Stokes hydrodynamics of a rigid body vibrating in the vicinity of solid surfaces is presented in [35]. Tractable semianalytical expressions for the complex hydrodynamic function of thin beams for selected gap-to-width ratios are discussed in [36] to analyse vibrations of flexible cantilevers near solid surfaces. Validation of this formulation against three-dimensional fluid-structure simulations and experimental results are presented in [37], along with complete expressions for the hydrodynamic function, in terms of both the oscillatory Reynolds number and the gap-to-width ratio. A thorough analysis of the modification of the beam modal properties owing to the encompassing fluid is presented in [38], and additional results on three-dimensional effects, stochastic vibrations and large-amplitude oscillations are reported in [39–42].

While most of the research on beams vibrating in fluids is focused on isolated structures, fluid-induced coupling in oscillating beams or plates is ubiquitous across a multitude of engineering and natural systems. Examples of such systems include arrays of piezoelectric fans [43,44], systems of ionic polymer metal composites [45], the so-called piezoelectric grass for energy harvesting from turbulence-induced vibrations [46], comb-drives of micropumps for biomedical applications [47], viscosimeters based on piezoelectric cantilevers or clamped–clamped membranes [48] and nanomats of carbon nanotubes [49,50]. Examples of natural systems where multiple structures interact in a liquid include the antennulary setal hairs of copepods, which are mechanoreceptive sensors used to find food, avoid predators and track mates [51,52], and the lateral line of fish, which is used to detect the presence of fixed obstacles or other motile bodies [53,54].

In this paper, we study steady-state small amplitude vibrations of two identical thin beams that are immersed in an otherwise quiescent, unbounded, viscous and incompressible fluid at a constant gap. We model the fluid–structure interaction through a matrix of hydrodynamic functions that incorporates the hydrodynamic coupling between the two beams. On account of the symmetry of the problem, the matrix has two independent elements, which are controlled by the gap-to-width ratio and the oscillatory Reynolds number, similar to the problem of a beam vibrating close to a solid surface [35–37]. These two hydrodynamic functions are identified from the analysis of the two-dimensional flow physics generated by a rigid lamina, oscillating in the vicinity of a fixed identical blade, as the key non-dimensional variables, the gap-to-width ratio and the frequency parameter, are parametrically changed. In this context, we adapt the boundary integral formulation of [35] to generate semianalytical expressions for the hydrodynamic functions. Such tractable expressions are, in turn, used to formulate a reduced order modal model for steady-state vibrations of the two coupled beams. To validate the proposed theoretical approach, we perform experiments on centimetre-size compliant beams undergoing underwater base excitation. Experimental results are compared with model predictions for the following cases: (i) both beams are subjected to the same base-excitation; (ii) the base excitations of the two beams have identical amplitude but are out of phase; and (iii) only one beam undergoes base excitation while the other is passive.

The paper is organized as follows. In §2, we present the problem statement and the governing equations. In §3, we introduce the two-dimensional fluid problem, illustrate the boundary integral formulation and establish semianalytical formulae for the hydrodynamic functions. In §4, we formulate a reduced order modal model for the analysis of the coupled beams. In §§5 and 6, we report on the validation of the proposed theoretical model through comparison with experimental results. Specifically, the experimental set-up and the results are described in §5, while a detailed discussion is presented in §6. Conclusion of this study is summarized in §7. Verification of the numerical procedure by comparison with computational fluid dynamics (CFD) simulations is offered in appendix A. In appendix B, we briefly present the mechanism that is designed to elicit out of phase base excitation for the coupled beams.

## 2. Problem statement and governing equations

We consider small amplitude oscillations of two identical cantilevered thin beams of rectangular cross section, spaced by a constant gap *g*, and immersed in an otherwise quiescent viscous fluid. The beams are modelled by using the classical linear Euler–Bernoulli beam theory [55] and the fluid is assumed to be Newtonian and incompressible. The beams have uniform cross section and are composed of the same homogeneous and isotropic material. We focus on the general case in which the cantilevers vibrate under harmonic base excitations of different amplitude and phase.

The geometry of the problem is displayed in figure 1 along with the pertinent nomenclature. In particular, we label the beams through the index *j*=1,2. We identify with *x* the common coordinate along the beams' axes, with *y* and *z* the coordinates along the width and thickness directions, respectively, and with *t* the time variable. The origin of the coordinate system is set at the centre of the left end of beam 1 in the rest position. We denote with *b* and *h* the beams' width and thickness, respectively, with *E* the Young's modulus of the beams' material, and with *ρ*_{s} the mass density of the beams' material per unit volume. The bending stiffness and mass per unit length are consequently computed as *K*=*Ebh*^{3}/12 and *ϱ*_{s}=*ρ*_{s}*bh*. (Note that *E* should be considered as the modulus for plane strain conditions, that is, the Young's modulus divided by (1−*ν*^{2}), where *ν* is Poisson's ratio because the beam cross section is very thin [56,57].) We use **w**(*x*,*t*) and **W**(*t*) to indicate the vector of the beams' elastic deflections [*w*_{1}(*x*,*t*),*w*_{2}(*x*,*t*)]^{T} and the vector of the base excitations [*W*_{1}(*t*),*W*_{2}(*t*)]^{T}, where T indicates matrix transposition. Thus, the beams' vibrations are governed by the following system of two partial differential equations
2.1where the vectors **H**(*x*,*t*) and **s**(*x*,*t*) represent the load exerted by the encompassing fluid on the beams and the structural damping, respectively.

The boundary conditions are
2.2aand
2.2bwhere *L* is the common length of the beams. Equations (2.2a) are used to block the displacement and rotation at the base, whereas equations (2.2b) are used to set to zero the bending moment and shear force at the cantilvers’ tip.

We hypothesize that in steady-state conditions, each of the two beams vibrates harmonically under the effect of the imposed base excitation , where *ω* is the radian frequency of vibration and *A*_{Wj}(*ω*) and *α*_{Wj}(*ω*) are the amplitude and the phase, respectively. By denoting phasors with a superimposed hat, we have and , with , Im[⋅] being the imaginary part (a similar notation is used for the real part, that is, Re[⋅]) and *j*=1,2.

We assume that *ω* is sufficiently small so that the beams vibrate primarily along the lowest mode shapes, which are characterized by a spatial wavelength considerably larger than the width. Moreover, we hypothesize that *h*≪*b* so that each cross section can be proxied by a lamina. In this case, three-dimensional phenomena associated with variations of the flow physics along the beams' axes can be neglected, and the hydrodynamic forcing vector **H**(*x*,*t*) is estimated by studying the two-dimensional flow induced by the motion of two identical rigid laminae; see also the detailed analysis on isolated beams in [41].

Through a standard non-dimensional analysis of the incompressible Navier–Stokes equations in the absence of bulk loads [58] and recalling that we are here interested in small amplitude vibrations, we find that **H**(*x*,*t*) can be written in the frequency domain as
2.3where *ρ*_{f} is the fluid mass per unit volume, is the phasor of the total deflection at *x* of the *j*th beam, and is a 2×2 matrix, which we term the hydrodynamic matrix. Such matrix describes the fluid–structure interaction and depends on two non-dimensional parameters, the former associated with the oscillatory Reynolds number [35–37] and the latter related to the geometry of the two-dimensional problem. These parameters are the gap-to-width ratio *γ* and the frequency parameter *β*, defined by
2.4aand
2.4bwhere *μ* is the dynamic viscosity of the fluid. Note that the oscillatory Reynolds number can be computed from the frequency parameter through *Re*=(*π*/2)*β* [35–37].

Different from related studies on the dynamics of cilia [17], we consider the effect of fluid inertia, while neglecting the presence of a mean flow and non-local fluid coupling. In a similar vein, our work differs from recent studies on shear flow of elastic rods on a substrate [49,50], whereby we focus on small structural vibrations in quiescent fluids, for which the flow field in a generic plane orthogonal to the *x*-axis is controlled by the in-plane motion of the local cross sections and long-range three-dimensional hydrodynamic coupling is negligible. While these three-dimensional phenomena are naturally treated within slender-body theories [17,49,50], the hypothesis of local fluid coupling is central to our approach, which posits a decomposition of the original three-dimensional problem into a two-dimensional flow analysis, in the plane of the cross sections, and a one-dimensional structural problem, along the beams' span. This decomposition allows for the independent analysis of structural vibration problems, without the need of further flow physics simulations, through the mere implementation of the hydrodynamic functions that are developed as part of this study. Based on the analysis in [42], we expect equation (2.3) to be accurate for oscillations as large as the width in the case of *β* of the order of one. As *β* increases, the advection term in the Navier–Stokes equation increases and the range of validity of the approximation is expected to drastically decrease, whereby only oscillations of the order of 10^{−2} *b* should be resolved for *β* approaching the limit value of 10^{4} considered in this study.

We select a classical frequency-independent hysteretic model for the structural damping *s*_{j}(*x*,*t*), so that the bending stiffness of each beam is replaced by *K*(1+i*η*) in the frequency domain, where *η* is the damping coefficient [55]. By rewriting equation (2.1) in the frequency domain and using equation (2.3), we obtain the following system of two linear partial differential equations, describing the beams' harmonic vibrations:
2.5where is the so-called mass ratio and *I*_{2} is the 2×2 identity matrix. Note that the mass ratio corresponds to the ratio between the mass per unit length of a fluid cylinder of diameter equal to *b* and the mass per unit length of the beam.

## 3. Analysis of hydrodynamic interactions

### (a) Properties of the hydrodynamic matrix

The *jj*th entry of the hydrodynamic matrix describes the hydrodynamic loading on the *j*th lamina as the other lamina is held fixed. Similarly, the *ij*th entry of , with *i*≠*j*, identifies the hydrodynamic loading on the *i*th lamina corresponding to the motion of the *j*th lamina. As the two laminae are identical, we have
3.1aand
3.1bwhere *Γ*(*β*,*γ*) and *Γ*_{c}(*β*,*γ*) are two independent complex hydrodynamic functions. The hydrodynamic matrix can be consequently rewritten as
3.2

### (b) Boundary integral formulation

We calculate the two hydrodynamic functions *Γ* and *Γ*_{c} using a boundary integral approach similar to the one proposed in [27]. We specifically adapt the formulation developed in [35] for studying the hydrodynamic loading exerted on an infinitely long rigid cylinder of arbitrary cross section in the presence of a solid surface. A similar approach is followed in [59], where hydrodynamic coupling between vibrating cantilevers is studied for a different geometry, in which the beams are placed side by side.

As mentioned earlier, the fluid is Newtonian and incompressible and we consider the body forces to be negligible. The moving laminae are regarded as rigid bodies, harmonically oscillating in the *z*-direction at the radian frequency *ω*, with complex amplitudes and . We focus on unsteady Stokes flow, for which both these amplitudes are considerably smaller than any other geometric scale of the problem, that is, *b* and *g*. Therefore, the fluid flow is the solution of [58]
3.3aand
3.3bwhere **u**(*y*,*z*,*t*)=[*u*_{y}(*y*,*z*,*t*),*u*_{z}(*y*,*z*,*t*)]^{T} is the velocity field, *p*(*y*,*z*,*t*) is the fluid pressure and ∇ and *Δ* are the Nabla and Laplace operators, respectively. No-slip boundary conditions on the velocity field are imposed at the two fluid–lamina interfaces and the fluid is assumed to be at rest at infinity, where the reference pressure is set to zero.

Following [27,35], we introduce the stream function *ψ*(*y*,*z*,*t*), such that ∂*ψ*/∂*z*=*u*_{y} and −∂*ψ*/∂*y*=*u*_{z}, to identically satisfy equation (3.3b). The vorticity component along the *x*-axis is consequently computed through
3.4

Computing the curl of equation (3.3a), we find the following equation for the vorticity field:
3.5which, accounting for equation (3.4), yields a fourth-order partial differential equation for *ψ*(*y*,*z*,*t*), namely
3.6In the frequency domain, equations (3.5) and (3.6) are written as
3.7aand
3.7bwhere *σ*^{2}=i*ωρ*_{f}/*μ*. By evaluating the divergence of equation (3.3a), we also find the following condition for the pressure field:
3.8

Application of Green's integral theorem allows for obtaining the following integral representation of [27]
3.9where *Ω*(*y*,*z*,*y*′,*z*′)=−(2*π*)^{−1}*K*_{0}(*σR*), *Ψ*(*y*,*z*,*y*′,*z*′) and are Green's functions for equation (3.7a), equation (3.7b) and equation (3.8), respectively. Here, *Ψ*(*y*,*z*,*y*′,*z*′)=*σ*^{−2}(*Ω*(*y*,*z*,*y*′,*z*′)−*G*(*y*,*z*,*y*′,*z*′)), and *K*_{0}(⋅) is the modified Bessel function of the third kind [60]. Moreover, ∂(⋅)/∂*l* and ∂(⋅)/∂*n* denote the directional derivatives along the tangent and normal vectors to the fluid boundary , respectively (figure 1); the latter being directed outward from the fluid region.

The fluid boundary depicted in figure 1 comprises (i) the boundaries and of the two laminae and (ii) the closing contour at infinity. The latter does not contribute to equation (3.9) owing to the fact that the fluid is at rest at infinity. In addition, the contribution from the boundaries of the cross sections corresponds to a simple integration along the *y*-axis as the beams are assumed to be infinitely thin. Therefore, equation (3.9) reduces to
3.10where and are the vorticity and pressure jumps across the *j*th lamina, respectively. We remark that the first two terms in the integral of equation (3.9) are zero, as and are continuous across both the rigid laminae [27,35].

By differentiating equation (3.10) with respect to *z*′ and *y*′, we calculate *u*_{y}(*y*′,*z*′,*ω*) and *u*_{z}(*y*′,*z*′,*ω*). By evaluating the velocity fields at *z*′=0 and *z*′=*g* and imposing no-slip boundary conditions, we obtain the following system of four coupled integral equations:
3.11a
3.11b
3.11c
and
3.11dNote that in the derivation of equation (3.11), we have used the fact that ∂^{2}*Ψ*(*y*,0,*y*′,0)/∂*z*′∂*y*, ∂^{2}*Ψ*(*y*,0,*y*′,0)/∂*y*′∂*z*, ∂^{2}*Ψ*(*y*,*g*,*y*′,*g*)/∂*z*′∂*y* and ∂^{2}*Ψ*(*y*,*g*,*y*′,*g*)/∂*y*′∂*z* are all equal to zero [35].

We solve equation (3.11) by following the numerical implementation reported in [27,35,37]. Briefly, the system of integral equations is transformed into a linear system for the nodal unknown jumps of the pressure and vorticity across the two laminae through a suitable discretization of the integrals. Here, we use a total of 500 elements per lamina and we compute the requisite integrals by using closed form formulae, when available [35], or 10 Gaussian quadrature points. We note that the coefficient matrix of the linear system, related to Green's function *Ψ* in the right-hand sides of equation (3.11), depends only on *ω* and *g*. On the other hand, the amplitude and phase of the two imposed vibrations influence only the vector of constants of the linear system. Thus, the hydrodynamic matrix is expected to share similarities with *Ψ* in its frequency dependence.

For brevity, we define the amplitude ratio as
3.12varying in the range [0,1]. With respect to the phase difference, we consistently define
3.13spanning the range [0,2*π*]. In figures 2 and 3, we report the differences in the pressure and vorticity across the two laminae for , *α*_{2}=0, *γ*=0.3, *β*=100, and different values of the amplitude ratio (*ξ*=0, 1, and 0.5) and phase difference (*δα*=0 and *π*). The lamina width is set to 10 mm and the encompassing fluid is considered as common water at room temperature, for which *ρ*_{f}=998.2 kg m^{−3} and *μ*=1.003×10^{−3} Pa s (the same water properties are used throughout the study).

For *ξ*=1, we find that and are identical in the case of in phase vibration (figure 2), while they are opposite in value when the oscillations are out of phase (figure 3). In the latter case, the pressure field is symmetric with respect to *z*=*g*/2, in contrast with the former case, where it is skew-symmetric. Notably, the vorticity jumps mirror this behaviour, whereby and are identical for out of phase vibrations (figure 3) and, vice versa, have opposite values for in phase oscillations (figure 2). Such symmetries are lost for *ξ*=0 and 0.5. We note that for a single lamina vibrating in a viscous fluid, the vorticity jump would be identically zero [27] in agreement with our observations for *δα*=0 for two laminae moving as one.

Further details on the two-dimensional flow physics can be found in appendix A, where the results obtained through the proposed boundary integral formulation are verified against CFD simulations.

### (c) Hydrodynamic functions

Once the pressure jumps across the two laminae are calculated, the hydrodynamic forces per unit beam length are evaluated as 3.14aand 3.14bSuch hydrodynamic loads are expressed in terms of the displacements of the two laminae by accounting for the hydrodynamic functions in equations (2.3) and (3.2), that is, we write 3.15aand 3.15b

The hydrodynamic functions *Γ* and *Γ*_{c} are conveniently estimated by solving equation (3.11) with and a selected, possibly unitary, amplitude of vibration for the second lamina through
3.16aand
3.16bFor each pair of gap-to-width ratio *γ* and frequency parameter *β*, we thus calculate the corresponding values of the two hydrodynamic functions *Γ* and *Γ*_{c}. We execute a total of 336 simulations, by considering 21 values for *β* spanning the range from 1 to 10 000 and 16 values of *γ* in the range from 0.01 to 10, as reported in figures 4 and 5.

Semianalytical formulae are obtained by fitting numerical results following the methodology presented in [37]. For example, to estimate *Γ*, we use
3.17where , , and *a*_{qr} are the complex valued fit coefficients. A similar notation is used for *Γ*_{c}. The results of the fits are listed in table 1 for and , which offer a good compromise between computational complexity and accuracy. The coefficients of determination are equal to 0.99994 and 0.99985 for the real and imaginary parts of *Γ* and 0.99941 and 0.99994 for the real and imaginary parts of *Γ*_{c}, respectively.

Figures 4 and 5 display the semianalytical expressions of *Γ* and *Γ*_{c} against numerical results from boundary element analysis. As expected, for increasing values of the separation, each lamina behaves as if it were isolated and unaffected by the presence of the other one. Therefore, , as and , where is the semianalytical formula for a single oscillating lamina in an unbounded fluid proposed in [30]. On the other hand, for and small values of *β*, the hydrodynamic damping approaches the classical lubrication theory expression, whereby Im[*Γ*(*β*,*γ*)]=Im[*Γ*_{c}(*β*,*γ*)]=−2/(*π*^{2}*βγ*^{3}). Thus, reducing the gap between the beams drastically increases the hydrodynamic damping along with the added mass effect, as displayed in figures 4 and 5.

For simplicity, we split *Γ* into three contributions
3.18where the last two terms vanish for increasing separations. The hydrodynamic loadings on the two laminae can be rewritten as
3.19aand
3.19bEquation (3.19) demonstrates that the force on each lamina is the sum of two terms, which are controlled by the absolute motion of the lamina with respect to the fluid and the relative motion with respect to the other lamina. Thus, the correction *δΓ* quantifies the effect of the hydrodynamic coupling from the other lamina on the individual added mass and viscous damping. Specifically, it quantifies the effect of the other lamina if it were moving at the same amplitude and phase, that is, if the two laminae were moving as one. Notably, such correction depends on both *γ* and *β*, displaying interesting behaviours similar to those identified in [35] for a single lamina vibrating close to a solid surface. For example, the real part of *δΓ* can change its sign as a function of *γ* and *β* indicating the possibility of increasing or decreasing the added mass effect.

We finally comment that, different from the case of an isolated beam, the slope of the imaginary parts of the hydrodynamic functions with respect to variations in *β* varies from half a decade per decade to one decade per decade, in agreement with [35].

## 4. Analysis of beams' vibration

We solve the system of two coupled equations in equation (2.5) through the Galerkin method, whereby we first diagonalize equation (2.5) to obtain two decoupled boundary value problems, and then we project the corresponding displacement fields following [29,30,42,61]. For convenience, we rewrite equation (2.5) in the following compact form
4.1where and is given in equation (3.2). Next, we solve the associated eigenvalue problem
4.2for the eigenvalue *λ* and the eigenvector **v**. By inspection of the matrix , we find that [1,1]^{T} and [1,−1]^{T} are eigenvectors with eigenvalues equal to
4.3aand
4.3bNotably, [1,1]^{T} corresponds to the symmetric vibration mode, in which both beams vibrate in phase with the same magnitude, and [1,−1]^{T} is associated with the skew-symmetric vibration mode, in which the beams vibrate out of phase with the same magnitude. Thus, we introduce the modal matrix
4.4which allows for transforming the beam deflection fields to the corresponding modal coordinates , that is,
4.5

In terms of the modal coordinates, equation (4.1) becomes
4.6where and are the *j*th component of and , respectively. Note that the boundary conditions for the modal coordinates are identical to those of the original displacement fields in equation (2.2).

To solve equation (4.6), we project the modal displacement fields on the in-vacuum flexural modes of an undamped cantilever beam of length *L*, see for example [55]. Such projection is implemented by writing the deflection field in modal coordinates as
4.7where are the complex modal coefficients for the *j*th beam. By substituting this expansion into equation (4.6), multiplying by *ψ*_{i}(*x*), integrating from 0 to *L*, and accounting for the orthogonality of the mode shapes, we find
4.8where *i*=1,2,…, is the *i*th in-vacuum resonance radian frequency, and *χ*_{i} is the *i*th eigenvalue of the characteristic equation . By replacing the modal coefficients in equation (4.7) and using equation (4.5), we compute the displacement field of each beam. In practice, we truncate the series in equation (4.7) to a finite number of *N*=10 summands.

## 5. Experimental set-up and procedures

We validate the proposed theoretical model by comparing predictions with experimental results, obtained by performing three different types of laboratory experiments on two centimetre-size compliant beams. Such beams vibrate underwater at room temperature under the effect of base excitation and their spacing is parametrically varied. Experiments are performed on Mylar cantilever beams of free length of vibration *L*=80 mm, width *b*=10 mm, thickness *h*=150 μm and mass per unit length *ϱ*_{s}=1.75 g m^{−1}.

### (a) Experimental set-up

The experimental set-up, displayed in figure 6, is composed of a Bruel & Kjaer 4810 mini shaker with a 2718 power amplifier, an Agilent 33210 A function generator and a Dantec Dynamics Phantom Speed Sense 9040 high-speed camera.

We consider three different experimental conditions (figure 6): (i) ‘in phase’ (IP), in which both beams are clamped to the shaker through a rigid fixture; (ii) ‘out of phase’ (OP), where the beams are excited by a dedicated mechanism described in appendix B which allows to obtain out of phase base excitations with identical amplitude; and (iii) ‘active passive’ (AP), where only one beam is activated through base excitation while the other is fixed to a stationary clamp and passively vibrates owing to the fluid coupling. In each condition, the beams are submerged in a water tank of approximately 40 l capacity. Three different spacing distances are considered in each experiment described by the gap-to-width ratios *γ*=0.3, 1 and 3. Thus, we consider a total of nine independent experiments.

To maintain a uniform sampling density, the camera frame rate of acquisition is set to 50 Hz for experiments with base excitation of up to 2 Hz. The frame rate is increased to 100 Hz for frequencies of excitation between 2 and 5 Hz, to 200 Hz for excitations between 5 and 10 Hz, and eventually to 400 Hz for frequencies larger than 10 Hz. Xcitex ProAnalyst software (Xcitex, Inc., http://www.xcitex.com) is used for video postprocessing to extract the position of the two beam bases *W*_{j}(*t*) and tips *W*_{j}(*t*)+*w*_{j}(*L*,*t*), *j*=1,2, at each frame. Such quantities are in turn used to compute the corresponding transfer functions ; further technical details on this procedure can be found in [30]. Note that the base excitation of the second beam is used to scale the elastic deflection of both beams and standardize the whole set of experimental conditions, which include a case where the first beam is held stationary.

### (b) In-air experiments

In-air experiments are initially carried out to identify Mylar properties. In this case, each beam is separately tested using the same set-up presented in figure 6 along with the acquisition procedure detailed above. In particular, the bending stiffness *K* and the damping coefficient *η* of the two cantilevers are obtained by comparing the experimentally measured transfer functions with model predictions obtained by replacing **H**=**0** in equation (2.1). In this case, we obtain a system of two uncoupled identical equations, each describing the in-vacuum vibrations of the corresponding beam.

Figure 7 displays experimental results against model predictions with *K*=10.07×10^{−6} *N* m^{2} and *η*=0.04. The frequency responses reported therein refer to the ratio between the phasor of the elastic deflection at the tip and the corresponding base excitation, that is, . We note that the resulting Young's modulus is *E*=3.58 G Pa. These experiments are executed by considering a nominal base excitation of *A*_{Wj}/*b*=0.01, *j*=1,2, obtained by selecting peak-to-peak shaker input voltages equal to 0.4 V.

### (c) Underwater experiments

We perform nine independent experiments in a rather broad range of variation of the frequency parameter *β*∈[20,1000] to elicit both the symmetric and skew-symmetric vibration modes of the coupled beams. For the IP experiment (*ξ*(0,*ω*)=1 and *δα*(0,*ω*)=0), only symmetric modes are expected to be excited; by contrast, only skew-symmetric modes are anticipated to be excited in the OP experiment (*ξ*(0,*ω*)=1,*δα*(0,*ω*)=*π*). In the AP experiment, both the mode types should contribute to the system response. IP experiments are performed by considering a nominal base excitation of *A*_{W1}/*b*=*A*_{W2}/*b*=0.01, corresponding to peak-to-peak input voltages to the shaker of 0.4 V. For AP experiments, we consistently set a nominal base excitation of *A*_{W2}/*b*=0.01, which corresponds to the same peak-to-peak input to the shaker. On the other hand, for the OP experiment a peak-to-peak voltage input of 2.0 V is needed to obtain the same level of base excitation, owing to the presence of the mechanism interposed between the shaker and the beams.

## 6. Results and discussion

A comparison between experimental results and theoretical predictions is synoptically presented in figures 8, 9 and 10, which, respectively, address *γ*=0.3,1 and 3. The first inspection of the results depicted therein indicates that predictions of the proposed modelling framework are in close agreement with experimental results across the entire range of variation of *γ* and *β* and for the three selected experimental conditions, that is, IP, OP and AP. Notably, the model is successful in accurately predicting both the resonance frequencies of the structures as well as the associated quality factors. We comment that, in agreement with figure 5, the presence of the second beam produces an observable difference between figures 8 and 9, while modest variations are evidenced by comparing figures 9 and 10.

To offer a quantitative assessment of the model accuracy, we adopt the normalized root mean square error *E*_{NRMS} defined in [30], namely,
6.1where *n*_{p} is the number of the data points and and are the experimentally measured and the theoretically estimated amplitude of each response at the *i*th data point, respectively. As indicated in the captions of figures 8–10, the discrepancy between experimental and theoretical results is generally of the order of 10%, with smaller errors in the condition IP and larger deviations for condition AP. We anticipate the source of such discrepancy to relate with the smaller motions that are elicited in the AP condition, and especially for larger gaps. For OP experiments, we also offer that a possible source of uncertainty should be ascribed to the frictions in the mechanism along with imperfections that cause the separation between the two beams to slightly vary along their span. The agreement between theoretical predictions and experimental findings is expected to be reduced for larger levels of base excitation and higher frequencies, for which advective and three-dimensional effects become significant [32]. Wall effects are also expected to play a role if the structures were excited in the vicinity of the tank walls, especially at larger amplitude and frequencies [42].

Comparing the frequency responses of the beams across the three selected gap-to-width ratios, we find that increasing *γ* differentially affects the response of the beams for IP and OP conditions. In the former case, increasing *γ* reduces both the resonance frequencies and the quality factors of the beams. On the other hand, for the OP condition, the opposite holds, whereby increasing *γ* increases both the resonance frequencies and quality factors of the fluid-coupled beams. To further elucidate this phenomenon, in figure 11, we report the resonance frequencies and quality factors computed from the proposed theoretical framework for conditions IP and OP, in which symmetric or skew-symmetric modes are individually excited. These quantities are estimated from the complex poles *ω*_{11},*ω*_{21},*ω*_{12},*ω*_{22},…, which are calculated from equation (6.2), that is,
6.2Specifically, the resonance frequency and the quality factor *Q*_{ji} for *i*=1,2,… and *j*=1,2 are computed as [55,62]
6.3aand
6.3bwhere *ω*_{ji} is given by equation (6.2).

## 7. Conclusion

In this paper, we studied the steady-state small amplitude vibrations of two identical thin beams that are immersed in an otherwise quiescent, unbounded, viscous and incompressible fluid. A matrix of hydrodynamic functions was introduced to model the fluid–structure interaction. Such matrix was identified from the analysis of the two-dimensional flow physics generated by two rigid laminae, as the key non-dimensional variables, that is, the gap-to-width ratio and the frequency parameter, are parametrically varied. We proposed a boundary integral formulation to describe the two-dimensional flow generated by the motion of the laminae and establishing, in turn, tractable semianalytical formulae for the hydrodynamic matrix.

The analysis showed that the distributed hydrodynamic loading experienced by each beam is not only related to its absolute motion, as in the case of isolated structures, but is also affected by the relative motion with respect to the other beam. Specifically, the added mass effect is magnified for decreasing gaps and reduced as the oscillatory Reynolds number increases. Similarly, hydrodynamic damping decreases as the gap and the oscillatory Reynolds number increase.

The obtained semianalytical expressions were used to formulate a reduced order modal model for the steady-state vibrations of the two coupled beams. To validate the proposed theoretical approach, we performed experiments on centimetre-size compliant beams undergoing base excitation underwater. Three different experiments were considered: (i) both beams were subjected to the same base excitation; (ii) the base excitations of the two beams had identical amplitude but were out of phase; and (iii) only one beam was subjected to base excitation while the other was passive. The gap was parametrically varied in each experiment to offer a comprehensive set of conditions for model validation. Theoretical predictions were found to be in good agreement with experimental findings in a wide frequency range, where the fundamental symmetric and skew-symmetric vibration modes were excited.

The formulation presented in this paper is expected to find application in the analysis of arrays of sharp-edged slender beams vibrating in viscous fluids. The proposed framework is anticipated to find use in the design of MEMS-based devices operating in liquids, including novel viscosimeters relying on multiple sensing microelements [48], arrays of piezoelectric fans [43,44], systems of ionic polymer metal composites [45] and energy harvesters based on coupled active beams [46]. Such methodology may also contribute to improve our understanding of natural systems comprising arrays of coupled slender structures, such as mechanoreceptive setal hairs of copepods [51,52] or fish lateral lines [53,54].

## Funding statement

This material is based upon work partially supported by the National Science Foundation under grant no. 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 (H2CU) through a scholarship to Carmela Intartaglia.

## Acknowledgements

Views expressed herein are those of the authors and not of the funding agencies.

## Appendix A. Computational fluid dynamics simulations

Here, we analyse the two-dimensional flow induced by the harmonic oscillations of two fluid-coupled identical laminae through CFD simulations with the twofold aim of investigating the accuracy of unsteady Stokes hydrodynamics and verifying the proposed boundary integral formulation.

Non-dimensional parameters explored in the CFD analysis target the range of moderate to high oscillatory Reynolds numbers, that is, *β*∈[30,3000]. Moreover, the amplitude ratio *ξ* spans the whole range [0,1] along with the phase difference *δα* that varies in [0,2*π*]. For brevity, we focus on the case *γ*=0.3 and we consider 45 simulations with , *ξ*={0,0.5,1} and *δα*={0,*π*/2,*π*,3*π*/2}. In each simulation, the laminae start from rest, respectively, at *z*=0 and *z*=*g*, , and simulations are executed for 3.3 complete cycles to achieve steady state. The physical parameters are identical to those considered in the boundary integral formulation, that is, *b*=10 mm, *ρ*_{f}=998.2 kg m^{−3} and *μ*=1.003×10^{−3} Pa s.

We retain the computational model, domain and solution strategy used in [29,30,42,61], by simply replacing the single lamina studied therein with two laminae. CFD analysis provides pressure and velocity fields in the computational domain and the time histories of the forces exerted by the fluid on the two laminae. In particular, the components along the *z*-direction are used to extract the hydrodynamic functions *Γ*(*β*,*γ*) and *Γ*_{c}(*β*,*γ*). In the case of small amplitude oscillations considered herein, these time histories are harmonic signals at the same frequency *ω* of the imposed displacement. In our analyses, we always discard the initial transient corresponding to the first cycle to extract the requisite phasors.

As an illustrative example of the flow fields obtained through CFD, in figure 12, we display the velocity magnitude, pressure and vorticity for the cases *β*=100, *γ*=0.3 and *δα*=0 and *π*, at different time instants in the oscillation period T.

**(a) Details on mesh refinement and implementation**

We model the ideal case of an infinitely thin lamina by selecting the thickness *h* to be equal to 10^{−2} *b*. Accounting for the problem symmetry with respect to the *z*-axis, only one half of the whole domain is simulated and the selected dimensions are 10 *b*×(20+*γ*) *b* to minimize boundary effects. No-slip boundary conditions are enforced at the fluid–lamina interfaces and at the right edge of the computational domain, a symmetry condition is set along the left edge, while outflow boundary conditions are imposed at the top and bottom edges, as depicted in figure 13.

To obtain a good compromise between computational effort and the accuracy required to simulate the complex flow physics in proximity of the two vibrating laminae, the computational domain is split into six regions with varying mesh size. The regions around the laminae, called CV-1, and the seven dynamic mesh buffer zones, denoted with DMB, are discretized by using structured quad-map regions, with cell sizes of 2×10^{−3} *b*. The other zones, called CV-2, CV-3, CV-4 and Bulk, are discretized by adopting a cell size of, respectively, 9×10^{−3} *b*, 25×10^{−3} *b*, 50×10^{−3} *b* and 80×10^{−3} *b*. To impose the continuity of field variables across these regions, the so-called interior and interface boundary conditions are set along each of the edges. The motion of the laminae in the fluid is simulated by implementing a dynamic meshing technique in the CV-1 regions. This allows for a mesh updating depending on the current position of each of the two laminae. The technique relies on the usage of an ad hoc-developed user-defined function (udf) that is based on a macro obtainable in the Ansys Fluent library (Ansys, Inc., Ansys Fluent 13.0. http://www.ansys.com).

The pressure-based transient laminar solver is used in this study. Pressure and velocity are coupled by using the SIMPLE scheme, while pressure, momentum and spatial derivatives of field variables are calculated by using PRESTO, QUICK and least squares cell-based schemes, respectively. The time evolution is computed by using a first-order implicit transient formulation. We select decreasing values of the time step for decreasing oscillation frequencies to guarantee the presence of at least 1000 samples per cycle. Furthermore, each time step value is chosen so that the ratio between the displacement of nodes on the laminae within a time step and the cell height is always smaller than 0.2 to avoid multiple splitting/collapsing of the cell during mesh updating (table 2).

**(b) Comparison with boundary integral approach**

In figures 14 and 15, the fluid loadings acting on the two laminae during the oscillation are reported for in phase and out of phase vibrations, respectively. Specifically, the dimensionless quantities , *j*=1,2, are displayed in the rows (*a*,*b*) of each figure for *ξ*=0, 0.5 and 1. Circular markers are obtained by using the proposed boundary integral formulation. Squared markers indicate results from CFD simulations. Continuous lines show the prediction of semianalytical expressions in equation (3.17).

We find an excellent agreement between values calculated with different numerical techniques. Specifically, the maximum difference between boundary integral formulation and CFD results, calculated as the relative error with respect to the latter ones, is lower than 4% in all the cases.

## Appendix B. Design of the mechanism for out of phase excitation

In figure 16, the schematic drawing of the planar mechanism designed for performing OP experiments is illustrated. Through this mechanical device, we impose to the two beams oscillations of equal amplitude with a phase shift of *π*.

The shaker sinusoidal motion is transmitted to an input slider to which it is connected by a bracket joint. The slider, moving on a linear guide, is connected by a pin to two rods of each 200 mm length, symmetrically positioned with respect to the shaker axis, forming an angle of *π*/4 with its direction. The other ends of the rods are connected by pins to two sliders, both coupled with a second linear guide that is perpendicular to the other one. The motions of these two sliders have the same amplitude, but opposite direction. As the beams have to be placed in close proximity, these sliders are connected to other two output sliders through two other links. Briefly, when the input slider oscillates along the *y*-axis in figure 16, the two output sliders move along the *z*-axis, in phase opposition. The set-up is located on a rigid plate, which is positioned on the tank through a frame. The shaker is fixed to a rigid bracket, anchored to the floor. The distortion of the output signal, evaluated as the relative error with respect to the input, is theoretically lower than 0.1% during the oscillation.

- Received June 14, 2013.
- Accepted November 5, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.