# Vibration Analysis of Euler-Bernoulli Beams Partially Immersed in a Viscous Fluid

Authors: Wafik Abassi; Adil El Baroudi; Fulgence Razafimahery
Journal: Physics Research International
 (2016)

---

## Abstract

The vibrational characteristics of a microbeam are well known to strongly depend on the fluid in which the beam is immersed. In this paper, we present a detailed theoretical study of the modal analysis of microbeams partially immersed in a viscous fluid. A fixed-free microbeam vibrating in a viscous fluid is modeled using the Euler-Bernoulli equation for the beams. The unsteady Stokes equations are solved using a Helmholtz decomposition technique in a two-dimensional plane containing the microbeams cross sections. The symbolic software Mathematica is used in order to find the coupled vibration frequencies of beams with two portions. The frequency equation is deduced and analytically solved. The finite element method using Comsol Multiphysics software results is compared with present method for validation and an acceptable match between them was obtained. In the eigenanalysis, the frequency equation is generated by satisfying all boundary conditions. It is shown that the present formulation is an appropriate and new approach to tackle the problem with good accuracy.


---

## Body

## 1. Introduction

The objective of this paper is to provide an analytical method to calculate the coupled frequencies of vibration of microbeams partially immersed in a viscous fluid. The microbeams are clamped on one edge while the other edge is free.The motivation of this work is to provide a theoretical model that can be used in the design and interpretation of density and viscosity sensors.Due to their size and potential for highly sensitive and low cost compact device applications, microstructures are becoming increasingly attractive for sensing applications and have been studied extensively in recent years. Microstructures are commonly used in atomic force microscopy (AFM) to probe surface properties and to measure interfacial forces [1–7], in biological and chemical sensors [8–10]. A precise modeling of the solid-fluid interaction and the determination of the frequency response enable the measurement of the density and the rheological behavior of fluids [11–16]. Reference [16] uses finite element analysis (FEA) method in order to predict the dynamic response of the cantilever beam. This method can be easily applied to the measurement of the fluid viscosity.For microstructures, fluid viscosity can greatly affect their frequency response. Reference [7] presented a rigorous theoretical model for the frequency response of cantilever beams that are undergoing flexural vibrations and immersed in viscous fluids, which is of particular relevance to applications of the AFM. The knowledge and understanding of the frequency analysis of microbeams are of fundamental practical importance in application to the AFM.The frequency analysis of a microbeam can be dramatically affected by the properties of the fluid in which it is immersed. Whereas calculation of the natural frequencies in vacuum can be performed routinely, analysis of the effects of immersion in fluid poses a formidable challenge. The modal response of an immersed microbeam can be considerably affected by the properties of fluid. The added mass effect due to the fluid structure interaction can, however, cause considerable variations in natural frequencies. The knowledge and understanding of this viscous fluid-structure coupling are lacking at present.In this contribution, we investigate the vibrational behavior of microbeams partially immersed in a viscous fluid, which describes the interrelation between the fluid’s density and viscosity. For a viscous fluid problem, the analytical formulation is based upon a convenient decomposition of the velocity field into two contributions, one being related to the scalar potential and the other being the vector potential. The solutions of the differential equations of motion turn out to be complex and can be conveniently treated with the aid of the symbolic software Mathematica. Furthermore, this work investigates the influence of the fluid’s viscosity on the vibrational behavior of the microbeams.


## 2. Modal Analysis of Beams and Frequency Equation

Modal analysis of elastic immersed structures is needed in every modern construction and should have wide engineering application. In this study, modal analysis is important to predict the dynamic behavior of the submerged beams. It is well known that the natural frequencies of the submerged elastic structures are different from those in vacuum. The effect of fluid forces on the immersed beam decreases the natural frequencies from those that would be measured in the vacuum. This decrease in the natural frequencies is caused by the increase of the kinetic energy of the fluid-beams system without a corresponding increase in strain energy. The Euler-Bernoulli beam is partially immersed inside rectangular fluid domain (Figure1). Consider a beam of length 3.06·10−2 [m], width 4.6·10−3 [m], and thickness 1.27·10−4 [m] as shown in Figure 1, which corresponds to a model developed in [16]. The interaction between the fluid and the Euler-Bernoulli beams is taken into account to calculate the natural frequencies and mode shapes of the coupled system. The dynamics of each beam portion are treated separately. It is assumed that the beam has aligned neutral axis.Figure 1
Sketch of beam composed of 2 uniform beam segments (denoted by(
1
) and (
2
)) partially immersed in a fluid-filled rectangular geometry.


### 2.1. Mathematical Formulation

