Three-Dimensional Investigation of the Stokes Eigenmodes in Hollow Circular Cylinder

Abstract

This paper studies the influence of boundary conditions on a fluid medium of finite depth. We determine the frequencies and the modal shapes of the fluid. The fluid is assumed to be incompressible and viscous. A potential technique is used to obtain in three-dimensional cylindrical coordinates a general solution for a problem. The method consists in solving analytically partial differential equations obtained from the linearized Navier-Stokes equation. A finite element analysis is also used to check the validity of the present method. The results from the proposed method are in good agreement with numerical solutions. The effect of the fluid thickness on the Stokes eigenmodes is also investigated. It is found that frequencies are strongly influenced.

1. Introduction

Flow modeling in confined domains leads mostly to Stokes models [1]. The use of this model in the field of microfluidics and nanofluidics [2, 3] and in the understanding the movement of microorganisms [4] currently in full swing. This recent years, several authors have focused on the modal analysis of this model. We can cite in particular the aspects related to problems arising from geophysics. These studies highlight the existence of slow wave or Stoneley waves [5–7]. Other authors [8, 9] have focused on the dynamic aspects in the context of fluid-structure interaction. The numeric aspects for coupled (or not) modal problem were discussed in [10–12].

The theory of potential flow of viscous fluid was introduced by [13]. All of his work on this topic is framed in terms of the effects of viscosity on the attenuation of small amplitude waves on a liquid-gas surface. The problem treated by Stokes was solved exactly using the linearized Navier-Stokes equations, without assuming potential flow, and was solved exactly by [14]. Reference [15] has identified the main events in the history of thought about potential flow of viscous fluids. The problem of Stokes flow in cylindrical domain has been investigated by several authors. References [16, 17] studied Stokes flow in a cylindrical container by an eigenfunction expansion procedure without the compressibility effect.

Potential flows through different kinds of geometry have been studied by many investigators for several applications. For example in vascular fluid dynamics, [18] presented the role of curvature in the wave propagation and in the development of a secondary flow. Reference [19] has studied the flow and pressure dynamics of the cerebrospinal fluid flow. Three-dimensional CSF flow studies have also been reported [20].

The knowledge of Stokes eigenmodes in a three-dimensional confined domain, in a square/cube [10], or in any bounded domain could provide some insight into the understanding or analysis of turbulent instantaneous flow field in a geometry as simple as, for instance, the driven cavity.

This paper deals with analytical and numerical modeling of three-dimensional incompressible viscous fluid using the potential technique and the finite element method. The analytical formulation is based upon a convenient decomposition of the velocity field into two contributions [9–21]: one is related to the scalar potential and the other is the vector potential.

The rest of the paper is organized as follows. In Section 2, we write the linearized Navier-Stokes equation, a simple representation for the velocity field in terms of the potential. The implementation of analytical solutions is discussed in detail. Section 3 describes boundary condition treatment for modeling Stokes eigenmodes. Section 4 investigates the case of fluid-solid interaction. Section 5 is completely devoted to the analytical and numerical results. Finally, the paper is ended by some conclusions.

2. Governing Equations

From conservations of mass and momentum and assuming that the fluid is viscous and incompressible, the motion of the fluid flowing in the cylindrical waveguide in the absence of body force is governed by 𝜌𝜕𝐯(𝑟,𝜃,𝑧)()𝜕𝑡=−∇𝑝𝑟,𝜃,𝑧+𝜂∇²𝐯(),()𝑟,𝜃,𝑧∇⋅𝐯𝑟,𝜃,𝑧=0, (1)(2) in which 𝐯(𝑟,𝜃,𝑧)={𝑣_𝑟,𝑣_𝜃,𝑣_𝑧}^𝑇(𝑟,𝜃,𝑧) is the fluid velocity vector, 𝑝(𝑟,𝜃,𝑧) is the fluid pressure, 𝜌 is the density of the fluid, and 𝜈 and 𝜂=𝜌𝜈 are the kinematic and dynamic fluid viscosity, respectively.

The obtained equations of motion are highly complex and coupled. However, a simpler set of equations can be obtained by introducing scalar potentials 𝜙, 𝜓, and 𝜒, known as the Helmholtz decomposition. The Helmholtz decomposition theorem [22] states that any vector 𝐯 can be written as the sum of two parts: one is curl-free and the other is solenoidal. In flow fields, the velocity is thereby decomposed into a potential flow and a viscous flow. In other words, the velocity 𝐯 can be decomposed into the following form: 𝐯=∇𝜙+∇×𝜓𝐞_𝐳+∇×∇×𝜒𝐞_𝐳, (3) where 𝐞_𝐳 is the unit vector along the 𝑧 direction. Substituting the above resolutions into (1) and (2), after some manipulations, the equation for the conservation of mass (2) becomes the Laplace equation: ∇²𝜙(𝑟,𝜃,𝑧)=0, (4) and the equation for the conservation of momentum (1) becomes the Helmholtz equations ∇²𝜓()−1𝑟,𝜃,𝑧𝜈𝜕𝜓(𝑟,𝜃,𝑧)∇𝜕𝑡=0,₂𝜒()−1𝑟,𝜃,𝑧𝜈𝜕𝜒(𝑟,𝜃,𝑧)𝜕𝑡=0, (5) where ∇²=(𝜕²/𝜕𝑟²)+(1/𝑟)(𝜕/𝜕𝑟)+(1/𝑟²)(𝜕²/𝜕𝜃²)+(𝜕²/𝜕𝑧²) is the Laplacian operator in polar coordinates, and time dependence has the form exp(𝑗𝜔𝑡). The pressure can be presented by 𝑝=−𝜌𝜕𝜙.𝜕𝑡 (6) Thus, the Navier-Stokes equations are reduced to formulations (4), (5), and (6). Of course, it does not mean that in any case such a decomposition of fluid’s velocity field gives considerable simplifications in solving the problem because boundary conditions for stresses are not separated. The vibrations are harmonic, with constant angular frequency 𝜔; as a result all potentials components and pressure will depend on time only through the factor exp(𝑗𝜔𝑡).

