# 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

Yves-Marie Scolan, Stéphane Etienne

► To cite this version:

Yves-Marie Scolan, Stéphane Etienne. 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. Journal of Engineering Mathematics, 2008, 61 (1), pp.17-34. 10.1007/s10665-007-9159-1. hal-00442681

HAL Id: hal-00442681

https://hal.science/hal-00442681v1

Submitted on 17 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# 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

Y.-M. Scolan

Ecole Centrale de Marseille, Technopôle Château Gombert, 13451, Marseille Cedex 20, France e-mail: ymscolan@ec-marseille.fr

S. Etienne

Département de Génie Mécanique, Ecole Polytechnique de Montréal, 2500, chemin Polytechnique, Montréal, QC, Canada H3T1J4
e-mail: stephane.etienne@polymtl.ca

**Abstract** Two-dimensional viscous flows around obstacles are considered in an unbounded liquid. The basic idea developed in Part 1 is further extended from a single body to multi-body configurations. This idea follows from the formulation by Quartapelle and Napolitano (AIAA J. 21:991–913, 1983) who proposed computation of the force and moment in incompressible viscous flow without explicitly calculating the pressure. The principle is the projection of Navier–Stokes equations on a set of functions. Surprisingly, these functions have a precise meaning in potential theory. They are the solutions which lead to the added masses and added moment of inertia for the potential flow around the studied arrangement of obstacles. By revisiting this problem, a general identity of the full coupled matrix of added masses and added moment of inertia is formulated. To this end conformal mappings for multi-body configurations are used. Robustness of the proposed algorithms is tested and illustrated. The obtained potential solution is merely a mathematical solution and it does not allow to describe the actual potential flow since the circulation is not accounted for. However, its interest is crucial for implementing the projection technique developed by Quartapelle and Napolitano. The interest of such a method is two-fold. Firstly, it provides a way of computing the force without explicitly calculating the pressure. Consequently and secondly, it offers an alternate way to validate the computation of the loads. In effect, these loads are always available from a direct integration of the Cauchy stress tensor (pressure plus friction). It is worth mentioning that the present technique allows an a posteriori computation of the pressure.

**Keywords** Conformal mappings · Multi-body configuration · Two-dimensional flows
# 1 Introduction

On the basis of the theoretical formulation developed by Quartapelle and Napolitano [1], Part I demonstrated how potential theory for two-dimensional flows provides useful mathematical solutions for the computation of the hydrodynamic viscous loads. Simply connected bodies were considered in Part I.

Part II shows the extension of this technique to multi-body configurations of arbitrary shape. Conformal transformations for $N$ bodies are implemented in order to provide a convenient way to describe the computational domain. This domain is suitable for the simulation of Navier-Stokes equations by using optimized algorithms, basically a Lagrangian transportation of the properties of the fluid: vorticity and, possibly, the turbulent properties as well, as in algebraic turbulent models.

As a first step, the boundary-value problems corresponding to the added masses and added moment of inertia calculations are formulated and then solved via a multipole technique. The corresponding velocity potential is precisely required in the approach proposed in [1]. By revisiting this potential problem, we provide a new formulation (to our knowledge) of the hydrodynamic-load coefficients in potential theory for multi-body configurations of arbitrary shape. Illustrative cases are computed for different shapes such as an array of square cylinders. Robustness and convergence of the algorithms are tested as well. It should be noted that the obtained potential solution is not physical in the sense that no circulation is accounted for. However, it is demonstrated here how useful this potential solution is as a basis of functions on which the Navier-Stokes equations can be projected in order to get the hydrodynamic loads in incompressible viscous flows.

The main feature of Quartapelle and Napolitano's approach is that the pressure is not explicitly required to calculate the forces and moment. However, it is shown how the pressure can be reconstructed a posteriori.

The Navier-Stokes solver for a multi-body configuration has not yet been extended to shapes other than circular sections. In future work, we plan to study configurations of more practical interest by incorporating the present formulations into the computational code developed by Etienne [2].

# 2 A mathematical solution in potential theory

We first focus on the flow around $N$ bodies in potential theory. This preliminary step is essential for the sequel but it should be noted that the corresponding potential flow is unphysical. The reason is that no circulation is introduced. The computation of more physical flows for $N$ bodies in potential theory is fully described in [3].

In the most general case (with no circulation), where $N$ bodies move along their three possible degrees of freedom in a fluid at rest, the boundary-value problem to be solved can be written as follows