In this section we present the general theory for the dynamic deflection of beams partially submerged in a viscous fluid. A schematic depiction of beams partially submerged in viscous fluid is displayed in Figure1. We begin by discussing some general assumptions and approximations taken into account in the present theoretical model. It is assumed that the cross section of the beams is uniform over its entire length and the length of the beams greatly exceeds its width. Also the beams are an isotropic linearly elastic solid and internal frictional effects are negligible. The amplitude of the vibrations of the beams is small. In addition, we shall neglect all torsional effects in the beams and only consider the flexural modes of vibration and we shall consider modes whose motion is strict in the x-direction. For beams vibrating at small amplitudes (small compared to the beams dimensions), the governing dynamical equations ignoring shear deformation and rotary inertia effects for the transverse deflection u
i
(
y
,
t
) (i
=
1 for submerged portion and i
=
2 for portion in vacuum) of uniform elastic beams can be written in the form(1)
EI
∂
4
u
1
y
,
t
∂
y
4
+
ρ
s
S
∂
2
u
1
y
,
t
∂
t
2
=
f
f
l
u
i
d
,
L
2
⩽
y
⩽
L
,
(2)
EI
∂
4
u
2
y
,
t
∂
y
4
+
ρ
s
S
∂
2
u
2
y
,
t
∂
t
2
=
0
,
L
⩽
y
⩽
2
L
,where u
1
(
y
,
t
) and u
2
(
y
,
t
) are the lateral deflections at distance y (spatial coordinate) along the length of the beams and t is time; EI, ρ
s, and S are the flexural rigidity, the mass per unit volume, and the cross-sectional area of the beam, respectively. f
f
l
u
i
d is the external force per unit length acting on the beam in the direction of the flexural displacement, which is caused by the viscous fluid in the beams. The dynamics of each beam portion are treated separately.For a beam moving in a viscous fluid, the applied loadf
f
l
u
i
d can be obtained by integrating the normal component of the total force exerted by the fluid over the beam section. To proceed with the analysis, the general form of f
f
l
u
i
d is required. We therefore examine the equations of motion for the fluid. For this, the viscous fluid mode shapes are first computed for the square domain fluid assuming the boundary conditions given in Figure 1. From conservations of mass and momentum, the motion of the fluid is governed by (see [17])(3)
∇
·
v
=
0
,
(4)
ρ
∂
v
∂
t
+
ρ
v
·
∇
v
=
∇
·
σ
,in which v
=
v
x
,
v
y
T is the fluid velocity vector, ρ is the density of the fluid, and σ is the total fluid stress tensor (pressure and viscous forces). Since the amplitude of vibration of each beam is small compared to its cross-sectional dimensions, it then follows that all nonlinear convective inertial effects in the fluid can be neglected, and the hydrodynamic loading on the beam will be a linear function of its displacement. This implies that the fluid dynamics can be modeled as an unsteady linear Stokes flow. Assuming that the fluid is Newtonian, its constitutive equation is given by(5)
σ
=
τ
-
p
I
with
τ
=
η
∇
v
+
∇
v
T
,where I is a unit tensor, η is the dynamic viscosity, and p is the fluid pressure. The boundary conditions which define the fluid domain are(6)
v
·
n
=
0
on
Γ
,
(7)
v
=
0
on
Γ
0
,
(8)
σ
·
n
=
0
on
Γ
s
.The equations of motion (3) and (4) are highly complex and coupled. However, a simpler set of equations can be obtained by introducing scalar potentials ϕ and ψ, known as the Helmholtz decomposition [18] in a way which permits easily transferring the vector problem formulation ((3) and (4)) to scalar problem formulation. In flow fields, the velocity is thereby decomposed into a potential flow and a viscous flow. In other words, the velocity v can be expressed as a sum of the gradient of a scalar potential ϕ and the curl of a vector potential Ψ as follows [18]:(9)
v
=
∇
ϕ
+
∇
×
Ψ
,where Ψ is a vector stream function. Using the problem symmetry, the vector potential Ψ reduces to a scalar equation; that is, Ψ
=
(
0,0
,
ψ
), using the condition ∇
·
Ψ
=
0, and substituting the above resolutions into (3) and (4), after some manipulations, the equation for the conservation of mass (3) becomes the Laplace equation(10)
∇
2
ϕ
=
0and the equation for the conservation of momentum (4) becomes(11)
∇
2
ψ
-
1
ν
∂
ψ
∂
t
=
0
,
p
+
ρ
∂
ϕ
∂
t
=
0
,where ν is the kinematic fluid viscosity and ∇
2
=
∂
2
/
∂
x
2
+
∂
2
/
∂
y
2 is the Laplacian operator. Thus, the Stokes equation is reduced to formulation (10) and (11). Moreover, adopting the Cartesian coordinate system (
x
,
y
) and using expression (9), the velocity components, v
=
(
v
x
,
v
y
), may be expressed as simple functions of compressional and shear wave potentials in the form(12)
v
x
=
∂
ϕ
∂
x
+
∂
ψ
∂
y
,
v
y
=
∂
ϕ
∂
y
-
∂
ψ
∂
xand the pertinent stress-velocity relations are(13)
σ
x
x
=
ρ
∂
ϕ
∂
t
+
2
η
∂
v
x
∂
x
,
σ
x
y
=
η
∂
v
x
∂
y
+
∂
v
y
∂
x
,
σ
y
y
=
ρ
∂
ϕ
∂
t
+
2
η
∂
v
y
∂
y
.


#### 2.1.1. Field Expansions