Applying the method of separation of variables, the solution of the equations for potentials, associated with an axial wave number 𝑘_𝑧, radial wave number 𝑘_𝜓, and circumferential mode parameter 𝑛, after considerable algebraic manipulations, can be shown to be 𝜙=𝐴𝐼_𝑛𝑘_𝑧𝑟+𝐵𝐾_𝑛𝑘_𝑧𝑟()()×𝑘sin𝑛𝜃cos𝑛𝜃cos_𝑧𝑧sin𝑘_𝑧𝑧,𝜓=𝐶𝐽_𝑛𝑘_𝜓𝑟+𝐷𝑌_𝑛𝑘_𝜓𝑟()()×𝑘cos𝑛𝜃sin𝑛𝜃cos_𝑧𝑧sin𝑘_𝑧𝑧,𝜒=𝐸𝐽_𝑛𝑘_𝜓𝑟+𝐹𝑌_𝑛𝑘_𝜓𝑟()()×𝑘sin𝑛𝜃cos𝑛𝜃sin_𝑧𝑧cos𝑘_𝑧𝑧. (7)(8)(9)𝐽_𝑛 and 𝑌_𝑛 are Bessel functions of the first and second kind of order 𝑛. 𝐼_𝑛 and 𝐾_𝑛 are modified Bessel functions of the first and second kinds of order 𝑛. 𝑛 and 𝑘_𝑧 are the azimuthal and axial wavenumbers. 𝑛=0,1,2,…, whereas 𝑘_𝑧 is found by satisfying the symmetry boundary condition on the 𝑧=0 and 𝑧=𝑙. Note that the velocity field 𝐯 has components that are symmetric or antisymmetric in 𝜃 and 𝑧. Following standard practice, the solutions with symmetric (antisymmetric) axial velocities are called the antisymmetric (symmetric) axial modes, respectively, with 𝑘^𝑎_𝑧 and 𝑘^𝑠_𝑧 denoting the corresponding eigenvalues. Thus, we have 𝑘^𝑎_𝑧=(𝑚−1)𝜋/𝑙 and 𝑘^𝑠_𝑧=𝑚𝜋/𝑙 (𝑚=0,1,2,…). However, the azimuthal modes corresponding to cos(𝑛𝜃) and sin(𝑛𝜃) are really the same, due to periodicity in the azimuthal direction; that is, there is no distinction in the values of 𝑛 for the two families. 𝐴, 𝐵, 𝐶, 𝐷, 𝐸, and 𝐹 are unknown coefficients which will be determined later by imposing the appropriate boundary conditions. The radial wave number 𝑘_𝜓 is related to the axial wave number 𝑘_𝑧 by 𝑘²_𝜓=𝛼²−𝑘²_𝑧,𝜔=𝑗𝜈𝛼², (10) where 𝑗=√−1. In the end, the pressure can be obtained by insertion of (7) into (6), 𝑝=𝜂𝛼²𝐴𝐼_𝑛𝑘_𝑧𝑟+𝐵𝐾_𝑛𝑘_𝑧𝑟cos(𝑛𝜃)cos𝑘_𝑧𝑧. (11)

Using (3) and taking into account (7), (8), and (9), the fluid particle velocities in terms of the Bessel functions can be written as 𝑣_𝑟=𝐴𝐼^′_𝑛𝑘_𝑧𝑟+𝐵𝐾^′_𝑛𝑘_𝑧𝑟−𝑛𝑟𝐶𝐽_𝑛𝑘_𝜓𝑟−𝑛𝑟𝐷𝑌_𝑛𝑘_𝜓𝑟+𝑘_𝑧𝐸𝐽^′_𝑛𝑘_𝜓𝑟+𝑘_𝑧𝐹𝑌^′_𝑛𝑘_𝜓𝑟×sin(𝑛𝜃)cos𝑘_𝑧𝑧,𝑣_𝜃=𝑛𝑟𝐴𝐼_𝑛𝑘_𝑧𝑟+𝑛𝑟𝐵𝐾_𝑛𝑘_𝑧𝑟−𝐶𝐽^′_𝑛𝑘_𝜓𝑟−𝐷𝑌^′_𝑛𝑘_𝜓𝑟+𝑛𝑘_𝑧𝑟𝐸𝐽_𝑛𝑘_𝜓𝑟+𝐹𝑌_𝑛𝑘_𝜓𝑟×cos(𝑛𝜃)cos𝑘_𝑧𝑧,𝑣_𝑧=𝑘²_𝜓𝐸𝐽_𝑛𝑘_𝜓𝑟+𝑘²_𝜓𝐹𝑌_𝑛𝑘_𝜓𝑟−𝑘_𝑧𝐴𝐼_𝑛𝑘_𝑧𝑟−𝑘_𝑧𝐵𝐾_𝑛𝑘_𝑧𝑟sin(𝑛𝜃)sin𝑘_𝑧𝑧. (12)

In the following section we consider a fluid with different boundary conditions.

3. Boundary Conditions and Frequency Equation

First we define Γ₁ (𝑟=𝑅₁), as the inner boundary of the fluid region, and Γ₂ (𝑟=𝑅₂), as the outer boundary. We begin the analysis of frequency equation for imposing the no-slip boundary conditions.

3.1. Stokes Eigenmodes in the Case of “No-Slip-No-Slip” Boundary Conditions

The no-slip hypothesis of fluid mechanics states that liquid velocity at a solid surface is equal to the velocity of the solid surface. Hence the no-slip boundary condition can be written as(1)the velocity is equal to zero at Γ₁: 𝑣_𝑟=𝑣_𝜃=𝑣_𝑧=0, (13)(2)the velocity is equal to zero at Γ₂: 𝑣_𝑟=𝑣_𝜃=𝑣_𝑧=0. (14)Combining these boundary conditions with (12) yields for each mode number (𝑛,𝑚) the following linear system: 𝐌₁{𝑥}={0},{𝑥}={𝐴𝐵𝐶𝐷𝐸𝐹}^𝑇, (15) where [𝐌₁] is a 6×6 matrix whose components 𝑎_𝑖𝑗 are given in Appendix A. For a non-trivial solution, the determinant of the matrix [𝐌₁] must be equal to zero det𝐌₁=0. (16)

The roots of (16) give the natural frequencies 𝜔 of the cylindrical oscillations. Figure 16 in Appendix D shows analytical calculation of eigenvalues 𝛼 and using (10) we obtain the natural frequencies 𝜔 [rd/s].

3.2. Stokes Eigenmodes in the Case of "No-Slip-Normal Stress-Free" Boundary Conditions

The pressure and velocity boundary conditions on the free surface are both formulated from the dynamic constraint of continuity of normal momentum flux across the free surface. The component of the stress tensor in the outward normal direction is therefore 𝜎⋅𝐧=−𝑝𝐈+𝜂∇𝐯+(∇𝐯)^𝑇,𝐧=0, (17) in which 𝐈 is a unit tensor. By applying this condition on Γ₁ we obtain 𝜎_𝑟𝑟=𝜎_𝑟𝜃=𝜎_𝑟𝑧=0. (18)

Combining these boundary conditions and (14) yields for each mode number (𝑛,𝑚) the following linear system: 𝐌₂{𝐱}={0},{𝐱}={𝐴𝐵𝐶𝐷𝐸𝐹}^𝑇, (19) where [𝐌₂] is a 6×6 matrix whose components 𝑏_𝑖𝑗 are given in Appendix B. For a nontrivial solution, the determinant of the matrix [𝐌₂] must be equal to zero det𝐌₂=0. (20)

