limit $\kappa \to 0$ and a plane-wave solution of the form

$$u = A \cos(k_c x) \qquad (59)$$

the $\Gamma_1$-term in Eq. (43) exerts an average force per wavelength $\lambda = 2\pi/|k_c|$ of

$$\left\langle \Gamma_1 \left[ (\nabla u)^2 + 2u \Delta u \right] \right\rangle_{\lambda} = -\frac{1}{2} \Gamma_1 A^2 k_c^2 \qquad (60)$$

where

$$\langle f(x) \rangle_{\lambda} \equiv \frac{1}{\lambda} \int_{0}^{\lambda} dx f(x) \qquad (61)$$

Comparing Eq. (60) with the average force exerted by a quadratic force $f = au^2$ for the wave solution (59),

$$\langle au^2 \rangle_{\lambda} = \frac{aA^2}{2} \qquad (62)$$

we can approximate the $\Gamma_1$-term by an 'equivalent' average force term of the form

$$\Gamma_1 \left[ (\nabla u)^2 + 2u \Delta u \right] \approx -\Gamma_1 k_c^2 u^2 \qquad (63a)$$

Similarly, the average force due to the $\Gamma_2$-term can be approximated by a cubic force. Since the $\Gamma_2$-term is antisymmetric in $u$, the corresponding mean force is obtained by averaging over the interval $[\lambda/4, 3\lambda/4]$, yielding

$$\Gamma_2 \left[ u (\nabla u)^2 + u^2 \Delta u \right] \approx -\frac{1}{2} \Gamma_2 k_c^2 u^3 \qquad (63b)$$

With these approximations, Eq. (43) reduces to the SH equation

$$\begin{aligned} \partial_t u = & \gamma_0 \Delta u - \gamma_2 \Delta^2 u - \\ & au - (b + \Gamma_1 k_c^2) u^2 - \left(c + \frac{\Gamma_2 k_c^2}{2}\right) u^3 \end{aligned} \qquad (64)$$

To facilitate direct comparison with results in the literature [14], it is convenient to rewrite Eq. (64) in the rescaled normal form

$$\partial_T \phi = -2 \Delta_X \phi - \Delta_X^2 \phi - A \phi - B \phi^2 - \phi^3 \qquad (65)$$

where

$$\begin{aligned} T &= \frac{4\gamma_2 t}{\gamma_0^2} = \frac{t}{3\gamma_0^2} \\ \Delta_X &= \frac{\Delta}{k_c^2} = \frac{\Delta}{6|\gamma_0|} \\ \phi &= \frac{u}{u_*} \end{aligned} \qquad (66)$$

$$\begin{aligned} u_* &= \sqrt{\frac{\gamma_0^2}{4\gamma_2(c + \Gamma_2 k_c^2/2)}} = \sqrt{\frac{\gamma_0^2}{(c/3) + \Gamma_2|\gamma_0|}} \\ A &= \frac{a}{a_c} = \frac{4a\gamma_2}{\gamma_0^2} = \frac{3a}{\gamma_0^2} \\ B &= \frac{2(b + \Gamma_1 k_c^2)}{\sqrt{(\gamma_0^2/\gamma_2)(c + \Gamma_2 k_c^2/2)}} = u_* \frac{(b/3) + 2|\gamma_0|\Gamma_1}{\gamma_0^2} \end{aligned} \qquad (67)$$

with $a_c$ given by Eq. (54).

**Excess film stress.** Our model accounts for the excess film stress $\Sigma_e = (\sigma/\sigma_c) - 1$ through the contribution $\tilde{a}_2\Sigma_e$ that appears in the coefficient $a$ of the linear force, see Eq. (44). This specific functional relationship between $a$ and $\Sigma_e$ follows from the requirement that, in the planar limit $\kappa \to 0$, our model must produce the same amplitude-stress relation as classical wrinkling theory [5]. More precisely, classical planar wrinkling theory predicts that the amplitude $A_{el}$ of the wrinkling solution $u = A_{el} \cos(kx)$, measured in units of the film thickness $h$, is equal to square of the excess film stress [9],

$$A_{el} = \sqrt{\Sigma_e}. \qquad (67)$$

This must be equal to the amplitude $A$ of the solution $u = A \cos(k_c x)$ of the standard SH equation (64), which is given by [14]

$$A = \frac{2}{\sqrt{3}} u_* \sqrt{1 - A} \qquad (68)$$

Equating the amplitudes yields

$$A = \frac{a}{a_c} = 1 - \frac{3\Sigma_e}{4u_*^2} \qquad (69)$$

Finally, considering the definition of $a$ in Eq. (44) in the planar case $\kappa = 0$ and inserting $a_c$ from Eq. (54), we find

$$\tilde{a}_2 = -\frac{3\gamma_2}{4u_*^2} \left( \frac{3E_s}{E_f} \right)^{4/3} \qquad (70)$$

We have thus determined all parameters of the effective theory with the exception of $\tilde{c}$, which enters through both $c$ in Eq. (44) and $u_*$ in (70), see Eq. (66). As shown in the next section, $\tilde{c}$ can be estimated by matching our model predictions with the experiments. Table I summarizes the results for our model parameters, retaining curvature terms up to $O(\kappa^3)$, and is also reproduced in Table 1 of the Main Text.

**Remark.** In the above derivation, we have matched the amplitude-stress relationship using the plane-wave wrinkling solution. Alternatively, one can match the amplitude-stress relationship of wrinkling patterns using the hexagonal solution of the SH equation. The general procedure is identical to the one above. The hexagonal solution $u_H$ of the standard SH equation is given by [14]

$$u_H = A \left[ \cos(k_c x) + 2 \cos \left( \frac{k_c x}{2} \right) \cos \left( \frac{\sqrt{3} k_c y}{2} \right) \right] \qquad (71)$$

with

$$A = \frac{2}{15} \left[ B u_* + \sqrt{(B u_*)^2 + 15 u_*^2 (1 - A)} \right] \qquad (72)$$

$B$ vanishes in the planar case, where $\kappa = 0$. Classical planar wrinkling theory predicts the same solution, Eq. (71), with amplitude [9]

$$A_{el} = \frac{2}{\sqrt{11 + 6\nu - 5\nu^2}} \sqrt{\Sigma_e} \qquad (73)$$
Comparing the two amplitudes, we obtain

$$A = 1 - \frac{15\Sigma_e}{(11 + 6\nu - 5\nu^2)u_*^2} \quad (74)$$

from which a and $\tilde{a}_2$ follow again by Eqs. (44) and (54).

Hence, matching amplitudes in the hexagonal phase again leads to a linear relationship between A and $\Sigma_e$, but the prefactor differs compared with the plane-wave solution considered above.

However, irrespective of the chosen matching procedure, we obtain a linear relation between $\tilde{a}_2$ and $\tilde{c}$, leaving only one fitting parameter, $\tilde{c}$. Due to the linear dependence between $\tilde{a}_2$ and $\tilde{c}$ for both hexagonal and labyrinth patterns, one can use either matching relation to fit the experiments. Since the plane-wave (labyrinth) solution (67) leads to slightly simpler formulas, we choose the plane-wave matching in the following; that is, we assume that A and $\Sigma_e$ are related by Eq. (69), and $\tilde{a}_2$ is given by Eq. (70). This convention means that we must determine $\tilde{c}$ by matching the theoretically predicted phase transition curve between the bistable and labyrinth states to the corresponding experimentally measured phase transition curve. The fitting value $\tilde{c}$ then automatically fixes the theoretical prediction for the second transition curve between hexagonal and bistable phase.

### Phase diagram of the wrinkling morphology

The reduction of our effective theory to a standard Swift-Hohenberg equation allows us to make approximative predictions regarding the wrinkling morphologies. For the normal form, Eq. (65), three different wrinkling phases emerge depending on the choice of parameters A and B [14]:

$$\text{Unwrinkled:} \qquad 1 - A \le 0 \quad (75a)$$

$$\text{Hexagons:} \qquad -\frac{1}{15}B^2 < 1 - A < \frac{4}{3}B^2 \quad (75b)$$

$$\text{Bistable phase:} \qquad \frac{4}{3}B^2 < 1 - A < \frac{16}{3}B^2 \quad (75c)$$

$$\text{Labyrinths:} \qquad \frac{16}{3}B^2 < 1 - A \quad (75d)$$

We recall that the coefficient B depends on $\kappa$ via $\Gamma_1$, $\Gamma_2$, b and c, see Eq. (66), while A depends on the excess film stress and $\kappa$ via Eq. (69). Substituting Eqs. (69) and (66) for A and B in Eq. (75), we obtain to leading order in $\kappa$ the following stability criteria:

$$\text{Unwrinkled:} \qquad \Sigma_e \le 0 \quad (76a)$$

$$\text{Hexagons:} \qquad -\frac{\rho}{20}\kappa^2 < \Sigma_e < \rho\kappa^2 \quad (76b)$$

$$\text{Bistable phase:} \qquad \rho\kappa^2 < \Sigma_e < 4\rho\kappa^2 \quad (76c)$$

$$\text{Labyrinths:} \qquad 4\rho\kappa^2 < \Sigma_e \quad (76d)$$

where

$$\rho = \frac{1}{[\eta^{1/3}(1 - \nu)\tilde{c}]^2} \quad (76e)$$

The free parameter $\tilde{c}$ is a function of $\eta = 3E_s/E_f$ and $\nu$. The functional form of $\tilde{c} = \tilde{c}(\eta, \nu)$ determines how the phase transition lines depend on those material properties. Comparison with our experimental data suggests, however, that the phase transition lines are in fact independent of $\eta$ and $\nu$ (Fig. 3 of Main Text). In this case, $\rho$ must be independent of $\eta$ and $\nu$, which means that

$$\tilde{c} = \frac{c_1}{(1 - \nu)\eta^{1/3}}, \quad (77)$$

with fit parameter $c_1$. Adopting the ansatz (77), we find that $c_1 = 0.0188 \pm 0.0002$ gives the best fit to the experimental data (Fig. 3 of Main Text). The resulting good agreement with the data suggests strongly that the critical curves, which separate the different wrinkling phases in our experiments, are independent of material properties.

### Hysteresis

We explain how the critical curves in Fig. 4 of Main Text are obtained. To this end, recall that the amplitude of the hexagonal solutions is given by Eq. (72). Inserting Eq. (69) for A gives

$$A = \frac{2}{15} \left[ B u_* + \sqrt{(B u_*)^2 + \frac{45 \Sigma_e}{4}} \right] \quad (78)$$

To compare with the labyrinth solutions, we consider the measurable difference $U_H$ of the maximal and minimal values of the hexagonal displacement field $u_H$ (measured in units of h). Equation (71) implies that

$$U_H = \max u_H - \min u_H = \frac{9}{2}A \quad (79)$$

The prefactor in the last expression is due to the fact that $u_H$ assumes its minimum at $-3A$ and its maximum at $+3A/2$.

To calculate the corresponding quantity $U_L$ for the labyrinth solution $u_L$, we assume that labyrinths are locally described by the plane-wave $u_L = A\cos(k_c x)$. We then obtain with Eqs. (68) and (69)

$$U_L = \max u_L - \min u_L = 2\sqrt{\Sigma_e} \quad (80)$$

Note that $U_L$ is independent of curvature, whereas the square-root law for $U_H$ is shifted horizontally and vertically by the $Bu_*$-terms (see Fig. 4b of Main Text). To first order in curvature, we have

$$Bu_* \approx \frac{3\kappa}{4c_1} \quad (81)$$
<figure><img src="image_51.png" /><figcaption>Fig. S1: Comparison of wrinkling morphologies in the hexagonal phase: (a) Planar analytical solution from Eq. (71), (b) numerical solution of the generalized Swift-Hohenberg theory from Eq. (43) (simulation parameters: $\gamma_0 = -0.08$, $a = 0.0151$, $c = 0.0095$, $R/h = 80$), and (c) 3D surface scan from experiments ($R = 20$mm, $h = 0.375$mm, $E_f = 2100$kPa, $E_s = 230$kPa). $\theta_1$ and $\theta_2$ denote polar angles.</figcaption></figure>

yielding for the horizontal and vertical shifts (see Fig. 4b in the Main Text)

$$ \delta\Sigma_e = \frac{\kappa^2}{20c_1^2}, \quad \delta U_H = \frac{9\kappa}{10c_1} \qquad (82) $$

Based on the results for $U_L$ and $U_H$, we expect two hysteresis cycles (see Fig. 4 of Main Text). The first cycle is due to the shift $\delta\Sigma_e$, which creates a subcritical stability zone for the hexagonal solutions [14]. Increasing $\Sigma_e$ from the unwrinkled phase, hexons emerge at the onset $\Sigma_e = 0$. However, once they have formed, hexons remain stable even if the excess stress $\Sigma_e$ is subsequently reduced below the wrinkling onset $\Sigma_e = 0$. Only for $\Sigma_e < -\delta\Sigma_e$, hexons lose stability and the unwrinkled solution remains as the only stable state. The second hysteresis cycle involves larger excess stresses (see Fig. 4a of Main Text). Starting from the hexagonal phase, hexons remain stable if the excess stress is increased into the bistable phase [14]. Only if $\Sigma_e$ is increased beyond the bistable-to-labyrinth transition line, hexons lose their stability and the system jumps to a labyrinth state. Upon decreasing $\Sigma_e$ from the labyrinth phase, the system adopts a different path as labyrinthins remain stable throughout the bistable phase; a transition to hexons happens when $\Sigma_e$ is decreased below the critical value that separates the hexagonal from the bistable phase.

## Validation of hexagonal patterns

To test our effective theory further, we compare the planar hexagon solution, Eq. (71), with the numerical solution of Eq. (43) on a spherical geometry (Fig. S1a,b). The good agreement between analytical and numerical solution confirms that the effective theory is indeed well approximated by the standard planar SH equation (65). These results corroborate that once hexons are selected, curvature has negligible influence on their morphology.

Moreover, analytical and numerical solutions compare well with the experimentally determined 3D surface scans (Fig. S1c), demonstrating that the presented effective theory is able to reproduce the morphological details of thin film buckling patterns.

## TOROIDAL GEOMETRIES

We consider a torus with major radius $R_1$ and minor radius $R_2$ measured in units of $h$. Using Eq. (1) and the standard surface parametrization $\mathbf{S} = [(R_1 + R_2 \cos \theta_2) \cos \theta_1, (R_1 + R_2 \cos \theta_2) \sin \theta_1, R_2 \sin \theta_2]$ with coordinates $(\theta_1, \theta_2) \in [0, 2\pi) \times [0, 2\pi)$ we obtain the metric tensor

$$ (a_{\alpha\beta}) = \begin{pmatrix} (R_1 + R_2 \cos \theta_2)^2 & 0 \\ 0 & R_2^2 \end{pmatrix} \qquad (83) $$

Equation (5a) yields for the curvature tensor

$$ (b_{\alpha\beta}) = \begin{pmatrix} -\cos \theta_2(R_1 + R_2 \cos \theta_2) & 0 \\ 0 & -R_2 \end{pmatrix} \qquad (84) $$

and, from Eq. (26), one finds mean and Gaussian curvature as

$$ \mathcal{H} = -\frac{1}{2} \left( \frac{\cos \theta_2}{R_1 + R_2 \cos \theta_2} + \frac{1}{R_2} \right) \qquad (85a) $$

$$ \mathcal{K} = \frac{\cos \theta_2}{R_1 R_2 + R_2^2 \cos \theta_2} \qquad (85b) $$

Recall that the coefficients $\tilde{a}$, $\tilde{\gamma}$, and $\tilde{a}_2$ were determined in the previous section using asymptotic comparison with the flat case. Hence, the expressions for these coefficients remain valid for toroidal geometries. Moreover, we observe that only the symmetry-breaking term in Eq. (34) contains contractions of the curvature tensor of first order, whereas the coefficients (33c) only depend on curvature at second or higher order. For instance, for
<figure><img src="image_52.png" /><figcaption>Fig. S2: Wrinkling morphologies on a toroidal geometry with $R_1$ = 80 and $R_2$ = 16 for increasing excess film stress: (a) $\Sigma_e$ = 0.25, (b) $\Sigma_e$ = 0.5 and (c) $\Sigma_e$ = 2.0. As in spherical geometries, one observes a transition from hexagonal to labyrinth-like patterns with increasing excess stress $\Sigma_e$. However, in contrast to the spherical case, the non-constant curvature on the torus can lead to local symmetry-breaking, i.e., at intermediate values of $\Sigma_e$ labyrinth patterns are more likely to emerge at the inner saddle-like regions of the torus, see (b), whereas hexagons remain stable in the outer regions of the torus, where the two principal curvatures have the same sign resulting in a larger mean curvature. Simulation parameters are $\gamma_0$ = -0.079, $c$ = 0.075, (a) $a$ = 0.017, (b) $a$ = 0.016, and (c) $a$ = 0.007, using an unstructured surface triangulation with > 26,000 nodes.</figcaption></figure>

a toroidal geometry with $R_1$ = 80, $R_2$ = 16, and $\eta$ = 0.33 as shown in Fig. S2, we have $\gamma_0$ = -0.079, so that $|\gamma_0|$ is much larger than the maximum curvature correction $\max_{\theta_1, \theta_2} |\mathcal{R}/6| \approx 0.0007$. We can therefore neglect curvature corrections in the coefficients (33c).

Simulations of Eq. (34) confirm that a symmetry-breaking transition from hexagonal to labyrinth-like structures can also be observed on toroidal geometries when the excess stress is increased (Fig. S2). However, as curvature is now spatially varying, we find that the transition occurs first in the inner regions of the torus, where the principle curvatures have opposite sign (Fig. S2b). In the outside regions, where both principal components have the same sign, hexagons remain stable for relatively larger overstresses, until they become also unstable eventually (Fig. S2c).

In summary, Fig. S2 illustrates that the generalized theory derived above can be applied to arbitrarily curved surfaces.

[1] P. G. Ciarlet. *Mathematical elasticity, vol. III: Theory of shells*. North Holland, 2000.

[2] D. Breid and A. J. Crosby. Curvature-controlled wrinkle morphologies. *Soft Matter*, 9:3624–3630, 2013.