Consider time harmonic motion throughout with angular frequencyω and with the e
x
p
(
i
ω
t
) dependence suppressed for simplicity. Applying the classical technique of separation of variables in the Cartesian coordinates, the solution of (10) and (11) after some manipulations, and taking into account boundary condition (6) can be shown to be(14)
ψ
=
s
i
n
β
x
A
cos
⁡
δ
y
+
B
sin
⁡
δ
y
,
(15)
ϕ
=
cos
⁡
β
x
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(16)
p
=
η
α
2
cos
⁡
β
x
C
cos
⁡
λ
y
+
D
sin
⁡
λ
ywith the complex-valued coefficients (17)
λ
=
j
β
,
δ
=
j
β
2
-
α
2
,
α
=
1
-
j
2
ω
ν
,where β
=
n
π
/
L and n is an integer. A, B, C, and D are unknown coefficients which will be determined later by imposing the appropriate boundary conditions. Hence, the fluid particle velocity field can be determined from (12) in the form(18)
v
x
=
sin
⁡
β
x
δ
B
cos
⁡
δ
y
-
A
sin
⁡
δ
y
-
β
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(19)
v
y
=
cos
⁡
β
x
λ
D
cos
⁡
λ
y
-
C
sin
⁡
λ
y
-
β
A
cos
⁡
δ
y
+
B
sin
⁡
δ
y
.Consequently, direct substitution of the expansions (14)–(19) into the stress-velocity relations (13), after some manipulations, leads to(20)
σ
x
x
=
η
cos
⁡
β
x
2
β
δ
B
cos
⁡
δ
y
-
A
sin
⁡
δ
y
-
2
β
2
+
α
2
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(21)
σ
y
y
=
η
cos
⁡
β
x
2
β
δ
A
sin
⁡
δ
y
-
B
cos
⁡
δ
y
-
2
λ
2
+
α
2
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(22)
σ
x
y
=
η
sin
⁡
β
x
2
β
2
-
α
2
A
cos
⁡
δ
y
+
B
sin
⁡
δ
y
+
2
β
λ
C
sin
⁡
λ
y
-
D
cos
⁡
λ
y
.In order to determine the unknown constants, the appropriate boundary conditions (7) and (8) must be explicitly used. Substituting the total fluid stress tensor components (20)–(22) into boundary conditions equations (8) and velocity components (18)-(19) into boundary conditions equations (7), and after some considerable algebraic manipulation, the scalar potentials ϕ and ψ become(23)
ψ
=
A
sin
⁡
β
x
cos
⁡
δ
y
+
β
2
δ
Φ
sin
⁡
δ
y
,
(24)
ϕ
=
A
cos
⁡
β
x
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,where (25)
Φ
=
2
δ
sin
⁡
δ
L
+
2
β
2
-
α
2
sin
⁡
λ
L
/
λ
2
β
2
cos
⁡
δ
L
-
2
β
2
-
α
2
cos
⁡
λ
L
.Using the expressions of ψ and ϕ from (23) and (24), in (12), we have the following expressions for the velocity components:(26)
v
x
=
A
sin
⁡
β
x
δ
β
2
δ
Φ
cos
⁡
δ
y
-
sin
⁡
δ
y
-
β
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,
(27)
v
y
=
A
cos
⁡
β
x
λ
β
λ
cos
⁡
λ
y
-
β
Φ
sin
⁡
λ
y
-
β
cos
⁡
δ
y
+
β
2
δ
Φ
sin
⁡
δ
y
.Components of the stress tensor can be obtained by substituting the expressions of velocity components v
x; v
y, from (26) and (27), and the expression of scalar potential ϕ from (24) in (13) which comes out as (28)
σ
x
x
=
A
η
cos
⁡
β
x
2
β
δ
β
2
δ
Φ
cos
⁡
δ
y
-
sin
⁡
δ
y
-
2
β
2
+
α
2
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,
(29)
σ
y
y
=
A
η
cos
⁡
β
x
2
β
δ
sin
⁡
δ
y
-
β
2
δ
Φ
cos
⁡
δ
y
-
2
λ
2
+
α
2
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,
(30)
σ
x
y
=
A
η
sin
⁡
β
x
2
β
2
-
α
2
cos
⁡
δ
y
+
β
2
δ
Φ
sin
⁡
δ
y
+
2
β
λ
β
Φ
sin
⁡
λ
y
-
β
λ
cos
⁡
λ
y
.Now we can integrate the normal component of the total force exerted by the fluid over the beam section in order to obtain the external force per unit length acting on the beam in the direction of the flexural displacement f
f
l
u
i
d. Hence, the deflection of the immersed part of the beam (which is equal to u
1
(
y
,
t
)
=
u
1
y
exp
⁡
(
j
ω
t
)) can be determined from (1) and (28) in the form(31)
u
1
′
′
′
′
y
-
ρ
s
S
ω
2
EI
u
1
y
=
2
A
β
η
cos
⁡
β
L
/
2
EI
β
2
Φ
cos
⁡
δ
y
-
δ
sin
⁡
δ
y
-
β
2
+
α
2
2
Φ
cos
⁡
λ
y
+
sin
⁡
λ
y
λ
,where primes denote differentiation with respect to the position variable y.


### 2.2. Method of Solution and Frequency Equation