The roots of (20) give the natural frequencies 𝜔 of the cylindrical oscillations.

3.3. Torsional, Flexural, and Breathing Stokes Eigenmode

The results presented in (16) and (20) are a general natural frequencies equation. For some simpler modes, the abovementioned method can be simplified. For example, we have the following.

3.3.1. Torsional Stokes Eigenmode

The torsion mode vibration is such a mode in which the scalar components of the velocity {𝑣_𝑟,𝑣_𝑧} are zeros and only the circumferential velocity 𝑣_𝜃 is independent of 𝜃. This condition is achieved if 𝜙=0 and 𝜒=0. Through (3) this gives for the nonvanishing components of displacement and stresses: 𝑣_𝜃=−𝜕𝜓.𝜕𝑟 (21)

Thus, the general solution for 𝜓 must be constructed from the set 𝜓(𝑟,𝑧)=𝐶𝐽₀𝑘_𝜓𝑟+𝐷𝑌₀𝑘_𝜓𝑟sin𝑘_𝑧𝑧. (22)

In this case the boundary conditions equations (13) and (14) become 𝑣_𝜃=0at𝑟=𝑅₁,𝑣_𝜃=0at𝑟=𝑅₂. (23)

Then, (15) becomes [𝐓]{𝐱}={𝟎},{𝐱}={𝐶𝐷}^𝑇. (24)[𝐓] is a 2×2 matrix whose components are calculated using Appendices A and B. Solving det[𝐓]=0 gives the torsional modes.

3.4. Longitudinal Stokes Eigenmode

Another simpler mode vibration is called longitudinal mode vibration in which 𝑣_𝜃=0 and {𝑣_𝑟,𝑣_𝑧} are independent of 𝜃. This means that the motion is confined to planes perpendicular to the 𝑧-axis, which can move, expand, and contract in their planes. The solution for the displacement field and stress vector follows from (3) and (17): 𝑣_𝑟=𝜕𝜙𝜕𝑟+𝑘_𝑧𝜕𝜒𝑣𝜕𝑟_𝑧=−𝑘_𝑧𝜙+𝑘²_𝜓𝜎𝜒,_𝑟𝑟𝜕=2𝜂²𝜙𝜕𝑟²−𝛼²2𝜙+𝑘_𝑧𝜕²𝜒𝜕𝑟²,𝜎_𝑟𝑧=𝜂−2𝑘_𝑧𝜕𝜙+𝑘𝜕𝑟^2𝜓−𝑘²_𝑧𝜕𝜒.𝜕𝑟 (25)

Thus, the general solution for 𝜙 and 𝜒 must be constructed from the set: 𝜙=𝐴𝐽₀𝑘_𝑧𝑟+𝐵𝑌₀𝑘_𝑧𝑟sin𝑘_𝑧𝑧,𝜒=𝐸𝐽₀𝑘_𝜓𝑟+𝐹𝑌₀𝑘_𝜓𝑟cos𝑘_𝑧𝑧. (26)

In this case the boundary conditions equations (14) and (18) become 𝜎_𝑟𝑟=𝜎_𝑟𝑧=0at𝑟=𝑅₁,𝑣_𝑟=𝑣_𝑧=0at𝑟=𝑅₂. (27)

Then, (19) becomes [𝐋]{𝐱}={0},{𝐱}={𝐴𝐵𝐸𝐹}^𝑇. (28)[𝐋] is a 4×4 matrix whose components are calculated using the Appendices A and B. Solving det[𝐋]=0 gives the longitudinal modes.

3.4.1. Flexural and Breathing Stokes Eigenmodes

The mode shape 𝑛=1 is called flexural mode vibration in which all components of the displacement are nonvanishing and depend on 𝑟, 𝜃, and 𝑧. The mode shape 𝑛⩾2 is called breathing mode vibration in which all components of the displacement are non-vanishing and depend on 𝑟, 𝜃, and 𝑧.

In the following we will introduce the interest of the fluid-structure interaction and we will study the influence of viscosity on the eigenmodes of an elastic solid. For this, the inner boundary of the fluid region Γ₁ is represented by an elastic wall.

4. Fluid-Structure Interaction

Fluid-structure interaction problems have long since attracted the attention of engineers and applied mathematics. The most important applications of this theory, are probably, structural acoustics [23], vibrations of fluid-conveying pipes [24, 25], and biomechanics. As these problems are rather complicated, some simplifications are typically adopted to facilitate their solving. In particular, it is quite typical to ignore viscosity effects (especially in structural acoustics) or to use local theories of interaction, such as the one referred to as thin layer or plane wave approximation.

4.1. Governing Equations of Elastic Media

The wave motion in an isotropic elastic medium is governed by the classical Navier’s equation: −𝜌_𝑠𝜔²𝐮=𝜇∇²𝐮+𝜆+𝜇∇∇⋅𝐮, (29) where 𝜌_𝑠 is the density, 𝜆, 𝜇 are the Lamé constants, and 𝐮(𝑟,𝜃,𝑧)={𝑢_𝑟,𝑢_𝜃,𝑢_𝑧}^𝑇(𝑟,𝜃,𝑧) is the vector displacement of particles.

The obtained equations of motion are highly complex and coupled. However, a simpler set of equations can be obtained by introducing scalar potentials Φ, Ψ, and Θ, known as the Helmholtz decomposition such that 𝐮=∇Φ+∇×Ψ𝐞_𝐳+∇×∇×Θ𝐞_𝐳. (30)

Substituting (30) into (29) leads to three sets of differential equations ∇²Φ()−𝜈𝑟,𝜃,𝑧²𝛼⁴𝑐²_𝐿Φ()∇𝑟,𝜃,𝑧=0,₂Ψ()−𝜈𝑟,𝜃,𝑧²𝛼⁴𝑐²_𝑇Ψ()∇𝑟,𝜃,𝑧=0,₂Θ()−𝜈𝑟,𝜃,𝑧²𝛼⁴𝑐²_𝑇Θ(𝑟,𝜃,𝑧)=0, (31) where 𝑐_𝐿=(𝜆+2𝜇)/𝜌_𝑠 and 𝑐_𝑇=𝜇/𝜌_𝑠 are the compressional and shear wave velocities in the solids, respectively. Applying the method of separation of variables, the solution of the equations for potentials, associated with an axial wave number 𝑘_𝑧, radial wave number (𝑘_Φ,𝑘_Ψ), and circumferential mode parameter 𝑛, after considerable algebraic manipulations, can be shown to be Φ=𝑎𝐼_𝑛𝑘_Φ𝑟+𝑏𝐾_𝑛𝑘_Φ𝑟sin(𝑛𝜃)cos𝑘_𝑧𝑧,Ψ=𝑐𝐼_𝑛𝑘_Ψ𝑟+𝑑𝐾_𝑛𝑘_Ψ𝑟cos(𝑛𝜃)cos𝑘_𝑧𝑧,Θ=𝑒𝐼_𝑛𝑘_Ψ𝑟+𝑓𝐾_𝑛𝑘_Ψ𝑟sin(𝑛𝜃)sin𝑘_𝑧𝑧. (32)