[3] M. M. Müller, M. Deserno, and J. Guven. Interface-mediated interactions between particles: A geometrical approach. *Phys. Rev. E*, 72(6):061407, 2005.

[4] W. Helfrich. Elastic properties of lipid bilayers: theory and possible experiments. *Z. Naturforsch.*, 28(11):693–703, 1973.

[5] Z. Y. Huang, W. Hong, and Z. Suo. Nonlinear analyses of wrinkles in a film bonded to a compliant substrate. *J. Mech. Phys. Solids*, 53:2101–2118, 2005.

[6] J. Swift and P. C. Hohenberg. Hydrodynamic fluctuations at the convective instability. *Phys. Rev. A*, 15(1):319–328, 1977.

[7] M. C. Cross and P. C. Hohenberg. Pattern formation outside of equilibrium. *Rev. Mod. Phys.*, 65(3):851–1112, 1993.

[8] H. G. Allen. *Analysis and Design of Structural Sandwich Panels*. Pergamon, New York, 1969.

[9] S. Cai, D. Breid, A. J. Crosby, Z. Suo, and J. W. Hutchinson. Periodic patterns and energy states of buckled films on compliant substrates. *J. Mech. Phys. Solids*, 59:1094–1114, 2011.

[10] L. D. Landau and E. M. Lifshitz. *Elasticity theory*. Pergamon Press, New York, 1975.

[11] Y. Jeong, Y.-C. Chen, M. K. Turksoy, S. Rana, G. Y. Tonga, B. Creran, A. Sanyal, A. J. Crosby, and V. M. Rotello. Tunable elastic modulus of nanoparticle monolayer films by host–guest chemistry. *Advanced Materials*, in press, doi: 10.1002/adma.201401226, 2014.

[12] G. Cao, X. Chen, C. Li, A. Ji, and Z. Cao. Self-assembled triangular and labyrinth buckling patterns of thin films on spherical substrates. *Phys. Rev. Lett.*, 100(3):036102, 2008.

[13] J. Yin, Z. Cao, C. Li, I. Sheinman, and X. Chen. Stress-driven buckling patterns in spheroidal core/shell structures. *Proc. Natl Acad. Sci.*, 105(49):19132–19135, 2008.

[14] A. A. Golovin and A. A. Nepomnyashchy. *Self-assembly, pattern formation and growth phenomena in nano-systems*. Springer, Dordrecht, 2006.
C. INTERACTION FLUIDE-STRUCTURE
ENTRE DEUX CYLINDRES
CONCENTRIQUES

Dans cette annexe, la publication parue dans la revue *Journal of Applied Mechanics*, sous le titre "Viscous theory for the vibrations of coaxial Cylinders : analytical formulas for the fluid forces and the modal added coefficients" est reproduite<fnref n="1" />.

<footnote n="1">¹ Nomination par Journal of Applied Mechanics pour le prix du meilleur article 2023.</footnote>
# Viscous Theory for the Vibrations of Coaxial Cylinders: Analytical Formulas for the Fluid Forces and the Modal Added Coefficients

This article addresses the small-amplitude forced beam vibrations of two coaxial finite-length cylinders separated by a viscous Newtonian fluid. A new theoretical approach based on an Helmholtz expansion of the fluid velocity vector is carried out, leading to a full analytical expression of the fluid forces and subsequently of the modal added mass and damping coefficients. Our theory shows that the fluid forces are linear combinations of the Fourier harmonics of the vibration modes. The coefficients of the linear combinations are shown to depend on the aspect ratio of the cylinders, on the separation distance, and on the Stokes number. As a consequence, the linear fluid forces do not have, in general, the same shape as the forced vibration mode, so that the fluid makes it possible to couple vibration modes with different wave numbers. Compared to the previous works, the present theory includes the viscous effects of the fluid, accounts for the finite length of the cylinders, does not rely on the assumption of a narrow annulus, and covers in a unique formulation all types of classical boundary conditions for an Euler-Bernoulli beam. The theoretical predictions for the modal added mass and damping coefficients (self and cross) are corroborated numerically, considering rigid, pinned-pinned, and clamped-free vibrations.
[DOI: 10.1115/1.4056910]

Keywords: vibration, fluid-structure interaction, fluid forces, coaxial cylinders, viscosity

**Romain Lagrange<fnref n="1" />**

Université Paris-Saclay, CEA,
Service d'Etudes Mécaniques et Thermiques,
F-91191 Gif-sur-Yvette, France
e-mail: romain.lagrange@cea.fr

**Maria Adela Puscas**

Université Paris-Saclay, CEA,
Service de Thermo-hydraulique et de
Mécanique des Fluides,
F-91191 Gif-sur-Yvette, France
e-mail: maria-adela.puscas@cea.fr

## 1 Introduction

This article addresses the hydrodynamic interaction between two elongated cylinders separated by a viscous quiescent fluid. This classical fluid-structure problem has been the topic of considerable works; see the book of Refs. [1,2] and the survey of Ref. [3] for an extensive bibliography; among which significant analytical contributions were provided by [4–14]. Usually, none of the theoretical approaches is a panacea for a large spectrum of boundary conditions and geometric complexities. Therefore, the analytical calculations mainly differ from each other by the structural models introduced to account for the free or forced vibrations of the inner and outer structures (rigid, beams, shells) and the properties of the coupling fluid (inviscid, viscous). For example, in Ref. [4], the case of flexible concentric cylinders imposed harmonic breathing shell modes (modes with no longitudinal wave number) was considered. A classical acoustic fluid equation for the pressure field was solved, leading to a complex solution expressed in terms of Bessel functions. In Ref. [7], a similar study was carried out on the small-amplitude vibrations (breathing mode) of flexible coaxial shells of infinite length, leading to exact expressions for the added mass coefficients (self and cross). In the same spirit, Ref. [8] considered the small oscillations of a flexible shell (regarded as a cantilever beam with a flexible support) surrounded by a rigid shell (representative of a reactor vessel). The slender-body formulation of Ref. [15] (two-dimensional hydrodynamic mass formulas for rigid fluid-coupled bodies) was used to perform a parametric study on the effects of the boundary stiffness and the confinement (narrow gap) on the modal added mass coefficients and the natural frequencies of the coupled system. This study

was extended by Ref. [16], who replaced the slender-body assumption with a wave equation for compressible and incompressible fluid flow.

One of the reasons for the high interest in the fluid-coupled cylinders problem is that the added mass of the fluid, especially for a narrow annulus, can be several times larger than the mass of the structures that constitute the boundaries of the annulus. Hence, the modal characteristics of the structures are not just modified by the presence of the fluid but are strongly affected by it. Viscous effects, which mainly occur near the interface of the structure, in a thin layer of the rotational flow, are also of great importance in the dynamic response of coupled systems. Actually, in industrial applications, like the one that motivated this study and presented later, the added mass and damping effects may be so strong that they are sometimes the only effects taken into account by engineers in the calculation of the vibrational characteristics of the system, even when the fluid is flowing.

## 2 Industrial Context and Motivation

The present study finds its motivations in the context of the vibrations of internal components of the Jules Horowitz Reactor. This test reactor is mainly dedicated to the research on materials and fuels for the nuclear power industry and nuclear medicine, see Ref. [17] for more details. Essentially, the internal structure of the Jules Horowitz reactor is a system of finite, coaxial cylinders immersed in a fluid. A simplified schematic diagram of the reactor is shown in Fig. 1. The reactor core is made of cylindrical fuel elements of length $L=0.7$ m. Aluminum cylindrical guide tubes, of length 2.1 m and outer diameter $2R_1 = 40$ mm, are inserted in the free space between fuel elements to guide experimental devices to be tested. The diameter of this free space is $2R_2 = 44$ mm so that the middle part of the guide tube is relatively confined (2 mm clearance). Upstream and downstream, the guide tube is less confined, and the clearance varies between 4 and 6 mm. The fuel

<footnote n="1">¹Corresponding author.</footnote>

<footnote>Contributed by the Applied Mechanics Division of ASME for publication in the JOURNAL OF APPLIED MECHANICS. Manuscript received December 15, 2022; final manuscript received February 8, 2023; published online March 6, 2023. Tech. Editor: Pradeep Sharma.</footnote>
<figure><img src="image_53.png" /><figcaption>Fig. 1 Schematic diagram of the Jules Horowitz reactor, from Ref. [18]: (a) 3D sketch of the Jules Horowitz reactor and (b) axial cross section of the Jules Horowitz reactor</figcaption></figure>

elements and the guide tubes are flexible and are subjected to random vibration induced by a turbulent axial flow from the primary circuit. The slender-body approximation is irrelevant here since, for instance, the middle part of the guide tube is not sufficiently slender, i.e., $L/R_1 = 35$.

This article reconsiders the analytical calculation of the hydrodynamics forces acting on coaxial elongated cylinders separated by a viscous stagnant Newtonian fluid. The cylinders are imposed small-amplitude harmonic beam vibrations.

The dynamics of the two cylinders (length L) is usually modeled by a system of forced oscillator equations. For a given couple of structural vibration modes ($W_1, W_2$) of cylinders ($C_1, C_2$), equations for the modal amplitudes ($A_1, A_2$) are given as follows:

$$M_{s_1} \frac{d^2 A_1}{dT^2} + C_{s_1} \frac{dA_1}{dT} + K_{s_1} A_1 = \langle W_1, \overline{F_{T_1}} + \overline{F_{T_1}} + \overline{F_{T_2}} \rangle_L \quad (1a)$$

$$M_{s_2} \frac{d^2 A_2}{dT^2} + C_{s_2} \frac{dA_2}{dT} + K_{s_2} A_2 = \langle W_2, \overline{F_{T_2}} + \overline{F_{T_2}} + \overline{F_{T_3}} \rangle_L \quad (1b)$$

where $T, M_{s_i}, C_{s_i}$, and $K_{s_i}$ are the dimensional time, modal structural mass, damping, and rigidity coefficients, respectively. The right-hand side forcing terms in Eq. (1) represent modal forces, i.e., linear fluid forces $\overline{F_{T_i}} + \overline{F_{ij}}$ projected onto vibration modes $W_i$, according to the scalar product $\langle W_i, \overline{F_{ij}} \rangle_L = \int_0^L W_i \overline{F_{ij}} dL$. The term $\overline{F_{T_i}}$ refers to a turbulent (random) linear fluid force, considered as independent of cylinder vibrations. The term $\overline{F_{ij}}$ refers to the linear fluid force acting on $C_i$ due to the vibration of $C_j$. This latter is assumed to

be a linear combination of the modal acceleration $d^2 A_j/dT^2$, velocity $dA_j/dT$, and position $A_j$, i.e.,

$$\langle W_i, \overline{F_{ij}} \rangle_L = -M_{ij} \frac{d^2 A_j}{dT^2} - C_{ij} \frac{dA_j}{dT} - K_{ij} A_j \quad (2)$$

where $M_{ij}, C_{ij}$, and $K_{ij}$ the modal added mass, damping, and rigidity coefficients, respectively (self coefficients if $i=j$, cross coefficients otherwise). As most often assumed, it is postulated that the modal added mass $M_{ij}$ does not depend on flow velocity, whereas $C_{ij}$ and $K_{ij}$ are expanded as follows:

$$C_{ij} = C_{ij}^{(0)} + \hat{C}_{ij}, \quad K_{ij} = K_{ij}^{(0)} + \hat{K}_{ij} \quad (3)$$

The term $C_{ij}^{(0)}$ (resp. $K_{ij}^{(0)}$) is the modal added damping (resp. stiffness) in a still fluid, and $\hat{C}_{ij}$ (resp. $\hat{K}_{ij}$) is a complementary damping (resp. stiffness) due to fluid flow velocity. Usually, surface stiffness effects in still fluid are neglected, i.e., $K_{ij}^{(0)} = 0$, so that Eq. (1) is rewritten as follows:

$$(M_{s_1} + M_{11}) \frac{d^2 A_1}{dT^2} + M_{12} \frac{d^2 A_2}{dT^2} + (C_{s_1} + C_{11}^{(0)} + \hat{C}_{11}) \frac{dA_1}{dT} + (C_{12}^{(0)} + \hat{C}_{12}) \frac{dA_2}{dT} + (K_{s_1} + \hat{K}_{11}) A_1 + \hat{K}_{12} A_2 = \langle W_1, \overline{F_{T_1}} \rangle_L \quad (4a)$$
$$ (M_{s2} + M_{22}) \frac{d^2 A_2}{dT^2} + M_{21} \frac{d^2 A_1}{dT^2} + (C_{s2} + C_{22}^{(0)} + \hat{C}_{22}) \frac{dA_2}{dT} \\ + (C_{21}^{(0)} + \hat{C}_{21}) \frac{dA_1}{dT} + (K_{s2} + \hat{K}_{22}) A_2 + \hat{K}_{21} A_1 = \langle W_2, \overline{F_{T2}} \rangle_L \quad (4b) $$

The modal coefficients in Eq. (4) determine the linear stability properties of the coupled system subjected to the axial flow. As a first step toward this goal, in what follows, we derive a new three-dimensional theoretical approach leading to a full analytical expression of coefficients $M_{ij}$ and $C_{ij}^{(0)}$. These expressions are presented in dimensionless forms, see Eq. (51), introducing rescaled quantities as follows:

$$ C_m^{(ij)} = \frac{M_{ij}}{L \rho R_i R_j \pi N_i N_j}, \quad C_v^{(ij)} = \frac{C_{ij}^{(0)}}{L \rho R_i R_j \pi \Omega_j N_i N_j} \quad (5) $$

where $N_i$ is the infinite norm of $W_i$, $R_i$ is the radius of $C_i$, $\Omega_i$ is the angular frequency of motion, and $\rho$ is the fluid density; see Sec. 3 for precise definitions.

This new theory (i) includes the viscous effects of the fluid so that added damping is accounted for, (ii) does not rely on the assumption of a narrow annulus, contrarily to [19–22], so that it applies to all sizes of the fluid gap (some internal components of the Jules Horowitz reactor, such as the middle part of the guide tube, are highly confined; others, such as the upstream and downstream parts of the guide tube, are less so), (iii) does not rely on the slender-body approximation, contrarily to Refs. [9,10,23], so that it applies to elongated cylinders with finite aspect ratios, and (iv) considers arbitrary forced beam vibrations so that it does not only apply to a special case of boundary conditions (the boundary conditions in the Jules Horowitz project being multiple).

This article is organized as follows. In Sec. 3, we introduce the problem and the governing equations. In Sec. 4, we present the theoretical formulation to estimate the modal fluid forces and the dimensionless modal fluid added coefficients. A comparison between theoretical predictions and numerical results is presented in Sec. 5, considering rigid, pinned-pinned, and clamped-free vibrations. Finally, some conclusions are drawn in Sec. 6.

## 3 Definition of the Problem and Governing Equations

We consider the hydrodynamic interaction of two coaxial flexible and circular cylinders $C_i$, ($i=1, 2$), with axis along $e_x$, length $L$, radii $R_i$, and boundaries $\partial C_i$, separated by a viscous fluid layer, as illustrated in Fig. 2. The cylinder $C_i$ is imposed a harmonic displacement $\mathcal{R}\{A_i(T)W_i(X)\}$, where $\mathcal{R}$ is the real part operator, $A_i(T) = Q_i e^{i\Omega_i T}$ is the complex amplitude, $i$ is the unit imaginary number, $W_i(X) = \cos(K_i X + \delta_i)e_y$ is the vibration mode, aligned with the $e_y$ direction orthogonal to $e_x$, $K_i$ its wave number, $\delta_i$ its phase, and $X \in [0, L]$ the axial coordinate.

**3.1 Euler Equations and Boundary Conditions.** The fluid is Newtonian, homogeneous, of volume mass density $\rho$ and kinematic

<figure><img src="image_54.png" /><figcaption>Fig. 2 Two coaxial cylinders separated by a viscous fluid layer</figcaption></figure>

viscosity $\nu$. Assuming that the amplitude of the displacement remains small compared to the thickness of the fluid layer, i.e., $Q_i \ll R_2 - R_1$, the nonlinear term of the Navier–Stokes equations is negligible. It follows that the flow expands as a linear combination of the form $\mathcal{R}\{\sum_{i=1}^2 e^{i\Omega_i T}(\mathbf{V}_i, P_i)\}$ where $(\mathbf{V}_i, P_i)$ is the solution of the Euler equations and boundary conditions

$$ \nabla \cdot \mathbf{V}_i = 0 \quad (6a) $$

$$ i\Omega_i \mathbf{V}_i + \frac{1}{\rho} \nabla P_i - \nu \Delta \mathbf{V}_i = \mathbf{0} \quad (6b) $$

$$ \mathbf{V}_i - i\Omega_i \mathbf{Q}_i \mathbf{W}_i = \mathbf{0} \quad \text{on } \partial C_i \quad (6c) $$

$$ \mathbf{V}_i = \mathbf{0} \quad \text{on } \partial \overline{C_i} \quad (6d) $$

$$ P_i = 0 \quad \text{at } X = 0 \quad \text{and} \quad X = L \quad (6e) $$

The third (resp. fourth) equation expresses the continuity of velocities at the boundary of the moving (resp. stationary) cylinder $C_i$ (resp. $\overline{C_i}$). The fluid force on $C_i$ due to the motion of $C_j$ is $\mathcal{R}\{\sum_{j=1}^2 e^{i\Omega_j T} \mathbf{F}_{ij}\}$ and is the sum of a pressure and a viscous term, i.e.,

$$ \mathbf{F}_{ij} = - \int_{\partial C_i} P_j \mathbf{n}_i \, dS_i + \rho \nu \int_{\partial C_i} [\nabla \mathbf{V}_j + (\nabla \mathbf{V}_j)^T] \cdot \mathbf{n}_i \, dS_i \quad (7) $$

In this equation, $\mathbf{n}_i$ is the outward normal unit vector to $\partial C_i$, $(\nabla \mathbf{V}_j)^T$ is the transposate tensor of $\nabla \mathbf{V}_j$, and $dS_i$ is an infinitesimal surface element of $\partial C_i$. The associated linear fluid force, noted $\overline{\mathbf{F}}_{ij}$, is defined as $\mathbf{F}_{ij} = \int_0^L \overline{\mathbf{F}}_{ij} \, dX$.