The general solutions of the ordinary differential equations (31) and (2) for the beams system, as shown in Figure 1, can be written in different segments in terms of trigonometric functions and represent propagating waves and hyperbolic functions represent evanescent waves as(32)
u
1
y
=
A
1
cos
⁡
Ω
y
+
B
1
sin
⁡
Ω
y
+
C
1
cosh
⁡
Ω
y
+
D
1
sinh
⁡
Ω
y
+
2
A
β
η
cos
⁡
β
L
/
2
EI
δ
sin
⁡
δ
y
-
β
2
Φ
cos
⁡
δ
y
Ω
4
-
δ
4
+
β
2
+
α
2
/
2
Ω
4
-
λ
4
Φ
cos
⁡
λ
y
+
sin
⁡
λ
y
λ
,
(33)
u
2
y
=
A
2
cos
⁡
Ω
y
+
B
2
sin
⁡
Ω
y
+
C
2
cosh
⁡
Ω
y
+
D
2
sinh
⁡
Ω
y
,where Ω is the flexural wavenumber and is given by(34)
Ω
=
ω
ρ
s
S
EI
1
/
4
.Equation (32) shows that the parameters (35)
ω
c
=
n
π
L
2
EI
ρ
s
S
,
ω
c
=
n
π
/
L
2
ρ
s
S
/
EI
+
1
/
νare resonance frequency parameters and can be plotted for different values of n.The nine constantsA
1, B
1, C
1, D
1, A
2, B
2, C
2, D
2, and A can be found by imposing the following boundary conditions. Equations (32) and (33) are the general solution for the vibration modes of beams partially immersed in fluid. In the case of a coupled system, the effect of fluid viscosity on the flexible beams must be considered. On the fluid-beams interface, the normal velocity must be continuous. Therefore, the fluid velocity and the deflection of a beam u
1 satisfy the relation (36)
v
x
x
=
L
/
2
=
j
ω
u
1
y
,
L
2
⩽
y
⩽
L
.Substituting (26) and (32) into (36), and after some manipulations, leads to the following equation:(37)
A
1
+
2
λ
ω
η
EI
Ω
4
-
δ
4
tan
⁡
β
L
/
2
β
2
Φ
cos
⁡
δ
y
-
δ
sin
⁡
δ
y
-
A
β
2
+
2
λ
ω
η
β
2
+
α
2
/
2
EI
Ω
4
-
λ
4
tan
⁡
β
L
/
2
Φ
cos
⁡
λ
y
+
sin
⁡
λ
y
λ
=
j
ω
sin
⁡
β
L
/
2
A
1
cos
⁡
Ω
y
+
B
1
sin
⁡
Ω
y
+
j
ω
sin
⁡
β
L
/
2
C
1
cosh
⁡
Ω
y
+
D
1
sinh
⁡
Ω
y
.Both sides of (37) are integrated over L
/
2
<
y
<
L to yield the following equation:(38)
A
=
j
ω
Υ
sin
⁡
β
L
/
2
A
1
I
1
+
B
1
I
2
+
C
1
I
3
+
D
1
I
4
,where Υ and the coefficients I
1, I
2, I
3, and I
4 are written down explicitly in Appendix E. Now the expression of the lateral movement of the submerged beams in contact with the fluid u
1 can be formulated using (32) and taking into account (38) as(39)
u
1
y
=
A
1
cos
⁡
Ω
y
+
I
1
U
y
+
B
1
sin
⁡
Ω
y
+
I
2
U
y
+
C
1
cosh
⁡
Ω
y
+
I
3
U
y
+
D
1
sinh
⁡
Ω
y
+
I
4
U
y
,where (40)
U
y
=
2
λ
η
ω
EI
Υ
tan
⁡
β
L
/
2
δ
sin
⁡
δ
y
-
β
2
Φ
cos
⁡
δ
y
Ω
4
-
δ
4
+
β
2
+
α
2
/
2
Ω
4
-
λ
4
Φ
cos
⁡
λ
y
+
sin
⁡
λ
y
λ
.Now, to derive the frequency equation of beams partially immersed in fluid, one assumes that, on the interface between each portion of the beam at the position y
=
L, the deflection, the rotation angle, the internal shear force, and bending moment of the beam must be continuous. This is satisfied when(41)
u
1
L
=
u
2
L
,
u
1
′
L
=
u
2
′
L
,
u
1
′
′
L
=
u
2
′
′
L
,
u
1
′
′
′
L
=
u
2
′
′
′
L
.To complete the formulation of the boundary-value problem, the four boundary conditions for the beam ends (y
=
L
/
2
=
a and y
=
2
L
=
b) considered in this work are specified as follows:(42)
u
1
′
′
a
=
0
,
u
1
′
′
′
a
=
0
,
u
2
b
=
0
,
u
2
′
b
=
0
.Combining the boundary conditions ((41) and (42)) with (33) and (39) yields the following linear homogeneous system of eight equations:(43)
M
X
=
0
,
X
=
A
1
,
B
1
,
C
1
,
D
1
,
A
2
,
B
2
,
C
2
,
D
2
T
,where M is a 8
×
8 matrix whose elements are designated as M
i
j. This system can have nontrivial solutions only if the determinant of the matrix M is zero, leading to the frequency equation. The matrix M is written down explicitly in Appendix F.


## 2.1. Mathematical Formulation

In this section we present the general theory for the dynamic deflection of beams partially submerged in a viscous fluid. A schematic depiction of beams partially submerged in viscous fluid is displayed in Figure1. We begin by discussing some general assumptions and approximations taken into account in the present theoretical model. It is assumed that the cross section of the beams is uniform over its entire length and the length of the beams greatly exceeds its width. Also the beams are an isotropic linearly elastic solid and internal frictional effects are negligible. The amplitude of the vibrations of the beams is small. In addition, we shall neglect all torsional effects in the beams and only consider the flexural modes of vibration and we shall consider modes whose motion is strict in the x-direction. For beams vibrating at small amplitudes (small compared to the beams dimensions), the governing dynamical equations ignoring shear deformation and rotary inertia effects for the transverse deflection u
i
(
y
,
t
) (i
=
1 for submerged portion and i
=
2 for portion in vacuum) of uniform elastic beams can be written in the form(1)
EI
∂
4
u
1
y
,
t
∂
y
4
+
ρ
s
S
∂
2
u
1
y
,
t
∂
t
2
=
f
f
l
u
i
d
,
L
2
⩽
y
⩽
L
,
(2)
EI
∂
4
u
2
y
,
t
∂
y
4
+
ρ
s
S
∂
2
u
2
y
,
t
∂
t
2
=
0
,
L
⩽
y
⩽
2
L
,where u
1
(
y
,
t
) and u
2
(
y
,
t
) are the lateral deflections at distance y (spatial coordinate) along the length of the beams and t is time; EI, ρ
s, and S are the flexural rigidity, the mass per unit volume, and the cross-sectional area of the beam, respectively. f
f
l
u
i
d is the external force per unit length acting on the beam in the direction of the flexural displacement, which is caused by the viscous fluid in the beams. The dynamics of each beam portion are treated separately.For a beam moving in a viscous fluid, the applied loadf
f
l
u
i
d can be obtained by integrating the normal component of the total force exerted by the fluid over the beam section. To proceed with the analysis, the general form of f
f
l
u
i
d is required. We therefore examine the equations of motion for the fluid. For this, the viscous fluid mode shapes are first computed for the square domain fluid assuming the boundary conditions given in Figure 1. From conservations of mass and momentum, the motion of the fluid is governed by (see [17])(3)
∇
·
v
=
0
,
(4)
ρ
∂
v
∂
t
+
ρ
v
·
∇
v
=
∇
·
σ
,in which v
=
v
x
,
v
y
T is the fluid velocity vector, ρ is the density of the fluid, and σ is the total fluid stress tensor (pressure and viscous forces). Since the amplitude of vibration of each beam is small compared to its cross-sectional dimensions, it then follows that all nonlinear convective inertial effects in the fluid can be neglected, and the hydrodynamic loading on the beam will be a linear function of its displacement. This implies that the fluid dynamics can be modeled as an unsteady linear Stokes flow. Assuming that the fluid is Newtonian, its constitutive equation is given by(5)
σ
=
τ
-
p
I
with
τ
=
η
∇
v
+
∇
v
T
,where I is a unit tensor, η is the dynamic viscosity, and p is the fluid pressure. The boundary conditions which define the fluid domain are(6)
v
·
n
=
0
on
Γ
,
(7)
v
=
0
on
Γ
0
,
(8)
σ
·
n
=
0
on
Γ
s
.The equations of motion (3) and (4) are highly complex and coupled. However, a simpler set of equations can be obtained by introducing scalar potentials ϕ and ψ, known as the Helmholtz decomposition [18] in a way which permits easily transferring the vector problem formulation ((3) and (4)) to scalar problem formulation. In flow fields, the velocity is thereby decomposed into a potential flow and a viscous flow. In other words, the velocity v can be expressed as a sum of the gradient of a scalar potential ϕ and the curl of a vector potential Ψ as follows [18]:(9)
v
=
∇
ϕ
+
∇
×
Ψ
,where Ψ is a vector stream function. Using the problem symmetry, the vector potential Ψ reduces to a scalar equation; that is, Ψ
=
(
0,0
,
ψ
), using the condition ∇
·
Ψ
=
0, and substituting the above resolutions into (3) and (4), after some manipulations, the equation for the conservation of mass (3) becomes the Laplace equation(10)
∇
2
ϕ
=
0and the equation for the conservation of momentum (4) becomes(11)
∇
2
ψ
-
1
ν
∂
ψ
∂
t
=
0
,
p
+
ρ
∂
ϕ
∂
t
=
0
,where ν is the kinematic fluid viscosity and ∇
2
=
∂
2
/
∂
x
2
+
∂
2
/
∂
y
2 is the Laplacian operator. Thus, the Stokes equation is reduced to formulation (10) and (11). Moreover, adopting the Cartesian coordinate system (
x
,
y
) and using expression (9), the velocity components, v
=
(
v
x
,
v
y
), may be expressed as simple functions of compressional and shear wave potentials in the form(12)
v
x
=
∂
ϕ
∂
x
+
∂
ψ
∂
y
,
v
y
=
∂
ϕ
∂
y
-
∂
ψ
∂
xand the pertinent stress-velocity relations are(13)
σ
x
x
=
ρ
∂
ϕ
∂
t
+
2
η
∂
v
x
∂
x
,
σ
x
y
=
η
∂
v
x
∂
y
+
∂
v
y
∂
x
,
σ
y
y
=
ρ
∂
ϕ
∂
t
+
2
η
∂
v
y
∂
y
.