The radial wave number (𝑘_Φ,𝑘_Ψ) is related to the axial wave number 𝑘_𝑧 by 𝑘²_Φ=𝜈²𝛼⁴𝑐²_𝐿+𝑘²_𝑧,𝑘²_Ψ=𝜈²𝛼⁴𝑐²_𝑇+𝑘²_𝑧 (33) and 𝑎, 𝑏, 𝑐, 𝑑, 𝑒, and 𝑓 are unknown coefficients which will be determined later by imposing the appropriate boundary conditions.

Using (30) the scalar components of the displacement vector 𝐮 in cylindrical coordinates can be expressed by 𝑢_𝑟=𝜕Φ−𝑛𝜕𝑟𝑟Ψ+𝑘_𝑧𝜕Θ,𝑢𝜕𝑟_𝜃=𝑛𝑟Φ−𝜕Ψ+𝜕𝑟𝑛𝑘_𝑧𝑟𝑢Θ,_𝑧=−𝑘_𝑧Φ−𝑘²_ΨΘ, (34) and the radial and tangential stresses are given by Hooke’s law as Σ_𝑟𝑟𝜈=𝜆²𝛼⁴𝑐²_𝐿Φ+2𝜇𝜕𝑢_𝑟,Σ𝜕𝑟_𝑟𝜃=𝜇𝜕𝑢_𝜃−𝑢𝜕𝑟𝑟^𝜃+1𝑟𝜕𝑢_𝑟,Σ𝜕𝜃_𝑟𝑧=𝜇𝜕𝑢_𝑟+𝜕𝑧𝜕𝑢_𝑧.𝜕𝑟 (35)

4.2. Boundary Condition and Frequency Equation

We define Γ₁ (𝑟=𝑅₁) as the boundary contact between the fluid region and the solid region and Γ₂ (𝑟=𝑅₂) as the outer boundary. The relevant boundary conditions can be taken as follows. (1)The normal components of the solid stresses must be zero on the interface Γ₀: Σ_𝑟𝑟=Σ_𝑟𝜃=Σ_𝑟𝑧=0. (36)(2)Kinematic boundary condition (velocity must be continuous) on the interface Γ₁: 𝑣_𝑟=𝑗𝜔𝑢_𝑟,𝑣_𝜃=𝑗𝜔𝑢_𝜃,𝑣_𝑧=𝑗𝜔𝑢_𝑧. (37)(3)Dynamic boundary condition (normal stresses must be continuous) on the interface Γ₁: 𝜎_𝑟𝑟=Σ_𝑟𝑟,𝜎_𝑟𝜃=Σ_𝑟𝜃,𝜎_𝑟𝑧=Σ_𝑟𝑧. (38)(4)On the interface Γ₂ at the outer cylinder the normal stress is equal to zero: 𝜎_𝑟𝑟=𝜎_𝑟𝜃=𝜎_𝑟𝑧=0. (39)

Combining these boundary conditions with (34)-(35) and taking into account (7)–(9) and (32) yields for each mode number (𝑛,𝑚) the following linear system: [𝐌]{𝐲}={0}, (40) where {𝐲}={𝐴𝐵𝐶𝐷𝐸𝐹𝑎𝑏𝑐𝑑𝑒𝑓}^𝑇 and [𝐌] is an twelfth-order operator matrix, which is given in Appendix C. For nontrivial solution, the determinant of the matrix 𝐌 must be equal to zero: det[𝐌]=0. (41)

This equation indicates a relationship between the dynamic fluid viscosity 𝜂, density of the fluid 𝜌, and the elastic constants. The roots of (41) give the infinite natural frequencies 𝜔.

5. Analytical, Numerical Results and Validation

Numerical calculations were performed on the example of the hollow cylinder (Figures 1, 2, and 3) with 𝑙=0.15 [m], and different dimensions such as the inner radius 𝑅₁=0.07 [m] and the outer radius 𝑅₂=0.09  [m] of fluid domain are used. The fluid used in the hollow cylinder for which density of 920 [kg·m⁻³] and dynamic’s viscosity of 0.1 [Pa·s] is assumed. The following values of parameters of the the elastic solid are in contact with a viscous liquid were assumed: 𝜌_𝑠=1150 [kg/m³], 𝜈_𝑠=0.48, 𝐸=3⋅10⁵ [Pa], and the inner radius of elastic solid 𝑅₀=0.068 [m].

Figure 1 

Configuration of the viscous oscillations of a cylindrical incompressible fluid of length 𝑙 in circular cylindrical coordinate system (𝑟,𝜃,𝑧). 𝑅1 and 𝑅2 are the inner and outer radius of fluid domain, respectively.

Figure 2 

Geometry of fluid domain.

Figure 3 

Geometry of fluid-solid interation model. 𝑅0 and 𝑅1 are the inner and outer radius of solid domain, respectively.

With the derived eigenfrequency equations, natural frequencies 𝜔_𝑛𝑚=Im(𝜈𝛼²_𝑛𝑚𝑗) for each pair of (𝑛,𝑚) are calculated in the software Mathematica. 𝑚 and 𝑛 denote the mode in axial and azimuthal (propagating clockwise around the vortex) direction, respectively. 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 (24) for torsion vibration. For longitudinal vibration, (28) can be used to determine the corresponding natural frequencies.

Tables 1 and 2 show the comparison of the first 18 natural frequencies and the corresponding mode shapes of viscous fluid by FEM and the present method (see (16) and (20)). For example in the case of “no-slip-normal stress-free” boundary conditions, in the first 18 natural frequencies, four correspond to flexural vibration (𝑛=1), nine to breathing vibration (𝑛⩾2), four to torsional vibration, and one to longitudinal vibration. The very good agreement is observed between the results of the present method and those of FEM and the relative difference ((FEM-Present)/Present) is ⩽1%. This shows that the algorithm implemented in Comsol Multiphysics [26, 27] software for numerical computation is highly reliable and accurate. This algorithm is based on the UMFPACK method [28]. Is the attention to use the numerical formulation in future for more general geometries.

Table 1 

Natural frequencies 𝜔 [rd/s] for various mode shapes in the case of “no-slip-no-slip” boundary conditions.

Table 2 

Natural nfrequencies 𝜔 [rd/s] for various mode shapes in the case of “no-slip-normal stress-free” boundary conditions.

Tables 1 and 2 and Figures 4, 5, and 6 show that natural frequencies are very sensitive to the nature of the boundary conditions. It is seen that the effect of the free surface is very interesting and decrease the frequency of a fluid confined in a rigid cylinder.

Figure 4 