$$ \left\{ \begin{array}{ll} \Delta\phi = 0, & \text{in the fluid domain } \Omega, \\ \vec{\nabla}\phi \cdot \vec{n} = \vec{n} \cdot \vec{V}(P) & \text{on each body } (B_k), \\ \vec{V}\phi \to 0, & \text{at infinity,} \end{array} \right. \qquad (1) $$

where $\phi$ is a perturbed velocity potential which is of mathematical interest mainly in the present study. Figure 1 illustrates the original configuration and the notations used are indicated as well. $\vec{V}(P)$ is the velocity of a point $P$ belonging to the contour $B_k$. The kinematics of each body is determined in a fixed coordinate system. This kinematics take into account

- a translational motion $\vec{W}_k$; this is the velocity of a point $O_k$ attached to the body $B_k$

- a rotational motion $\vec{\Omega}_k = \Omega_k \vec{z}$ directed in the normal direction $\vec{z}$ of the plane flow about the point $O_k$ (which hence moves with the body).

The Neumann condition at point $P$ is written

$$ \vec{n} \cdot \vec{V}(P) = \vec{n} \cdot \vec{W}_k + \vec{n} \cdot (\vec{\Omega}_k \wedge O_k P) = \vec{n} \cdot \vec{W}_k + \vec{\Omega}_k \cdot (O_k P \wedge \vec{n}). \qquad (2) $$

Due to the linearity (with respect to the data of the kinematics) of the boundary conditions on each body, we can decompose $\phi$ into a sum of several contributions, each describing a given flow. The superposition of all of these
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Physical configuration and notations. Current point on the body contour: *P*, translational velocity of body *k*: *W*<sub>*k*</sub>, angular velocity: *Ω*<sub>*k*</sub>, center of rotation: *O*<sub>*k*</sub></figcaption></figure>

provides the general flow in (1). The boundary-value problem (1) is hence decomposed as

$$ 
\left\{ 
\begin{array}{ll} 
\Delta\phi = 0, & \text{in } \Omega, \\ 
\vec{\nabla}\phi \cdot \vec{n} = \vec{n} \cdot \vec{V}(P) & \text{on } (B_k), \\ 
\vec{\nabla}\phi \to 0, & \text{at } \infty 
\end{array} 
\right\} 
\equiv \sum_{k=1}^{N} \sum_{\ell=1}^{3} 
\left\{ 
\begin{array}{ll} 
\Delta\phi_{k\ell} = 0, & \text{in } \Omega, \\ 
\vec{\nabla}\phi_{k\ell} \cdot \vec{n} = n_{k\ell} & \text{on } (B_k), \\ 
\vec{\nabla}\phi_{k\ell} \cdot \vec{n} = 0 & \text{sur } (B_j)_{j \neq k}, \\ 
\vec{\nabla}\phi_{k\ell} \to 0, & \text{at } \infty, 
\end{array} 
\right. 
\quad (3) 
$$

where the Neumann condition on $B_k$ is defined as

$$ n_{k\ell} = \begin{cases} \vec{n} \cdot \vec{x} & \ell = 1, \\ \vec{n} \cdot \vec{y} & \ell = 2, \\ \vec{n} \cdot (\vec{z} \wedge \vec{r}) & \ell = 3. \end{cases} \qquad (4) $$

By adding the elementary solutions, we obtain the total velocity potential of the flow,

$$ \phi = \sum_{k=1}^{N} \sum_{\ell=1}^{3} U_{k\ell} \phi_{k\ell}, \qquad (5) $$

where $U_{k\ell}$ denotes the generalized amplitudes of the rotation–translation of body $k$ along its degree of freedom $\ell$.

### 3 Method of solution for an elementary problem

The method of solution for an elementary problem consists in searching for an analytical complex function $F$ which describes the two-dimensional flow around several bodies, with one of them, the body $k$, moving along its degree of freedom $\ell$ and the other bodies being at rest. For the sake of simplicity we omit the indices $k$ and $\ell$ to formulate the corresponding boundary-value problem which reads

$$ 
\left\{ 
\begin{array}{ll} 
\Delta\phi = 0, & \text{in the fluid domain } \Omega, \\ 
\vec{\nabla}\phi \cdot \vec{n} = \begin{cases} \vec{n} \cdot \vec{x}, \\ \vec{n} \cdot \vec{y}, \\ \vec{n} \cdot (\vec{z} \wedge \vec{r}) \end{cases}, & \text{on the body } (B_k), \\ 
\vec{\nabla}\phi \cdot \vec{n} = 0 & \text{on the bodies } (B_j)_{j \neq k}, \\ 
\vec{\nabla}\phi \to 0, & \text{at infinity.} 
\end{array} 
\right. 
\qquad (6) 
$$
On the basis of the theoretical developments in [4, pp. 328–329] it is often advantageous to formulate the same problem for the stream function $\psi$. The main advantage is to get Dirichlet conditions on the body surface.

$$ egin{cases} 
\Delta\psi = 0, & \text{in the fluid domain } \Omega, \\
\psi = \begin{cases} -y, & \text{on the body } (B_k), \\ \frac{r^2}{2}, & \text{on the bodies } (B_j)_{j \neq k}, \end{cases} \\
\nabla\psi \to 0, & \text{at infinity.}
\end{cases} 
\quad (7) $$

The additive constant in $\psi$ on body $B_k$ is not required, it can be introduced in the constants $q_j$ attached to the other bodies.

The complex function $F$ combines the solutions of both problems (6) and (7), say $F(z) = \phi + i\psi$.

The right-hand sides of the Dirichlet conditions contain geometric information of a given body. However, there is a complex function $z = f(\zeta)$ which connects the complex coordinate $z$ in the physical plane to the complex coordinate $\zeta$ defined in a transformed plane where the contours of bodies $B_k$ are all circular with radius $\rho_k$ and centered at $\zeta = d_k$. This means that there is a parametrization of $z$ describing the contour $B_k$ with an argument $\alpha$ ranging in the interval $[-\pi : \pi]$.

The function $F$ denotes the complex potential of the flow in both planes (physical and transformed) as soon as we can provide a conformal transformation $z = f(\zeta)$ linking one point $z$ to the point $\zeta$ as described above.

The main problem is to prove that $f$ is effectively conformal. The theorem (see [5], p. 110) stipulates that $f$ is conformal if and only if (i) $f$ is a one-to-one function, (ii) $f$ is analytical and (iii) $f'(\zeta) \neq 0$ everywhere in the domain, that is to say, at all points not interior to the contours $B_k$. We can easily assert that $f$ will be analytical since it will be constructed for that purpose: all the poles are inside $B_k$.

Condition (3) is more difficult to prove as a constraint on function $f$ which could be approximately suitable as explained in Venkatesan [6]. Hence $f$ must be constructed without having to check the conformality a posteriori.

The classical way to proceed is to elaborate an iterative procedure by using a succession of conformal mappings. At convergence circular shapes are obtained. The main problem is to invert the used conformal mapping at each step. In practice we do not expect any difficulties, since only Laurent series are used. Then a Newton algorithm suffices.

As a result, a continuous function is obtained which provides the coordinates of a point on a body contour in the physical plane as a function of an azimuthal coordinate on a circle in the transformed plane.

The parametrization is then well defined in order to calculate the right-hand sides in the Dirichlet conditions. As a consequence the complex potential $F$ is deduced. As a first result we can proceed to the computation of the added masses and added moment of inertia. Then the solutions are introduced in Quartapelle and Napolitano’s formulation.

## 4 Computation of the conformal mapping $f$

The extension of conformal mappings for a single body to $N$ bodies has received much attention in the past. Lagally [7] already showed how to perform the transformation around two cylinders into an annulus. This offered an easier way of determining the corresponding potential flow. Ives [8], Halsey [9] and Suddhoo [10] used the same approach to treat non-circular shapes like airfoils. The book by Prosnak [11] provides a large range of applications as well. Scolan and Faltinsen [12] introduced these transformations to compute viscous flows around two circular cylinders. On the basis of the theoretical developments by Basset [13], Venkatesan elaborated the multipole approach to compute potential flow around two bodies of arbitrary shape. Venkatesan underlined that conformality of the transformation must be checked. For more than two bodies Wegmann [14] and Luchini and Manzo [15] proposed some improvements of the transformations. As much work has already been done in the past, the main outlines of the technique is described in the sequel. However, several applications illustrate the elaborated technique.
For N bodies, a straightforward technique is, in fact, well known. It consists in an iterative procedure involving a succession of conformal mappings like the transformation by Theodorsen and Garrick [16]. First a specific transformation is applied to each successively considered shape. For example, an ellipse is transformed by using the Laurent series

$$z_0 = d_j + \frac{1}{2} e^{i\gamma_j} \left( (a_j + b_j)z_1 + \frac{a_j - b_j}{z_1} \right), \qquad (8)$$

where $(z_0, z_1)$ are the complex coordinates of a point in the initial plane and the transformed plane, respectively. The parameters $(a_j, b_j, \gamma_j, d_j)$ are the semi-axes, orientation and position in the physical plane of the jth ellipse. This mapping function transforms the ellipse into a unit circle. This transformation is inverted to get the image point $z_1$ of $z_0$ belonging to other shapes.

For flat plates, degenerated transformations like (8) with a zero semi-axis are used. However, difficulties appear as soon as the first plate has been transformed since the image of the other plates become curved cuts in the plane. When the curved cuts are simple circular arcs, we find the mapping function in [5, p. 149, Art. 34]. Applications can then be found in [17]. More generally, the curved cuts are not necessarily circular arcs. However, it can be shown that the same transformation is well suited,

$$z_0 = \frac{1}{2} \left( z_1 + \frac{a^2}{z_1} \right), \qquad (9)$$

provided that we center the cut at the origin of the coordinate system and we set the tips of the cut along the real axis. Special attention must be paid to the determination of the argument of complex variables, i.e., the location of the branch cut.

For rectangles, it is getting more complicated but the transformations are also well known. The so-called Karmann-Trefftz transformations remove the four corners of each shape. The mapping function reads

$$z_1 = s \frac{1 + A^{1/\beta}}{1 - A^{1/\beta}}, \quad \beta = 2 - \frac{\alpha}{\pi}, \quad A = \frac{z_0 - Z_t}{z_0 - Z_c}, \quad s = \frac{Z_t - Z_c}{2\beta}, \qquad (10)$$

where $\alpha$ is the inner angle at a given corner. The complex variables $(Z_t, Z_c)$ are complex coordinates of the corner and the geometrical center of the shape, respectively. This transformation is otherwise easily inverted as

$$z_0 = \frac{Z_t - Z_c B^\beta}{1 - B}, \quad B = \frac{z_1 - s}{z_1 + s}. \qquad (11)$$

Other (and more irregular) shapes can be considered, such as circles with fins (see [18] or [19]), for which specific transformations lead to circles.

More generally, the first step must end up with shapes that are almost circular. Then the Theodorsen-Garrick transformation is used within an iterative procedure. The corresponding mapping function is

$$z_0 = z_b + z_1 \cdot \exp \left[ \sum_{n=0}^{\infty} \beta_n z_1^{-n} \right], \qquad (12)$$

where $z_b$ is the geometrical center of the considered body. The complex coefficients $\beta_n$ are computed via a fixed-point algorithm. In practice an accuracy of $10^{-12}$ is imposed. At convergence, perfect circles are obtained. The more irregular the initial physical shapes, the longer the iterative scheme will take. Usually the conditions established by Warchawski [20] inform us on the speed of the scheme.

We illustrate the procedure for two rectangles. Figure 2 shows the original shapes, the transformed shapes by using the Karmann-Trefftz transformations for each corner of each shape and the reconstructed shapes from the transformed plane where the body contours are circular. Here seven iterations are required to converge towards circles with a machine accuracy of $10^{-13}$.

The procedure is also illustrated for four flat plates. Figure 3 shows the successive transformations between the physical plane to almost circular shapes, before the Theodorsen-Garrick transformations start. A specific transformation is elaborated to transform a curved cut into a circle. Actually the (curiously) simple transformation (9) is
<figure><img src="image_3.png" /><figcaption>Fig. 2 (a) Original shapes. (b) Transformed shapes by using the Karmann–Trefftz transformations for each corner of each shape. (c) Obtained circular contours and reconstructed shapes from mapping inverse transformation</figcaption></figure>

suitable whether the curvature of the cut changes sign or not along the cut. To our knowledge, this is not mentioned in textbooks devoted to conformal mapping such as those of Kober [21] or Nehari [22].

Another application is considered illustrating a large number of bodies. Here 25 identical squares are studied. Figure 4b shows the shapes obtained before the Theodorsen–Garrick transformations start. Then the obtained array of circular sections is plotted (see Fig. 4a). The rate of convergence of the successive Theodorsen–Garrick transformations is also plotted (see Fig. 4c). At convergence 220 iterations have been carried out. Note that the slight change of orientation of the array is due to the successive combination of the conformal transformations. However, as expected, the relative positions of each cylinder in the array are not affected by the transformations.

The conformal mapping function provides a parametrization of each contour by means of an argument within the range $[-\pi : \pi]$. Due to the conformality of the transformation $f$, we make sure that the complex potential $F$ in both planes $z$ and $\zeta$ is identical, provided that $z = f(\zeta)$. In the next section, we show how to calculate $F$.

## 5 Computation of the complex potential $F$

The solution of the boundary-value problem (6) is sought. In practice, the general expression of the corresponding complex potential $F = \phi + i\psi$ is written as

$$F(\zeta) = \sum_{j=1}^{N} \sum_{n=1}^{\infty} a_{jn} \left( \frac{\rho_j}{\zeta - d_j} \right)^n. \qquad (13)$$
<figure><img src="image_4.png" /><figcaption>Fig. 3 Example of four plates: successive transformations between the physical plane (e) to almost circular shapes (a) before the iterative procedure of the Theodorsen–Garrick transformations starts</figcaption></figure>

The numerical implementation needs to truncate the infinite series to M terms. It is then checked that the series is convergent; this convergence is actually a function of the proximity of the bodies and their smoothness as well. In order to perform the identification $Im(F) = \psi$ with the Dirichlet conditions in (7), the complex potential F is transported in a coordinate system centered on body k, that is to say, at $\zeta = d_k + \rho_k e^{i\alpha}$. The following series is used

$$ \frac{1}{(1+x)^n} = \sum_{p=0}^{\infty} G(p, n)x^p, \qquad (14) $$

where $G(p, n)$ denotes the binomial coefficients

$$ G(p, n) = (-1)^p \frac{(p + n - 1)!}{(n - 1)!p!}. \qquad (15) $$

We obtain

$$ F(\zeta) = \alpha_{k0} + \sum_{n=1}^{\infty} \alpha_{kn} e^{in\alpha} + a_{kn} e^{-in\alpha}, \quad \text{at } \zeta = d_k + \rho_k e^{i\alpha}, \qquad (16) $$

where the unknown coefficients $\alpha_{kn}$ are given by

$$ \alpha_{kn} = \sum_{j=1, j \neq k}^{N} \sum_{p=1}^{\infty} a_{jp} \left[ \frac{\rho_j}{D_{kj}} \right]^p \left[ \frac{\rho_k}{D_{kj}} \right]^n G(n, p). \qquad (17) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4 Example of 25 identical squares in array. (a) Final and initial shapes. (b) Shapes obtained before the iterative procedure of Theodorsen–Garrick transformations. (c) Rate of convergence of the successive Theodorsen–Garrick transformations</figcaption></figure>

Next we proceed to the Fourier decomposition of $-y(\alpha)$, $x(\alpha)$ and $\frac{1}{2}r^2(\alpha)$, calculated in the coordinate system centered on body k

$$ \left. \begin{array}{l} -y(\alpha) \\ x(\alpha) \\ \frac{1}{2}r^2(\alpha) \end{array} \right\} = \sum_{n=0}^{\infty} A_n \cos n\alpha + B_n \sin n\alpha, \qquad (18) $$

where the Fourier coefficients ($A_n$, $B_n$) are now perfectly known from Sect. 4. To illustrate the convergence of the series, we successively consider two configurations: 2 ellipses and 2 rectangles. Figure 5 shows the considered couple of physical shapes, the variations of $-y(\alpha)$, $x(\alpha)$ and $\frac{1}{2}r^2(\alpha)$ in terms of $\alpha$ for one of the shapes and the corresponding variations of $\sqrt{A_n^2 + B_n^2}$ in terms of the mode number $n$. The series indeed converge at a reasonable rate. The identifications between (18) and (16) provide

$$ \mathfrak{J}\mathrm{m}(\alpha_{kn}) + \mathfrak{J}\mathrm{m}(a_{kn}) = A_n, \qquad (19) $$

$$ \mathfrak{R}\mathrm{e}(\alpha_{kn}) - \mathfrak{R}\mathrm{e}(a_{kn}) = B_n. \qquad (20) $$

Note that the constant $A_0$ will not be required for the purpose of force calculation. Alternatively, in a more compact form, by noting $C_n = A_n + iB_n$, we obtain

$$ iC_n = a_{kn} - \overline{\alpha_{kn}}, \quad n \ge 1, \qquad (21) $$
<figure><img src="image_6.png" /><figcaption>Fig. 5 Illustration for two ellipses (a) and two rectangles (b). From left to right: couple of physical shapes, variations with α of −y(α), x(α) and 1/2 r²(α), corresponding variations of √A²n + B²n with the mode number n</figcaption></figure>

where the overline denotes the complex conjugate. In order to close the problem, it remains to check the Dirichlet conditions on the other body contours, (B_j)_{j≠k}. Written successively in the coordinate system attached to each body, we get

$$0 = a_{\ell n} - \bar{\alpha}_{\ell n}, \quad n \ge 1, \ell \ne k. \qquad (22)$$

Numerically a linear system of order N × 2M is solved. Once this is done, the complex potential is known throughout the fluid domain. This is precisely the variable required in the formulation of [1] to perform the computation of the loads in viscous flows. As a preliminary we revisit the formulation of the added masses and added moment of inertia.

## 6 Hydrodynamic loads in potential theory

Having arrived at this stage, we must clarify some things. In effect, we are considering to derive the hydrodynamic loads computed from the solution of the BVP (6). However this BVP can only provide the added masses and added moment of inertia and nothing else. Hence the obtained loads are part of the actual solution as shown in the following developments. The hydrodynamic loads acting on body B_i are written as follows

$$\vec{F}_i = - \int_{B_i} p \vec{n} \, ds, \quad \vec{M}_i = - \int_{B_i} p (\vec{r} \wedge \vec{n}) \, ds, \qquad (23)$$

where $\vec{n}$ is the normal pointing into the fluid and p denotes the pressure. The moment is calculated with respect to a point O_i attached to the body B_i. The pressure is given by Bernoulli's formula

$$p = -\rho\phi_{,t} - \frac{1}{2}\rho(\vec{\nabla}\phi)^2, \qquad (24)$$

where ρ is the density of the fluid. It should be noted that a direct application of Blasius's theorem is hardly possible at that stage. It would mean that integrals in (23) are performed by using the residue theorem. However, since complicated and intricate conformal mappings are implemented before ending up with circles, the determination of the coefficients of the poles (depending on a precise calculation of the Jacobian of the whole transformation) becomes a tricky problem.
Focusing on the unsteady components of the loads only, we hence proceed as in [18] and arrive at expressions of the forces and moment containing the following integrals

$$ - \rho \int_{B_i} \phi \left\{ \begin{matrix} \vec{n} \\ \vec{r} \end{matrix} \wedge \vec{n} \right\} ds = - \rho \int_{B_i} \phi \phi_{,n} ds. \qquad (25) $$

These coefficients form a tensor of order 4 and, by definition, we define the added masses and added moment of inertia by the quantities

$$ m_{ijkl} = -\rho \int_{B_i} \phi_{k\ell} \phi_{ij,n} ds \qquad (26) $$

which provide the added force or the added moment of inertia acting on body i and induced by the acceleration of the fluid originating from the motion of body j along a given degree of freedom. The potential $\phi_{k\ell}$ is the solution of problem (6) for body k and its degree of freedom $\ell$ and is therefore given by

$$ \phi_{k\ell} = \Re(e(F_{k\ell})). \qquad (27) $$

However, $\phi_{,n}$ corresponds to the solution of the boundary-value problem (6) for body i and its degree of freedom $N^o j$,

$$ \phi_{ij,n} = \Re\left(\frac{dF_{ij}}{dz}n_z\right), \qquad (28) $$

where $n_z$ is the complex normal on $B_i$

$$ n_z(\alpha) = e^{i\alpha} \frac{J(\alpha)}{|J(\alpha)|}, \quad J = \frac{df}{d\zeta} \qquad (29) $$

and the complex velocity (its conjugate in fact) $\frac{dF_{ij}}{dz}$ is alternatively written as

$$ \frac{dF_{ij}}{dz} = \frac{1}{J(\alpha)} \frac{dF_{ij}}{d\zeta}. \qquad (30) $$

Integration is obviously performed in the transformed plane where the body contour $B_i$ is a circle, yielding

$$ m_{ijkl} = \rho \rho_i \int_0^{2\pi} \Re(e(F_{k\ell})\Re(e(\frac{dF_{ij}}{d\zeta}e^{i\alpha}))d\alpha, \qquad (31) $$

where it should be realized that the elemental arc length is $ds = \rho_i|J(\alpha)|d\alpha$. It remains to calculate $F_{k\ell}$ and $F_{ij}$ by transporting their expressions in the coordinate system attached to body $N^o i$, that is to say, at $\zeta = d_i + \rho_i e^{i\alpha}$. The developments (16) are used to express $F_{k\ell}$ in this coordinate system, namely

$$ F_{k\ell} = \alpha_{i0} + \sum_{n=1}^{\infty} \alpha_{in} e^{in\alpha} + a_{in} e^{-in\alpha}. \qquad (32) $$

Then $\frac{dF_{ij}}{d\zeta}e^{i\alpha}$ is turned into

$$ \frac{dF_{ij}}{d\zeta}e^{i\alpha} = \frac{1}{\rho_i} \sum_{n=1}^{\infty} n \left( \alpha_{in} e^{in\alpha} - a_{in} e^{-in\alpha} \right) \qquad (33) $$

by using the following identity

$$ G(p-1, n+1) = -\frac{p}{n}G(p, n). \qquad (34) $$

At this stage, for the sake of simplicity, the coefficients $a_{in}$ and $\alpha_{in}$ have no indices $(k\ell$ or $ij)$ as they appear in $F_{k\ell}$ or $\frac{dF_{ij}}{d\zeta}e^{i\alpha}$ respectively. This is not possible any longer. From the results (19) and (20), the coefficient $m_{ijkl}$ finally read

$$ m_{ijkl} = \pi \rho \sum_{n=1}^{\infty} n \left( 2 \Re(m(\bar{a}_{in}^{(k\ell)} C_n^{(ij)}) + |C_n^{(ij)}|^2 \delta_{ik} \delta_{j\ell} \right). \qquad (35) $$
Clearly, for bodies that are far apart, Eq. 35 reduces to the expressions elaborated in [19] since, in that case $\Re(a_{in}^{(k\ell)}) \approx -B_{n}^{(k\ell)}$ and $\Im(a_{in}^{(k\ell)}) \approx A_{n}^{(k\ell)}$.

In order to check identity (35), a standard integral method is implemented. The number of panels is approximately twice the number of modes in Fourier series (18); in the present case the number of panels is about $2^8$ on each body. It should be also noted that the integral equation is discretized by constant functions on each panel. Table 1 summarizes the comparisons for the configuration with two rectangles already studied in Sect. 4. It is shown that the discrepancy between the two approaches is less than 1%. In order to better compare the present multipole technique with the integral-method results, the case of a single square is considered. The relative error is calculated with respect to the theoretical results in [18]. Figure 6 shows that the convergence is slow (roughly $O(\delta^{1.2})$ where $\delta$ is the average size of a panel), while the multipole technique provides a reasonable convergence at an incomparably lower computational cost.

The summation of the infinite series (35) must be truncated to, say, $N_m$ number of modes. This truncation mainly depends on the proximity of the bodies but also on the rate of convergence of Fourier series (18). In order to check that $N_m$ is well chosen, convergence tests are performed. A limiting case is considered for two very close identical squares. The gap is $g/2a = 0.015$ where $2a$ is the length of one side. The convergence of the successive Theodorsen–Garrick transformations is obtained after 106 iterations. Figure 7 shows the variation of the coefficients $m_{111}, m_{121}, m_{131}, m_{112}$ and $m_{122}$, with the retained number of modes $N_m$. These coefficients have been made non-dimensional with their respective value for $N_m = 100$. It is shown that 100 terms practically suffice to get converged values. It should be noted that the first three coefficients are the diagonal terms of the auto-influence matrix. The last two coefficients determine the influence between the 2 squares.

The next parametric study evaluates the gap between twin squares at which their mutual influence vanishes. First comparisons are made with some known results; for example, the coefficients of the auto-influence matrix (auto-induced forces and moment). Figure 8 shows the variation of $m_{111}, m_{121}$ and $m_{131}$ in terms of the gap (rendered non-dimensional with the length of one side of the square). The relative errors in comparison with semi-analytical results by Newman [23] are plotted in Fig. 9. It is shown that hardly any influence is felt from the presence of a square as soon as this square is located at four times the side length. On the other hand, the present approach agrees with semi-analytical results [18]; the relative error is within 5 digits for the added mass and within 4 digits for the added moment of inertia. Regarding the force-moment induced on one cylinder due the motion of the other, Fig. 10 sums up the variations of the significant coefficients in terms of the gap. The coefficients are made non-dimensional with either the coefficients $m_{11}$ or the coefficients $m_{66}$ corresponding to an isolated square.

The variations of the force coefficients are now illustrated for two similar ellipses as their aspect ratio varies from zero to unity. This means that we progressively shift from two parallel flat plates up to circles. The spacing is such that the shapes touch each other when we end up with circles. The orientation of the plate is vertical. The added moment of inertia must be compared to $\frac{1}{8}\pi\rho a^4$, corresponding to that of an isolated flat plate; this coefficient is used to make the moment non-dimensional. The added masses are made non-dimensional with $\pi\rho a^2$. In all cases $a$ is the highest semi-axis. Figure 11 shows the variations of the coefficients. It is worth noticing the non-monotonical variations of some coefficients, for example, $m_{1323}$ which corresponds to the moment induced on one body due to the accelerated rotation of the other one. The minimum occurs at about the aspect ratio 0.55. This is the case for $m_{1213}, m_{1322}$ and $m_{1223}$ but for different values of the aspect ratios. Note that the latter two coefficients are identical, but of opposite sign.

# 7 Hydrodynamic loads in viscous flow

There are two main results obtained in the developments described above. One is the derivation of a new formula for the added masses and added moment of inertia in potential flow. The second result is that the obtained potential solution makes it possible to implement the formulation given in [1].
<table><caption>**Table 1** Comparison of $m_{ijkl}$ for the two rectangles studied in Sect. 4</caption><thead><tr><th>i</th><th>k</th><th>j</th><th>ℓ</th><th>m<sub>ijkl</sub> (semi-ana)</th><th>m<sub>ijkl</sub> (numer.)</th><th>relative error</th></tr></thead><tbody><tr><td>1</td><td>1</td><td>1</td><td>1</td><td>8.655072E + 00</td><td>8.670987E + 00</td><td>1.835453E - 03</td></tr><tr><td>1</td><td>1</td><td>1</td><td>2</td><td>1.942108E + 00</td><td>1.944632E + 00</td><td>1.297853E - 03</td></tr><tr><td>1</td><td>1</td><td>1</td><td>3</td><td>-3.717540E - 03</td><td>-3.829630E - 03</td><td>2.926902E - 02</td></tr><tr><td>1</td><td>1</td><td>2</td><td>1</td><td>1.942108E + 00</td><td>1.944632E + 00</td><td>1.297874E - 03</td></tr><tr><td>1</td><td>1</td><td>2</td><td>2</td><td>1.952990E + 00</td><td>1.960361E + 00</td><td>3.759925E - 03</td></tr><tr><td>1</td><td>1</td><td>2</td><td>3</td><td>-3.000904E - 03</td><td>-3.113974E - 03</td><td>3.631055E - 02</td></tr><tr><td>1</td><td>1</td><td>3</td><td>1</td><td>-3.717540E - 03</td><td>-3.831116E - 03</td><td>2.964567E - 02</td></tr><tr><td>1</td><td>1</td><td>3</td><td>2</td><td>-3.000904E - 03</td><td>-3.115026E - 03</td><td>3.663602E - 02</td></tr><tr><td>1</td><td>1</td><td>3</td><td>3</td><td>2.358275E + 00</td><td>2.372969E + 00</td><td>6.192253E - 03</td></tr><tr><td>1</td><td>2</td><td>1</td><td>1</td><td>-4.880324E - 01</td><td>-4.894823E - 01</td><td>2.962201E - 03</td></tr><tr><td>1</td><td>2</td><td>1</td><td>2</td><td>-8.636832E - 01</td><td>-8.656954E - 01</td><td>2.324456E - 03</td></tr><tr><td>1</td><td>2</td><td>1</td><td>3</td><td>-3.908437E - 01</td><td>-3.921240E - 01</td><td>3.265083E - 03</td></tr><tr><td>1</td><td>2</td><td>2</td><td>1</td><td>-5.330805E - 01</td><td>-5.343968E - 01</td><td>2.463144E - 03</td></tr><tr><td>1</td><td>2</td><td>2</td><td>2</td><td>-3.178861E - 01</td><td>-3.186361E - 01</td><td>2.353680E - 03</td></tr><tr><td>1</td><td>2</td><td>2</td><td>3</td><td>-1.223604E - 02</td><td>-1.226253E - 02</td><td>2.159737E - 03</td></tr><tr><td>1</td><td>2</td><td>3</td><td>1</td><td>1.292289E - 01</td><td>1.296690E - 01</td><td>3.393922E - 03</td></tr><tr><td>1</td><td>2</td><td>3</td><td>2</td><td>7.124675E - 02</td><td>7.148246E - 02</td><td>3.297370E - 03</td></tr><tr><td>1</td><td>2</td><td>3</td><td>3</td><td>2.536999E - 03</td><td>2.562158E - 03</td><td>9.819273E - 03</td></tr><tr><td>2</td><td>1</td><td>1</td><td>1</td><td>-4.880324E - 01</td><td>-4.894915E - 01</td><td>2.980849E - 03</td></tr><tr><td>2</td><td>1</td><td>1</td><td>2</td><td>-5.330805E - 01</td><td>-5.343723E - 01</td><td>2.417473E - 03</td></tr><tr><td>2</td><td>1</td><td>1</td><td>3</td><td>1.292289E - 01</td><td>1.296207E - 01</td><td>3.022663E - 03</td></tr><tr><td>2</td><td>1</td><td>2</td><td>1</td><td>-8.636832E - 01</td><td>-8.656911E - 01</td><td>2.319500E - 03</td></tr><tr><td>2</td><td>1</td><td>2</td><td>2</td><td>-3.178861E - 01</td><td>-3.186285E - 01</td><td>2.330048E - 03</td></tr><tr><td>2</td><td>1</td><td>2</td><td>3</td><td>7.124675E - 02</td><td>7.145205E - 02</td><td>2.873307E - 03</td></tr><tr><td>2</td><td>1</td><td>3</td><td>1</td><td>-3.908437E - 01</td><td>-3.922900E - 01</td><td>3.686924E - 03</td></tr><tr><td>2</td><td>1</td><td>3</td><td>2</td><td>-1.223604E - 02</td><td>-1.228303E - 02</td><td>3.825193E - 03</td></tr><tr><td>2</td><td>1</td><td>3</td><td>3</td><td>2.536999E - 03</td><td>2.562342E - 03</td><td>9.890406E - 03</td></tr><tr><td>2</td><td>2</td><td>1</td><td>1</td><td>1.931193E + 01</td><td>1.935089E + 01</td><td>2.013638E - 03</td></tr><tr><td>2</td><td>2</td><td>1</td><td>2</td><td>9.659566E + 00</td><td>9.671692E + 00</td><td>1.253671E - 03</td></tr><tr><td>2</td><td>2</td><td>1</td><td>3</td><td>1.098452E - 02</td><td>1.131665E - 02</td><td>2.934871E - 02</td></tr><tr><td>2</td><td>2</td><td>2</td><td>1</td><td>9.659566E + 00</td><td>9.671690E + 00</td><td>1.253532E - 03</td></tr><tr><td>2</td><td>2</td><td>2</td><td>2</td><td>1.232494E + 01</td><td>1.235513E + 01</td><td>2.443852E - 03</td></tr><tr><td>2</td><td>2</td><td>2</td><td>3</td><td>2.791748E - 02</td><td>2.843748E - 02</td><td>1.828546E - 02</td></tr><tr><td>2</td><td>2</td><td>3</td><td>1</td><td>1.098452E - 02</td><td>1.132160E - 02</td><td>2.977343E - 02</td></tr><tr><td>2</td><td>2</td><td>3</td><td>2</td><td>2.791748E - 02</td><td>2.845000E - 02</td><td>1.871753E - 02</td></tr><tr><td>2</td><td>2</td><td>3</td><td>3</td><td>1.810144E + 01</td><td>1.822819E + 01</td><td>6.953258E - 03</td></tr></tbody></table>

The indices *i* and *k* correspond the body $B_i$ and $B_k$, respectively. The index *j* corresponds to the force ($j = 1$ or $j = 2$) and moment ($j = 3$) acting on the body $B_i$. The index *ℓ* corresponds to the translation ($ℓ = 1$ or $ℓ = 2$) and rotation ($ℓ = 3$) of the body $B_k$. Column 5: present semi-analytical results. Column 6: results from a standard integral method. Column 7: relative error between the two computations
<figure><img src="image_7.png" /><figcaption>Fig. 6 Comparison of results obtained either by a panel method (constant function over each panel) or the present multipole approach. The relative error between the two approaches is plotted in terms of the number of panels. The references are the theoretical results by Newman [18]: $m_{11} \approx 4.75376$ and $m_{66} \approx 0.72457574$; the multipole technique provides: $m_{11} \approx 4.75376166$ and $m_{66} \approx 0.72457944$</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 Left: physical configuration; the gap between the two faces of the squares is $g/2a = 0.015$ where $2a$ is the length of one face. Right: variation of coefficients $m_{1111}$, $m_{1212}$, $m_{1313}$, $m_{1121}$ and $m_{1222}$, with the computed number of modes $N_m$. The coefficients are made nondimensional with their respective value for $N_m = 100$</figcaption></figure>

As mentioned previously, this is an alternative to the integration of the pressure following its explicit calculation. This can be done for a two-dimensional configuration by integrating the tangential pressure gradient resulting from the Navier–Stokes equation written on the contour (supposed fixed in time). For three-dimensional configurations, it is more complicated since a Poisson equation must be solved.

The principle follows from the inner product in the Hilbert space $L^2(\Omega)$ (space of square integrable functions) of the pressure gradient with the gradient of a function which is precisely the velocity potential solution of (6).

The pressure gradient is given by the Navier–Stokes equation (conservation of momentum). After some algebra, the final identity reads

$$ \frac{1}{\rho} \int_{\partial \Omega} p \phi_{,n} ds = \int_{\partial \Omega} \phi(\vec{u}_{,t} \cdot \vec{n}) ds + \frac{1}{2} \int_{\partial \Omega} \vec{u}_{,n}^2 \phi_{,n} ds + \int_{\Omega} \vec{\nabla} \phi \cdot (\vec{u} \wedge \text{rot} \vec{u}) d\omega + \nu \int_{\partial \Omega} \vec{\nabla} \phi \cdot (\vec{n} \wedge \text{rot} \vec{u}) ds, \quad (36) $$

where $\partial \Omega$ denotes the entire contour of the fluid domain, say $\Sigma \cup \sum_{j=1}^N B_j$, and $\Sigma$ is a control surface where the flow is matched with an undisturbed flow (there is no vorticity).

In the left-hand side of (36), we have to extract the integration of the pressure on a single body. For that purpose, $\phi$ must be a solution of a boundary-value of type (6). In the right-hand side, the different contributions must be calculated and some of them may disappear, depending on the considered flow conditions.
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Two squares of length 2a (same configuration than in Fig. 7). Variation of $m_{1111}$, $m_{1212}$ and $m_{1313}$ with the gap g. Nondimensionalization: coefficients with $ρa^2$ or $ρa^4$, gap with 2a. Comparison with theoretical results in [18]</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Fig. 9** Two squares of length 2a (same configuration than in Fig. 7). Variation with the gap g/2a of the relative errors between the semi-analytical results for single square in [18] and the added masses $m_{1212}$ and $m_{2121}$ and added moment of inertia $m_{1313}$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10** Two squares of length 2a (same configuration than in Fig. 7). Variations of some significant coefficients with the gap g/2a. The coefficients are made nondimensional with coefficients for a single body (see [18]): $m_{11} = 4.75376$ or $m_{66} = 0.72457574$</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Fig. 11 Top: physical configuration of two identical flat plates, ellipses and circles. The semi axes of the ellipses are (a, b), the largest semi-axis (a) is oriented vertically. The aspect ratio b/a varies from zero to unity: flat plate to a circle. Variation of some coefficients with the aspect ratio b/a. Nondimensionalization: added masses with ρa² and added moment of inertia with ρa⁴</figcaption></figure>

As an application we formulate the force acting on fixed bodies placed in a uniform flow. Then identity (36) simplifies since $\vec{u}$ and its time derivative vanish on the solid contour except on $\Sigma$. Protas et al. [24] and Pan and Chew [25] examined in detail the present formulation for a single body. Part I presented some complementary calculations as well. Hence, only the main outlines are presented here. In particular, Eq. 36 can be decomposed in three contributions corresponding to the following three integrals:

- on the control surface $\Sigma$, they are denoted $I_{acc}$ and called inertia force since acceleration of the fluid at infinity appears,

- on the fluid domain, they are denoted $I_{\omega}$ and called vorticity force,

- on the body contour, they are denoted $I_{v}$ and called viscous force since they contribute to the friction force.

## 7.1 Inertia forces $I_{acc}$

The corresponding contributions read

$$I_{acc} = \int_{\Sigma} \phi_{k\ell}(\vec{u}_{,t} \cdot \vec{n})ds, \qquad (37)$$
where $\phi_{k\ell}$ is the solution of the boundary-value problem (6) relative to a given body (k being the one on which the force is calculated) and for the degree of freedom $\ell$, thus providing the force along this degree of freedom. The complex potential is expressed as follows

$$F_{k\ell}(\zeta) = \sum_{j=1}^{N} \sum_{n=1}^{\infty} a_{jn}^{(k\ell)} \left( \frac{\rho_j}{\zeta - d_j} \right)^n . \qquad (38)$$

Its real part gives $\phi_{k\ell}$. The bodies are assumed to be placed in an unsteady flow $\vec{v}_{\infty} = u_{\infty}\vec{x} + v_{\infty}\vec{y}$ whose Cartesian components are functions of time only. Asymptotically far from the bodies, the flow is considered as a potential flow; therefore the components of the velocity will only depend on the conformal transformation. It is easily checked that this Jacobian is a constant when $|z| \to \infty$. We denote this Jacobian by $J^\infty = (dz/d\zeta)_\infty$. This means that we can determine the modification of the orientation of the normal on $\Sigma$ due to the transformation. Hence the inner product $(\vec{u}_{,t} \cdot \vec{n})$ reduces to

$$(\vec{u}_{,t} \cdot \vec{n})_z = \dot{u}_{\infty} \cos(\alpha_{\infty} + \alpha) + \dot{v}_{\infty} \sin(\alpha_{\infty} + \alpha), \qquad (39)$$

where $\alpha_{\infty}$ denotes the argument of the complex variable $J^\infty$. The integration is advantageously performed in the transformed plane. The surface $\Sigma$ is a circle of large radius $R_\Sigma$ and a point on this surface is $\zeta = R_\Sigma e^{i\alpha}$, yielding finally, as $R_\Sigma \to \infty$,

$$I_{\text{acc}} = \pi \sum_{j=1}^{N} \rho_j \Re \left[ J^\infty a_{j1}^{(k\ell)} (\dot{u}_\infty - i\dot{v}_\infty) \right] \qquad (40)$$

## 7.2 Viscous forces $I_v$

The total friction force has two contributions. The integral $I_v$ must now be calculated. Also, we must integrate the friction appearing in the Cauchy stress tensor. On body $B_k$ the latter contribution reads

$$\vec{F}_f = \rho v \int_{B_k} \omega \vec{t} \, ds, \qquad (41)$$

where the system $(\vec{n}, \vec{t}, \vec{z})$ is right-handed. The former contribution $I_v$ is turned into

$$I_v = \rho v \sum_{j=1}^{N} \int_{B_j} \omega \phi_{k\ell,s} \, ds, \qquad (42)$$

where $\phi_{k\ell,s}$ denotes the tangential gradient of the potential. The index $\ell$ refers to the direction of the force (1 or 2 corresponding to the $\vec{x}$ or $\vec{y}$-direction, respectively). The tangential derivative of $\phi$ is conveniently expressed by using the corresponding complex potential

$$\phi_{k\ell,s} = -\frac{1}{|J|} \Re m \left( \frac{dF_{k\ell}}{d\zeta} e^{i\alpha} \right). \qquad (43)$$

By using (33), we have that $I_v$ can finally be written as

$$I_v = -\rho v \sum_{j=1}^{N} \int_{B_j} \frac{\omega_{\zeta}}{|J|^2} \Re m \left[ \sum_{n=1}^{\infty} n \left( \alpha_{jn}^{(k\ell)} e^{in\alpha} - a_{jn}^{(k\ell)} e^{-in\alpha} \right) \right], \qquad (44)$$

where the vorticity $\omega_\zeta$ in the transformed plane is introduced.<fnref n="1" /> The expression (44) is suitable for computation since all terms in the integrand (including $\omega_\zeta/|J|^2$) can be expressed as Fourier series.

<footnote n="1">¹ The vorticity calculated in the $\zeta$-plane differs from the vorticity calculated in the $z$-plane by a factor $|J|^2$, say $\omega(\zeta) = |J|^2\omega(z)$ since locally the circulation of a discrete vortex would be the same in the two planes.</footnote>
## 7.3 Vorticity forces $I_\omega$

This contribution is examined within the framework of a Lagrangian approach for the vorticity transportation. The vorticity field is represented by means of a finite sum of $N_v$ discrete vortices:

$$ \omega(x, y, t) = \sum_{q=1}^{N_v} \Gamma_q \delta(x - x_q) \delta(y - y_q). \qquad (45) $$

Each discrete voxtex carries a small amount of circulation $\Gamma_q$. The force $I_\omega$ then reduces to

$$ I_\omega = -\rho \sum_{q=1}^{N_v} \Gamma_q (v_q \phi_{,x} - u_q \phi_{,y}), \qquad (46) $$

where $(u_q, v_q)$ are the Cartesian coordinates of the physical velocity at a point vortex $q$. The gradient of the velocity potential is derived from the complex potential:

$$ \phi_{,x} - i\phi_{,y} = -\frac{1}{J} \sum_{j=1}^{N} \sum_{n=1}^{\infty} na_{jn} \frac{\rho_j^n}{(\zeta - d_j)^{n+1}}. \qquad (47) $$

These expressions are used for the computation of the force and moment by choosing the proper potential $\phi$. It should be noted that we end up with a triple summation which might be time-consuming. However, two of them (over the modes $n$ and over the bodies $j$) can be achieved once and for all at the nodes of a fictitious mesh covering the support of the vorticity pattern only, namely the same which is used to compute the velocity $(u_q, v_q)$ at a point vortex $q$.

## 7.4 A posteriori pressure formulation

More interesting is the computation of the pressure. It is decomposed as a Fourier series on a given body contour $(B_k)$. Consistent with the previous developments, the pressure is parametrized with the angular coordinate $\alpha$ defined on the transformed plane

$$ p(\alpha) = \sum_{n=0}^{\infty} A_n \cos n\alpha + B_n \sin n\alpha, \qquad (48) $$

where the Fourier coefficients $(A_n, B_n)$ must be calculated. It is shown that solving the following boundary-value problem,

$$ \left\{ \begin{array}{ll} \Delta\phi = 0, & \text{in the fluid domain } \Omega, \\ \vec{\nabla}\phi \cdot \vec{n} = \frac{1}{|J|} \begin{cases} \cos m\alpha, & \text{on the body surface } (B_k) (m > 0), \\ \sin m\alpha, & \text{on the body surface } (B_j) (j \neq k), \end{cases} \\ \vec{\nabla}\phi \cdot \vec{n} = 0, & \text{on the body surface } (B_j) (j \neq k), \\ \vec{\nabla}\phi \to 0, & \text{at infinity,} \end{array} \right. \qquad (49) $$

provides a set of functions needed to compute the coefficients $(A_n, B_n)$ in Eq. 18. The method of solution developed in Sect. 5 is used here. It leads to linear systems to be solved:

$$ \alpha_{kn} - \bar{a}_{kn} = -\frac{\rho_k}{m} \delta_{nm}, \quad \forall n > 0, \qquad (50) $$

$$ \alpha_{jn} - \bar{a}_{jn} = 0 \quad j \neq k, \quad \forall n > 0. \qquad (51) $$

These linear systems are similar to those that have already been solved but with different right-hand sides. The corresponding solutions can be calculated once and for all as long as the relative positions of the bodies do not change in time, which is the case in the present application.
## 8 Conclusion

The aim of this study has been to improve the computation of the hydrodynamic loads on multi-body configurations of arbitrary shape. This has been applied to two-dimensional configurations only. Two aspects of the problem have been investigated. One aspect has led to the revisiting of the formulation of the added masses and added moment of inertia. An original formulation was established based on a multipole approach. The other aspect offers the opportunity to perform computations of viscous flows in simpler domains as soon as complicated shapes are considered. Optimized algorithms exist already for arrays of circular sections. They can now be extended to arbitrary shapes via the implementation of conformal mappings that are well known. Robust algorithms have been developed here in order to tackle these conformal transformations. Rectangles and flat plates can be handled up to small gaps and rates of convergence are shown to be fast enough.

More applications of practical interest can be considered. As an example, the arrays of cylinders in exchangers of nuclear-power plants have typically non-circular sections. It is expected that the present approach will bring some new insights into the evaluation of hydrodynamic loads acting on such devices. However, it is doubtful whether or not the present approach can compete with other approaches, such as an integral-equation technique. Without doubt, computational tests have shown that the gain is significant. On the other hand, when zero-thickness bodies (like flat plates, bodies with fins,...), for which an integral-equation technique may suffer from numerical difficulties, the present approach seems to remain unchallenged as long as a suitable conformal-mapping transformation exists.

## References

1. Quartapelle L, Napolitano M (1983) Force and moment in incompressible flows. AIAA J 21:911-913

2. Etienne S (1999) Contribution a la modélisation de l'écoulement de fluide visqueux autour de faisceaux de cylindres circulaires. PhD Thesis, University of Aix-Marseille 2

3. Crowdy DG (2006) Calculating the lift on a finite stack of circular aerofoils. Proc Roy Soc A 462:1387-1407

4. Kochin NE, Kibel IA, Roze NV (1964) Theoretical hydromechanics. Intersciences Publishers, New-York

5. Lavrentyev M, Shabat B (1977) Méthodes de la théorie des fonctions d'une variable complexe. Edition du MIR (in French), Moscou

6. Venkatesan SK (1985) Added mass of two cylinders. J Ship Res 19(4):234-240

7. Lagally M (1929) Die reibungslose Strömung im Aussengebiet zweier Kreise (English translation: the frictionless flow in the region around two circles). ZAMM 9(4):299-305

8. Ives DC (1976) A modern look at conformal mapping including multiply connected regions. AIAA J 14(8):1006-1011

9. Halsey N (1979) Potential flow analysis of multielement airfoils using conformal mapping. AIAA J 17:1281-1288

10. Suddhoo A (1985) Inviscid compressible flow past multi-element airfoils. PhD Thesis, University of Manchester

11. Prosnak WJ (1987) Computation of fluid motions in multiply connected domains. Wissenschaft + Technik, Brau Karlsruhe

12. Scolan Y-M, Faltinsen O (1993) Numerical prediction of viscous flows around multi-bodies by a Vortex Method. Proc. 6th international conference on numerical ship hydrodynamics, Iowa, USA, pp 439-453

13. Basset AB (1888) A treatise on hydrodynamics. Cambridge University Press

14. Wegmann R (1986) An iterative method for the conformal mapping of doubly connected regions. J Comput Appl Math 14(1-2):79-98

15. Luchini P, Manzo F (1989) Flow around simply and multiply connected bodies-a new iterative scheme for conformal mapping. AIAA J 27(3):345-351

16. Theodorsen T, Garrick IE (1933) General potential theory of arbitrary wing sections. NACA Rept. 452:177-209

17. Chu JH, Choi DH, Kim MU (1999) Two-dimensional viscous flow past a circular arc. Fluid Dyn Res 25(5):275-291

18. Newman JN (1977) Marine hydrodynamics. MIT Press, Cambridge

19. Scolan Y-M (2005) Some aspects of the potential flow around rotating bodies. Comptes Rendus Mecanique 333(6):487-492

20. Warchawski SE (1945) On Theodorsen's method of conformal mapping of nearly circular regions. Quar J Appl Math 3(2):12-28

21. Kober H (1957) Dictionary of conformal representations. Dover Publications, Inc

22. Nehari Z (1975) Conformal mapping. Dover Publications, New York

23. Newman JN (1979) Added moment of inertia of two-dimensional cylinders. J Ship Res 23:1-8

24. Protas B, Styczek A, Nowakowski A (2000) An effective approach to computation of forces in viscous incompressible flows. J Comput Phys 159(2):231-245

25. Pan LS, Chew YT (2002) A general formula for calculating forces on a 2-d arbitrary body in incompressible flow. J Fluids Struct 16(1):71-82