### 2.1.1. Field Expansions

Consider time harmonic motion throughout with angular frequencyω and with the e
x
p
(
i
ω
t
) dependence suppressed for simplicity. Applying the classical technique of separation of variables in the Cartesian coordinates, the solution of (10) and (11) after some manipulations, and taking into account boundary condition (6) can be shown to be(14)
ψ
=
s
i
n
β
x
A
cos
⁡
δ
y
+
B
sin
⁡
δ
y
,
(15)
ϕ
=
cos
⁡
β
x
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(16)
p
=
η
α
2
cos
⁡
β
x
C
cos
⁡
λ
y
+
D
sin
⁡
λ
ywith the complex-valued coefficients (17)
λ
=
j
β
,
δ
=
j
β
2
-
α
2
,
α
=
1
-
j
2
ω
ν
,where β
=
n
π
/
L and n is an integer. A, B, C, and D are unknown coefficients which will be determined later by imposing the appropriate boundary conditions. Hence, the fluid particle velocity field can be determined from (12) in the form(18)
v
x
=
sin
⁡
β
x
δ
B
cos
⁡
δ
y
-
A
sin
⁡
δ
y
-
β
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(19)
v
y
=
cos
⁡
β
x
λ
D
cos
⁡
λ
y
-
C
sin
⁡
λ
y
-
β
A
cos
⁡
δ
y
+
B
sin
⁡
δ
y
.Consequently, direct substitution of the expansions (14)–(19) into the stress-velocity relations (13), after some manipulations, leads to(20)
σ
x
x
=
η
cos
⁡
β
x
2
β
δ
B
cos
⁡
δ
y
-
A
sin
⁡
δ
y
-
2
β
2
+
α
2
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(21)
σ
y
y
=
η
cos
⁡
β
x
2
β
δ
A
sin
⁡
δ
y
-
B
cos
⁡
δ
y
-
2
λ
2
+
α
2
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(22)
σ
x
y
=
η
sin
⁡
β
x
2
β
2
-
α
2
A
cos
⁡
δ
y
+
B
sin
⁡
δ
y
+
2
β
λ
C
sin
⁡
λ
y
-
D
cos
⁡
λ
y
.In order to determine the unknown constants, the appropriate boundary conditions (7) and (8) must be explicitly used. Substituting the total fluid stress tensor components (20)–(22) into boundary conditions equations (8) and velocity components (18)-(19) into boundary conditions equations (7), and after some considerable algebraic manipulation, the scalar potentials ϕ and ψ become(23)
ψ
=
A
sin
⁡
β
x
cos
⁡
δ
y
+
β
2
δ
Φ
sin
⁡
δ
y
,
(24)
ϕ
=
A
cos
⁡
β
x
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,where (25)
Φ
=
2
δ
sin
⁡
δ
L
+
2
β
2
-
α
2
sin
⁡
λ
L
/
λ
2
β
2
cos
⁡
δ
L
-
2
β
2
-
α
2
cos
⁡
λ
L
.Using the expressions of ψ and ϕ from (23) and (24), in (12), we have the following expressions for the velocity components:(26)
v
x
=
A
sin
⁡
β
x
δ
β
2
δ
Φ
cos
⁡
δ
y
-
sin
⁡
δ
y
-
β
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,
(27)
v
y
=
A
cos
⁡
β
x
λ
β
λ
cos
⁡
λ
y
-
β
Φ
sin
⁡
λ
y
-
β
cos
⁡
δ
y
+
β
2
δ
Φ
sin
⁡
δ
y
.Components of the stress tensor can be obtained by substituting the expressions of velocity components v
x; v
y, from (26) and (27), and the expression of scalar potential ϕ from (24) in (13) which comes out as (28)
σ
x
x
=
A
η
cos
⁡
β
x
2
β
δ
β
2
δ
Φ
cos
⁡
δ
y
-
sin
⁡
δ
y
-
2
β
2
+
α
2
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,
(29)
σ
y
y
=
A
η
cos
⁡
β
x
2
β
δ
sin
⁡
δ
y
-
β
2
δ
Φ
cos
⁡
δ
y
-
2
λ
2
+
α
2
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,
(30)
σ
x
y
=
A
η
sin
⁡
β
x
2
β
2
-
α
2
cos
⁡
δ
y
+
β
2
δ
Φ
sin
⁡
δ
y
+
2
β
λ
β
Φ
sin
⁡
λ
y
-
β
λ
cos
⁡
λ
y
.Now we can integrate the normal component of the total force exerted by the fluid over the beam section in order to obtain the external force per unit length acting on the beam in the direction of the flexural displacement f
f
l
u
i
d. Hence, the deflection of the immersed part of the beam (which is equal to u
1
(
y
,
t
)
=
u
1
y
exp
⁡
(
j
ω
t
)) can be determined from (1) and (28) in the form(31)
u
1
′
′
′
′
y
-
ρ
s
S
ω
2
EI
u
1
y
=
2
A
β
η
cos
⁡
β
L
/
2
EI
β
2
Φ
cos
⁡
δ
y
-
δ
sin
⁡
δ
y
-
β
2
+
α
2
2
Φ
cos
⁡
λ
y
+
sin
⁡
λ
y
λ
,where primes denote differentiation with respect to the position variable y.