Variations of flexural eigenmode for different boundary conditions.

Figure 5 

Variations of breathing eigenmode for different boundary conditions.

Figure 6 

Variations of torsional eigenmode for different boundary conditions.

Figures 7, 8, and 9 show, respectively, the two modal shapes of the torsional, flexural, and breathing vibrations. The modal shape can be regarded as the mode (𝑛,𝑚), where 𝑛 is the modal number in the circumferential direction and 𝑚 is the modal number in the axial direction. The modal shapes are not in order with the parameters 𝑛 and 𝑚. This feature of cylindrical vibration is different from that of beam vibration in which the order increases with the modal parameter. Therefore in the vibration of the cylinder, one should be careful to find the right mode of the vibration.

(a) (𝑛,𝑚)=(0,1)

(b) (𝑛,𝑚)=(0,2)

Figure 7 

The torsional modal shapes of (𝑛,𝑚): streamlines of components of velocity field.

(a)  (𝑛,𝑚)=(1,1)

(b) (𝑛,𝑚)=(1,2)

Figure 8 

The flexural modal shapes of (1,𝑚): streamlines of components of velocity field.

(a) (𝑛,𝑚)=(2,1)

(b) (𝑛,𝑚)=(2,2)

Figure 9 

The breathing modal shapes of (2,𝑚): streamlines of components of velocity field.

In this paper, the effects of boundary conditions and of cylindrical parameters on the natural frequencies of cylindrical viscous fluid are presented with the present method. In these studies, investigations are carried out to study the effects of circumferential mode 𝑛, axial mode 𝑚, and fluid thickness 𝑒=𝑅₂−𝑅₁ on the frequencies. Influence of viscosity on the natural frequencies of an elastic solid is also investigated.

First, one investigates how the natural frequencies vary with the axial mode 𝑚. Figure 10 shows that the natural frequencies increase as the axial mode 𝑚 increases except when 𝑚=0. This value of 𝑚 corresponds to 2D problem of viscous oscillations.

Figure 10 

Variations of the natural frequencies 𝜔 for different values 𝑚 with the mode 𝑚 with the mode 𝑛.

Secondly, one investigates how the frequencies vary with the fluid thickness. Figures 11, 12, and 13 show that the fluid thickness has a strong influence on the natural frequencies.

Figure 11 

Variations of the torsional eigenmode (𝑛=0) 𝜔 for different fluid’s thickness 𝑒 with the mode 𝑚.

Figure 12 

Variations of the flexural eigenmode (𝑛=1) 𝜔 for different fluid’s thickness 𝑒 with the mode 𝑚.

Figure 13 

Variations of the breathing eigenmode (𝑛=2) 𝜔 for different fluid’s thickness 𝑒 with the mode 𝑚.

Thirdly, one investigates how the dense fluid (added mass) affects the natural frequencies. Tables 3 and 4 show the coupled and uncoupled natural frequencies varying with circumferential and axial mode (𝑛,𝑚). As 𝑛 increases, the difference between the coupled and uncoupled natural frequencies increases. Figures 14 and 15 show that the presence of fluid has no influence on the modal shapes of an elastic solid.

Table 3 

Uncoupled natural frequencies 𝜔 [rd/s] for various mode shapes of elastic solid.

Table 4 

Coupled natural nfrequencies 𝜔 [rd/s] for various mode shapes.

(a) 𝜔=0.952

(b) 𝜔=7.389

(c) 𝜔=11.456

(d) 𝜔=20.289

Figure 14 

The modal shapes (𝑛,𝑚) of elastic solid without fluid: the colours pertain to the displacement filed.

(a) 𝜔=0.702

(b) 𝜔=0.749

(c) 𝜔=0.645

(d) 𝜔=0.813

Figure 15 

The modal shapes (𝑛,𝑚) of elastic solid with fluid: the colours pertain to the displacement filed.

Figure 16 

Analytical calculation of eigenvalues 𝛼.

6. Conclusion

We have presented an analytic solution for the Stokes eigenmodes of a viscous incompressible cylindrical fluid. Using the Helmholtz decomposition for the velocity field, we obtain an eigenvalue problem for 𝜔. The analytical results are in very good agreement with FEM results. This analytic method clearly distinguishes between the potential and rotational components and the contributions of each to various flow variables can be analyzed if one wishes so. The present solution is for cylindrical viscous fluid and is not limited to a cylindrical shape. The same scalar potentials can be used to obtain stokes eigenmodes in the case of a spherical geometry. Finally, a knowledge of the stokes eigenmodes of viscous fluid is likely to be of use in performing a dynamic analysis by modal projection method.

Appendices

A.

The matrix 𝐌₁ in (16) is defined as follows: 𝐌₁=𝑎₁₁𝑎₁₂𝑎₁₃𝑎₁₄𝑎₁₅𝑎₁₆𝑎₂₁𝑎₂₂𝑎₂₃𝑎₂₄𝑎₂₅𝑎₂₆𝑎₃₁𝑎₃₂00𝑎₃₅𝑎₃₆𝑎₄₁𝑎₄₂𝑎₄₃𝑎₄₄𝑎₄₅𝑎₄₆𝑎₅₁𝑎₅₂𝑎₅₃𝑎₅₄𝑎₅₅𝑎₅₆𝑎₆₁𝑎₆₂00𝑎₆₅𝑎₆₆, (A.1) where 𝑎₁₁=𝐼^′_𝑛𝑘_𝑧𝑅₁,𝑎₁₂=𝐾^′_𝑛𝑘_𝑧𝑅₁,𝑎₁₃𝑛=−𝑅₁𝐽_𝑛𝑘_𝜓𝑅₁,𝑎₁₄𝑛=−𝑅₁𝑌_𝑛𝑘_𝜓𝑅₁,𝑎₁₅=𝑘_𝑧𝐽^′_𝑛𝑘_𝜓𝑅₁,𝑎₁₆=𝑘_𝑧𝑌^′_𝑛𝑘_𝜓𝑅₁,𝑎₂₁=𝑛𝑅₁𝐼_𝑛𝑘_𝑧𝑅₁,𝑎₂₂=𝑛𝑅₁𝐾_𝑛𝑘_𝑧𝑅₁,𝑎₂₃=−𝐽^′_𝑛𝑘_𝜓𝑅₁,𝑎₂₄=−𝑌^′_𝑛𝑘_𝜓𝑅₁,𝑎₂₅=𝑛𝑘_𝑧𝑅₁𝐽_𝑛𝑘_𝜓𝑅₁,𝑎₂₆=𝑛𝑘_𝑧𝑅₁𝑌_𝑛𝑘_𝜓𝑅₁,𝑎₃₁=−𝑘_𝑧𝐼_𝑛𝑘_𝑧𝑅₁,𝑎₃₂=−𝑘_𝑧𝐾_𝑛𝑘_𝑧𝑅₁,𝑎₃₅=𝑘²_𝜓𝐽_𝑛𝑘_𝜓𝑅₁,𝑎₃₆=𝑘²_𝜓𝑌_𝑛𝑘_𝜓𝑅₁,𝑎₄₁=𝐼^′_𝑛𝑘_𝑧𝑅₂,𝑎₄₂=𝐾^′_𝑛𝑘_𝑧𝑅₂,𝑎₄₃𝑛=−𝑅₂𝐽_𝑛𝑘_𝜓𝑅₂,𝑎₄₄𝑛=−𝑅₂𝑌_𝑛𝑘_𝜓𝑅₂,𝑎₄₅=𝑘_𝑧𝐽^′_𝑛𝑘_𝜓𝑅₂,𝑎₄₆=𝑘_𝑧𝑌^′_𝑛𝑘_𝜓𝑅₂,𝑎₅₁=𝑛𝑅₂𝐼_𝑛𝑘_𝑧𝑅₂,𝑎₅₂=𝑛𝑅₂𝐾_𝑛𝑘_𝑧𝑅₂,𝑎₅₃=−𝐽^′_𝑛𝑘_𝜓𝑅₂,𝑎₅₄=−𝑌^′_𝑛𝑘_𝜓𝑅₂,𝑎₅₅=𝑛𝑘_𝑧𝑅₂𝐽_𝑛𝑘_𝜓𝑅₂,𝑎₅₆=𝑛𝑘_𝑧𝑅₂𝑌_𝑛𝑘_𝜓𝑅₂,𝑎₆₁=−𝑘_𝑧𝐼_𝑛𝑘_𝑧𝑅₂,𝑎₆₂=−𝑘_𝑧𝐾_𝑛𝑘_𝑧𝑅₂,𝑎₆₅=𝑘²_𝜓𝐽_𝑛𝑘_𝜓𝑅₂,𝑎₆₆=𝑘²_𝜓𝑌_𝑛𝑘_𝜓𝑅₂. (A.2)