**3.2 Dimensionless Euler Equations and Boundary Conditions.** Let $R_2 - R_1$ and $\Omega_i^{-1}$ be the characteristic length and time. We define $w_i$, $v_i$, $p_i$, and $f_{ij}$ as follows:

$$ \begin{gathered} W_i(X) = w_i \left(\frac{X}{L}\right) e_y, \quad \mathbf{V}_i = Q_i \Omega_i \mathbf{v}_i, \quad P_i = \rho (R_2 - R_1) Q_i \Omega_i^2 p_i, \\ \mathbf{F}_{ij} = \rho \pi R_i L R_j Q_i \Omega_j^2 \mathbf{f}_{ij} \end{gathered} \quad (8) $$

which are the dimensionless displacement, fluid velocity, fluid pressure, and fluid force, respectively. The dimensionless vibration mode is given as follows:

$$ w_i(x) = \cos(k_i x + \delta_i) \quad (9) $$

where $k_i = K_i L$ is the dimensionless wave number and $x = X/L$ is the dimensionless axial coordinate. To reduce the number of parameters of the problem, we also introduce the rescaled quantities

$$ l = \frac{L}{R_1}, \quad \epsilon = \frac{R_2}{R_1}, \quad \text{Sk}_i = \frac{(R_2 - R_1)^2 \Omega_i}{\nu} \quad (10) $$

as the aspect ratio, radius ratio, and Stokes number, respectively. Introducing the dimensionless operators $\nabla^* = (R_2 - R_1)\nabla$ and $\Delta^* = (R_2 - R_1)^2\Delta$, the dimensionless Euler equations and boundary conditions are given as follows:

$$ \nabla^* \cdot \mathbf{v}_i = 0 \quad (11a) $$

$$ i\mathbf{v}_i + \nabla^* p_i - \frac{1}{\text{Sk}_i} \Delta^* \mathbf{v}_i = \mathbf{0} \quad (11b) $$

$$ \mathbf{v}_i = i w_i \mathbf{e}_y \quad \text{on } \partial C_i \quad (11c) $$

$$ \mathbf{v}_i = \mathbf{0} \quad \text{on } \partial \overline{C_i} \quad (11d) $$
$$p_i = 0 \quad \text{at } x=0 \quad \text{and} \quad x=1 \qquad (11e)$$

The dimensionless fluid force is expressed as follows:

$$\mathbf{f}_{ij} = (\varepsilon - 1)\varepsilon^{1-j} \frac{1}{\pi} \left[ \int_{\partial C_i} -p_j \mathbf{n}_i \, d\mathbf{s}_i + \frac{1}{\text{Sk}_j} \int_{\partial C_i} \boldsymbol{\tau}_j \, d\mathbf{s}_i \right] \quad (12)$$

where $d\mathbf{s}_i = dS_i/(R_i L)$ and $\boldsymbol{\tau}_j = [\nabla^* \mathbf{v}_j + (\nabla^* \mathbf{v}_j)^T] \cdot \mathbf{n}_i$ is the dimensionless viscous stress vector. The associated dimensionless linear fluid force, noted $\bar{\mathbf{f}}_{ij}$, is defined as $\mathbf{f}_{ij} = \int_0^1 \bar{\mathbf{f}}_{ij} \, dx$.

# 4 Theoretical Solution

## 4.1 Helmholtz Decomposition.
We seek a solution of Eq. (11) as a superposition of an irrotational and a divergence-free flow (Helmholtz decomposition):

$$\mathbf{v}_i = \nabla^* \phi_i + \nabla^* \times \boldsymbol{\psi}_i \qquad (13)$$

where $\phi_i$ and $\boldsymbol{\psi}_i$ some unknown potential and stream functions, respectively. Introducing this decomposition in Eq. (11) yields the dimensionless Euler equations:

$$\Delta^* \phi_i = 0 \qquad (14a)$$

$$\nabla^* \times (\Delta^* \boldsymbol{\psi}_i - i \text{Sk}_i \boldsymbol{\psi}_i) - \text{Sk}_i \nabla^* (i \phi_i + p_i) = \mathbf{0} \qquad (14b)$$

Taking the divergence and the curl of Eq. (14b) yields two equations:

$$p_i = -i\phi_i \qquad (15a)$$

$$\Delta^* \boldsymbol{\psi}_i + \beta_i^2 \boldsymbol{\psi}_i = \mathbf{0} \quad \text{with} \quad \beta_i = \sqrt{-i \text{Sk}_i} \qquad (15b)$$

from which the pressure and the stream functions can be determined.

## 4.2 Fluid Functions.
To solve the Laplace Eq. (14a) and the Helmholtz Eq. (15b), along with the boundary conditions (11c)-(11e), we introduce the dimensionless cylindrical coordinates $(r, \theta, x)$ and the associated physical basis $\mathcal{B} = (\mathbf{e}_r, \mathbf{e}_\theta, \mathbf{e}_x)$. In this basis, the displacement $w_i e_y$ expands as $w_i(\cos(\theta), -\sin(\theta), 0)$ so that $\phi_i$ and $\boldsymbol{\psi}_i$ shall be functions of cos $(\theta)$ and sin $(\theta)$. Also, as the gradient and the curl operators $\nabla^*$ and $\nabla^* \times$ are first-order derivatives in $x$, the functions $\phi_i$ and $\boldsymbol{\psi}_i$ shall be linear combinations of $w_i$ and $dw_i/dx$. The displacement function $w_i(x) = \cos(k_i x + \delta_i)$ is defined for $x \in [0, 1]$ and in general does not satisfy $w_i(0) = w_i(1) = 0$. In order to fulfill the pressure boundary conditions, $p_i(0) = p_i(1) = 0$, it is thus necessary to extend $w_i$ to an odd function of period 2 whose Fourier series converges to $w_i$ for $x \in [0, 1]$ i.e.,

$$\forall x \in [0,1], \quad w_i(x) = \sum_{n=1}^{\infty} w_{in}(x) \qquad (16)$$

where

$$w_{in}(x) = b_n(k_i, \delta_i) \sin(n\pi x) \quad \text{and} \\ b_n(k_i, \delta_i) = 2 \int_0^1 \cos(kx + \delta) \sin(n\pi x) \, dx \qquad (17)$$

Finally, to uncouple the $\mathbf{e}_r$ and $\mathbf{e}_\theta$ projections of the Helmholtz Eq. (15b), one should seek a stream function $\boldsymbol{\psi}_i$, whose radial and orthodial components have the same variations in the polar coordinate $r$ but are in phase quadrature for the angle $\theta$. Accounting for all the aforementioned conditions, we seek fluid functions of the form $(\phi_i, \boldsymbol{\psi}_i) = \sum_{n=1}^{\infty} (\phi_{in}, \boldsymbol{\psi}_{in})$ with

$$\phi_{in}(r, \theta, x) = \Phi_{in}(r) \cos(\theta) w_{in}(x) \qquad (18a)$$

$$\boldsymbol{\psi}_{in}(r, \theta, x) = \begin{pmatrix} 0 \\ 0 \\ \Psi_{in}(r) \sin(\theta) \end{pmatrix} w_{in}(x) + \frac{\varepsilon - 1}{l} \begin{pmatrix} A_{in}(r) \sin(\theta) \\ A_{in}(r) \cos(\theta) \\ 0 \end{pmatrix} \frac{dw_{in}}{dx}(x) \qquad (18b)$$

where the matrix notation in Eq. (18b) is used to express the components of $\boldsymbol{\psi}_{in}$ in the basis $\mathcal{B}$. Plugging Eqs. (17) and (18a) into Eq. (15a) yields a $n$ th Fourier harmonics for the pressure field of the form

$$p_{in}(r, \theta, x) = -i \Phi_{in}(r) \cos(\theta) b_n(k_i, \delta_i) \sin(n\pi x) \qquad (19)$$

so that the pressure boundary conditions (11e) are automatically satisfied. The dimensionless differential operators applied to Eq. (18) are given as follows:

$$\nabla^* \phi_{in} = \begin{pmatrix} \frac{d\Phi_{in}}{dr} \cos(\theta) \\ -\frac{1}{r} \Phi_{in} \sin(\theta) \\ 0 \end{pmatrix} w_{in} + \frac{\varepsilon - 1}{l} \begin{pmatrix} 0 \\ 0 \\ \Phi_{in} \cos(\theta) \end{pmatrix} \frac{dw_{in}}{dx} \qquad (20a)$$

$$\Delta^* \phi_{in} = \left[ \frac{d^2 \Phi_{in}}{dr^2} + \frac{1}{r} \frac{d\Phi_{in}}{dr} + \left( \gamma_n^2 - \frac{1}{r^2} \right) \Phi_{in} \right] \cos(\theta) w_{in} \qquad (20b)$$

$$\nabla^* \times \boldsymbol{\psi}_{in} = \begin{pmatrix} \left(\frac{1}{r}\Psi_{in} - \gamma_n^2 A_{in}\right) \cos(\theta) \\ -\left(\frac{d\Psi_{in}}{dr} - \gamma_n^2 A_{in}\right) \sin(\theta) \\ 0 \end{pmatrix} w_{in} + \frac{\varepsilon - 1}{l} \begin{pmatrix} 0 \\ 0 \\ \frac{dA_{in}}{dr} \cos(\theta) \end{pmatrix} \frac{dw_{in}}{dx} \qquad (20c)$$

$$\Delta^* \boldsymbol{\psi}_{in} = \begin{pmatrix} 0 \\ 0 \\ \left[ \frac{d^2 \boldsymbol{\psi}_{in}}{dr^2} + \frac{1}{r} \frac{d \boldsymbol{\psi}_{in}}{dr} + \left( \gamma_n^2 - \frac{1}{r^2} \right) \boldsymbol{\psi}_{in} \right] \sin(\theta) \end{pmatrix} w_{in} \qquad (20d)$$

$$+ \frac{\varepsilon - 1}{l} \begin{pmatrix} \left( \frac{d^2 A_{in}}{dr^2} + \frac{1}{r} \frac{dA_{in}}{dr} + \gamma_n^2 A_{in} \right) \sin(\theta) \\ \left( \frac{d^2 A_{in}}{dr^2} + \frac{1}{r} \frac{dA_{in}}{dr} + \gamma_n^2 A_{in} \right) \cos(\theta) \\ 0 \end{pmatrix} \frac{dw_{in}}{dx} \qquad (20e)$$

where

$$\gamma_n = i \frac{n\pi}{l} (\varepsilon - 1) \qquad (21)$$

The magnitude of $\gamma_n$ is $(n\pi/L)/(R_2 - R_1)$ so that it expresses the ratio between the wave-number of the $n$ th Fourier component of the displacement $w_i$ and the thickness $R_2 - R_1$ of the duct.

From Eq. (20b), the Laplace Eq. (14a) is written as follows:

$$r^2 \frac{d^2 \Phi_{in}}{dr^2} + r \frac{d \Phi_{in}}{dr} + \left[ (\gamma_n r)^2 - 1 \right] \Phi_{in} = 0 \qquad (22)$$

From Eq. (20e), the Helmholtz Eq. (15b) projected onto the basis $\mathcal{B}$ has the same component along the $\mathbf{e}_r$ and $\mathbf{e}_\theta$ directions:

$$r^2 \frac{d^2 A_{in}}{dr^2} + r \frac{dA_{in}}{dr} + (\alpha_{in} r)^2 A_{in} = 0 \qquad (23)$$

whereas the component along $\mathbf{e}_x$ is given as follows:

$$r^2 \frac{d^2 \boldsymbol{\psi}_{in}}{dr^2} + r \frac{d \boldsymbol{\psi}_{in}}{dr} + \left[ (\alpha_{in} r)^2 - 1 \right] \boldsymbol{\psi}_{in} = 0 \qquad (24)$$
where

$$ \alpha_{in}^2 = \gamma_n^2 + \beta_i^2 = \gamma_n^2 - iSk_i \quad (25) $$

The solutions of (22)–(24) are linear combinations of Hankel functions of the first and second kind, so that the fluid functions are written as follows:

$$ \begin{bmatrix} \Phi_{in}(r) \\ \Psi_{in}(r) \\ A_{in}(r) \end{bmatrix} = [B(\gamma_n r, \alpha_{in} r)][C(r_i, \gamma_n, \alpha_{in})] \quad (26) $$

where $r_i = R_i/(R_2 - R_1)$ is the dimensionless radial coordinate of the points lying on the boundary $\partial C_i$ of the moving cylinder and

$$ [B(\gamma_n r, \alpha_{in} r)] = \begin{bmatrix} H_1^{(1)}(\gamma_n r) & H_1^{(2)}(\gamma_n r) & 0 & 0 & 0 & 0 \\ 0 & 0 & H_1^{(1)}(\alpha_{in} r) & H_1^{(2)}(\alpha_{in} r) & 0 & 0 \\ 0 & 0 & 0 & 0 & H_0^{(1)}(\alpha_{in} r) & H_0^{(2)}(\alpha_{in} r) \end{bmatrix} \quad (27) $$

For a given set of parameters $r_i, \gamma_n$, and $\alpha_{in}$, the vector $[C(r_i, \gamma_n, \alpha_{in})]$ is composed of six constants determined by the conditions (11c) and (11d) at the cylinders boundaries. In terms of the fluid functions, the projection of Eqs. (11c) and (11d) onto the basis $\mathcal{B}$

yields two sets of boundary conditions:

$$ \left\{ \begin{aligned} & \left| \frac{d\Phi_{in}}{dr} + \frac{1}{r}\Psi_{in} - \gamma_n^2 A_{in} \right| = i \\ & \left| \frac{1}{r}\Phi_{in} + \frac{d\Psi_{in}}{dr} - \gamma_n^2 A_{in} \right| = i \text{ at } r=r_i, \\ & \left| \Phi_{in} + \frac{dA_{in}}{dr} \right| = 0 \end{aligned} \right. \quad (28) $$

where $\bar{r}=(r_1+r_2)-r_i$ is the radial coordinate of the points lying on the boundary $\partial \bar{C}_i$ of the stationary cylinder. The solution of the linear system (28) is

$$ [C(r_i, \gamma_n, \alpha_{in})] = \begin{bmatrix} [M(r_i, \gamma_n, \alpha_{in})] \\ [M(\bar{r}_i, \gamma_n, \alpha_{in})] \end{bmatrix}^{-1} [Q] \quad (29) $$

where $[Q]=i[1 \ 1 \ 0 \ 0 \ 0 \ 0]^T$ and

$$ [M(r, \gamma, \alpha)] = \begin{bmatrix} \gamma H_1^{(1)\prime}(\gamma r) & \gamma H_1^{(2)\prime}(\gamma r) & \frac{H_1^{(1)}(\alpha r)}{r} & \frac{H_1^{(2)}(\alpha r)}{r} & -\gamma^2 H_0^{(1)}(\alpha r) & -\gamma^2 H_0^{(2)}(\alpha r) \\ \frac{H_1^{(1)}(\gamma r)}{r} & \frac{H_1^{(2)}(\gamma r)}{r} & \alpha H_1^{(1)\prime}(\alpha r) & \alpha H_1^{(2)\prime}(\alpha r) & -\gamma^2 H_0^{(1)}(\alpha r) & -\gamma^2 H_0^{(2)}(\alpha r) \\ H_1^{(1)}(\gamma r) & H_1^{(2)}(\gamma r) & 0 & 0 & \alpha H_0^{(1)\prime}(\alpha r) & \alpha H_0^{(2)\prime}(\alpha r) \end{bmatrix} \quad (30) $$

In Eq. (30), the prime subscript notation on the Hankel functions stands for the derivative with the $r$ variable. Finally, plugging Eq. (29) into Eq. (26) yields an expression of the fluid functions

$$ \begin{bmatrix} \Phi_{in}(r) \\ \Psi_{in}(r) \\ A_{in}(r) \end{bmatrix} = [B(\gamma_n r, \alpha_{in} r)] \begin{bmatrix} [M(r_i, \gamma_n, \alpha_{in})] \\ [M(\bar{r}_i, \gamma_n, \alpha_{in})] \end{bmatrix}^{-1} [Q] \quad (31) $$

### 4.2.1 Inviscid Limit, $Sk_i \to \infty$, i.e., $|\alpha_{in}| \to \infty$
The particular case of an inviscid fluid is recovered from the present viscous theory by taking the limit $Sk_i \to \infty$ so that $|\alpha_{in}| \to \infty$. In such a case, the stream function vanishes, i.e., $\psi_i = 0$, whereas the potential function is given as follows:

$$ \phi_{in}^{(\text{POT})}(r, \theta, x) = \Phi_{in}^{(\text{POT})}(r) \cos(\theta) w_{in}(x) \quad (32) $$

where $\Phi_{in}^{(\text{POT})}$ is the solution of the Laplace equation:

$$ r^2 \frac{d^2 \Phi_{in}^{(\text{POT})}}{dr^2} + r \frac{d \Phi_{in}^{(\text{POT})}}{dr} \left[ (\gamma_n r)^2 - 1 \right] \Phi_{in}^{(\text{POT})} = 0 \quad (33) $$

For the inviscid problem, the boundary conditions (28) reduce to

$$ \frac{d\Phi_{in}^{(\text{POT})}}{dr} = i \text{ at } r = r_i, \quad \frac{d\Phi_{in}^{(\text{POT})}}{dr} = 0 \text{ at } r = \bar{r}_i \quad (34) $$

so that the inviscid potential function is given as follows:

$$ \Phi_{in}^{(\text{POT})}(r) = [B^{(\text{POT})}(\gamma_n r)] \begin{bmatrix} [M^{(\text{POT})}(r_i, \gamma_n)] \\ [M^{(\text{POT})}(\bar{r}_i, \gamma_n)] \end{bmatrix}^{-1} [Q^{(\text{POT})}] \quad (35) $$

where $[Q^{(\text{POT})}] = i[1 \ 0]^T$ and

$$ [B^{(\text{POT})}(\gamma_n r)] = [H_1^{(1)}(\gamma_n r) H_1^{(2)}(\gamma_n r)] \quad (36) $$

$$ [M^{(\text{POT})}(r, \gamma)] = [\gamma H_1^{(1)\prime}(\gamma r) \gamma H_1^{(2)\prime}(\gamma r)] \quad (37) $$

### 4.2.2 Slender-Body Limit, $l \to \infty$, i.e., $|\gamma_n| \to 0$
The particular case of infinitely long cylinders has been covered by Ref. [10] who solved the two-dimensional viscous fluid problem by the mean of a stream function solution of a Laplace equation. The solution of Ref. [10] is recovered from the present three-dimensional theory by taking the limit $l \to \infty$ so that $|\gamma_n| \to 0$. In such a case, one has to disregard the pressure boundary conditions (11e) and consequently the Fourier expansion of $w_i$ introduced in Eq. (16). In such a case, the fluid functions reduce to

$$ \phi_i^{(2D)}(r, \theta) = \Phi_i^{(2D)}(r) \cos(\theta) \quad (38a) $$

$$ \psi_i^{(2D)}(r, \theta) = \Psi_i^{(2D)}(r) \sin(\theta) e_x \quad (38b) $$

where $\Psi_i^{(2D)}$ and $\psi_i^{(2D)}$ are solutions of the Laplace and the Helmholtz equations:

$$ r^2 \frac{d^2 \Phi_i^{(2D)}}{dr^2} + r \frac{d \Phi_i^{(2D)}}{dr} - \Phi_i^{(2D)} = 0 \quad (39a) $$

$$ r^2 \frac{d^2 \Psi_i^{(2D)}}{dr^2} + r \frac{d \Psi_i^{(2D)}}{dr} + [ (\beta_i r)^2 - 1 ] \Psi_i^{(2D)} = 0 \quad (39b) $$
<figure><img src="image_55.png" /><figcaption>Fig. 3 Magnitude and angle of $-i\Phi_{in}$ in the thickness of the duct, for $i=1$ (black), $i=2$ (red) and $n \in \{1, 2, 3\}$. The solid lines correspond to a Stokes number $Sk_i=250$ and an aspect ratio $l=1$. The dashed lines correspond to the inviscid limit $Sk_i \to \infty$. The dotted lines correspond to the limit of infinitely long cylinders, $l \to \infty$, for $Sk_i=250$. On the right picture, the solid and dotted lines are indistinguishable. The radius ratio is $\epsilon = 1.1$. (Color version online.)</figcaption></figure>

For the two-dimensional fluid problem, the boundary conditions (28) reduce to

$$ 
\left\{ 
\begin{array}{l}
\displaystyle \left| \frac{d\Phi_i^{(2D)}}{dr} + \frac{1}{r} \Psi_i^{(2D)} \right| = i \\
\\
\displaystyle \frac{1}{r} \Phi_i^{(2D)} + \frac{d\Psi_i^{(2D)}}{dr} = i \\
\\
\text{at } r = r_i, \quad \left| \frac{d\Phi_i^{(2D)}}{dr} + \frac{1}{r} \Psi_i^{(2D)} \right| = 0 \\
\\
\displaystyle \frac{1}{r} \Phi_i^{(2D)} + \frac{d\Psi_i^{(2D)}}{dr} = 0
\end{array}
\right.
\quad (40) 
$$

so that the fluid functions are given as follows:

$$ \begin{bmatrix} \Phi_i^{(2D)}(r) \\ \Psi_i^{(2D)}(r) \end{bmatrix} = [B^{(2D)}(r, \beta_i r)] \begin{bmatrix} [M^{(2D)}(r_i, \beta_i)]^{-1} \\ [M^{(2D)}(\bar{r}_i, \beta_i)] \end{bmatrix} [Q^{(2D)}] \quad (41) $$

<figure><img src="image_56.png" /></figure>

<figure><img src="image_57.png" /><figcaption>Fig. 4 Color map of the dimensionless pressure field $\mathfrak{R}\{p_i\}$ generated by the motion of $C_i$, which vibrates in the first mode of a clamped-free beam. The aspect ratio is $l=4$. The Stokes number is number $Sk_i=250$. The radius ratio is $\epsilon = 1.1$.</figcaption></figure>

where $[Q^{(2D)}] = i[\begin{matrix} 1 & 1 & 0 & 0 \end{matrix}]^T$ and

$$ [B^{(2D)}(r, \beta_i r)] = \begin{bmatrix} r & \frac{1}{r} & 0 & 0 \\ 0 & 0 & H_1^{(1)}(\beta_i r) & H_1^{(2)}(\beta_i r) \end{bmatrix} \quad (42a) $$

$$ [M^{(2D)}(r, \beta)] = \begin{bmatrix} 1 & -\frac{1}{r^2} & \frac{H_1^{(1)}(\beta r)}{r} & \frac{H_1^{(2)}(\beta r)}{r} \\ 1 & \frac{1}{r^2} & \beta H_1^{(1)'}(\beta r) & \beta H_1^{(2)'}(\beta r) \end{bmatrix} \quad (42b) $$

By having determined the fluid functions, we now proceed with analyzing the variations of $-i\Phi_{in}$, which is, according to Eq. (19), directly related to the pressure field $p_i$ in the duct. In Fig. 3, we show that the magnitude and the angle of $-i\Phi_{in}$ do not vary that much with the radial coordinate $r$, so that $p_i$ is at the first-order constant in the thickness of the duct. The amplitude of $p_i$ increases with the aspect ratio $l$ of the cylinders, whereas the Stokes number $Sk_i$ mainly affects the phase shift between $p_i$ and $w_i$. Finally, for two cylinders vibrating in the same mode, i.e., $w_1 = w_2$, with the same frequency, i.e., $Sk_1 = Sk_2$, the pressures $p_1$ and $p_2$ have similar amplitudes and are in phase opposition. As a consequence, $p_1$ and $p_2$ tend to cancel each other out so that the total pressure in the duct is close zero. In Fig. 4, we show the pressure field in the half-plane $(e_x, e_y)$, $y>0$, when $C_j$ vibrates in the first mode of a clamped-free configuration. As expected from the variations of $-i\Phi_{in}$, we recover that the pressure mostly varies along the $x$-direction and that $\mathfrak{R}\{p_1\} \approx -\mathfrak{R}\{p_2\}$.

### 4.3 Linear Fluid Forces.
The viscous stress vector at the cylinders boundaries is obtained by substituting the Helmholtz decomposition (13) in $\tau_j = [\nabla^* v_j + (\nabla^* v_j)^T] \cdot e_r$. By using the boundary conditions (28), the components of $\tau_{jn}$ in the basis $\mathcal{B}$ simplify to

$$ \tau_{jn} = \begin{pmatrix} 0 \\ (\alpha_{jn}^2 \Psi_{jn} - \gamma_n^2 \Phi_{jn}) \sin(\theta) \\ 0 \end{pmatrix} w_{jn} + \frac{\epsilon - 1}{l} \times \begin{pmatrix} 0 \\ (2 \frac{d\Phi_{jn}}{dr} + \frac{\Phi_{jn} + \Psi_{jn}}{r} - (\alpha_{jn}^2 + \gamma_n^2)A_{jn}) \cos(\theta) \\ 0 \end{pmatrix} \frac{dw_{jn}}{dx} \quad (43) $$
<figure><img src="image_58.png" /><figcaption>Fig. 5 Magnitude and angle of $h_n^{(ij)}$ as functions of $|\gamma_n|$, for $i=j=1$ (black), $i=j=2$ (red) and $i \neq j$ (blue). The solid lines correspond to a Stokes number $Sk_j = 250$. The dashed lines (black and red are indistinguishable) correspond to the inviscid limit $Sk_j \to \infty$. The dotted lines correspond to the limit of infinitely long cylinders, $l \to \infty$, for $Sk_j = 250$. The dashed-dotted line corresponds to the limit of infinitely short cylinders, $l \to 0$. On the right picture, the solid, dotted, and dash-dotted lines are indistinguishable. The radius ratio is $\epsilon = 1.1$. (Color version online.)</figcaption></figure>

Plugging Eq. (43) and the pressure equation (15a) into Eq. (12) and assuming small displacement and deformation of the cylinders, i.e., $\mathbf{n}_i \approx (-1)^{i+1}\mathbf{e}_r$, show that the linear fluid force is aligned with the direction of motion, i.e., $\overline{f}_{ij} = \overline{f}_{ij}\mathbf{e}_y$, with

$$ \overline{f}_{ij} = (-1)^{i+1}(\epsilon - 1)\epsilon^{1-j} \sum_{n=1}^{\infty} i \left[ \Phi_{jn}(r_i) + \frac{\gamma_n^2 \Phi_{jn}(r_i) - \alpha_{jn}^2 \Psi_{jn}(r_i)}{\gamma_n^2 - \alpha_{jn}^2} \right] w_{jn} \quad (44) $$

The term $i\Phi_{jn}(r_i)$ arises from the integration of the pressure field on the circumference of $\partial C_i$, whereas the second term entering in the sum sign arises from the integration of the $\mathbf{e}_y$ component of the viscous stress vector $\boldsymbol{\tau}_{jn}$. From Eq. (44), the linear fluid force $\overline{f}_i = \overline{f}_{i1} + \overline{f}_{i2}$ acting on the cylinder $C_i$ is written as follows:

$$ \left[ \frac{\overline{f}_1}{\overline{f}_2} \right] = (\epsilon - 1) \sum_{n=1}^{\infty} \left[ \begin{matrix} h_n^{(11)} & h_n^{(12)} \\ h_n^{(21)} & h_n^{(22)} \end{matrix} \right] \left[ \begin{matrix} b_{1n} \\ b_{2n} \end{matrix} \right] \sin(n\pi x) \quad (45) $$

where $b_{jn} = b_n(k_j, \delta_j)$ is the n-th Fourier coefficient of $w_j = \cos(k_j x + \delta_j)$, see Eq. (17) and

$$ h_n^{(ij)} = h_n^{(ji)} = (-1)^{i+1} \epsilon^{1-j} \left[ \Phi_{jn}(r_i) + \frac{\gamma_n^2 \Phi_{jn}(r_i) - \alpha_{jn}^2 \Psi_{jn}(r_i)}{\gamma_n^2 - \alpha_{jn}^2} \right] \quad (46) $$

In Eq. (45), $h_n^{(ij)}$ acts as a transfer function, which relates the fluid force on $C_i$ to the n-th Fourier harmonics of $w_j$. Since $\alpha_{jn}$ depends on the Stokes number, $h_n^{(ij)}$ is actually a function of $(\epsilon, \gamma_n, Sk_j)$.

In Fig. 5, we show the evolution of the magnitude and the angle of $h_n^{(ij)}$ as functions of $|\gamma_n|$. We observe that the magnitude decreases monotonically from the asymptotic limit of an infinitely long cylinder ($|\gamma_n| \to 0$, i.e., $l \to \infty$) to the limit of an infinitely short cylinder ($l \gamma_n \to \infty$, i.e., $l \to 0$). This decrease indicates that the n-th Fourier harmonics of the vibration mode is all the more amplified that the aspect ratio $l$ of the cylinders is large. For a given aspect ratio $l$ of the cylinders, this decrease also indicates that the fluid has a different effect on the Fourier harmonics of the vibration mode, amplifying its low-order harmonics more than its large ones. As a direct consequence, the shape of the linear fluid force $\overline{f}_{ij}$ is generally different from that of the vibration mode $w_j$. The proximity of the curves drawn for $Sk_j=250$ and $Sk_j \to \infty$ (inviscid limit) shows that the Stokes number has a very weak effect on the magnitude of $h_n^{(ij)}$. Actually, the Stokes number mainly affects the angle of $h_n^{(ij)}$, which is poorly sensitive on $|\gamma_n|$. Finally, we shall note that the magnitude of $h_n^{(ij)}$ does not depend that much on the values of the indices $i$ and $j$. It follows that the linear fluid force acting on

$C_i$, due to the motion of $C_j$, has the same magnitude as the linear fluid force acting on $C_j$, due to the motion of $C_i$, providing that the two cylinders vibrate in the same mode. Still, the angle of $h_n^{(ij)}$ indicates that $\overline{f}_{ii}$ and $\overline{f}_{ij}$, $i \neq j$, are in phase opposition.

### 4.4 Modal Fluid Forces and Matrix of the Dimensionless Modal Fluid Added Coefficients.
The modal fluid force $\langle \overline{w}_i, \overline{f}_i \rangle$ is defined as the projection of $\overline{f}_i$ onto a vibration mode $\overline{w}_i(x)$, i.e.,

$$ \langle \overline{w}_i, \overline{f}_i \rangle = \int_0^1 \overline{w}_i(x) \overline{f}_i(x) dx \quad (47) $$

where

$$ \overline{w}_i(x) = \cos(k_i x + \delta_i) \quad (48) $$

From Eq. (45), it follows that the modal fluid forces are written as follows:

$$ \left[ \begin{matrix} \langle \overline{w}_1, \overline{f}_1 \rangle \\ \langle \overline{w}_2, \overline{f}_2 \rangle \end{matrix} \right] = (\epsilon - 1) \sum_{n=1}^{\infty} \int_0^1 \left[ \begin{matrix} h_n^{(11)} \overline{w}_1(x) & h_n^{(12)} \overline{w}_1(x) \\ h_n^{(21)} \overline{w}_2(x) & h_n^{(22)} \overline{w}_2(x) \end{matrix} \right] \times \begin{matrix} b_{1n} \\ b_{2n} \end{matrix} \sin(n\pi x) dx \quad (49) $$

Introducing the n-th Fourier coefficient $\overline{b}_{in}$ of $\overline{w}_i$, Eq. (49) simplifies to

$$ \left[ \begin{matrix} \langle \overline{w}_1, \overline{f}_1 \rangle \\ \langle \overline{w}_2, \overline{f}_2 \rangle \end{matrix} \right] = \frac{(\epsilon - 1)}{2} \sum_{n=1}^{\infty} \left[ \begin{matrix} h_n^{(11)} \overline{b}_{1n} & h_n^{(12)} \overline{b}_{1n} \\ h_n^{(21)} \overline{b}_{2n} & h_n^{(22)} \overline{b}_{2n} \end{matrix} \right] \left[ \begin{matrix} b_{1n} \\ b_{2n} \end{matrix} \right] \quad (50) $$

By splitting $\langle \overline{w}_i, \overline{f}_i \rangle$ as $\langle \overline{w}_i, \overline{f}_{i1} \rangle + \langle \overline{w}_i, \overline{f}_{i2} \rangle$, we introduce the dimensionless normalized modal mass and damping coefficients, defined as $C_m^{(ij)} - iC_v^{(ij)} = \langle \overline{w}_i, \overline{f}_{ij} \rangle / (N_i N_j)$, where $N_j = \sup_{x \in [0,1]} (|w_j|)$ is the infinite norm of the mode $w_j$. From Eq. (50), the modal matrix
eventually is written as follows:

$$
\begin{bmatrix} C_m^{(11)} - iC_v^{(11)} & C_m^{(12)} - iC_v^{(12)} \\ C_m^{(21)} - iC_v^{(21)} & C_m^{(22)} - iC_v^{(22)} \end{bmatrix} = \frac{(\epsilon - 1)}{2} \\
\times \sum_{n=1}^{\infty} \begin{bmatrix} h_n^{(11)} \frac{b_{1n}b_{1n}}{N_1N_1} & h_n^{(12)} \frac{b_{1n}b_{2n}}{N_1N_2} \\ h_n^{(21)} \frac{b_{2n}b_{1n}}{N_2N_1} & h_n^{(22)} \frac{b_{2n}b_{2n}}{N_2N_2} \end{bmatrix} \quad (51)
$$

**4.5 Euler-Bernoulli Beam Vibration Modes.** So far the theory has been derived for a vibration mode $w_j$ of the form $w_j = \cos(k_j x + \delta_j)$. Considering the classical boundary conditions of an Euler-Bernoulli beam, the shape of $w_j$ is expressed as follows:

$$ w_j(x) = \chi_j^{(1)} \cosh(\lambda_j x) + \chi_j^{(2)} \cos(\lambda_j x) + \chi_j^{(3)} \sinh(\lambda_j x) + \chi_j^{(4)} \sin(\lambda_j x) \quad (52) $$

so that $b_{jn}$ is obtained from the superposition of trigonometric modes with different wave numbers, leading to

$$ b_{jn} = \chi_j^{(1)} b_n(i\lambda_j, 0) + \chi_j^{(2)} b_n(\lambda_j, 0) + \chi_j^{(3)}(-i)b_n(i\lambda_j, -\frac{\pi}{2}) \\ + \chi_j^{(4)} b_n(\lambda_j, -\frac{\pi}{2}) \quad (53) $$

The expansion of Eq. (53) then reduces to

$$
\begin{align} b_{jn} = & -2 \frac{\pi ((\pi^2 n^2 - \lambda_j^2) \cosh(\lambda_j) (-1)^n - \pi^2 n^2 + \lambda_j^2) n_j^{(1)}}{\pi^4 n^4 - \lambda_j^4} \chi_j^{(1)} \nonumber \\ & -2 \frac{\pi ((\pi^2 n^2 + \lambda_j^2) \cos(\lambda_j) (-1)^n - \pi^2 n^2 - \lambda_j^2) n_j^{(2)}}{\pi^4 n^4 - \lambda_j^4} \chi_j^{(2)} \nonumber \\ & -2 \frac{\pi (\pi^2 n^2 - \lambda_j^2) \sinh(\lambda_j) (-1)^n n_j^{(3)}}{\pi^4 n^4 - \lambda_j^4} \chi_j^{(3)} \nonumber \\ & -2 \frac{\pi (\pi^2 n^2 + \lambda_j^2) \sin(\lambda_j) (-1)^n n_j^{(4)}}{\pi^4 n^4 - \lambda_j^4} \chi_j^{(4)} \tag{54} \end{align}
$$

for $\lambda_j \neq n\pi$ and

$$ b_{jn} = \frac{-\cosh(n\pi)(-1)^n+1}{n\pi} \chi_j^{(1)} - \frac{(-1)^n\sinh(n\pi)\chi_j^{(3)} + \chi_j^{(4)}}{n\pi} \quad (55) $$

for $\lambda_j = n\pi$. Obviously, the mode used for the projection is written as follows:

$$ \overline{w_i}(x) = \chi_i^{(1)} \cosh(\overline{\lambda_i}x) + \chi_i^{(2)} \cos(\overline{\lambda_i}x) + \chi_i^{(3)} \sinh(\overline{\lambda_i}x) + \chi_i^{(4)} \sin(\overline{\lambda_i}x) \quad (56) $$

so that the Fourier coefficient $\overline{b_{in}}$ appearing in Eq. (51) is simply given by either Eqs. (54) or (55), replacing $\lambda_j$ by $\overline{\lambda_i}$, $\chi_j^{(3)}$ by $\overline{\chi_i^{(3)}}$, and $\chi_j^{(4)}$ by $\overline{\chi_i^{(4)}}$.

<figure><img src="image_59.png" /><figcaption>**Fig. 6** Modal self and cross-added mass and damping coefficients as functions of the aspect ratio, $l$, for $i=j=1$ (black), $i=j=2$ (red), and $i \neq j$ (blue and purple). The cylinder $C_j$ vibrates in a rigid mode. The solid lines correspond to a Stokes number $Sk_j = 250$. The dashed lines (black and red are indistinguishable) correspond to the inviscid limit $Sk_j \to \infty$. The dotted lines correspond to the limit of infinitely long cylinders, $l \to \infty$, for $Sk_j = 250$. The radius ratio is $\epsilon = 1.1$. (Color version online.)</figcaption></figure>
<figure><img src="image_60.png" /><figcaption>Fig. 7 Modal self and cross-added mass and damping coefficients as functions of the aspect ratio, *l*, for *i*=*j*=1 (black), *i*=*j*=2 (red), and *i*≠*j* (blue and purple). The cylinder $C_j$ vibrates in the first mode of a pinned-pinned beam. The solid lines correspond to a Stokes number $Sk_j$ = 250. The dashed lines (black and red are indistinguishable) correspond to the inviscid limit $Sk_j \to \infty$. The dotted lines correspond to the limit of infinitely long cylinders, *l* → ∞, for $Sk_j$ = 250. The radius ratio is $\epsilon$ = 1.1. (Color version online.)</figcaption></figure>

Equation (51) along with Eqs. (54) and (55) constitute a full ana-
lytical expression of the modal added coefficients. This expression
shows that the added coefficients depend on the characteristics of
the imposed vibrations, through the Fourier coefficients $b_{jn}$, on
the aspect ratio *l*, on the radius ratio $\epsilon$, and on the Stokes number,
through the transfer function term $h_n^{(j)}$, given by Eq. (46). In
Figs. 6–8, we show the evolution of the modal coefficients
(taking $\overline{w_i} = w_i$), as functions of the aspect ratio *l*, when $C_j$ ($j$=1
or $j$=2) vibrates in a rigid mode (Fig. 6), i.e., $w_j=1$, in the first
mode of a pinned-pinned beam (Fig. 7), i.e., $w_j=\sin(\pi x)$, and in
the first mode of a clamped-free beam (Fig. 8), i.e., $w_j$ given by
Eq. (52), where $\lambda_j=1.8751$, $\chi_j^{(1)}=1$, $\chi_j^{(2)}=-1$, $\chi_j^{(3)}=-\sigma$,
$\chi_j^{(4)}=\sigma$, and $\sigma=0.73410$. In all cases, we observe that the self
LTE (resp. cross) added coefficients are positive (resp. negative) mono-
tonic increasing (resp. decreasing) functions of the aspect ratio, with
an horizontal asymptote corresponding to the slender-body limit,
see Sec. 4.2.2, as $l \to \infty$. More importantly, the comparison of the
curves plotted for $Sk=250$ and $Sk \to \infty$ (inviscid fluid) shows
that the added mass coefficients (self and cross) are weakly influ-
enced by the viscosity, contrary to the damping coefficients
(which obviously vanish as $Sk \to \infty$).

Finally, in Fig. 9, we show that the linear fluid forces $\overline{f}_{ij}$ do not
have the same shape as the forced vibration mode (here, the first
mode of a clamped-free beam), especially for low values of the
aspect ratio $l$. This is of course related to the fact that the pressure

field is imposed to vanish at x=0 and x=1, which is not necessarily
the case of the imposed displacement (for a clamped-free beam,
$w_j(1) \neq 0$). In such a case, the transfer function $h_n^{(j)}$, see Eq. (46)
and the discussion in Sec. 4.3, affects differently the Fourier har-
monics of the imposed vibration mode, so that the fluid makes it
possible to couple modes with different wave numbers.

After presenting the theoretical approach, we now proceed to its
evaluation by comparing the predictions with numerical
simulations.

5 Numerical Assessment

The numerical simulations of the Navier–Stokes equations are per-
formed with the open-source code TRIOCFD, a computational fluid
dynamics code for incompressible or quasi-incompressible fluids,
coupled with an arbitrary Lagrange–Euler module. In our simulations,
the Navier–Stokes equations are time discretized using a backward
Euler scheme of the first order. A hybrid finite element volume
method, see Refs. [24,25], for unstructured grids is used for the
space discretization. The nonlinear convective term is approximated
using a second-order MUSCL scheme upwind. Discrete linear
systems are solved using iterative solvers from the PETSc library.
The discretization of the domain is based on MGCADSurf-MGTetra
mesh generator of the SALOME platform. The reader is referred to
Refs. [23,24,26,27] for additional details on the code TRIOCFD.
<figure><img src="image_61.png" /><figcaption>**Fig. 8** Modal self and cross-added mass and damping coefficients as functions of the aspect ratio, *l*, for *i*=*j*=1 (black), *i*=*j*=2 (red), and *i*≠*j* (blue and purple). The cylinder *C*<sub>j</sub> vibrates in the first mode of a clamped-free beam. The solid lines correspond to a Stokes number *Sk*<sub>j</sub>=250. The dashed lines (black and red are indistinguishable) correspond to the inviscid limit *Sk*<sub>j</sub> → ∞. The dotted lines correspond to the limit of infinitely long cylinders, *l* → ∞, for *Sk*<sub>j</sub>=250. The radius ratio is *ε* = 1.1. (Color version online.)</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>**Fig. 9** Magnitude of $\bar{f}_{ij}$ as function of the dimensionless axial coordinate $x$, for an aspect ratio $l \in \{1, 4, 16, 64\}$, for $i=j=1$ (black), $i=j=2$ (red), and $i \neq j$ (blue). The cylinder $C_j$ vibrates in the first mode of a clamped-free beam. The solid green line indicates the shape of the vibration mode. The solid lines correspond to a Stokes number $Sk_j=250$. The dashed lines (black and red are indistinguishable) correspond to the inviscid limit $Sk_j \to \infty$. The radius ratio is $\epsilon = 1.1$. (Color version online.)</figcaption></figure>

**5.1 Numerical Setup.** The numerical domain consists of two coaxial cylinders $C_j$ with radii $R_j$ and lengths $L_j$, see Fig. 2. Either the inner or the outer cylinder is imposed a displacement of the form $W(X)\sin(\Omega T)$, whereas a zero pressure field is set on $X=0$ and $X=L$. The fluid domain is closed with wall boundaries (zero velocity) in the $e_y$ and $-e_y$ directions. Simulations are performed under an imposed displacement of amplitude $W=5 \cdot 10^{-5}$ [uol] (units of length) and a forcing frequency $\Omega/(2\pi)=10$ [uof] (units of frequency). The inner radius is $R_1=0.02$ [uol], and the radius of the outer cylinder is $R_2=0.022$ [uol]. The fluid properties are $\nu = 1.007 \times 10^{-6}$ [uol]$^2$/ [uot] and $\rho = 1000$ [uom]/[uol]$^3$, where [uom] stands for units of mass. Given this set of parameters, the Stokes number is sufficiently low, i.e., $Sk_j=250$, to be representative of a viscous fluid.

**5.2 Fluid Force Matrix Structure.** As a first evaluation of the theoretical predictions, we consider the case in which either the inner or the outer cylinder is imposed a rigid, a pinned-pinned, or a clamped-free mode of vibration (first mode). The aspect ratio is fixed to $l=15$, setting $L_1=0.3$ [uol] and $L_2=0.31$ [uol] in the numerical simulations. For such an aspect ratio, the mesh is composed of 9 million tetrahedrons; the time-step is $1 \times 10^{-4}$ [uot] (units of times) and the CPU time is approximately 5 h. The modal added coefficients are obtained from the projection of the
<table><caption>Table 1 Cylinders vibrating in a rigid mode</caption><thead><tr><th>Rigid</th><th>Present theory</th><th>Numerics TRIOCFD</th><th>Relative deviation (%)</th></tr></thead><tbody><tr><td></td><td>( 9.93 - 1.07i -9.80 + 0.968i )</td><td>( 9.92 - 1.12i -9.84 + 1.03i )</td><td>( 0.10 + 4.5i 0.41 + 6.0i )</td></tr><tr><td></td><td>( -9.80 + 0.968i 9.78 - 0.880i )</td><td>( -9.79 + 1.08i 9.83 - 1.01i )</td><td>( 0.10 + 10i 0.99 + 13i )</td></tr></tbody></table>

Note: Matrix of the dimensionless fluid added coefficients. The aspect ratio is $l=15$, the radius ratio is $\epsilon = 1.1$, and the Stokes number is $Sk = 250$.

<table><caption>Table 2 Pinned-pinned cylinders vibrating in the first mode</caption><thead><tr><th>Pinned-pinned</th><th>Present theory</th><th>Numerics TRIOCFD</th><th>Relative deviation (%)</th></tr></thead><tbody><tr><td></td><td>( 5.50 - 0.590i -5.44 + 0.540i )</td><td>( 5.41 - 0.610i -5.44 + 0.560i )</td><td>( 1.7 + 3.3i 0.00 + 5.3i )</td></tr><tr><td></td><td>( -5.44 + 0.540i 5.42 - 0.490i )</td><td>( -5.35 + 0.560i 5.43 - 0.520i )</td><td>( 1.2 + 5.3i 0.20 + 7.7i )</td></tr></tbody></table>

Note: Matrix of the dimensionless modal fluid added coefficients. The aspect ratio is $l=15$, the radius ratio is $\epsilon = 1.1$, and the Stokes number is $Sk = 250$.

<table><caption>Table 3 Clamped-free cylinders vibrating in the first mode</caption><thead><tr><th>Clamped-free</th><th>Present theory</th><th>Numerics TRIOCFD</th><th>Relative deviation (%)</th></tr></thead><tbody><tr><td></td><td>( 2.17 - 0.230i -2.14 + 0.210i )</td><td>( 2.19 - 0.240i -2.15 + 0.220i )</td><td>( 0.90 + 4.2i 0.50 + 4.5i )</td></tr><tr><td></td><td>( -2.14 + 0.210i 2.13 - 0.190i )</td><td>( -2.13 + 0.240i 2.12 - 0.220i )</td><td>( 0.50 + 12i 0.50 + 13i )</td></tr></tbody></table>

Note: Matrix of the dimensionless modal fluid added coefficients. The aspect ratio is $l=15$, the radius ratio is $\epsilon = 1.1$, and the Stokes number is $Sk = 250$.

<figure><img src="image_63.png" /><figcaption>Fig. 10 Modal self and cross-added mass and damping coefficients as functions of the aspect ratio, $l$. The cylinder $C_1$ vibrates in the first mode of a clamped-free beam. The Stokes number is $Sk_1 = 250$. The radius ratio is $\epsilon = 1.1$.</figcaption></figure>
numerical fluid forces onto the same vibration mode as that imposed, i.e., $w_i = \overline{w_i}$. The modal coefficients of the first (resp. second) column of the matrix (51) are extracted from the numerical fluid forces acting on the two cylinders, when the vibration mode is imposed onto the inner (resp. outer) cylinder, the outer (resp. inner) cylinder being fixed.

The results of the comparative study between numerical and theoretical predictions are given in Table 1 for the rigid vibration, in Table 2 for the pinned-pinned vibration, and in Table 3 for the clamped-free vibration. Whatever the type of vibration, we observe an excellent agreement between numerical and theoretical predictions, whether for the added mass or for the added damping coefficients. The relative deviation between theoretical and numerical predictions actually arises from two causes. First, the theoretical approach is fully linear since the convective term of the Navier-Stokes equations is neglected. In our numerical simulations, this term is nonzero even if it is retained as small as possible by imposing a small-amplitude vibration. Second, some relative deviations given in our tables are computed from small values of added coefficients, especially for the damping terms. Tiny differences between numerical and theoretical predictions on such small values automatically yield important relative deviations that should be considered with a critical eye.

## 5.3 Clamped-Free Cylinders Vibrating in the First Mode.

As a second evaluation of the theoretical predictions, we consider the case in which the inner cylinder vibrates in the first mode of a clamped-free beam, considering various values of the aspect ratio $l$ and a low Stokes number, i.e., $Sk_j=250$, representative of a viscous fluid. The numerical simulations are performed setting $L_1 \in \{0.1, 0.14, 0.2, 0.3, 0.4, 0.7, 1.4\}$ [uol] and $L_2 \in \{0.11, 0.15, 0.21, 0.31, 0.41, 0.71, 1.41\}$ [uol] so that the aspect ratio is $l \in \{5, 7, 10, 15, 20, 35, 70\}$. The modal added coefficients are obtained from the projection of the numerical fluid forces onto the same vibration mode as that imposed, i.e., $w_i = \overline{w_i}$.

The results of the comparative study between numerical and theoretical predictions are shown in Fig. 10. There again, an excellent agreement is observed between the two approaches, whether for the added mass or for the added damping coefficients. More specifically, the variations of the coefficients with the aspect ratio are correctly reproduced, thus validating our analytical formulation (51) for the modal added coefficients.

## 6 Conclusion

In the present work, we have considered the fluid-structure interaction problem of two coaxial cylinders, separated by a viscous Newtonian fluid and imposed small-amplitude harmonic beam vibrations. A new theoretical approach based on the Helmholtz expansion of the fluid velocity vector is carried out, leading to a full analytical expression of the modal fluid forces along with the associated mass and damping coefficients. This new theory, which includes the viscous effects of the fluid, does not rely on the assumption of a narrow gap, accounts for the aspect ratio of the cylinders, and applies to all types of classical boundary conditions, thus extending considerably the results of Refs. [21,22]. The new theoretical approach shows that the fluid forces are linear combinations of the Fourier harmonics of the vibration modes. The coefficients of the linear combination depend on the aspect ratio of the cylinders, on the separation distance, and on the Stokes number. The theoretical predictions for the modal self and cross-added mass and damping coefficients are corroborated numerically, considering rigid, pinned-pinned, and clamped-free vibrations. This study shall foster further analytical developments to include the effect of axial fluid flow on the added coefficients and determine subsequent linear stability properties of the coupled system.

## Conflict of Interest

There are no conflicts of interest.

# Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## References

[1] Paidoussis, M. P., 2014, *Fluid-Structure Interactions: Slender Structures and Axial Flow*, Vol. 1, 2nd ed., Elsevier, London.

[2] Paidoussis, M. P., 2016, *Fluid-Structure Interactions. Volume 2: Slender Structures and Axial Flow*, 2nd ed., Academic Press, London, UK.

[3] Brown, S. J., 1982, "A Survey of Studies Into the Hydrodynamic Response of Fluid-Coupled Circular Cylinders," *J. Pressure. Vessel. Technol.*, **104**(1), pp. 2-19.

[4] Magrab, E. B., and Burroughs, C., 1971, "Forced Harmonic and Random Vibrations of Concentric Cylindrical Shells Immersed in Acoustic Fluids," *J. Acoust. Soc. Am.*, **52**(3), pp. 858-864.

[5] Chen, S. S., and Wambsganss, M. W., 1972, "Parallel-Flow-Induced Vibration of Fuel Rods," *Nucl. Eng. Des.*, **188**(2), pp. 253-278.

[6] Levin, L., and Milan, D., 1973, "Coupled Breathing Vibrations of Two Thin Cylindrical Coaxial Shells in Fluid," Vibration Problems in Industry, International Symposium, Keswick, UK.

[7] Krajcinovic, D., 1974, "Vibrations of Two Coaxial Cylindrical Shells Containing Fluid," *Nucl. Eng. Des.*, **30**(2), pp. 242-248.

[8] Bowers, G., and Horvay, G., 1974, "Beam Modes of Vibration of a Thin Cylindrical Shell Flexibly Supported and Immersed in Water Inside of a Coaxial Cylindrical Container of Slightly Larger Radius," *Nucl. Eng. Des.*, **26**(2), pp. 291-298.

[9] Chen, S. S., Wambsganss, M. W., and Jendrzejczyk, J. A., 1976, "Added Mass and Damping of a Vibrating Rod in Confined Viscous Fluids," *ASME J. Appl. Mech.*, **43**, pp. 325-329.

[10] Yeh, T. T., and Chen, S. S., 1978, "The Effect of Fluid Viscosity on Coupled Tube/fluid Vibrations," *J. Sound. Vib.*, **59**(3), pp. 453-467.

[11] Au-Yang, M. K., 1976, "Free Vibration of Fluid-Coupled Coaxial Cylindrical Shells of Different Lengths," *ASME J. Appl. Mech.*, **43**(3), pp. 480-484.

[12] Au-Yang, M. K., 1977, "Generalized Hydrodynamic Mass for Beam Mode Vibration of Cylinders Coupled by Fluid Gap," *ASME J. Appl. Mech.*, **44**(1), pp. 172-173.

[13] Au-Yang, M. K., 1986, "Dynamics of Coupled Fluid-Shells," *J. Vib. Acoust. Stress Reliab. Des.*, **108**(3), pp. 339-347.

[14] Paidoussis, M., and Ostoja-Starzewski, M., 1981, "Dynamics of a Flexible Cylinder in Subsonic Axial Flow," *AIAA. J.*, **19**(11), pp. 1467-1475.

[15] Fritz, R. J., 1972, "The Effect of Liquids on the Dynamic Motions of Immersed Solids," *J. Eng. Ind.*, **94**(1), pp. 167-173.

[16] Horvay, G., and Bowers, G., 1975, "Influence of Entrained Water Mass on the Vibration Modes of a Shell," *J. Fluid. Eng.*, **97**(2), pp. 211-216.

[17] Bergamaschi, Y., Bouilloux, Y., Chantoin, P., Guigon, B., Bravo, X., Germain, C., Rommens, M., and Tremodeux, P., 2002, "Jules Horowitz Reactor, Basic Design," *Proceedings of ENC, San Carlos de Bariloche, Argentina*, Nov. 3-8.

[18] Laurens, M., "Présentation IFS. Bloc Pile RJH, 2020," TA-6515496 Ind. A.

[19] Mateescu, D., and Paidoussis, M. P., 1984, "Annular-Flow-Induced Vibrations of an Axially Variable Body of Revolution in a Duct of Variable Crosssection," *Proceedings of the ASME Symposium on Flow-Induced Vibrations*, ASME Winter Annual Meeting, M. P. Paidoussis and M. K. Au-Yang, eds., New Orleans, LA, ASME, New York, Vol. 4, pp. 53-69.

[20] Paidoussis, M. P., Mateescu, D., and Sim, W.-G., 1990, "Dynamics and Stability of a Flexible Cylinder in a Narrow Coaxial Cylindrical Duct Subjected to Annular Flow," *ASME J. Appl. Mech.*, **57**(1), pp. 232-240.

[21] Lagrange, R., and Puscas, M. A., 2022, "Hydrodynamic Interaction Between Two Flexible Finite Length Coaxial Cylinders: New Theoretical Formulation and Numerical Validation," *ASME J. Appl. Mech.*, **89**(8), p. 081006.

[22] Lagrange, R., Puscas, M. A., Piteau, P., Delaune, X., and Antunes, J., 2022, "Modal Added-Mass Matrix of an Elongated Flexible Cylinder Immersed in a Narrow Annular Fluid, Considering Various Boundary Conditions. New Theoretical Results and Numerical Validation," *J. Fluids Structures*, **114**, p. 103754.

[23] Panunzio, D., Puscas, M.-A., and Lagrange, R., 2022, "FSI-Vibrations of Immersed Cylinders. Simulations With the Engineering Open-Source Code TrioCFD. Test Cases and Experimental Comparisons," *Comptes Rendus. Mécanique*, **350**(G3), pp. 451-476.

[24] Angeli, P. E., Puscas, M. A., Fauchet, G., and Cartalade, A., 2017, "FVCA8 Benchmark for the Stokes and Navier-Stokes Equations with the TrioCFD Code-Benchmark Session," Finite Volumes for Complex Applications VIII—Methods and Theoretical Aspects, Lille, France, June 12-16, pp. 181-202.

[25] Fiorini, C., Després, B., and Puscas, M. A., 2020, "Sensitivity Equation Method for the Navier-Stokes Equations Applied to Uncertainty Propagation," *Int. J. Numerical Methods Fluids*, **93**, pp. 71-92.

[26] Puscas, M. A., Monasse, L., Em, A., Tenaud, C., and Mariotti, C., 2015, "A Conservative Embedded Boundary Method for an Inviscid Compressible Flow Coupled With a Fragmenting Structure," *Int. J. Numerical Methods Eng.*, **103**, pp. 970-995.

[27] Donea, J., Huerta, A., Ponthot, J. P., and Rodriguez-Ferran, A., 2004, "Arbitrary Lagrangian-Eulerian Methods," *Encyclopedia of Computational Mechanics*, E. Stein, R. Borst, and T. J. R. Hughes, eds., American Cancer Society.
D. INTERACTION FLUIDE-STRUCTURE
ENTRE DEUX CYLINDRES
PARALLÈLES

Dans cette annexe, la publication parue dans la revue *Journal of Fluids and Structures*, sous le titre "A new analytical approach for modeling the added mass and hydrodynamic interaction of two cylinders subjected to large motions in a potential stagnant fluid" est reproduite.
# A new analytical approach for modeling the added mass and hydrodynamic interaction of two cylinders subjected to large motions in a potential stagnant fluid

Check for
updates

Romain Lagrange<sup>a,*</sup>, Xavier Delaune<sup>a</sup>, Philippe Piteau<sup>a</sup>, Laurent Borsoi<sup>a</sup>, José Antunes<sup>b</sup>

## ARTICLE INFO

**Article history:**

Received 25 April 2017

Received in revised form 2 October 2017

Accepted 1 December 2017

**Keywords:**

Potential flow

Added mass

Fluid force

Cylinders interaction

## ABSTRACT

A potential theory is presented for the problem of two moving circular cylinders, with possibly different radii, large motions, immersed in an perfect stagnant fluid. We show that the fluid force is the superposition of an added mass term, related to the time variations of the potential, and a quadratic term related to its spatial variations. We provide new simple and exact analytical expressions for the fluid added mass coefficients, in which the effect of the confinement is made explicit. The self-added mass (resp. cross-added mass) is shown to decrease (resp. increase) with the separation distance and increase (resp. decreases) with the radius ratio. We then consider the case in which one cylinder translates along the line joining the centers with a constant speed. We show that the two cylinders are repelled from each other, with a force that diverges to infinity at impact. We extend our approach to the case in which one cylinder is imposed a sinusoidal vibration. We show that the force on the stationary cylinder and the vibration displacement have opposite (resp. identical) axial (resp. transverse) directions. For large vibration amplitudes, this force is strongly altered by the nonlinear effects induced by the spatial variations of the potential. The force on the vibrating cylinder is in phase with the imposed displacement and is mainly driven by the added mass term. The results of this paper are of particular interest for engineers who need to understand the essential features associated with the vibration of a solid body in a still fluid.

© 2017 Elsevier Ltd. All rights reserved.

## 1. Introduction

The interaction between a fluid and moving bodies is a fundamental problem which finds many applications, for example in turbomachinery (Furber and Ffowcswilliams, 1979), fish schooling (Nair and Kanso, 2007), heat exchangers tube banks (Chen, 1975, 1977), vibration of flexible risers (Le Cunff et al., 2002), biomechanics of plants (De Langre, 2008), or energy harvesting (Doare and Michelin, 2011; Singh et al., 2012; Michelin and Doare, 2013; Virot et al., 2016) from flapping flags (Eloy et al., 2008)... The reader should refer to the book of Païdoussis et al. (2014) for a complete bibliography of important works in this field. A classical approach to understand the fluid dynamics surrounding interacting bodies is to consider a potential flow model, in which viscous, rotational vortex and wake effects are disregarded. For large Reynolds numbers and

<footnote><sup>a</sup> Den-Service d'études mécaniques et thermiques (SEMT), CEA, Université Paris-Saclay, F-91191, Gif-sur-Yvette, France</footnote>

<footnote><sup>b</sup> Centro de Ciências e Tecnologias Nucleares, Instituto Superior Técnico, Universidade de Lisboa, Estrada Nacional 10, Km 139.7, 2695-066 Bobadela LRS, Portugal</footnote>

<footnote>* Corresponding author.
E-mail address: romain.lagrange@cea.fr (R. Lagrange).</footnote>
<figure><img src="image_64.png" /><figcaption>Fig. 1. Schematic diagram of the system: two moving circular cylinders $C_i$ with radii $R_i$, centers $O_i$, velocities $\dot{\mathbf{U}}_i(T)$, are immersed in an inviscid fluid $\Omega$. The motions of $C_i$ generate an incompressible and irrotational flow. The instantaneous center-to-center distance is $D$. The bipolar coordinates $\sigma \in [0, 2\pi]$ and $\tau = \ln(L_1/L_2)$ are used to track the position of a fluid particle $P$. The points $(B, F_1, F_2)$ are defined in such a way that points lying on the cylinder boundaries $\partial\Omega_i$ have bipolar coordinates $\tau = \tau_i$.</figcaption></figure>

outside the boundary layers, the potential flow theory is expected to provide a good approximation of the solution, or some of its related characteristics as the added-mass, see Chen (1987). The very first discussions of potential flows around two circular cylinders probably originate from Hicks (1879), who studied the motion of a cylindrical pendulum inside another cylinder filled with fluid. To solve this problem, Hicks considered the potential due to flow singularities distributed over the fluid domain and the cylinder surfaces. The strengths and locations of the singularities were expressed as a set of iterative equations and the potential in an integral form, suitable for numerical computation (Lamb, 1945). Many investigators (Greenhill, 1882; Basset, 1888; Carpenter, 1958; Birkhoff, 1960; Gibert and Sagner, 1980; Landweber and Shahshahan, 1991) then tried to clarify and simplify the method of singularities by Hicks to derive a more tractable expression for the fluid potential. Still, this approximate method remains difficult to apply and as the gap between the cylinders becomes small, a high density of singularities is necessary for the solution to converge. More recently, analytical approaches based on complex analysis and conformal mapping methods (Wang, 2004; Burton et al., 2004; Tchieu et al., 2010) have overcome this problem and general theories to determine the flow through multiple bodies have been proposed (Scolan and Etienne, 2008; Crowdy, 2006, 2010). The two cylinders problem can be solved in this framework, but here we provide a more flexible method which yields new analytical and simplified expressions for the added mass coefficients. We also show that the large motions of the cylinders generate some strong spatial variations of the potential which produce nonlinear inertial effects of the fluid forces.

This paper is organized as follows. Section 2 presents the problem of two parallel circular cylinders immersed in an inviscid fluid and introduces the governing dimensionless numbers. In Section 3, we solve the potential flow problem, provide analytical expressions of the added mass terms and determine the fluid force on the cylinders. In Section 4 we test our predictions against published results and consider the case of a vibrating cylinder. Finally, some conclusions are drawn in Section 5.

## 2. Definition of the problem

Let $C_i$, ($i = 1, 2$), be two circular cylinders with radii $R_i$, immersed in an inviscid fluid of volume mass density $\rho$. We aim to determine the 2D flow generated by the arbitrary motions of $C_i$, (Fig. 1). For convenience we introduce rescaled quantities to reduce the number of parameters of the problem. We use $R_2$ to normalize lengths, some characteristic speed $V_0$ to normalize velocities and $\rho V_0^2$ to normalize pressures. We define

$$r = \frac{R_1}{R_2}, d = \frac{D}{R_2}, \mathbf{u}_i = \frac{\mathbf{U}_i}{R_2}, t = \frac{TV_0}{R_2}, p = \frac{P}{\rho V_0^2}, \mathbf{f}_i = \frac{\mathbf{F}_i}{\rho V_0^2 R_2}, \quad (1)$$

as the radius ratio, the dimensionless center-to-center distance, cylinder displacement, time, fluid pressure and fluid force, respectively. We also introduce the dimensionless separation distance, $\varepsilon = d - (r + 1)$ such that $\varepsilon = 0$ corresponds to cylinders in contact.

## 3. Theoretical model

### 3.1. Fluid equations

It is assumed that the fluid is inviscid and the flow generated by the motion of the circular cylinders is incompressible and irrotational. We thus model the flow with a potential function $\varphi$ which satisfies the Laplace equation

$$\Delta\varphi = 0 \quad \text{in } \Omega, \quad (2)$$
<figure><img src="image_65.png" /><figcaption>Fig. 2. Bipolar coordinate system. The dashed circles are iso-values of σ while the solid circles are iso-values of τ. The two foci points F₁ are defined such that ∂Ωᵢ has a coordinate τᵢ. The points on the x-axis have σ = π if they lie in between F₁ and F₂, and σ = 0 otherwise.</figcaption></figure>

and the inviscid boundary conditions

$$ (\nabla\varphi - \dot{\mathbf{u}}_i) \cdot \mathbf{n}_i = 0 \quad \text{on } \partial\Omega_i. \qquad (3) $$

In the above eqs., Ω is the fluid domain and **n**ᵢ is the outward normal unit vector to ∂Ωᵢ (Fig. 1). To account for an unperturbed fluid at infinity we also require that

$$ \varphi \to 0 \text{ at infinity.} \qquad (4) $$

The fluid pressure p is derived from the unsteady Bernoulli equation

$$ p - p_{\infty} = - \left( \frac{\partial \varphi}{\partial t} + \frac{(\nabla \varphi)^2}{2} \right), \qquad (5) $$

where p_∞ is the pressure at infinity. The fluid force on Cᵢ is obtained by integration of (5) on ∂Ωᵢ

$$ \mathbf{f}_i = \int_{\partial\Omega_i} \left( \frac{\partial\varphi}{\partial t} + \frac{(\nabla\varphi)^2}{2} \right) \mathbf{n}_i ds_i, \qquad (6) $$

where dsᵢ is an infinitesimal element of integration.

The fluid equations must be solved on the instantaneous configuration to account for the possible large displacements of Cᵢ. This leads us to define a coordinate system based on the kinematics of Cᵢ and in which the fluid eqs. can be solved analytically.

## 3.2. Bipolar coordinates

Let (B, F₁, F₂) be three points defined as **OB** = x_B**e**_x and **BF**₁ = -a**e**_x and **BF**₂ = a**e**_x, a > 0. From these three points, we introduce the bipolar coordinates σ ∈ [0, 2π] and τ ∈ ℝ to track the position of a fluid particle, as shown in Fig. 1. Since the iso-τ contours are circles with centers on the x-axis (Fig. 2), one can align ∂Ωᵢ with an iso-contour τᵢ by setting the position of B and Fᵢ. In Appendix A, we show that

$$ x_B = (r^2 - 1) / (2d), \qquad (7a) $$

$$ a = \sqrt{d^2 - (1+r)^2} \sqrt{d^2 - (1-r)^2} / (2d), \qquad (7b) $$

yield iso-contours for the two cylinders

$$ \tau_1 = -\sinh^{-1}(a/r) \quad \text{and} \quad \tau_2 = \sinh^{-1}(a). \qquad (7c) $$

Introducing (e_σ, e_τ) as the coordinate system basis, see Fig. 2, Appendix A, the outward normal unit vector **n**ᵢ to ∂Ωᵢ is

$$ \mathbf{n}_i = (-1)^{i-1} \mathbf{e}_\tau (\sigma, \tau_i), \qquad (8) $$
and the infinitesimal measure $ds_i$ of $\partial\Omega_i$ is

$$ds_i = \kappa_{\sigma\tau_i} d\sigma, \qquad (9)$$

with $\kappa_{\sigma\tau} = a/(\cosh(\tau) - \cos(\sigma))$.

Finally, the gradient and the Laplace operators in bipolar coordinates are

$$\nabla\varphi = \frac{1}{\kappa_{\sigma\tau}} \left( \frac{\partial\varphi}{\partial\sigma} \mathbf{e}_{\sigma} + \frac{\partial\varphi}{\partial\tau} \mathbf{e}_{\tau} \right), \qquad (10a)$$

$$\Delta\varphi = \left(\frac{1}{\kappa_{\sigma\tau}}\right)^2 \left(\frac{\partial^2\varphi}{\partial\sigma^2} + \frac{\partial^2\varphi}{\partial\tau^2}\right). \qquad (10b)$$

### 3.3. Solution of the fluid equations

In bipolar coordinates the fluid Eqs. (2)–(4) read

$$\frac{\partial^2\varphi}{\partial\sigma^2} + \frac{\partial^2\varphi}{\partial\tau^2} = 0 \quad \text{in } \Omega, \qquad (11a)$$

$$\frac{\partial\varphi}{\partial\tau} = g_{ix}\dot{u}_{ix} + g_{iy}\dot{u}_{iy} \quad \text{on } \tau = \tau_i, i = 1, 2, \qquad (11b)$$

$$\varphi \to 0 \quad \text{as } (\sigma, \tau) \to (0, 0), \qquad (11c)$$

with $\dot{u}_{ix}$, $\dot{u}_{iy}$ the $\mathbf{e}_x$ and $\mathbf{e}_y$ components of $\dot{\mathbf{u}}_i$, respectively. In (11), $g_{ix}$ and $g_{iy}$ are the $2\pi$ periodic functions with Fourier components $g_{in} = -2nae^{-n|\tau_i|}$

$$g_{ix}(\sigma) = \kappa_{\sigma\tau_i} \mathbf{e}_x \cdot \mathbf{e}_\tau(\sigma, \tau_i) = -a \frac{\cos(\sigma) \cosh(\tau_i) - 1}{(\cosh(\tau_i) - \cos(\sigma))^2} = \sum_{n=1}^{\infty} g_{in} \cos(n\sigma), \qquad (12a)$$

$$g_{iy}(\sigma) = \kappa_{\sigma\tau_i} \mathbf{e}_y \cdot \mathbf{e}_\tau(\sigma, \tau_i) = -a \frac{\sin(\sigma) \sinh(\tau_i)}{(\cosh(\tau_i) - \cos(\sigma))^2} = \sum_{n=1}^{\infty} g_{in} \sin(\tau_i) \sin(n\sigma). \qquad (12b)$$

The fluid problem being linear in $\dot{u}_{ix}$ and $\dot{u}_{iy}$, the potential $\varphi$ is a superposition of two functions $\varphi_x$ and $\varphi_y$, solutions of (11) for $\dot{u}_{iy} = 0$ and $\dot{u}_{ix} = 0$, respectively. The calculus of $\varphi_x$ and $\varphi_y$ is reported in Appendix B for the sake of clarity. We find that the potential is

$$\begin{aligned} \varphi &= \dot{u}_{1x} \sum_{n=1}^{\infty} \varphi_{1n} [\cos(n\sigma) \cosh(n(\tau - \tau_2)) - \cosh(n\tau_2)] \\ &\quad + \dot{u}_{1y} \sum_{n=1}^{\infty} (-\varphi_{1n}) \sin(n\sigma) \cosh(n(\tau - \tau_2)) \\ &\quad + \dot{u}_{2x} \sum_{n=1}^{\infty} \varphi_{2n} [\cos(n\sigma) \cosh(n(\tau - \tau_1)) - \cosh(n\tau_1)] \\ &\quad + \dot{u}_{2y} \sum_{n=1}^{\infty} \varphi_{2n} \sin(n\sigma) \cosh(n(\tau - \tau_1)) \\ &= \dot{\mathbf{u}}_1 \cdot \mathbf{h}_1 + \dot{\mathbf{u}}_2 \cdot \mathbf{h}_2, \end{aligned} \qquad (13)$$

with $\varphi_{in}$ given by (B.4).

The calculation of the fluid pressure follows from the Bernoulli equation (5), leading to

$$p - p_{\infty} = -(\ddot{\mathbf{u}}_1 \cdot \mathbf{h}_1 + \ddot{\mathbf{u}}_2 \cdot \mathbf{h}_2 + p_q), \qquad (14)$$

with

$$\begin{aligned} p_q &= \dot{\mathbf{u}}_1 \cdot \frac{\partial \mathbf{h}_1}{\partial t} + \dot{\mathbf{u}}_2 \cdot \frac{\partial \mathbf{h}_2}{\partial t} + \frac{[\nabla (\dot{\mathbf{u}}_1 \cdot \mathbf{h}_1 + \dot{\mathbf{u}}_2 \cdot \mathbf{h}_2)]^2}{2}, \\ &= -(\dot{\mathbf{u}}_1 \cdot \nabla \mathbf{h}_1 + \dot{\mathbf{u}}_2 \cdot \nabla \mathbf{h}_2) \cdot \dot{\mathbf{u}}_0 + \frac{[\nabla (\dot{\mathbf{u}}_1 \cdot \mathbf{h}_1 + \dot{\mathbf{u}}_2 \cdot \mathbf{h}_2)]^2}{2}, \end{aligned} \qquad (15)$$

and $\dot{\mathbf{u}}_0 = (\dot{\mathbf{u}}_1 + \dot{\mathbf{u}}_2)/2$. The terms $\ddot{\mathbf{u}}_i \cdot \mathbf{h}_i$ come from the time variation of the fluid potential due to the motion of the cylinders. The term $p_q$ is quadratic in $(\dot{\mathbf{u}}_1, \dot{\mathbf{u}}_2)$ and comes from the spatial variations of the fluid potential. This variation is due to the
motion of the cylinders (first right hand side term of (15)) and to the convective term of the Navier-Stokes equation (last RHS term).

The integration of (14) on $\partial\Omega_i$ yields the fluid force on $C_i$

$$ \begin{pmatrix} f_{1x0} \\ f_{1y0} \\ f_{2x0} \\ f_{2y0} \end{pmatrix} = -\pi [C_a] \begin{pmatrix} \ddot{u}_{1x0} \\ \ddot{u}_{1y0} \\ \ddot{u}_{2x0} \\ \ddot{u}_{2y0} \end{pmatrix} + \begin{pmatrix} f_{1qx0} \\ f_{1qy0} \\ f_{2qx0} \\ f_{2qy0} \end{pmatrix}, \qquad (16) $$

with $[C_a]$ the added mass matrix

$$ [C_a] = \begin{pmatrix} m_{11} & 0 & m_{12} & 0 \\ 0 & m_{11} & 0 & -m_{12} \\ m_{12} & 0 & m_{11} & 0 \\ 0 & -m_{12} & 0 & m_{11} \end{pmatrix}, \qquad (17a) $$

$$ m_{11} = \sum_{n=1}^{\infty} \frac{4na^2 e^{-2n\tau_2}}{\tanh[n(\tau_2 - \tau_1)]}, \qquad (17b) $$

$$ m_{12} = \sum_{n=1}^{\infty} \frac{-4na^2 e^{-n(\tau_2 - \tau_1)}}{\sinh[n(\tau_2 - \tau_1)]}. \qquad (17c) $$

These are new analytical expressions for the self and cross-added mass terms $m_{11}$ and $m_{12}$. The self-added mass $m_{11}$ relates the fluid force on $C_i$ to its own acceleration. The cross-added mass $m_{12}$ relates the fluid force on $C_i$ to the acceleration of $C_j$, $j \neq i$. These new expressions have the advantage of making explicit the effect of the confinement on the added mass terms: $m_{11}$ and $m_{12}$ depend on $a$ and $\tau_i$, which are functions of the separation distance, according to (7).

The terms $\mathbf{f}_{iq}$ are quadratic in $(\dot{\mathbf{u}}_1, \dot{\mathbf{u}}_2)$ and come from the integration of $p_q$ on $\partial\Omega_i$

$$ \mathbf{f}_{iq} = \int_{0}^{2\pi} p_q(\sigma, \tau_i) \mathbf{n}_i \kappa_{\sigma \tau_i} d\sigma. \qquad (18) $$

These quadratic terms (computed numerically, see Appendix C) represent the part of the fluid force due to the spatial variations of the fluid potential.

Finally, note that if the cylinders are subject to an incident unsteady flow $\mathbf{v}_{\infty}(t)$, the fluid potential modifies to

$$ \varphi = \mathbf{v}_{\infty} \cdot \mathbf{OM} + (\dot{\mathbf{u}}_1 - \mathbf{v}_{\infty}) \cdot \mathbf{h}_1 + (\dot{\mathbf{u}}_2 - \mathbf{v}_{\infty}) \cdot \mathbf{h}_2, \qquad (19) $$

where $\mathbf{OM}$ is the position vector of the fluid particle $P$. The first right hand side term ensures that $\nabla\varphi \to \mathbf{v}_{\infty}$ far from the cylinders. The other terms follow from (13) where $\dot{\mathbf{u}}_i - \mathbf{v}_{\infty}$ is the velocity of $C_i$ in the reference frame attached to the incident flow. Noting $[I]$ the identity matrix, the fluid force is

$$ \begin{pmatrix} f_{1x0} \\ f_{1y0} \\ f_{2x0} \\ f_{2y0} \end{pmatrix} = \pi ([C_a] + [I]) \begin{pmatrix} \dot{v}_{\infty x0} \\ \dot{v}_{\infty y0} \\ \dot{v}_{\infty x0} \\ \dot{v}_{\infty y0} \end{pmatrix} - \pi [C_a] \begin{pmatrix} \ddot{u}_{1x0} \\ \ddot{u}_{1y0} \\ \ddot{u}_{2x0} \\ \ddot{u}_{2y0} \end{pmatrix} + \begin{pmatrix} f_{1qx0} \\ f_{1qy0} \\ f_{2qx0} \\ f_{2qy0} \end{pmatrix}, \qquad (20) $$

where $p_q$ and $\mathbf{f}_{iq}$ are given by (15)–(18) under the change $\dot{\mathbf{u}}_i \to \dot{\mathbf{u}}_i - \mathbf{v}_{\infty}$. The Eq. (20) is similar to the Morison equation (Morison et al., 1950), which gives the inline force on a body subject to an unsteady flow. From an experimental stand point, (20) shows that $[C_a]$ can be extracted either from the force on cylinders vibrating in a fluid at rest or from the force on fixed cylinders subject to an unsteady fluid flow.

# 4. Results

## 4.1. Added mass

In the previous section, we have obtained new simple analytical expressions for the added mass coefficients $m_{11}$ and $m_{12}$. Their evolution with the separation distance $\varepsilon$ and the radius ratio $r$ is shown in Fig. 3. We find that $m_{11}$ (resp. $m_{12}$) decreases (resp. increases) monotonically with the dimensionless separation distance, and increases (resp. decreases) with the radius ratio. When the cylinders are in close proximity, i.e. $\varepsilon \to 0$, the confinement is maximum and the added mass coefficients become unbounded, as expected. When the two cylinders are far apart, i.e. $\varepsilon \to \infty$, they both behave like an isolated cylinder in an infinite fluid domain, $m_{11} \to 1$ and $m_{12} \to 0$. To validate our observations, we have reported in Fig. 3 the theoretical predictions of the literature (Mazur, 1970; Landweber and Shahshahan, 1991), for $r = 1$. Unlike the current method, Mazur (1970) used a conformal mapping method to solve the potential problem and extracted the added mass coefficients from
<figure><img src="image_66.png" /><figcaption>Fig. 3. Self-added mass $m_{11}$ (solid lines) and cross-added mass $m_{12}$, given by Eq. (17), as functions of the separation distance $\epsilon$ and the radius ratio $r$.</figcaption></figure>

<caption>**Table 1**</caption>

Added mass coefficients from various references. Note that Gibert's theory (Gibert and Sagner, 1980) is slightly off as it relies on an asymptotic expansion that is valid for large $\epsilon$. The table is given for $\epsilon = 1$ and $r = 1$.

<table><thead><tr><th></th><th>Present reference</th><th>Païdoussis et al. (1984)</th><th>Chen (1978)</th><th>Suss (1977)</th><th>Dalton and Helfinstine (1971)</th><th>Gibert and Sagner (1980)</th></tr></thead><tbody><tr><td>$m_{11}$</td><td>1.0319</td><td>1.0320</td><td>1.0325</td><td>1.0328</td><td>1.0313</td><td>1.0250</td></tr><tr><td>$m_{12}$</td><td>-0.2269</td><td>-0.2269</td><td>-0.2265</td><td>-0.2266</td><td>-0.2245</td><td>-0.2250</td></tr></tbody></table>

the kinetic energy of the fluid. On his side, Landweber and Shahshahan (1991) extended the method of images by Hicks (1879), Herman (1887) and extracted the added mass coefficients from the fluid force acting on the cylinders. We obtain an excellent agreement with those authors (and some others reproduced in Table 1), thereby validating our prediction for $m_{11}$ and $m_{12}$.

## 4.2. Quadratic fluid force

### 4.2.1. Uniform linear motion

The fluid force is the superposition of an added mass term and a quadratic term $\mathbf{f}_{iq}$. In the previous section, we have validated the prediction for the added mass term and we now proceed with assessing the validity of $\mathbf{f}_{iq}$. To do so we consider the case in which $C_1$ is stationary while $C_2$ translates along the line joining the centers with a constant characteristic speed $V_0$. The dimensionless displacement is written $\mathbf{u} = t\mathbf{e}_x$ with $t = TV_0/R_2$ the dimensionless time. The y-component of the fluid force is $f_{iy} = 0$ by symmetry.

In Fig. 4, we plot $f_{ix}$, given by Eq. (16), as a function of the separation distance $\epsilon$. We find that $f_{1x}$ (resp. $f_{2x}$) is negative (resp. positive) so that the cylinders are repelled from each other. The magnitude of this repelling force decreases monotonically with $\epsilon$, and increases with the radius ratio $r$. When the cylinders are near contact, the force $f_{ix}$ becomes singular and diverges to infinity. When the cylinders are far apart the fluid force on $C_i$ vanishes, a result known as the d'Alembert's paradox. In Fig. 4 the results by Bampalas and Graham (2008) and Landweber and Shahshahan (1991) are presented for comparison. In Bampalas and Graham (2008) the force is computed using a Möbius conformal mapping and a series of image singularities. Here again, we find an excellent agreement with those authors thereby validating our expression (18) for $\mathbf{f}_{iq}$.

### 4.2.2. Sinusoidal vibrations

Thus far, we have successfully compared the results of our model to those found in the literature. We now proceed with extending those results to the case in which $C_1$ is stationary while $C_2$ is imposed a sinusoidal vibration with amplitude $U_0$ and angular frequency $\Omega$. We choose $V_0 = \Omega R_2$ as the characteristic speed so that the dimensionless displacement writes $\mathbf{u}(t) = u_0 \sin(t) \mathbf{u}_0$, with $\mathbf{u}_0 = \cos(\gamma) \mathbf{e}_{x_0} + \cos(\gamma) \mathbf{e}_{y_0}$. Here $\gamma = (\mathbf{e}_{x_0}, \mathbf{u}_0)$ is the angle made by the line joining the centers in the initial configuration and the direction of vibration. The initial separation distance is $\epsilon_0 = 1$ and three vibration amplitudes are considered, $u_0 = \{0.1, 0.5, 0.9\}$. The time evolution of the fluid force is shown in Figs. 5 and 6 for $\gamma = 0$ and $\gamma = \pi/4$, respectively. Note that for $\gamma = 0$ the problem is symmetric about the x-axis, such that $f_{iy} = 0$ and we only show the evolution of $f_{ix}$. For $\gamma = \pi/4$ the fluid force is represented through its components ($f_{ix_0}, f_{iy_0}$) in the initial basis ($\mathbf{e}_{x_0}, \mathbf{e}_{y_0}$).
<figure><img src="image_67.png" /><figcaption>**Fig. 4.** Fluid force, computed from Eq. (16), as a function of the separation distance $\epsilon$ and the radius ratio $r$. The cylinder $c_1$ is stationary. The cylinder $c_2$ translates along the line joining the centers with a constant speed.</figcaption></figure>

For both cases ($\gamma = 0, \gamma = \pi/4$) we find that $f_{1x}$ and $u_{2x}$ have opposite signs, see Figs. 5(a) and 6(a). Thus, the fluid always pushes the stationary cylinder in an axial direction that is in phase opposition with the vibration displacement. On the contrary, $f_{1y}$ and $u_{2y}$ have the same sign, see Fig. 6(b), such that the fluid pushes the stationary cylinder in a transverse direction that is in phase with the imposed displacement. We also observe that $f_1$ has a high harmonic distortion due to the non-negligible effect of the quadratic term $f_{1q}$, see Fig. 5(b). It follows that the force on the stationary cylinder is due both to the time and spatial variations of the fluid potential.

The force $f_2$ on the moving cylinder has the same direction as the imposed displacement, see Figs. 5(c) and 6(c–d). The quadratic term $f_{2q}$ remains negligible, see Fig. 5(d), such that $f_2$ has a low harmonic distortion and approximates to $f_2 = -\pi m_{11} \dot{u}_2$. Thus, the force on the moving cylinder is mainly due to the time variations of the fluid potential. Finally, we note that the magnitude $|f_i|$ is maximum as the distance between the cylinders is minimum, and it increases with the radius ratio. This is of course related to a confinement effect which induces an increase of the added mass terms $|m_{11}|$ and $|m_{12}|$, as already pointed out in Section 4.1.

The online movies ($r = 0.5, u0 = 0.9$) show the fluid pressure distribution and its time evolution as the cylinder vibrates. As expected, a positive (resp. negative) pressure appears in the wake of the moving cylinder as it decelerates (resp. accelerates).

## 5. Conclusion

We have considered the problem of two moving circular cylinders immersed in an inviscid fluid. A potential theory based on a bipolar coordinate system shows that the fluid force on the cylinders is the superposition of an added mass term and a quadratic term. The added mass term comes from the time variations of the fluid potential while the quadratic term is related to its spatial variations. We provide new simple and exact analytical expressions for the fluid added mass coefficients, in which the effect of the confinement is made explicit. The self-added mass $m_{11}$ (resp. cross-added mass $m_{12}$) decreases (resp. increases) with the separation distance and increases (resp. decreases) with the radius ratio. When the separation distance tends to zero, the confinement is maximum and the added mass coefficients $m_{11}$ and $m_{12}$ diverge to infinity. When the two cylinders are far apart, the added mass coefficients are those of a single cylinder in a fluid of infinite extent, $m_{11} \to 1$ and $m_{12} \to 0$. We then consider the case in which one cylinder is stationary while the other is translating along the line of centers with a constant speed. Our results indicate that the fluid forces are repulsive, diverge to infinity when the cylinders are in contact, and decrease to zero as they separate. These observations are in excellent agreement with published results. Finally, we consider the case in which one cylinder is fixed while the other is imposed a sinusoidal vibration. We show that the force on the stationary cylinder and the vibration displacement have opposite axial directions but identical transverse directions. This force is strongly altered by the nonlinear effects due to the spatial variations of the fluid potential. On the other hand, the force on the vibrating cylinder is mainly due to the time variations of the potential.

## Acknowledgment

The authors would like to thank Prof. C. Eloy for his valuable comments and suggestions to improve the manuscript.
<figure><img src="image_68.png" /><figcaption>Fig. 5. Fluid force, computed from Eq. (16), as a function of the dimensionless time t. The cylinder c₁ is stationary. The cylinder c₂ vibrates along the line joining the centers with a dimensionless amplitude u₀. The radius ratios are: r = 0.5 (dashed lines), r = 1 (solid lines), r = 2 (dotted lines). The dimensionless vibration amplitudes are u₀ = 0.1 (black color), u₀ = 0.5 (blue color) and u₀ = 0.9 (red color). The initial separation distance is ε₀ = 1. The temporal and spatial evolutions of the fluid pressure and velocity fields for r = 0.5 and u₀ = 0.9 can be observed in the corresponding movies. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

## Appendix A. Relations between coordinates

The relations between the cartesian coordinates (x₀, y₀) (initial configuration), (x, y) (instantaneous configuration) and the bipolar coordinates (σ, τ) are, see Fig. A.1

$$ \begin{pmatrix} x_0 \\ y_0 \end{pmatrix} = \frac{1}{2} \begin{pmatrix} u_{1x_0} + u_{2x_0} \\ u_{1y_0} + u_{2y_0} \end{pmatrix} + \mathbf{P} \begin{pmatrix} x \\ y \end{pmatrix}, \qquad (A.1) $$

$$ \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} x_B + a \sinh(\tau) / (\cosh(\tau) - \cos(\sigma)) \\ a \sin(\sigma) / (\cosh(\tau) - \cos(\sigma)) \end{pmatrix} = \mathbf{g}(\sigma, \tau), \qquad (A.2) $$

$$ \begin{pmatrix} \sigma \\ \tau \end{pmatrix} = \left( \pi - 2\tan^{-1} \left( \frac{2ay}{a^2 - (x - x_B)^2 - y^2 + \sqrt{(a^2 - (x - x_B)^2 - y^2)^2 + 4a^2y^2}} \right) \right), \qquad (A.3) $$
<figure><img src="image_69.png" /></figure>

<figure><img src="image_70.png" /></figure>

<figure><img src="image_71.png" /></figure>

<figure><img src="image_72.png" /><figcaption>**Fig. 6.** Fluid force, computed from Eq. (16), as a function of the dimensionless time *t*. The cylinder *c*₁ is stationary. The cylinder *c*₂ vibrates along the direction ($e_{x_0}$, $u_0$) = π/4 with a dimensionless amplitude *u*₀. Same legend as Fig. 5.</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>**Fig. A.1.** Coordinate systems used for the two cylinders problem. The dashed lines represent the initial configuration whose frame of reference is ($O_0$, $e_{x_0}$, $e_{y_0}$). In this frame a fluid particle *P* has cartesian coordinates ($x_0$, $y_0$). The solid lines represent the instantaneous configuration whose frame of reference is ($O$, $e_x$, $e_y$). In this frame a fluid particle *P* has cartesian coordinates ($x$, $y$) and bipolar coordinates ($\sigma$, $\tau$).</figcaption></figure>
with

$$ \mathbf{P} = \frac{1}{\sqrt{(d_0 + u_{2x_0} - u_{1x_0})^2 + (u_{2y_0} - u_{1y_0})^2}} \begin{pmatrix} d_0 + u_{2x_0} - u_{1x_0} & -(u_{2y_0} - u_{1y_0}) \\ u_{2y_0} - u_{1y_0} & d_0 + u_{2x_0} - u_{1x_0} \end{pmatrix}, \qquad (A.4) $$

the passage matrix from the basis $(\mathbf{e}_{x_0}, \mathbf{e}_{y_0})$ to $(\mathbf{e}_x, \mathbf{e}_y)$.

From (A.2) it follows that

$$ y^2 + [x - (x_B + a \coth(\tau))]^2 = (a/\sinh(\tau))^2, \qquad (A.5) $$

so that iso-$\tau$ contours are circles of radius $a/\sinh(\tau)$ and center at $x = x_B + a \coth(\tau)$. To align $\partial C_i$ with an iso-contour $\tau_i$ we thus require

$$ x_B + a \coth(\tau_1) = -d/2, \qquad (A.6a) $$

$$ a/\sinh(\tau_1) = -r, \qquad (A.6b) $$

$$ x_B + a \coth(\tau_2) = d/2, \qquad (A.6c) $$

$$ a/\sinh(\tau_2) = 1, \qquad (A.6d) $$

which yields the solution (7) for $x_B$, $a$ and $\tau_i$.

The bipolar unit vectors $\mathbf{e}_\sigma$ and $\mathbf{e}_\tau$ shown in Fig. 2 are expressed in the basis $(\mathbf{e}_x, \mathbf{e}_y)$ as

$$ \mathbf{e}_{\sigma} = \frac{1}{|\partial \mathbf{g}/\partial \sigma|} \frac{\partial \mathbf{g}}{\partial \sigma} = \frac{1}{\cosh(\tau) - \cos(\sigma)} \begin{pmatrix} -\sin(\sigma) \sinh(\tau) \\ \cos(\sigma) \cosh(\tau) - 1 \end{pmatrix}, \qquad (A.7a) $$

$$ \mathbf{e}_{\tau} = \frac{1}{|\partial \mathbf{g}/\partial \sigma|} \frac{\partial \mathbf{g}}{\partial \tau} = \frac{1}{\cosh(\tau) - \cos(\sigma)} \begin{pmatrix} 1 - \cos(\sigma) \cosh(\tau) \\ -\sin(\sigma) \sinh(\tau) \end{pmatrix}. \qquad (A.7b) $$

## Appendix B. Determination of the fluid potential

### B.0.3. Determination of $\varphi_x$

The functions $\varphi_x$ and $\varphi_y$ are determined with a separation of variables method with a $2\pi$ periodic function of $\sigma$, Alassar and El-Gebeily (2009).

For $\dot{u}_{iy} = 0$, the fluid Eqs. (11a) and (11c) are automatically satisfied by the potential function

$$ \varphi_x = \dot{u}_{1x} \sum_{n=1}^{\infty} \varphi_{1xn} [\cos(n\sigma) \cosh(n(\tau - \tau_2)) - \cosh(n\tau_2)] \\ + \dot{u}_{2x} \sum_{n=1}^{\infty} \varphi_{2xn} [\cos(n\sigma) \cosh(n(\tau - \tau_1)) - \cosh(n\tau_1)]. \qquad (B.1) $$

The first terms in the brackets are harmonic, $2\pi$ periodic in $\sigma$ and yield a zero fluid velocity at one of the cylinder boundary. The second constant terms are added to ensure that $\varphi_x$ vanishes at infinity, i.e. as $(\sigma, \tau) \to (0, 0)$.

The boundary condition (11b) yields two equations for $\varphi_{ixn}$

$$ \sum_{n=1}^{\infty} \varphi_{1xn} n \cos(n\sigma) \sinh(n(\tau_1 - \tau_2)) = g_{1x}, \qquad (B.2a) $$

$$ \sum_{n=1}^{\infty} \varphi_{2xn} n \cos(n\sigma) \sinh(n(\tau_2 - \tau_1)) = g_{2x}, \qquad (B.2b) $$

whose solutions are obtained by taking the scalar product

$$ \langle f, h \rangle = \frac{1}{\pi} \int_{0}^{2\pi} f(\sigma) h(\sigma) d\sigma, \qquad (B.3) $$

of both sides by $\cos(n\sigma)$, leading to

$$ \varphi_{1xn} = \frac{\langle \cos(n\sigma), g_{1x}(\sigma) \rangle}{n \sinh(n(\tau_1 - \tau_2))} = \frac{-2nae^{n\tau_1}}{n \sinh(n(\tau_1 - \tau_2))} = \varphi_{1n}, \qquad (B.4a) $$

$$ \varphi_{2xn} = \frac{\langle \cos(n\sigma), g_{2x}(\sigma) \rangle}{n \sinh(n(\tau_2 - \tau_1))} = \frac{-2nae^{-n\tau_2}}{n \sinh(n(\tau_2 - \tau_1))} = \varphi_{2n}. \qquad (B.4b) $$
### B.0.4. Determination of $\varphi_y$

For $\dot{u}_{ix} = 0$, the fluid Eqs. (11a) and (11c) are automatically satisfied by the potential function

$$ \varphi_y = \dot{u}_{1y} \sum_{n=1}^{\infty} \varphi_{1yn} \sin(n\sigma) \cosh(n(\tau - \tau_2)) \\ + \dot{u}_{2y} \sum_{n=1}^{\infty} \varphi_{2yn} \sin(n\sigma) \cosh(n(\tau - \tau_1)). \qquad (B.5) $$

The boundary condition (11b) yields two equations for $\varphi_{1yn}$

$$ \sum_{n=1}^{\infty} \varphi_{1yn} n \sin(n\sigma) \sinh(n(\tau_1 - \tau_2)) = g_{1y}, \qquad (B.6a) $$

$$ \sum_{n=1}^{\infty} \varphi_{2yn} n \sin(n\sigma) \sinh(n(\tau_2 - \tau_1)) = g_{2y}, \qquad (B.6b) $$

whose solution are obtained by taking the dot product of both sides by $\sin(n\sigma)$

$$ \varphi_{1yn} = \frac{\langle \sin(n\sigma), g_{1y}(\sigma) \rangle}{n \sinh(n(\tau_1 - \tau_2))} = \frac{2nae^{n\tau_1}}{n \sinh(n(\tau_1 - \tau_2))} = -\varphi_{1n}, \qquad (B.7a) $$

$$ \varphi_{2yn} = \frac{\langle \sin(n\sigma), g_{2y}(\sigma) \rangle}{n \sinh(n(\tau_2 - \tau_1))} = \frac{-2nae^{-n\tau_2}}{n \sinh(n(\tau_2 - \tau_1))} = \varphi_{2n}. \qquad (B.7b) $$

## Appendix C. Strategy for coding

In this Appendix, we give the relevant equations and a strategy for coding the fluid force (16) on the cylinders

$$ \begin{pmatrix} f_{1x0} \\ f_{1y0} \\ f_{2x0} \\ f_{2y0} \end{pmatrix} = -\pi [C_a] \begin{pmatrix} \ddot{u}_{1x0} \\ \ddot{u}_{1y0} \\ \ddot{u}_{2x0} \\ \ddot{u}_{2y0} \end{pmatrix} + \begin{pmatrix} f_{1qx0} \\ f_{1qy0} \\ f_{2qx0} \\ f_{2qy0} \end{pmatrix}. \qquad (C.1) $$

The very first step is to define the initial configuration by setting some values to the radius ratio $r$ and the center-to-center distance $d_0$. Note that $d_0$ must be greater than $r + 1$ to avoid contact between the cylinders. Then, the cylinders are imposed a particular kinematics by setting some values to the $x_0$ and $y_0$ components of the displacement vectors. For example, in the case of Section 4.2.2, the cylinder $C_1$ is stationary, i.e. $u_{1x0} = u_{1y0} = 0$, while $C_2$ is imposed a sinusoidal vibration, i.e. $u_{2x0} = \sin(t)$, $u_{2y0} = 0$, with $t$ a free parameter. The velocity and the acceleration of the cylinders follow from the derivatives of $u_{1x0}$ and $u_{1y0}$. From the initial configuration and the kinematics of the cylinders, we compute the instantaneous center-to-center distance $d$

$$ d = \sqrt{(d_0 + u_{2x0} - u_{1x0})^2 + (u_{2y0} - u_{1y0})^2}, \qquad (C.2) $$

and the position $a$ of the foci points of the bipolar coordinate system

$$ a = \sqrt{d^2 - (1+r)^2} \sqrt{d^2 - (1-r)^2} / (2d). \qquad (C.3) $$

Having determined $a$, the instantaneous bipolar coordinates of the cylinders are given by (7c)

$$ \tau_1 = -\sinh^{-1}(a/r) \quad \text{and} \quad \tau_2 = \sinh^{-1}(a). \qquad (C.4) $$

The coefficients $m_{11}$ and $m_{12}$ of the added mass matrix $[C_a]$ follow from the analytical expressions (17)

$$ m_{11} = \sum_{n=1}^{\infty} \frac{4na^2 e^{-2n\tau_2}}{\tanh[n(\tau_2 - \tau_1)]}, \qquad (C.5a) $$

$$ m_{12} = \sum_{n=1}^{\infty} \frac{-4na^2 e^{-n(\tau_2 - \tau_1)}}{\sinh[n(\tau_2 - \tau_1)]}. \qquad (C.5b) $$
At this point, the added-mass term of the fluid force is determined and we are left with the computation of the quadratic term $f_{iq}$. This term is given by (18) and its components are

$$f_{iqx} = (-1)^{i-1} a \int_{0}^{2\pi} p_q \frac{1 - \cos(\sigma) \cosh(\tau_i)}{[\cosh(\tau_i) - \cos(\sigma)]^2} d\sigma, \quad i = 1, 2, \qquad (C.6a)$$

$$f_{iqy} = (-1)^{i-1} a \int_{0}^{2\pi} p_q \frac{-\sin(\sigma) \sinh(\tau_i)}{[\cosh(\tau_i) - \cos(\sigma)]^2} d\sigma, \quad i = 1, 2, \qquad (C.6b)$$

with, see (15)

$$p_q = -\frac{1}{2} [A] \begin{pmatrix} \dot{u}_{1x} + \dot{u}_{2x} \\ \dot{u}_{1y} + \dot{u}_{2y} \end{pmatrix} + \frac{1}{2} [A]^T [A], \qquad (C.7)$$

$$[A] = (\dot{u}_{1x} \quad \dot{u}_{1y}) \begin{pmatrix} \frac{\partial h_{1x}}{\partial x} & \frac{\partial h_{1x}}{\partial y} \\ \frac{\partial h_{1y}}{\partial x} & \frac{\partial h_{1y}}{\partial y} \end{pmatrix} + (\dot{u}_{2x} \quad \dot{u}_{2y}) \begin{pmatrix} \frac{\partial h_{2x}}{\partial x} & \frac{\partial h_{2x}}{\partial y} \\ \frac{\partial h_{2y}}{\partial x} & \frac{\partial h_{2y}}{\partial y} \end{pmatrix}, \qquad (C.8)$$

and, see (13)

$$\mathbf{h}_1 = \sum_{n=1}^{\infty} \varphi_{1n} [\cos(n\sigma) \cosh(n(\tau - \tau_2)) - \cosh(n\tau_2)] \mathbf{e}_x \\ + \sum_{n=1}^{\infty} (-\varphi_{1n}) \sin(n\sigma) \cosh(n(\tau - \tau_2)) \mathbf{e}_y, \qquad (C.9a)$$

$$\mathbf{h}_2 = \sum_{n=1}^{\infty} \varphi_{2n} [\cos(n\sigma) \cosh(n(\tau - \tau_1)) - \cosh(n\tau_1)] \mathbf{e}_x \\ + \sum_{n=1}^{\infty} \varphi_{2n} \sin(n\sigma) \cosh(n(\tau - \tau_1)) \mathbf{e}_y. \qquad (C.9b)$$

A chain rule between the bipolar coordinates $(\sigma, \tau)$ and the cartesian coordinates $(x, y)$, see (A.3), is used to compute the partial derivatives in $[A]$. By plugging (C.7) into (C.6), the integrals $f_{iqx}$ and $f_{iqy}$ are computed numerically and the components $f_{iqx_0}$ and $f_{iqy_0}$ are obtained from the relation

$$\begin{pmatrix} f_{iqx_0} \\ f_{iqy_0} \end{pmatrix} = \mathbf{P} \begin{pmatrix} f_{iqx} \\ f_{iqy} \end{pmatrix}, \qquad (C.10)$$

with $\mathbf{P}$ the passage matrix given by (A.4).

Finally, the sum of the added mass and quadratic terms yields the fluid force on the cylinders.

## Appendix D. Supplementary data

Supplementary material related to this article can be found online at https://doi.org/10.1016/j.jfluidstructs.2017.12.002.

## References

Alassar, R.S., El-Gebeily, M.A., 2009. Inviscid flow past two cylinders. ASME Trans. J. Fluid Eng. 131, 054501.

Bampalas, N., Graham, J.M.R., 2008. Flow-induced forces arising during the impact of two circular cylinders. J. Fluid Mech. 616, 205–234.

Basset, A.B., 1888. A Treatise on Hydrodynamics. Deighton, Bell and Co.

Birkhoff, G., 1960. Hydrodynamics, second ed. Princeton University Press, Princeton, New Jersey.

Burton, D.A., Gratus, J., Tucker, R.W., 2004. Hydrodynamic forces on two moving discs. Theoret. Appl. Mech. 31, 153–188.

Carpenter, L.H., 1958. On the motion of two cylinders in an ideal fluid. J. Res. Nat. Bur. Stand. 61, 83–87.

Chen, S.S., 1975. Vibration of nuclear fuel bundles. Nucl. Eng. Des. 35, 399–422.

Chen, S.S., 1977. Dynamics of heat exchanger tube banks. J. Fluid. Eng. 99, 462–469.

Chen, S.S., 1978. Crossflow-induced vibrations of heat exchanger tube banks. Nucl. Eng. Des. 47, 67–86.

Chen, S.S., 1987. Flow-induced Vibration of Circular Cylindrical Structures. Hemisphere Publishing.

Crowdy, D.G., 2006. Analytical solutions for uniform potential flow past multiple cylinders. Eur. J. Mech. B Fluids 25, 459–470.

Crowdy, D.G., 2010. A new calculus for two-dimensional vortex dynamics. Theor. Comput. Fluid Dyn. 24, 9.

Dalton, C., Helfinstine, R.A., 1971. Potential flow past a group of circular cylinders. Trans. ASME J. Basic Eng. 93, 636–642.

De Langre, E., 2008. Effects of wind on plants. Annu. Rev. Fluid Mech. 40, 141–168.

Doare, O., Michelin, S., 2011. Piezoelectric coupling in energy-harvesting fluttering flexible plates: Linear stability analysis and conversion efficiency. J. Fluids Struct. 27, 1357–1375.

Eloy, C., Lagrange, R., Souilliez, C., Schouveiler, L., 2008. Aeroelastic instability of cantilevered flexible plates in uniform flow. J. Fluid Mech. 611, 97–106.
Furber, S.B., Frowcswilliams, J.E., 1979. Is the Weis-Fogh principle exploitable in turbomachinery? J. Fluid Mech. 94, 519–540.

Gibert, R.J., Sagner, M., 1980. Vibration of structures in a static fluid medium. La Houille Blanche 1/2, 204–262.

Greenhill, A.G., 1882. Functional images in cartesianes. Quaterly J. Pure Appl. Math. 18, 231.

Herman, R.A., 1887. On the motion of two spheres in fluid and allied problem. Quaterly J. Pure Appl. Math. 22, 204–262.

Hicks, W.M., 1879. On the motion of two cylinders in a fluid. Quaterly J. Pure Appl. Math. 16, 113–140, 193–219.

Lamb, H., 1945. Hydrodynamics, sixth ed. Dover, New York.

Landweber, L., Shahshahan, A., 1991. Added masses and forces on two bodies approaching central impact in an inviscid fluid. In: Technical Rep. 346. Iowa Institute of Hydraulic Research.

Le Cunff, C., Biolley, F., Fontaine, E., Etienne, S., Facchinetti, M.L., 2002. Vortex-induced vibrations of risers: Theoretical, numerical and experimental investigation. Oil Gas Sci. Technol. 57, 59–69.

Mazur, V.Y., 1970. Motion of two circular cylinders in an ideal fluid. Izvestiya Akademii Nauk SSSR, Mekhanika Zhidkosti I Gaza 6, 80–84.

Michelin, S., Doare, O., 2013. Energy harvesting efficiency of piezoelectric flags in axial flows. J. Fluid Mech. 714, 489–504.

Morison, J.R., O'Brien, M.P., Johnson, J.W., Schaaf, S.A., 1950. The force exerted by surface waves on piles. Petroleum Trans. Am. Inst. Min. Eng. 189, 149–154.

Nair, S., Kanso, E., 2007. Hydrodynamically coupled rigid bodies. J. Fluid Mech. 592, 393–411.

Païdoussis, M.P., Mavriplis, D., Price, S.J., 1984. A potential-flow theory for the dynamics of cylinder arrays in cross-flow. J. Fluid Mech. 146, 227–252.

Païdoussis, M.P., Price, S.J., De Langre, E., 2014. Fluid-Structure Interactions: Cross-Flow-Induced Instabilities. Cambridge University Press.

Scolan, Y.M., Etienne, S., 2008. On the use of conformal mapping for the computation of hydrodynamic forces acting on bodies of arbitrary shape in viscous flow. Part 2: multi-body configuration. J. Eng. Math. 61, 17–34.

Singh, K., Michelin, S., De Langre, E., 2012. Energy harvesting from axial fluid-elastic instabilities of a cylinder. J. Fluids Struct. 30, 159–172.

Suss, S., 1977. Stability of a Cluster of Flexible Cylinders in Bounded Axial Flow, Ph.D. Thesis. McGill.

Tchieu, A.A., Crowdy, D., Leonard, A., 2010. Fluid-structure interaction of two bodies in an inviscid fluid. Phys. Fluids 22, 107101.

Virot, E., Amandolese, X., Hemon, P., 2016. Coupling between a flag and a spring-mass oscillator. J. Fluids Struct. 65, 447–454.

Wang, Q.X., 2004. Interaction of two circular cylinders in inviscid fluid. Phys. Fluids 16, 4412.