## 2.1.1. Field Expansions

Consider time harmonic motion throughout with angular frequencyω and with the e
x
p
(
i
ω
t
) dependence suppressed for simplicity. Applying the classical technique of separation of variables in the Cartesian coordinates, the solution of (10) and (11) after some manipulations, and taking into account boundary condition (6) can be shown to be(14)
ψ
=
s
i
n
β
x
A
cos
⁡
δ
y
+
B
sin
⁡
δ
y
,
(15)
ϕ
=
cos
⁡
β
x
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(16)
p
=
η
α
2
cos
⁡
β
x
C
cos
⁡
λ
y
+
D
sin
⁡
λ
ywith the complex-valued coefficients (17)
λ
=
j
β
,
δ
=
j
β
2
-
α
2
,
α
=
1
-
j
2
ω
ν
,where β
=
n
π
/
L and n is an integer. A, B, C, and D are unknown coefficients which will be determined later by imposing the appropriate boundary conditions. Hence, the fluid particle velocity field can be determined from (12) in the form(18)
v
x
=
sin
⁡
β
x
δ
B
cos
⁡
δ
y
-
A
sin
⁡
δ
y
-
β
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(19)
v
y
=
cos
⁡
β
x
λ
D
cos
⁡
λ
y
-
C
sin
⁡
λ
y
-
β
A
cos
⁡
δ
y
+
B
sin
⁡
δ
y
.Consequently, direct substitution of the expansions (14)–(19) into the stress-velocity relations (13), after some manipulations, leads to(20)
σ
x
x
=
η
cos
⁡
β
x
2
β
δ
B
cos
⁡
δ
y
-
A
sin
⁡
δ
y
-
2
β
2
+
α
2
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(21)
σ
y
y
=
η
cos
⁡
β
x
2
β
δ
A
sin
⁡
δ
y
-
B
cos
⁡
δ
y
-
2
λ
2
+
α
2
C
cos
⁡
λ
y
+
D
sin
⁡
λ
y
,
(22)
σ
x
y
=
η
sin
⁡
β
x
2
β
2
-
α
2
A
cos
⁡
δ
y
+
B
sin
⁡
δ
y
+
2
β
λ
C
sin
⁡
λ
y
-
D
cos
⁡
λ
y
.In order to determine the unknown constants, the appropriate boundary conditions (7) and (8) must be explicitly used. Substituting the total fluid stress tensor components (20)–(22) into boundary conditions equations (8) and velocity components (18)-(19) into boundary conditions equations (7), and after some considerable algebraic manipulation, the scalar potentials ϕ and ψ become(23)
ψ
=
A
sin
⁡
β
x
cos
⁡
δ
y
+
β
2
δ
Φ
sin
⁡
δ
y
,
(24)
ϕ
=
A
cos
⁡
β
x
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,where (25)
Φ
=
2
δ
sin
⁡
δ
L
+
2
β
2
-
α
2
sin
⁡
λ
L
/
λ
2
β
2
cos
⁡
δ
L
-
2
β
2
-
α
2
cos
⁡
λ
L
.Using the expressions of ψ and ϕ from (23) and (24), in (12), we have the following expressions for the velocity components:(26)
v
x
=
A
sin
⁡
β
x
δ
β
2
δ
Φ
cos
⁡
δ
y
-
sin
⁡
δ
y
-
β
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,
(27)
v
y
=
A
cos
⁡
β
x
λ
β
λ
cos
⁡
λ
y
-
β
Φ
sin
⁡
λ
y
-
β
cos
⁡
δ
y
+
β
2
δ
Φ
sin
⁡
δ
y
.Components of the stress tensor can be obtained by substituting the expressions of velocity components v
x; v
y, from (26) and (27), and the expression of scalar potential ϕ from (24) in (13) which comes out as (28)
σ
x
x
=
A
η
cos
⁡
β
x
2
β
δ
β
2
δ
Φ
cos
⁡
δ
y
-
sin
⁡
δ
y
-
2
β
2
+
α
2
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,
(29)
σ
y
y
=
A
η
cos
⁡
β
x
2
β
δ
sin
⁡
δ
y
-
β
2
δ
Φ
cos
⁡
δ
y
-
2
λ
2
+
α
2
β
Φ
cos
⁡
λ
y
+
β
λ
sin
⁡
λ
y
,
(30)
σ
x
y
=
A
η
sin
⁡
β
x
2
β
2
-
α
2
cos
⁡
δ
y
+
β
2
δ
Φ
sin
⁡
δ
y
+
2
β
λ
β
Φ
sin
⁡
λ
y
-
β
λ
cos
⁡
λ
y
.Now we can integrate the normal component of the total force exerted by the fluid over the beam section in order to obtain the external force per unit length acting on the beam in the direction of the flexural displacement f
f
l
u
i
d. Hence, the deflection of the immersed part of the beam (which is equal to u
1
(
y
,
t
)
=
u
1
y
exp
⁡
(
j
ω
t
)) can be determined from (1) and (28) in the form(31)
u
1
′
′
′
′
y
-
ρ
s
S
ω
2
EI
u
1
y
=
2
A
β
η
cos
⁡
β
L
/
2
EI
β
2
Φ
cos
⁡
δ
y
-
δ
sin
⁡
δ
y
-
β
2
+
α
2
2
Φ
cos
⁡
λ
y
+
sin
⁡
λ
y
λ
,where primes denote differentiation with respect to the position variable y.