B.

The matrix 𝐌₂ in (20) is defined as follows: 𝐌₂=𝑏₁₁𝑏₁₂𝑏₁₃𝑏₁₄𝑏₁₅𝑏₁₆𝑏₂₁𝑏₂₂𝑏₂₃𝑏₂₄𝑏₂₅𝑏₂₆𝑏₃₁𝑏₃₂𝑏₃₃𝑏₃₄𝑏₃₅𝑏₃₆𝑎₄₁𝑎₄₂𝑎₄₃𝑎₄₄𝑎₄₅𝑎₄₆𝑎₅₁𝑎₅₂𝑎₅₃𝑎₅₄𝑎₅₅𝑎₅₆𝑎₆₁𝑎₆₂00𝑎₆₅𝑎₆₆, (B.1) where 𝑏₁₁=𝜂2𝐼_𝑛^′′𝑘_𝑧𝑅₁−𝛼²𝐼_𝑛𝑘_𝑧𝑅₁,𝑏₁₂=𝜂2𝐾_𝑛^′′𝑘_𝑧𝑅₁−𝛼²𝐾_𝑛𝑘_𝑧𝑅₁,𝑏₁₃=2𝜂𝑛𝑅²₁𝐽_𝑛𝑘_𝜓𝑅₁−𝑅₁𝐽^′_𝑛𝑘_𝜓𝑅₁,𝑏₁₄=2𝜂𝑛𝑅²₁𝑌_𝑛𝑘_𝜓𝑅₁−𝑅₁𝑌^′_𝑛𝑘_𝜓𝑅₁,𝑏₁₅=2𝜂𝑘_𝑧𝐽_𝑛^′′𝑘_𝜓𝑅₁,𝑏₁₆=2𝜂𝑘_𝑧𝑌_𝑛^′′𝑘_𝜓𝑅₁,𝑏₂₁=2𝜂𝑛𝑅²₁𝑅₁𝐼^′_𝑛𝑘_𝑧𝑅₁−𝐼_𝑛𝑘_𝑧𝑅₁,𝑏₂₂=2𝜂𝑛𝑅²₁𝑅₁𝐾^′_𝑛𝑘_𝑧𝑅₁−𝐾_𝑛𝑘_𝑧𝑅₁,𝑏₂₃=−𝜂2𝐽_𝑛^′′𝑘_𝜓𝑅₁+𝑘²_𝜓𝐽_𝑛𝑘_𝜓𝑅₁,𝑏₂₄=−𝜂2𝑌_𝑛^′′𝑘_𝜓𝑅₁+𝑘²_𝜓𝑌_𝑛𝑘_𝜓𝑅₁,𝑏₂₅=2𝜂𝑛𝑘_𝑧𝑅²₁𝑅₁𝐼^′_𝑛𝑘_𝑧𝑅₁−𝐼_𝑛𝑘_𝑧𝑅₁,𝑏₂₆=2𝜂𝑛𝑘_𝑧𝑅²₁𝑅₁𝐾^′_𝑛𝑘_𝑧𝑅₁−𝐾_𝑛𝑘_𝑧𝑅₁,𝑏₃₁=−2𝜂𝑘_𝑧𝐼^′_𝑛𝑘_𝑧𝑅₁,𝑏₃₂=−2𝜂𝑘_𝑧𝐾^′_𝑛𝑘_𝑧𝑅₁,𝑏₃₃=𝜂𝑛𝑘_𝑧𝑅₁𝐽_𝑛𝑘_𝜓𝑅₁,𝑏₃₄=𝜂𝑛𝑘_𝑧𝑅₁𝑌_𝑛𝑘_𝜓𝑅₁,𝑏₃₅=𝜂𝑘²_𝜓−𝑘²_𝑧𝐽^′_𝑛𝑘_𝜓𝑅₁,𝑏₃₆=𝜂𝑘²_𝜓−𝑘²_𝑧𝑌^′_𝑛𝑘_𝜓𝑅₁. (B.2)

C.

The operator matrix 𝐌 in (41) is defined as follows: ,𝑏𝐌=𝐀𝐁𝐂𝐃𝐀=000000000000000000₁₁𝑏₁₂𝑏₁₃𝑏₁₄𝑏₁₅𝑏₁₆𝑏₂₁𝑏₂₂𝑏₂₃𝑏₂₄𝑏₂₅𝑏₂₆𝑏₃₁𝑏₃₂𝑏₃₃𝑏₃₄𝑏₃₅𝑏₃₆,𝑐𝐁=¹⁷𝑐₁₈𝑐₁₉𝑐₁₁₀𝑐₁₁₁𝑐₁₁₂𝑐₂₇𝑐₂₈𝑐₂₉𝑐₂₁₀𝑐₂₁₁𝑐₂₁₂𝑐₃₇𝑐₃₈𝑐₃₉𝑐₃₁₀𝑐₃₁₁𝑐₃₁₂𝑐₄₇𝑐₄₈𝑐₄₉𝑐₄₁₀𝑐₄₁₁𝑐₄₁₂𝑐₅₇𝑐₅₈𝑐₅₉𝑐₅₁₀𝑐₅₁₁𝑐₅₁₂𝑐₆₇𝑐₆₈𝑐₆₉𝑐₆₁₀𝑐₆₁₁𝑐₆₁₂,𝑎𝐂=¹¹𝑎₁₂𝑎₁₃𝑎₁₄𝑎₁₅𝑎₁₆𝑎₂₁𝑎₂₂𝑎₂₃𝑎₂₄𝑎₂₅𝑎₂₆𝑎₃₁𝑎₃₂00𝑎₃₅𝑎₃₆𝑐₁₀₁𝑐₁₀₂𝑐₁₀₃𝑐₁₀₄𝑐₁₀₅𝑐₁₀₆𝑐₁₁₁𝑐₁₁₂𝑐₁₁₃𝑐₁₁₄𝑐₁₁₅𝑐₁₁₆𝑐₁₂₁𝑐₁₂₂𝑐₁₂₃𝑐₁₂₄𝑐₁₂₅𝑐₁₂₆,𝑐𝐃=⁷⁷𝑐₇₈𝑐₇₉𝑐₇₁₀𝑐₇₁₁𝑐₇₁₂𝑐₈₇𝑐₈₈𝑐₈₉𝑐₈₁₀𝑐₈₁₁𝑐₈₁₂𝑐₉₇𝑐₉₈00𝑐₉₁₁𝑐₉₁₂,000000000000000000 (C.1) where 𝑐₁₇=2𝜇𝐼_𝑛^′′𝑘_Φ𝑅₀+𝜆𝜈²_𝑓𝛼⁴2𝜇𝑐²_𝐿𝐼_𝑛𝑘_Φ𝑅₀,𝑐₁₈=2𝜇𝐾_𝑛^′′𝑘_Φ𝑅₀+𝜆𝜈²_𝑓𝛼⁴2𝜇𝑐²_𝐿𝐾_𝑛𝑘_Φ𝑅₀,𝑐₁₉=2𝜇𝑛𝑅²₀𝐼_𝑛𝑘_Ψ𝑅₀−𝑅₀𝐼^′_𝑛𝑘_Ψ𝑅₀,𝑐₁₁₀=2𝜇𝑛𝑅²₀𝐾_𝑛𝑘_Ψ𝑅₀−𝑅₀𝐾^′_𝑛𝑘_Ψ𝑅₀,𝑐₁₁₁=2𝜇𝑘_𝑧𝐼_𝑛^′′𝑘_Ψ𝑅₀,𝑐₁₁₂=2𝜇𝑘_𝑧𝐾_𝑛^′′𝑘_Ψ𝑅₀,𝑐₂₇=2𝜇𝑛𝑅²₀𝑅₀𝐼^′_𝑛𝑘_Φ𝑅₀−𝐼_𝑛𝑘_Φ𝑅₀,𝑐₂₈=2𝜇𝑛𝑅²₀𝑅₀𝐾^′_𝑛𝑘_Φ𝑅₀−𝐾_𝑛𝑘_Φ𝑅₀,𝑐₂₉=𝜇𝑘²_Ψ𝐼_𝑛𝑘_Ψ𝑅₀−2𝐼_𝑛^′′𝑘_Ψ𝑅₀,𝑐₂₁₀=𝜇𝑘²_Ψ𝐾_𝑛𝑘_Ψ𝑅₀−2𝐾_𝑛^′′𝑘_Ψ𝑅₀,𝑐₂₁₁=2𝜇𝑛𝑘_𝑧𝑅²₀𝑅₀𝐼^′_𝑛𝑘_Ψ𝑅₀−𝐼_𝑛𝑘_Ψ𝑅₀,𝑐₂₁₂=2𝜇𝑛𝑘_𝑧𝑅²₀𝑅₀𝐾^′_𝑛𝑘_Ψ𝑅₀−𝐾_𝑛𝑘_Ψ𝑅₀,𝑐₃₇=−2𝜇𝑘_𝑧𝐼^′_𝑛𝑘_Φ𝑅₀,𝑐₃₈=−2𝜇𝑘_𝑧𝐾^′_𝑛𝑘_Φ𝑅₀,𝑐₃₉=𝜇𝑛𝑘_𝑧𝑅₀𝐼^′_𝑛𝑘_Ψ𝑅₀,𝑐₃₁₀=𝜇𝑛𝑘_𝑧𝑅₀𝐾^′_𝑛𝑘_Ψ𝑅₀,𝑐₃₁₁=−𝜇𝑘²_Ψ+𝑘²_𝑧𝐼^′_𝑛𝑘_Ψ𝑅₀,𝑐₃₁₂=−𝜇𝑘²_Ψ+𝑘²_𝑧𝐾^′_𝑛𝑘_Ψ𝑅₀,𝑐₄₇=2𝜇𝐼_𝑛^′′𝑘_Φ𝑅₁+𝜆𝜈²_𝑓𝛼⁴2𝜇𝑐²_𝐿𝐼_𝑛𝑘_Φ𝑅₁,𝑐₄₈=2𝜇𝐾_𝑛^′′𝑘_Φ𝑅₁+𝜆𝜈²_𝑓𝛼⁴2𝜇𝑐²_𝐿𝐾_𝑛𝑘_Φ𝑅₁,𝑐₄₉=2𝜇𝑛𝑅²₁𝐼_𝑛𝑘_Ψ𝑅₁−𝑅₁𝐼^′_𝑛𝑘_Ψ𝑅₁,𝑐₄₁₀=2𝜇𝑛𝑅²₁𝐾_𝑛𝑘_Ψ𝑅₁−𝑅₁𝐾^′_𝑛𝑘_Ψ𝑅₁,𝑐₄₁₁=2𝜇𝑘_𝑧𝐼_𝑛^′′𝑘_Ψ𝑅₁,𝑐₄₁₂=2𝜇𝑘_𝑧𝐾_𝑛^′′𝑘_Ψ𝑅₁,𝑐₅₇=2𝜇𝑛𝑅²₁𝑅₁𝐼^′_𝑛𝑘_Φ𝑅₁−𝐼_𝑛𝑘_Φ𝑅₁,𝑐₅₈=2𝜇𝑛𝑅²₁𝑅₁𝐾^′_𝑛𝑘_Φ𝑅₁−𝐾_𝑛𝑘_Φ𝑅₁,𝑐₅₉=𝜇𝑘²_Ψ𝐼_𝑛𝑘_Ψ𝑅₁−2𝐼_𝑛^′′𝑘_Ψ𝑅₁,𝑐₅₁₀=𝜇𝑘²_Ψ𝐾_𝑛𝑘_Ψ𝑅₁−2𝐾_𝑛^′′𝑘_Ψ𝑅₁,𝑐₅₁₁=2𝜇𝑛𝑘_𝑧𝑅²₁𝑅₁𝐼^′_𝑛𝑘_Ψ𝑅₁−𝐼_𝑛𝑘_Ψ𝑅₁,𝑐₅₁₂=2𝜇𝑛𝑘_𝑧𝑅²₁𝑅₁𝐾^′_𝑛𝑘_Ψ𝑅₁−𝐾_𝑛𝑘_Ψ𝑅₁,𝑐₆₇=−2𝜇𝑘_𝑧𝐼^′_𝑛𝑘_Φ𝑅₁,𝑐₆₈=−2𝜇𝑘_𝑧𝐾^′_𝑛𝑘_Φ𝑅₁,𝑐₆₉=𝜇𝑛𝑘_𝑧𝑅₁𝐼^′_𝑛𝑘_Ψ𝑅₁,𝑐₆₁₀=𝜇𝑛𝑘_𝑧𝑅₁𝐾^′_𝑛𝑘_Ψ𝑅₁,𝑐₆₁₁=−𝜇𝑘²_Ψ+𝑘²_𝑧𝐼^′_𝑛𝑘_Ψ𝑅₁,𝑐₆₁₂=−𝜇𝑘²_Ψ+𝑘²_𝑧𝐾^′_𝑛𝑘_Ψ𝑅₁,𝑐₇₇=𝜈𝛼²𝐼^′_𝑛𝑘_Φ𝑅₁,𝑐₇₈=𝜈𝛼²𝐾^′_𝑛𝑘_Φ𝑅₁,𝑐₇₉=−𝜈𝛼²𝑛𝑅₁𝐼^′_𝑛𝑘_Ψ𝑅₁,𝑐₇₁₀=−𝜈𝛼²𝑛𝑅₁𝐾^′_𝑛𝑘_Ψ𝑅₁,𝑐₇₁₁=𝜈𝛼²𝑘_𝑧𝐼^′_𝑛𝑘_Ψ𝑅₁,𝑐₇₁₂=𝜈𝛼²𝑘_𝑧𝐾^′_𝑛𝑘_Ψ𝑅₁,𝑐₈₇=𝜈𝛼²𝑛𝑅₁𝐼_𝑛𝑘_Φ𝑅₁,𝑐₈₈=𝜈𝛼²𝑛𝑅₁𝐾_𝑛𝑘_Φ𝑅₁,𝑐₈₉=−𝜈𝛼²𝐼^′_𝑛𝑘_Ψ𝑅₁,𝑐₈₁₀=−𝜈𝛼²𝐾^′_𝑛𝑘_Ψ𝑅₁,𝑐₈₁₁=𝜈𝛼²𝑛𝑘_𝑧𝑅₁𝐼_𝑛𝑘_Ψ𝑅₁,𝑐₈₁₂=𝜈𝛼²𝑛𝑘_𝑧𝑅₁𝐾_𝑛𝑘_Ψ𝑅₁,𝑐₉₇=−𝜈𝛼²𝑘_𝑧𝐼_𝑛𝑘_Φ𝑅₁,𝑐₉₈=−𝜈𝛼²𝑘_𝑧𝐾_𝑛𝑘_Φ𝑅₁,𝑐₉₁₁=−𝜈𝛼²𝑘²_Ψ𝐼_𝑛𝑘_Ψ𝑅₁,𝑐₉₁₂=−𝜈𝛼²𝑘²_Ψ𝐾_𝑛𝑘_Ψ𝑅₁,𝑐₁₀₁=𝜂2𝐼_𝑛^′′𝑘_𝑧𝑅₂−𝛼²𝐼_𝑛𝑘_𝑧𝑅₂,𝑐₁₀₂=𝜂2𝐾_𝑛^′′𝑘_𝑧𝑅₂−𝛼²𝐾_𝑛𝑘_𝑧𝑅₂,𝑐₁₀₃𝑛=2𝜂𝑅²₂𝐽_𝑛𝑘_𝜓𝑅₂−𝑅₂𝐽^′_𝑛𝑘_𝜓𝑅₂,𝑐₁₀₄𝑛=2𝜂𝑅²₂𝑌_𝑛𝑘_𝜓𝑅₂−𝑅₂𝑌^′_𝑛𝑘_𝜓𝑅₂,𝑐₁₀₅=2𝜂𝑘_𝑧𝐽_𝑛^′′𝑘_𝜓𝑅₂,𝑐₁₀₆=2𝜂𝑘_𝑧𝑌_𝑛^′′𝑘_𝜓𝑅₂,𝑐₁₁₁𝑛=2𝜂𝑅²₂𝑅₂𝐼^′_𝑛𝑘_𝑧𝑅₂−𝐼_𝑛𝑘_𝑧𝑅₂,𝑐₁₁₂𝑛=2𝜂𝑅²₂𝑅₂𝐾^′_𝑛𝑘_𝑧𝑅₂−𝐾_𝑛𝑘_𝑧𝑅₂,𝑐₁₁₃=−𝜂2𝐽_𝑛^′′𝑘_𝜓𝑅₂+𝑘²_𝜓𝐽_𝑛𝑘_𝜓𝑅₂,𝑐₁₁₄=−𝜂2𝑌_𝑛^′′𝑘_𝜓𝑅₂+𝑘²_𝜓𝑌_𝑛𝑘_𝜓𝑅₂,𝑐₁₁₅=2𝜂𝑛𝑘_𝑧𝑅²₂𝑅₂𝐼^′_𝑛𝑘_𝑧𝑅₂−𝐼_𝑛𝑘_𝑧𝑅₂,𝑐₁₁₆=2𝜂𝑛𝑘_𝑧𝑅²₂𝑅₂𝐾^′_𝑛𝑘_𝑧𝑅₂−𝐾_𝑛𝑘_𝑧𝑅₂,𝑐₁₂₁=−2𝜂𝑘_𝑧𝐼^′_𝑛𝑘_𝑧𝑅₂,𝑐₁₂₂=−2𝜂𝑘_𝑧𝐾^′_𝑛𝑘_𝑧𝑅₂,𝑐₁₂₃=𝜂𝑛𝑘_𝑧𝑅₂𝐽_𝑛𝑘_𝜓𝑅₂,𝑐₁₂₄=𝜂𝑛𝑘_𝑧𝑅₂𝑌_𝑛𝑘_𝜓𝑅₂,𝑐₁₂₅=𝜂𝑘²_𝜓−𝑘²_𝑧𝐽^′_𝑛𝑘_𝜓𝑅₂,𝑐₁₂₆=𝜂𝑘²_𝜓−𝑘²_𝑧𝑌^′_𝑛𝑘_𝜓𝑅₂. (C.2)

D.

See Figure 16.