## 2.2. Method of Solution and Frequency Equation

The general solutions of the ordinary differential equations (31) and (2) for the beams system, as shown in Figure 1, can be written in different segments in terms of trigonometric functions and represent propagating waves and hyperbolic functions represent evanescent waves as(32)
u
1
y
=
A
1
cos
⁡
Ω
y
+
B
1
sin
⁡
Ω
y
+
C
1
cosh
⁡
Ω
y
+
D
1
sinh
⁡
Ω
y
+
2
A
β
η
cos
⁡
β
L
/
2
EI
δ
sin
⁡
δ
y
-
β
2
Φ
cos
⁡
δ
y
Ω
4
-
δ
4
+
β
2
+
α
2
/
2
Ω
4
-
λ
4
Φ
cos
⁡
λ
y
+
sin
⁡
λ
y
λ
,
(33)
u
2
y
=
A
2
cos
⁡
Ω
y
+
B
2
sin
⁡
Ω
y
+
C
2
cosh
⁡
Ω
y
+
D
2
sinh
⁡
Ω
y
,where Ω is the flexural wavenumber and is given by(34)
Ω
=
ω
ρ
s
S
EI
1
/
4
.Equation (32) shows that the parameters (35)
ω
c
=
n
π
L
2
EI
ρ
s
S
,
ω
c
=
n
π
/
L
2
ρ
s
S
/
EI
+
1
/
νare resonance frequency parameters and can be plotted for different values of n.The nine constantsA
1, B
1, C
1, D
1, A
2, B
2, C
2, D
2, and A can be found by imposing the following boundary conditions. Equations (32) and (33) are the general solution for the vibration modes of beams partially immersed in fluid. In the case of a coupled system, the effect of fluid viscosity on the flexible beams must be considered. On the fluid-beams interface, the normal velocity must be continuous. Therefore, the fluid velocity and the deflection of a beam u
1 satisfy the relation (36)
v
x
x
=
L
/
2
=
j
ω
u
1
y
,
L
2
⩽
y
⩽
L
.Substituting (26) and (32) into (36), and after some manipulations, leads to the following equation:(37)
A
1
+
2
λ
ω
η
EI
Ω
4
-
δ
4
tan
⁡
β
L
/
2
β
2
Φ
cos
⁡
δ
y
-
δ
sin
⁡
δ
y
-
A
β
2
+
2
λ
ω
η
β
2
+
α
2
/
2
EI
Ω
4
-
λ
4
tan
⁡
β
L
/
2
Φ
cos
⁡
λ
y
+
sin
⁡
λ
y
λ
=
j
ω
sin
⁡
β
L
/
2
A
1
cos
⁡
Ω
y
+
B
1
sin
⁡
Ω
y
+
j
ω
sin
⁡
β
L
/
2
C
1
cosh
⁡
Ω
y
+
D
1
sinh
⁡
Ω
y
.Both sides of (37) are integrated over L
/
2
<
y
<
L to yield the following equation:(38)
A
=
j
ω
Υ
sin
⁡
β
L
/
2
A
1
I
1
+
B
1
I
2
+
C
1
I
3
+
D
1
I
4
,where Υ and the coefficients I
1, I
2, I
3, and I
4 are written down explicitly in Appendix E. Now the expression of the lateral movement of the submerged beams in contact with the fluid u
1 can be formulated using (32) and taking into account (38) as(39)
u
1
y
=
A
1
cos
⁡
Ω
y
+
I
1
U
y
+
B
1
sin
⁡
Ω
y
+
I
2
U
y
+
C
1
cosh
⁡
Ω
y
+
I
3
U
y
+
D
1
sinh
⁡
Ω
y
+
I
4
U
y
,where (40)
U
y
=
2
λ
η
ω
EI
Υ
tan
⁡
β
L
/
2
δ
sin
⁡
δ
y
-
β
2
Φ
cos
⁡
δ
y
Ω
4
-
δ
4
+
β
2
+
α
2
/
2
Ω
4
-
λ
4
Φ
cos
⁡
λ
y
+
sin
⁡
λ
y
λ
.Now, to derive the frequency equation of beams partially immersed in fluid, one assumes that, on the interface between each portion of the beam at the position y
=
L, the deflection, the rotation angle, the internal shear force, and bending moment of the beam must be continuous. This is satisfied when(41)
u
1
L
=
u
2
L
,
u
1
′
L
=
u
2
′
L
,
u
1
′
′
L
=
u
2
′
′
L
,
u
1
′
′
′
L
=
u
2
′
′
′
L
.To complete the formulation of the boundary-value problem, the four boundary conditions for the beam ends (y
=
L
/
2
=
a and y
=
2
L
=
b) considered in this work are specified as follows:(42)
u
1
′
′
a
=
0
,
u
1
′
′
′
a
=
0
,
u
2
b
=
0
,
u
2
′
b
=
0
.Combining the boundary conditions ((41) and (42)) with (33) and (39) yields the following linear homogeneous system of eight equations:(43)
M
X
=
0
,
X
=
A
1
,
B
1
,
C
1
,
D
1
,
A
2
,
B
2
,
C
2
,
D
2
T
,where M is a 8
×
8 matrix whose elements are designated as M
i
j. This system can have nontrivial solutions only if the determinant of the matrix M is zero, leading to the frequency equation. The matrix M is written down explicitly in Appendix F.


## 3. Analytical and Numerical Results and Validation

The model that was developed in Section2 can be validated by comparing the analytical results to model results calculated using Comsol Multiphysics FEM Simulation Software. This comparison used beams with Young’s modulus being E = 200·109 [Pa], Poisson’s ratio being ν
=
0.3, and density being ρ
=
9450
[
k
g
/
m
3
]. The fluid used in the rectangular container for which density of 965 [kg·m−3] and dynamic’s viscosity of 0.1 [Pa·s] is assumed.For coupled vibration analysis of beams partially immersed in a viscous fluid, the accuracy of the present method has been compared with the results obtained with Comsol Multiphysics FEM Simulation Software. The FEM used model had585039 number of degrees of freedom and was analyzed using the algorithm based on the UMFPACK method [19].In this paper the comparison of the values of the angular frequency parameterω is given in Table 1. As one can see from the comparison, very good agreement with those of FEM is obtained.Table 1
The first 24 coupled natural angular frequenciesω of the partially immersed beams in a viscous fluid. Comparison of frequency between FEM and present method.
Order
Present
FEM
Difference(
%
)
1
4.90
4.86
0.82
2
18.18
17.91
1.50
3
20.26
20.01
1.24
4
21.11
20.86
1.19
5
24.68
24.50
0.73
6
30.48
30.20
0.92
7
31.59
31.44
0.47
8
39.24
38.91
0.84
9
44.68
44.51
0.38
10
47.65
47.31
0.71
11
50.20
49.72
0.96
12
54.68
54.02
1.22
13
56.10
55.33
1.39
14
61.47
61.25
0.35
15
61.90
61.36
0.88
16
64.57
64.15
0.65
17
67.29
66.72
0.85
18
73.47
72.99
0.65
19
79.63
79.21
0.53
20
80.06
79.42
0.80
21
85.67
85.41
0.30
22
87.58
87.09
0.56
23
91.33
90.85
0.52
24
97.44
97.10
0.35With the derived eigenfrequency equations, natural frequenciesω are calculated in the software Mathematica. To validate the analytical results, the natural frequencies and mode shapes are also computed using Comsol Multiphysics FEM Simulation Software. The natural frequencies are computed directly from determinant of (43).Table1 shows the comparison of the first 24 natural frequencies and the corresponding mode shapes of viscous fluid by FEM and the present method (43). The good agreement is observed between the results of the present method and those of FEM and the relative difference (
100
×
(
P
r
e
s
e
n
t
-
F
E
M
)
/
F
E
M
) is ⩽2%. This shows that the algorithm implemented in Comsol Multiphysics [20] software for numerical computation is highly reliable and accurate. Mesh refinement can significantly decrease the relative difference.The two possible first coupled mode shapes are represented in Figure2. A very good agreement with Figure  3 of [16] confirms the adequate implementation of the method in the FEM software.Figure 2
Two possible first coupled mode shapes.Figures3 and 4 show the influence of fluid’s viscosity on the odd and even natural frequencies. It is seen that the effect of viscosity is very interesting and decreases the frequency of beams. These figures show also for small fluid’s viscosity that the model is not suitable to describe properly the vibrational behavior of beams immersed in a viscous fluid. In this case, the most appropriate model is the model developed in [21]. In other words, the viscosity terms can be neglected and the resulting model is called inertial coupling [22].Figure 3
The influence of fluid’s viscosity on the first five odd natural frequencies.Figure 4
The influence of fluid’s viscosity on the first five even natural frequencies.Figures5(a)–5(f)–8(a)–8(f) obtained by Comsol show the first 24 coupled modal shapes with the corresponding eigenvalues ω. In these figures, solid lines denote streamlines (i.e., Stokes eddies) caused by system boundaries and fluid structure interaction. In addition, the formation of Stokes eddies is not affected by the presence of an obstacle (beam) in the fluid domain. As expected these Stokes eddies increase when the angular frequency is raised like in the case without obstacle. Finally, these Stokes eddies are symmetric with respect to the obstacle.Figure 5
View of displacement filled of beams and Stokes eddies, showing various mode shapes for the coupled vibration.
(a)
Mode 1
(b)
Mode 2
(c)
Mode 3
(d)
Mode 4
(e)
Mode 5
(f)
Mode 6Figure 6
View of displacement filled of beams and Stokes eddies, showing various mode shapes for the coupled vibration.
(a)
Mode 7
(b)
Mode 8
(c)
Mode 9
(d)
Mode 10
(e)
Mode 11
(f)
Mode 12Figure 7
View of displacement filled of beams and Stokes eddies, showing various mode shapes for the coupled vibration.
(a)
Mode 13
(b)
Mode 14
(c)
Mode 15
(d)
Mode 16
(e)
Mode 17
(f)
Mode 18Figure 8
View of displacement filled of beams and Stokes eddies, showing various mode shapes for the coupled vibration.
(a)
Mode 19
(b)
Mode 20
(c)
Mode 21
(d)
Mode 22
(e)
Mode 23
(f)
Mode 24


## 4. Conclusion

This paper has presented an analytical method to understand the modal behavior of beams partially immersed in a viscous fluid. The validity of the present solution is solidly confirmed numerically. The main application of the proposed method is in characterization of rheological properties of viscous materials [16]. In addition the study can be also extended further for the optimization tool for complex engineering design problems.



---

*Source: 6761372-2016-02-21.xml*