2.2.1. Discrete surface and its immersion in the grid

We enforce the no-penetration condition at a finite number $N_s$ of discrete surface forcing points; the space of scalar data on these Lagrangian points is denoted by $\mathcal {S}^{N_s}$. In particular, let us define $\boldsymbol{r}_{x}, \boldsymbol{r}_{y} \in \mathcal {S}^{N_s}$ as the vectors of $x$ and $y$ coordinates of the surface points. Each surface point $p$ is associated with a small straight segment of length $\Delta S_p$. Some of the calculations will require information about the local surface orientation. For this purpose, we define vectors $\boldsymbol{\mathsf{n}}_{x}, \boldsymbol{\mathsf{n}}_{y} \in \mathcal {S}^{N_s}$ of components of the discrete surface unit normals. We also define the space $\mathcal {V}^{N_s}$ to hold vector-valued data, such as velocity, on the immersed surface points. For convenience, let us also define unit vectors $\boldsymbol{\mathsf{e}}_{p}$ on this space, equal to 1 at surface point $p$ ($1 \leqslant p \leqslant N_s$) and zero at every other point. For example, the $x$ coordinate of point $p$ is picked out of the vector $\boldsymbol{r}_{x}$ by projection onto the $p$th unit vector

(2.18)\begin{equation} \boldsymbol{\mathsf{e}}_{p}^{T}\boldsymbol{r}_{x}. \end{equation}

Each of the surface point spaces has an associated inner product that includes the surface length, e.g.

(2.19)\begin{equation} {\langle \boldsymbol{\mathsf{g}}_1,\boldsymbol{\mathsf{g}}_2 \rangle}_{\mathcal{S}^{N_s}} = \sum_p \Delta S_p g_{1,p} g_{2,p}, \end{equation}

for any $\boldsymbol{\mathsf{g}}_1, \boldsymbol{\mathsf{g}}_2 \in \mathcal {S}^{N_s}$, so that the inner product approximates a surface integral. Another vector we will make substantial use of in this paper is $\boldsymbol{\mathsf{1}}_{\mathcal{S}} \in \mathcal {S}^{N_s}$, a vector of ones on all surface points.

From any vector $\boldsymbol{\mathsf{g}} \in \mathcal {S}^{N_s}$, we can also form a diagonal $N_s\times N_s$ operator $\boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{g}}}$ with the entries of the vector along the diagonal. When this operator acts upon another vector $\boldsymbol{\mathsf{h}}\in \mathcal {S}^{N_s}$, it represents the Hadamard (i.e. element-by-element) product of the two vectors, $\boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{g}}} \boldsymbol{\mathsf{h}} = \boldsymbol{\mathsf{g}}\circ \boldsymbol{\mathsf{h}} \in \mathcal {S}^{N_s}$. Note that $\boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{g}}}\boldsymbol{\mathsf{h}} = \boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{h}}} {\boldsymbol{\mathsf{g}}}$, and that $\boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{g}}}\boldsymbol{\mathsf{1}}_{\mathcal{S}} = \boldsymbol{\mathsf{g}}$.

Similar to (2.13), surface data are immersed into the grid with the regularization operator $ {\mathcal {N}}{}{\boldsymbol{\mathsf{R}}_\mathcal {S}}: \mathcal {S}^{N_s} \mapsto \mathcal {N}$. Grid data are interpolated onto the surface points with the interpolation operator $ {\mathcal {S}}{}{\boldsymbol{\mathsf{E}}_\mathcal {N}}: \mathcal {N}\mapsto \mathcal {S}^{N_s}$; $ {\mathcal {N}}{}{\boldsymbol{\mathsf{R}}_\mathcal {S}}$ can be constructed (and we will assume it has) so that it is the transpose of the interpolation operator, $ {\mathcal {N}}{}{\boldsymbol{\mathsf{R}}_\mathcal {S}} = {\mathcal {S}}{}{\boldsymbol{\mathsf{E}}_\mathcal {N}}^{T}$, with respect to the grid and surface inner products defined in this paper. Furthermore, note that $ {\mathcal {N}}{}{\boldsymbol{\mathsf{R}}_\mathcal {S}}$ and $ {\mathcal {S}}{}{\boldsymbol{\mathsf{E}}_\mathcal {N}}$ can be constructed with a different choice for the discrete Dirac delta function than the one used for the vortex regularization. In this work, we use the smoothed three-point function from Yang et al. (Reference Yang, Zhang, Li and He2009) (figure 1) and we use a uniform spacing between the surface points.

2.2.2. The immersed surface potential flow problem

The surface's motion is specified by a velocity distribution $\boldsymbol {v}_{b}$, represented discretely by components $\boldsymbol{\mathsf{v}}_{b,x}, \boldsymbol{\mathsf{v}}_{b,y} \in \mathcal {S}^{N_s}$. For rigid bodies, this surface motion can be described by a streamfunction. Specifically, translation at velocity $(U,V)$ and rotation at angular velocity $\varOmega$ would be described equivalently by velocity components

(2.20a,b)\begin{equation} \boldsymbol{v}_{b,x} = U \boldsymbol{\mathsf{1}}_{\mathcal{S}} - \varOmega \boldsymbol{r}_{y}, \quad \boldsymbol{v}_{b,y} = V \boldsymbol{\mathsf{1}}_{\mathcal{S}} + \varOmega \boldsymbol{r}_{x} \end{equation}

or by a surface streamfunction $\boldsymbol{\mathsf{s}}_{b} \in \mathcal {S}^{N_s}$

(2.21)\begin{equation} \boldsymbol{\mathsf{s}}_{b} = U\boldsymbol{r}_{y} - V\boldsymbol{r}_{x} -\tfrac{1}{2}\varOmega (\boldsymbol{\mathsf{D}}_{\boldsymbol{r}_{x}}\boldsymbol{r}_{x} + \boldsymbol{\mathsf{D}}_{\boldsymbol{r}_{y}}\boldsymbol{r}_{y}). \end{equation}

The no-penetration condition can be imposed by setting the discrete streamfunction equal to that of the surface, up to a uniform value, $\boldsymbol{\mathsf{s}}_{0} \in \mathcal {S}^{N_s}$

(2.22)\begin{equation} {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} \boldsymbol{\mathsf{s}} = \boldsymbol{\mathsf{s}}_{b} - {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} \boldsymbol{\mathsf{s}}_{\infty} - \boldsymbol{\mathsf{s}}_{0}. \end{equation}

For later shorthand, we will denote the difference between the body motion streamfunction and interpolated uniform flow streamfunction by $\boldsymbol{\mathsf{s}}'_{b} \equiv \boldsymbol{\mathsf{s}}_{b} - {\mathcal {S}}{}{\boldsymbol{\mathsf{E}}_\mathcal {N}} \boldsymbol{\mathsf{s}}_{\infty }$. This modified streamfunction simply consists of subtracting the components $(U_\infty,V_\infty )$ of the uniform flow from $(U,V)$ in (2.21). The uniform value $\boldsymbol{\mathsf{s}}_{0}$ is left unspecified and will later serve the role of enforcing a constraint on circulation. For now, we will suppose that it can be set arbitrarily.

The no-penetration constraint is enforced in the basic potential flow problem (2.5) with the help of a vector of Lagrange multipliers, $\boldsymbol{\mathsf{f}} \in \mathcal {S}^{N_s}$, on the surface points. The modified potential flow problem is thus

(2.23)\begin{equation} \boldsymbol{\mathsf{L}} \boldsymbol{\mathsf{s}} ={-}(\boldsymbol{\mathsf{w}} + {\mathcal{N}}{}{\boldsymbol{\mathsf{R}}_\mathcal{S}} \boldsymbol{\mathsf{f}} ). \end{equation}

In fact, by simple comparison with the vorticity $\boldsymbol{\mathsf{w}}$, it is clear that the vector $\boldsymbol{\mathsf{f}}$ represents the strength of the discrete bound vortex sheet on the surface and serves as another source term of the Poisson equation. Suppose we consider the bound vortex sheet $\gamma (u)$ that emerges from the analogous continuous problem on the undiscretized surface, where $u$ is the arc-length parameter along the surface. At each point $p$, the discrete solution $\boldsymbol{\mathsf{f}}$ is approximately equal to this continuous solution

(2.24)\begin{equation} \boldsymbol{\mathsf{e}}_{p}^{T}\boldsymbol{\mathsf{f}} \simeq \gamma(u_p). \end{equation}

The vector of Lagrange multipliers $\boldsymbol{\mathsf{f}}$ is initially unknown. Thus, the potential flow problem in the presence of the impenetrable surface is

(2.25)\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{L}} & {\mathcal{N}}{}{\boldsymbol{\mathsf{R}}_\mathcal{S}} \\ {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} & 0 \end{bmatrix} \begin{pmatrix} \boldsymbol{\mathsf{s}} \\ \boldsymbol{\mathsf{\,f}} \end{pmatrix} = \begin{pmatrix} -\boldsymbol{\mathsf{w}} \\ \boldsymbol{\mathsf{s}}'_{b} - \boldsymbol{\mathsf{s}}_{0} \end{pmatrix}. \end{equation}

This problem (2.25) has the structure of a generic saddle-point problem (Benzi, Golub & Liesen Reference Benzi, Golub and Liesen2005). We will encounter many such systems in this work, so in Appendix A we summarize a solution approach, based on block-LU (LU, lower-upper) decomposition. The generated solution algorithm of (2.25) is

(2.26)\begin{gather} \boldsymbol{\mathsf{L}} \boldsymbol{\mathsf{s}}^{*} ={-}\boldsymbol{\mathsf{w}} , \end{gather}

(2.27)\begin{gather}\boldsymbol{\mathsf{S}} \boldsymbol{\mathsf{f}} = \boldsymbol{\mathsf{s}}'_{b} - \boldsymbol{\mathsf{s}}_{0} - {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} \boldsymbol{\mathsf{s}}^{*} , \end{gather}

(2.28)\begin{gather}\boldsymbol{\mathsf{s}} = \boldsymbol{\mathsf{s}}^{*} - \boldsymbol{\mathsf{L}}^{{-}1} {\mathcal{N}}{}{\boldsymbol{\mathsf{R}}_\mathcal{S}} \boldsymbol{\mathsf{f}}, \end{gather}

where the Schur complement $\boldsymbol{\mathsf{S}}$ is

(2.29)\begin{equation} \boldsymbol{\mathsf{S}} ={-}{\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} \boldsymbol{\mathsf{L}}^{{-}1} {\mathcal{N}}{}{\boldsymbol{\mathsf{R}}_\mathcal{S}}. \end{equation}

Based on the properties of the matrices comprising $\boldsymbol{\mathsf{S}}$, this operator is symmetric and negative definite, and therefore invertible. Its inverse $\boldsymbol{\mathsf{S}}^{-1}$, also symmetric, maps a surface distribution of streamfunction to a corresponding bound vortex sheet strength. Note that the computation of $\boldsymbol{\mathsf{f}}$ through $\boldsymbol{\mathsf{S}}^{-1}$ is sensitive to both the ratio of discrete surface spacing to the grid spacing $\Delta S/\Delta x$ and the choice of discrete Dirac delta function. On the one hand, small values for $\Delta S/\Delta x$ and discrete Dirac delta functions with small support generally both lead to more high-frequency noise in $\boldsymbol{\mathsf{f}}$. This is because the underlying continuous problem for $\boldsymbol{\mathsf{S}}^{-1}$ is a Fredholm integral equation of the first kind, which is ill posed. As a result, the discrete analogue problem is poorly conditioned (Goza et al. Reference Goza, Liska, Morley and Colonius2016). On the other hand, values for $\Delta S/\Delta x$ that are too high can lead to streamlines penetrating a surface. We found that values between one to four can provide a good balance between the smoothness of $\boldsymbol{\mathsf{f}}$ and the accuracy of the streamlines near a surface.

We can describe this algorithm in words: first, solve for the intermediate streamfunction field, associated with vorticity in the fluid, but without regard for the presence of the surface. Second, find the bound vortex sheet whose associated streamfunction cancels the difference between the specified streamfunction on the surface and the intermediate streamfunction evaluated on the surface. Finally, correct the intermediate streamfunction field for the influence of the bound vortex sheet. A version of the Julia code that implements this algorithm, as well as the algorithms in the following sections, are available in the authors’ Github repository (Beckers & Eldredge Reference Beckers and Eldredge2021).

We give two examples of the streamfunction with a body present and show the associated vortex sheet strength. Figure 3 shows a vortex near a circular cylinder and figure 4 shows a circular cylinder that translates horizontally. In both cases, the vortex sheet strength is in good agreement with the analytical solution from potential flow theory.

Figure 3.

(a) Contours (——) of the discrete streamfunction for a point vortex ($\bullet$) with strength $\varGamma _{v,1}$ at $(R_v,0)$ near a circular cylinder of radius $R_c$ and with a bound circulation $-\varGamma _{v,1}$. The inset shows a closeup of the nodes ($\circ$) and surface points ($\blacksquare$) in the boxed area (- - - -). (b) The scaled discrete vortex sheet strength (——) as a function of the angle $\theta$ measured counter-clockwise from the positive $x$-axis. Here, $k=1$ corresponds to the right-most point on the surface and increases counter-clockwise. Overlaid is the exact continuous solution (- - - -, red dashed). The simulation is performed with $R_v/R_c = 3/2$, $\Delta x/R_c=0.03$ and $\Delta S/\Delta x=2$.

Figure 4.

(a) Contours (——) of the discrete streamfunction for a horizontally translating circular cylinder with radius $R$, and (b) its scaled discrete vortex sheet strength (——). Overlaid is the exact continuous solution (- - - -, red dashed). The simulation is performed with $\Delta x/R=0.03$ and $\Delta S/\Delta x=2$.

For later use, we note that the solution of (2.25) can also be written in inverse form using (A7)

(2.30)\begin{equation} \begin{pmatrix} \boldsymbol{\mathsf{s}} \\ \,\boldsymbol{\mathsf{f}} \end{pmatrix} = \begin{bmatrix} \boldsymbol{\mathsf{L}}^{{-}1} + \boldsymbol{\mathsf{L}}^{{-}1} {\mathcal{N}}{}{\boldsymbol{\mathsf{R}}_\mathcal{S}}\boldsymbol{\mathsf{S}}^{{-}1} {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} \boldsymbol{\mathsf{L}}^{{-}1} & -\boldsymbol{\mathsf{L}}^{{-}1}{\mathcal{N}}{}{\boldsymbol{\mathsf{R}}_\mathcal{S}} \boldsymbol{\mathsf{S}}^{{-}1} \\ -\boldsymbol{\mathsf{S}}^{{-}1} {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}}\boldsymbol{\mathsf{L}}^{{-}1} & \boldsymbol{\mathsf{S}}^{{-}1}\end{bmatrix} \begin{pmatrix} -\boldsymbol{\mathsf{w}} \\ \boldsymbol{\mathsf{s}}'_{b} - \boldsymbol{\mathsf{s}}_{0} \end{pmatrix}. \end{equation}

The matrix operator in (2.30) is the inverse of the basic saddle-point system.

2.4.1. Using the Kutta condition in a steady-state problem

We will first take the steady-state approach to enforce the Kutta condition: allow the bound circulation to be set appropriately, with the implicit understanding that there is a starting vortex of equal and opposite circulation at infinity that preserves the Kelvin circulation theorem. The Lagrange multiplier for this constraint will not be $\varGamma _{b}$, but $s_{0}$, similar to (2.33). We also use the circulation removal operator to adjust the imposed surface streamfunction

(2.43)\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{L}} & {\mathcal{N}}{}{\tilde{\boldsymbol{\mathsf{R}}}_\mathcal{S}} & 0 \\ {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} & 0 & \boldsymbol{\mathsf{1}}_{\mathcal{S}} \\ 0 & \boldsymbol{\mathsf{e}}_{k}^{T} & 0 \end{bmatrix} \begin{pmatrix} \boldsymbol{\mathsf{s}} \\ \tilde{\boldsymbol{\mathsf{f}}} \\ s_{0} \end{pmatrix} = \begin{pmatrix} -\boldsymbol{\mathsf{w}} \\ \boldsymbol{\mathsf{P}}^{\varGamma}\boldsymbol{\mathsf{s}}'_{b} \\ 0 \end{pmatrix}. \end{equation}

This block system, like the earlier one in (2.25), has a saddle-point form, and we can reduce it by the same block-LU decomposition to develop a solution algorithm. We will interpret it in the general form (A3), with the upper-left $2\times 2$ block taking the role of $\mathcal {A}$, the solution vector $\boldsymbol{\mathsf{x}}$ and constraint force $\boldsymbol{\mathsf{y}}$ set, respectively, to

(2.44a,b)\begin{equation} \boldsymbol{\mathsf{x}} = \begin{pmatrix} \boldsymbol{\mathsf{s}} \\ \, \tilde{\boldsymbol{\mathsf{f}}} \end{pmatrix},\quad \boldsymbol{\mathsf{y}} = s_{0}, \end{equation}

the remaining operators set to

(2.45a–c)\begin{equation} \mathcal{B}_{2} = \big[0 \quad \boldsymbol{\mathsf{e}}_{k}^{T}\big], \quad \mathcal{B}_{1}^{T} = \begin{bmatrix} 0 \\ \boldsymbol{\mathsf{1}}_{\mathcal{S}} \end{bmatrix},\quad \mathcal{C} = 0, \end{equation}

and the right-hand side vectors set to

(2.46a,b)\begin{equation} \boldsymbol{\mathsf{r}}_{\boldsymbol{\mathsf{1}}_{\mathcal{N}}} = \begin{pmatrix} -\boldsymbol{\mathsf{w}} \\ \boldsymbol{\mathsf{P}}^{\varGamma}\boldsymbol{\mathsf{s}}'_{b} \end{pmatrix},\quad \boldsymbol{\mathsf{r}}_{\boldsymbol{\mathsf{2}}} = 0. \end{equation}

We note that block $\mathcal {A}$ has the original form of the system before the Kutta constraint (2.25), although with the slight modification of a re-scaled regularization operator, and we already have the inverse of $\mathcal {A}$ available from (2.30). The solution of this original system forms the intermediate solution of the full system endowed with the Kutta condition

(2.47a,b)\begin{equation} \tilde{\boldsymbol{\mathsf{f}}}^{*} = \tilde{\boldsymbol{\mathsf{S}}}^{{-}1} \left(\boldsymbol{\mathsf{P}}^{\varGamma}\boldsymbol{\mathsf{s}}'_{b} + {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}}\boldsymbol{\mathsf{L}}^{{-}1}\boldsymbol{\mathsf{w}} \right),\quad \boldsymbol{\mathsf{s}}^{*} ={-}\boldsymbol{\mathsf{L}}^{{-}1}\left(\boldsymbol{\mathsf{w}} + {\mathcal{N}}{}{\tilde{\boldsymbol{\mathsf{R}}}_\mathcal{S}} \tilde{\boldsymbol{\mathsf{f}}}^{*}\right). \end{equation}

Then, using the general procedure outlined in Appendix A, the solution of the full system (2.43) is easy to develop; its Schur complement is simply

(2.48)\begin{equation} \mathcal{S} ={-}1. \end{equation}

Applying the general solution equations, and using the property (2.41a,b) to simplify the resulting operators, it can be shown that the solution is

(2.49)\begin{gather} s_{0} = \boldsymbol{\mathsf{e}}_{k}^{T} \tilde{\boldsymbol{\mathsf{f}}}^{*} , \end{gather}

(2.50)\begin{gather}\begin{pmatrix} \boldsymbol{\mathsf{s}} \\ \,\tilde{\boldsymbol{\mathsf{f}}} \end{pmatrix} = \begin{pmatrix} \boldsymbol{\mathsf{s}}^{*} \\ \,\tilde{\boldsymbol{\mathsf{f}}}^{*} \end{pmatrix} - \begin{bmatrix} -\boldsymbol{\mathsf{L}}^{{-}1}{\mathcal{N}}{}{\tilde{\boldsymbol{\mathsf{R}}}_\mathcal{S}} \boldsymbol{\mathsf{1}}_{\mathcal{S}}\\ \boldsymbol{\mathsf{1}}_{\mathcal{S}}\end{bmatrix} \boldsymbol{\mathsf{e}}_{k}^{T} \tilde{\boldsymbol{\mathsf{f}}}^{*}. \end{gather}

The entire solution can be written more compactly as

(2.51)\begin{gather} \tilde{\boldsymbol{\mathsf{f}}} = \boldsymbol{\mathsf{P}}^{K}_{k} \tilde{\boldsymbol{\mathsf{S}}}^{{-}1} \left( \boldsymbol{\mathsf{P}}^{\varGamma}\boldsymbol{\mathsf{s}}'_{b} + {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}}\boldsymbol{\mathsf{L}}^{{-}1} \boldsymbol{\mathsf{w}} \right), \end{gather}

(2.52)\begin{gather}\boldsymbol{\mathsf{s}} ={-}\boldsymbol{\mathsf{L}}^{{-}1}\left(\boldsymbol{\mathsf{w}} + {\mathcal{N}}{}{\tilde{\boldsymbol{\mathsf{R}}}_\mathcal{S}} \tilde{\boldsymbol{\mathsf{f}}} \right), \end{gather}

where we have defined the Kutta projection operator,

(2.53)\begin{equation} \boldsymbol{\mathsf{P}}^{K}_{k} \equiv \boldsymbol{\mathsf{I}} - \boldsymbol{\mathsf{1}}_{\mathcal{S}} \boldsymbol{\mathsf{e}}_{k}^{T}, \end{equation}

which acts upon the (non-singular part of the) bound vortex sheet vector, subtracting the value at point $k$ from every point, including at $k$ itself.

Note that the Lagrange multiplier for the Kutta condition takes the simple value given by (2.49), revealing that the additional streamfunction on the surface is exactly the value of the intermediate bound vortex sheet at the Kutta point $k$.

The application of the Kutta condition to a steady-state problem is demonstrated in figure 7 on the flat plate problem that was introduced in the previous section. By constraining the trailing edge point of $\tilde {\boldsymbol{\mathsf{f}}}$, its whole distribution is shifted upward such that the trailing edge value equals zero. The resulting streamfunction indicates that the flow then indeed leaves the trailing edge smoothly.

Figure 7.

(a) Contours (——) of the steady, discrete streamfunction for a flat plate of length $c$ at $30^{\circ }$ in a uniform flow with enforcement of the Kutta condition at the trailing edge. (b) Zoomed-in image of the leading edge showing the streamlines for three values of the streamfunction, (c) the non-singular part of the associated discrete vortex sheet strength for $\Delta x/c=0.01$ with an inset enlarging the boxed area (— — —) and (d) the variation of its error with grid spacing for different values for $\Delta S/\Delta x$: $\cdot \cdot \cdot \cdot \cdot \,\cdot$, and $\diamond$, $\Delta S/\Delta x=1$; – – – (broken dashed blue) and $\square$ (blue), $\Delta S/\Delta x=2$; - - - - (green dashed) and $\circ$ (green), $\Delta S/\Delta x=3$; —— (red solid line) and $\triangle$ (red), $\Delta S/\Delta x=4$. Overlaid is an error (——) that scales as $\Delta x$.

Figure 7 also shows a spatial grid refinement analysis of the non-singular part of the vortex sheet strength. For this analysis, we can use a multiplicative decomposition for the continuous vortex sheet strength, $\gamma (\xi ) = \gamma _0(\xi ) \tilde {\gamma }(\xi )$, analogous to (2.38). The continuous counterpart of $\boldsymbol{\mathsf{f}}_{0}$ can be found as the bound vortex sheet strength of a circular cylinder with a point vortex of strength $\varGamma _{0}$ at its centre, conformally mapped to a flat plate. Due to the definition of $\boldsymbol{\mathsf{f}}_{0}$, the circulation $\varGamma _{0}$ depends on the grid spacing. Therefore, if we define $\gamma _0$ as the continuous counterpart of $\boldsymbol{\mathsf{f}}_{0}/\varGamma _{0}$, we obtain the grid-independent solution $\tilde {\gamma }(\xi ) = - {\rm \pi}c U_\infty \sin (\alpha ) (1-\xi )$ as the continuous counterpart of $\varGamma _{0}\,\tilde {\boldsymbol{\mathsf{f}}}$ for a flat plate with the Kutta condition enforced. We can then define the error for the non-singular part of the vortex sheet strength as $\epsilon _{\varGamma _{0}\,\tilde {\boldsymbol{\mathsf{f}}}} = \Vert \varGamma _{0}\,\tilde {\boldsymbol{\mathsf{f}}} - \tilde {\gamma }(\xi (\boldsymbol{r}_{x},\boldsymbol{r}_{y})\Vert _2 / \Vert \tilde {\gamma }(\xi (\boldsymbol{r}_{x},\boldsymbol{r}_{y})\Vert _2$. The refinement analysis verifies that the immersed boundary projection method using the vortex sheet strength decomposition is approximately first-order accurate in $\Delta x$. Furthermore, the figure shows the streamlines near the plate, the vortex sheet strength and the error for four different values of $\Delta S/\Delta x$. These confirm that a lower value increases the noise in the vortex sheet strength, but a higher value can lead to the streamlines penetrating the surface, as discussed previously.

2.4.2. Using the Kutta condition to set a new vortex element

In the previous section, we used the Kutta condition to set the bound circulation but did not explicitly create a new vortex element. This vortex element was assumed to lie at infinity so that its effect was negligible except insofar as it left equal but opposite circulation about the body.

In this section, we will create a new vortex element in the vicinity of the edge at which we are applying the Kutta condition. We will thus seek to establish the strength of this new element and to do so in such a manner that the overall circulation of the flow is conserved. Once the element is created, it will be allowed to advect with the local fluid velocity.

Let us assume that the new vortex element (which we label with the subscript 1) is introduced at some point in physical space, and that its immersion into the Cartesian grid is described by $\boldsymbol{\mathsf{d}}_{1}$ and that its strength (i.e. its circulation) is $\varGamma _{v,1}$. Thus, the fluid vorticity after this new element's introduction can be written as

(2.54)\begin{equation} \boldsymbol{\mathsf{w}} + \varGamma_{v,1}\boldsymbol{\mathsf{d}}_{1}. \end{equation}

The Kutta condition (2.42) is still to be enforced. We also seek to ensure that the total circulation is zero to satisfy Kelvin's circulation theorem. (We are assuming that the flow has started from rest.) Let us denote the circulation of the existing fluid vorticity $\boldsymbol{\mathsf{w}}$ by

(2.55)\begin{equation} \varGamma_{\boldsymbol{\mathsf{w}}} = {\langle \boldsymbol{\mathsf{1}}_{\mathcal{N}},\boldsymbol{\mathsf{w}} \rangle}_{\mathcal{N}}. \end{equation}

Then, the circulation constraint is

(2.56)\begin{equation} {\langle \boldsymbol{\mathsf{1}}_{\mathcal{S}},\boldsymbol{\mathsf{f}} \rangle}_{\mathcal{S}^{N_s}} + \varGamma_{v,1} + \varGamma_{\boldsymbol{\mathsf{w}}} = 0. \end{equation}

The circulation of the bound vortex sheet $\boldsymbol{\mathsf{f}}$ can be re-written in terms of the non-singular part of the sheet as ${\langle \boldsymbol{\mathsf{1}}_{\mathcal{S}},\boldsymbol{\mathsf{f}} \rangle }_{\mathcal {S}^{N_s}} = {\langle\, \boldsymbol{\mathsf{f}}_{0},\tilde {\boldsymbol{\mathsf{f}}} \rangle }_{\mathcal {S}^{N_s}} = \,\hat {{f}_{0}}^{T} \tilde {\boldsymbol{\mathsf{f}}}$.

With these two constraints, the overall saddle-point system of equations is

(2.57)\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{L}} & {\mathcal{N}}{}{\tilde{\boldsymbol{\mathsf{R}}}_\mathcal{S}} & 0 & \boldsymbol{\mathsf{d}}_{1} \\ {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} & 0 & \boldsymbol{\mathsf{1}}_{\mathcal{S}} & 0\\ 0 & \boldsymbol{\mathsf{e}}_{k}^{T} & 0 & 0 \\ 0 & \hat{{f}_{0}}^{T} & 0 & 1 \end{bmatrix} \begin{pmatrix} \boldsymbol{\mathsf{s}} \\ \tilde{\boldsymbol{\mathsf{f}}} \\ s_{0} \\ \varGamma_{v,1} \end{pmatrix} = \begin{pmatrix} -\boldsymbol{\mathsf{w}} \\ \boldsymbol{\mathsf{P}}^{\varGamma}\boldsymbol{\mathsf{s}}'_{b} \\ 0 \\ -\varGamma_{\boldsymbol{\mathsf{w}}} \end{pmatrix}. \end{equation}

Again, the basic saddle-point matrix constitutes the upper left $2\times 2$ block $\mathcal {A}$ and the solution vector $\boldsymbol{x}_{\mathcal{N}}$ is as before. The constraint force vector is

(2.58)\begin{equation} \boldsymbol{\mathsf{y}} = \begin{pmatrix} s_{0} \\ \varGamma_{v,1}\end{pmatrix}, \end{equation}

and the remaining vectors and operators are now

(2.59a–d)\begin{equation} \boldsymbol{\mathsf{r}}_{\boldsymbol{\mathsf{2}}} = \begin{pmatrix} 0 \\ -\varGamma_{\boldsymbol{\mathsf{w}}} \end{pmatrix},\quad \mathcal{B}_{2} = \begin{bmatrix} 0 & \boldsymbol{\mathsf{e}}_{k}^{T} \\ 0 & \hat{{f}_{0}}^{T} \end{bmatrix}, \quad \mathcal{B}_{1}^{T} = \begin{bmatrix} 0 & \boldsymbol{\mathsf{d}}_{1} \\ \boldsymbol{\mathsf{1}}_{\mathcal{S}} & 0\end{bmatrix},\quad \mathcal{C} ={-}\begin{bmatrix} 0 & 0 \\ 0 & 1\end{bmatrix}. \end{equation}

The solution algorithm follows, once again, from the equations in Appendix A. After carrying out the block matrix multiplications, it can be shown that the Schur complement (A2) is the $2\times 2$ matrix

(2.60)\begin{equation} \mathcal{S} = \begin{bmatrix} -1 & \boldsymbol{\mathsf{e}}_{k}^{T} \tilde{\boldsymbol{f}}_{1} \\ -\hat{{f}_{0}}^{T} \boldsymbol{\mathsf{1}}_{\mathcal{S}} & 1+\hat{{f}_{0}}^{T} \tilde{\boldsymbol{f}}_{1}, \end{bmatrix} \end{equation}

where, for convenience, we have defined

(2.61)\begin{equation} \tilde{\boldsymbol{f}}_{1} = \tilde{\boldsymbol{\mathsf{S}}}^{{-}1} {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} \boldsymbol{\mathsf{L}}^{{-}1} \boldsymbol{\mathsf{d}}_{1}, \end{equation}

which represents the (non-singular part of the) strength of the vortex sheet that ‘reacts’ to the presence of a unit-strength vortex $\boldsymbol{\mathsf{d}}_{1}$ immersed into the grid, cancelling that vortex's induced velocity on the surface. The term $\hat {{f}_{0}}^{T} \tilde {\boldsymbol{f}}_{1}$ represents this sheet's bound circulation and $\boldsymbol{\mathsf{e}}_{k}^{T} \tilde {\boldsymbol{f}}_{1}$ is its contribution to the Kutta condition at point $k$. The problem (A6) for the constraint forces $s_{0}$ and $\varGamma _{v,1}$ is then

(2.62)\begin{equation} \begin{bmatrix} -1 & \boldsymbol{\mathsf{e}}_{k}^{T} \tilde{\boldsymbol{f}}_{1} \\ -\hat{{f}_{0}}^{T} \boldsymbol{\mathsf{1}}_{\mathcal{S}} & 1+\hat{{f}_{0}}^{T} \tilde{\boldsymbol{f}}_{1}, \end{bmatrix} \begin{pmatrix} s_{0} \\ \varGamma_{v,1}\end{pmatrix} = \begin{pmatrix} -\boldsymbol{\mathsf{e}}_{k}^{T} \tilde{\boldsymbol{\mathsf{f}}}^{*} \\ -\varGamma_{\boldsymbol{\mathsf{w}}} - \hat{{f}_{0}}^{T} \tilde{\boldsymbol{\mathsf{f}}}^{*} \end{pmatrix}, \end{equation}

where the intermediate solution $\tilde {\boldsymbol{\mathsf{f}}}^{*}$ is available from (2.47a,b).

The determinant of this Schur complement matrix is $-1 - \hat {{f}_{0}}^{T} \boldsymbol{\mathsf{P}}^{K}_{k} \tilde {\boldsymbol{f}}_{1}$, which represents the negative of the circulation of the unit vortex and its associated vortex sheet, after the Kutta condition has been enforced on this sheet. It is straightforward then to calculate the strength of the new vortex $\varGamma _{v,1}$ and the additional uniform surface streamfunction, $s_{0}$

(2.63a,b)\begin{equation} \varGamma_{v,1} = \frac{-\varGamma_{\boldsymbol{\mathsf{w}}} - \hat{{f}_{0}}^{T} \boldsymbol{\mathsf{P}}^{K}_{k} \tilde{\boldsymbol{\mathsf{f}}}^{*}}{1 + \hat{{f}_{0}}^{T} \boldsymbol{\mathsf{P}}^{K}_{k} \tilde{\boldsymbol{f}}_{1}},\quad s_{0} = \boldsymbol{\mathsf{e}}_{k}^{T} (\tilde{\boldsymbol{\mathsf{\,f}}}^{*} + \varGamma_{v,1}\tilde{\boldsymbol{\,f}}_{1}). \end{equation}

From these, we can then obtain the vortex sheet strength and the fluid streamfunction,

(2.64a,b)\begin{equation} \tilde{\boldsymbol{\mathsf{f}}} = \boldsymbol{\mathsf{P}}^{K}_{k} (\,\tilde{\boldsymbol{\mathsf{f}}}^{*} + \varGamma_{v,1} \tilde{\boldsymbol{\,f}}_{1} ),\quad \boldsymbol{\mathsf{s}} ={-}\boldsymbol{\mathsf{L}}^{{-}1}(\boldsymbol{\mathsf{w}} +\varGamma_{v,1}\boldsymbol{\mathsf{d}}_{1} + {\mathcal{N}}{}{\tilde{\boldsymbol{\mathsf{R}}}_\mathcal{S}} \tilde{\boldsymbol{\mathsf{f}}}). \end{equation}

The intermediate solution, which corresponds to the flow associated with existing vorticity, is corrected here with the new vortex $1$ to enforce the Kutta condition at the point $k$.

We now first demonstrate the enforcement of the Kutta condition in unsteady flow on the flat plate problem with a point vortex near the trailing edge to enforce the Kutta condition at that edge. We position the point vortex at a distance $10 \Delta t U_{\infty }$ from the edge in the direction of the free stream, perpendicular to the plate. Figure 8 demonstrates that, because of the proximity of the point vortex to the flat plate, $\tilde {\boldsymbol{\mathsf{f}}}$ exhibits a quick variation at the surface points that lie closest to the point vortex. The value at the trailing edge point itself is still constrained to zero and the flow again leaves the edge smoothly. This situation corresponds to the flow right after impulsively starting a uniform flow around a flat plate and the point vortex now represents the starting vortex. The second demonstration of the method (figure 9) is the repeated enforcement of the Kutta condition on a NACA0012 airfoil while advancing the positions of the point vortices in time after impulsively starting the flow. This simulation, and all following simulations, use forward Euler time stepping, unless noted otherwise, and new point vortices are inserted at one third of the way from the edge to the last released vortex from that edge. The streamlines in the figure show that a strong initial vortex and weaker subsequent vortices were created and convected downstream. In turn, at each time step, the airfoil obtained a circulation that enforces the Kutta condition.

Figure 8.

(a) Contours (——) of the unsteady, discrete streamfunction for a flat plate of length $c$ at $30^{\circ }$ in a uniform flow with release of vorticity into a point vortex ($\bullet$) for enforcement of the Kutta condition at the trailing edge. (b) The non-singular part of the associated discrete vortex sheet strength. The simulation is performed with $\Delta x/c=0.01$ and $\Delta S/\Delta x=2$.

Figure 9.

Contours (——) of the unsteady, discrete streamfunction for a NACA0012 airfoil with chord length $c$ at $20^{\circ }$, three convective times after impulsively starting a uniform flow $U_{\infty }$. Vorticity is released into a stream of point vortices ($\bullet$) for enforcement of the Kutta condition at the trailing edge. The simulation is performed with $\Delta x/c=0.01$, $\Delta S/\Delta x=1.5$ and $\Delta t U_{\infty }/c=0.075$.

2.4.3. Applying more than one Kutta condition on a body

Suppose we wish to enforce the Kutta condition at two edges of the body – at points $k_{1}$ and $k_{2}$ – instead of one. Each such point has a constraint

(2.65)\begin{equation} \boldsymbol{\mathsf{e}}_{k_{j}}^{T}\, \tilde{\boldsymbol{\mathsf{f}}} = 0,\quad j = 1,2. \end{equation}

For two such constraints, we need two Lagrange multipliers: the strengths of two new vortices, $\varGamma _{v,1}$ and $\varGamma _{v,2}$, immersed into the grid with $\boldsymbol{\mathsf{d}}_{1}$ and $\boldsymbol{\mathsf{d}}_{2}$, respectively; and we still need the Lagrange multiplier $s_{0}$ to ensure that Kelvin's circulation theorem is also enforced. The system in the previous section is thus easily generalized to the following:

(2.66)\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{L}} & {\mathcal{N}}{}{\tilde{\boldsymbol{\mathsf{R}}}_\mathcal{S}} & 0 & \boldsymbol{\mathsf{d}}_{1} & \boldsymbol{\mathsf{d}}_{2} \\ {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} & 0 & \boldsymbol{\mathsf{1}}_{\mathcal{S}} & 0 & 0\\ 0 & \boldsymbol{\mathsf{e}}_{k_{1}}^{T} & 0 & 0 & 0 \\ 0 & \boldsymbol{\mathsf{e}}_{k_{2}}^{T} & 0 & 0 & 0 \\ 0 & \hat{{f}_{0}}^{T} & 0 & 1 & 1 \end{bmatrix} \begin{pmatrix} \boldsymbol{\mathsf{s}} \\ \tilde{\boldsymbol{\mathsf{f}}} \\ s_{0} \\ \varGamma_{v,1} \\ \varGamma_{v,2} \end{pmatrix} = \begin{pmatrix} -\boldsymbol{\mathsf{w}} \\ \boldsymbol{\mathsf{P}}^{\varGamma}\boldsymbol{\mathsf{s}}'_{b} \\ 0 \\ 0 \\ -\varGamma_{\boldsymbol{\mathsf{w}}} \end{pmatrix}. \end{equation}

The system is reduced in the same manner as before, with the same intermediate solution obtained from the basic system (2.25). Now, the Schur complement problem for the constraint forces takes the form

(2.67)\begin{equation} \begin{bmatrix} -1 & \boldsymbol{\mathsf{e}}_{k_{1}}^{T} \tilde{\boldsymbol{f}}_{1} & \boldsymbol{\mathsf{e}}_{k_{1}}^{T} \tilde{\boldsymbol{f}}_{2} \\ -1 & \boldsymbol{\mathsf{e}}_{k_{2}}^{T} \,\tilde{\boldsymbol{f}}_{1} & \boldsymbol{\mathsf{e}}_{k_{2}}^{T} \,\tilde{\boldsymbol{f}}_{2} \\ -\hat{{f}_{0}}^{T} \boldsymbol{\mathsf{1}}_{\mathcal{S}} & 1+\hat{{f}_{0}}^{T} \tilde{\boldsymbol{f}}_{1} & 1+\hat{{f}_{0}}^{T} \tilde{\boldsymbol{f}}_{2}, \end{bmatrix} \begin{pmatrix} s_{0} \\ \varGamma_{v,1}\\ \varGamma_{v,2}\end{pmatrix} = \begin{pmatrix} -\boldsymbol{\mathsf{e}}_{k_{1}}^{T}\, \tilde{\boldsymbol{\mathsf{f}}}^{*} \\ -\boldsymbol{\mathsf{e}}_{k_{2}}^{T}\, \tilde{\boldsymbol{\mathsf{f}}}^{*} \\ -\varGamma_{\boldsymbol{\mathsf{w}}} - \hat{{f}_{0}}^{T} \tilde{\boldsymbol{\mathsf{f}}}^{*} \end{pmatrix}, \end{equation}

where we have now defined bound vortex sheets associated with each of the two new vortices (with unit strengths)

(2.68)\begin{equation} \tilde{\boldsymbol{f}}_{j} = \tilde{\boldsymbol{\mathsf{S}}}^{{-}1} {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{N}} \boldsymbol{\mathsf{L}}^{{-}1} \boldsymbol{\mathsf{d}}_{j}, \end{equation}

for $j = 1,2$. It is interesting to note that, if we take the difference between the two Kutta constraints, we obtain

(2.69)\begin{equation} \left(\boldsymbol{\mathsf{e}}_{k_{1}}^{T}- \boldsymbol{\mathsf{e}}_{k_{2}}^{T}\right)(\,\tilde{\boldsymbol{\mathsf{f}}}^{*} + \varGamma_{v,1} \,\tilde{\boldsymbol{f}}_{1} + \varGamma_{v,2} \,\tilde{\boldsymbol{f}}_{2}) = 0. \end{equation}

It can be shown that this Schur complement problem can be split into

(2.70) \begin{equation} \begin{bmatrix} 1 + \,\hat{{f}_{0}}^{T}\boldsymbol{\mathsf{P}}^{K}_{k_{1}}\,\tilde{\boldsymbol{f}}_{1} & 1 + \,\hat{{f}_{0}}^{T}\boldsymbol{\mathsf{P}}^{K}_{k_{1}}\,\tilde{\boldsymbol{f}}_{2} \\ 1 + \,\hat{{f}_{0}}^{T}\boldsymbol{\mathsf{P}}^{K}_{k_{2}}\,\tilde{\boldsymbol{f}}_{1} & 1 + \,\hat{{f}_{0}}^{T}\boldsymbol{\mathsf{P}}^{K}_{k_{2}}\,\tilde{\boldsymbol{f}}_{2} \end{bmatrix} \begin{pmatrix} \varGamma_{v,1}\\ \varGamma_{v,2} \end{pmatrix} ={-}\begin{pmatrix} \varGamma_{\boldsymbol{\mathsf{w}}} + \,\hat{{f}_{0}}^{T} \boldsymbol{\mathsf{P}}^{K}_{k_{1}}\,\tilde{\boldsymbol{\mathsf{f}}}^{*}\\ \varGamma_{\boldsymbol{\mathsf{w}}} + \,\hat{{f}_{0}}^{T} \boldsymbol{\mathsf{P}}^{K}_{k_{2}}\,\tilde{\boldsymbol{\mathsf{f}}}^{*}\end{pmatrix} \end{equation}

and

(2.71) \begin{equation} s_{0} = \tfrac{1}{2} \left(\boldsymbol{\mathsf{e}}_{k_{1}}^{T}+\boldsymbol{\mathsf{e}}_{k_{2}}^{T}\right)(\,\tilde{\boldsymbol{\mathsf{f}}}^{*} + \varGamma_{v,1} \,\tilde{\boldsymbol{f}}_{1} + \varGamma_{v,2} \,\tilde{\boldsymbol{f}}_{2}). \end{equation}

The latter equation, when combined with (2.69), reveals that the value of the vortex sheet strength $\tilde {\boldsymbol{\mathsf{f}}}^{*} + \varGamma _{v,1}\,\tilde {\boldsymbol{f}}_{1} + \varGamma _{v,2}\,\tilde {\boldsymbol{f}}_{2}$ is the same at both Kutta points and equal to $s_{0}$.

Equation (2.70) can be solved easily for the strengths of the two new point vortices. Then, the solution for the vortex sheet strength and streamfunction are

(2.72a,b) \begin{align} \tilde{\boldsymbol{\mathsf{f}}} \!=\! \tfrac{1}{2}\left(\boldsymbol{\mathsf{P}}^{K}_{k_{1}} \!+\! \boldsymbol{\mathsf{P}}^{K}_{k_{2}}\right)(\,\tilde{\boldsymbol{\mathsf{f}}}^{*} \!+\! \varGamma_{v,1}\,\tilde{\boldsymbol{f}}_{1} \!+\! \varGamma_{v,2}\,\tilde{\boldsymbol{f}}_{2} ),\quad \boldsymbol{\mathsf{s}} ={-}\boldsymbol{\mathsf{L}}^{{-}1}(\boldsymbol{\mathsf{w}} +\varGamma_{v,1}\boldsymbol{\mathsf{d}}_{1} + \varGamma_{v,2}\boldsymbol{\mathsf{d}}_{2} + {\mathcal{N}}{}{\tilde{\boldsymbol{\mathsf{R}}}_\mathcal{S}} \,\tilde{\boldsymbol{\mathsf{f}}}). \end{align}

We now apply this method in figure 10 to enforce the Kutta condition at the leading and trailing edge of our flat plate problem. We position a point vortex close to each edge and observe again that $\tilde {\boldsymbol{\mathsf{f}}}$ shows strong variation at the surface points closest to the two point vortices. The contours of the streamfunction indicate that the flow indeed leaves the edges smoothly. Like the previous case, this solution corresponds to the flow right after impulsively starting a uniform flow around a flat plate, but unlike the previous case, the flow now separates at the leading edge.

Figure 10.

(a) Contours (——) of the unsteady, discrete streamfunction for a flat plate of length $c$ at $30^{\circ }$ in a uniform flow with release of vorticity into two point vortices ($\bullet$) for enforcement of the Kutta condition at both edges. (b) The non-singular part of the associated discrete vortex sheet strength. The simulation is performed with $\Delta x/c=0.01$ and $\Delta S/\Delta x=2$.

It should be observed that these solutions are posed in a manner easily extensible to an arbitrary number of edges.

2.6.1. Pressure

Here, we devise a means of computing the pressure, both throughout the flow field, $\boldsymbol{\mathsf{p}} \in \mathcal {C}$, as well as on the surface of a body. As we will show below, our immersed boundary treatment naturally gives rise to the jump in pressure across this surface, $\Delta \boldsymbol{\mathsf{p}} = \boldsymbol{\mathsf{p}}^{+} - \boldsymbol{\mathsf{p}}^{-} \in \mathcal {S}^{N_s}$ (with $+$ in the direction of the surface normal, $\boldsymbol{\mathsf{n}}_{} \in \mathcal {V}^{N_s}$). Thus, to distinguish the pressures on either side of the surface from one another, we use the fact that the interpolation of $\boldsymbol{\mathsf{p}}$ onto the surface produces the average of surface values, $ {\mathcal {S}}{}{\boldsymbol{\mathsf{E}}_\mathcal {C}} \boldsymbol{\mathsf{p}} = (\boldsymbol{\mathsf{p}}^{+} + \boldsymbol{\mathsf{p}}^{-})/2 \in \mathcal {S}$, where $ {\mathcal {S}}{}{\boldsymbol{\mathsf{E}}_\mathcal {C}}$ interpolates data from cell centres to the surface. It is thus easy to see that the pressure on either side is

(2.77)\begin{equation} \boldsymbol{\mathsf{p}}^{{\pm}} = {\mathcal{S}}{}{\boldsymbol{\mathsf{E}}_\mathcal{C}} \boldsymbol{\mathsf{p}} \pm \tfrac{1}{2} \Delta\boldsymbol{\mathsf{p}}. \end{equation}

Thus, we seek $\Delta \boldsymbol{\mathsf{p}}$ and $\boldsymbol{\mathsf{p}}$. It should not be a surprise that our starting point for these quantities is the Euler equations. However, our approach exploits the fact that the methodology we have presented thus far already solves the Euler equations in the fluid – in vorticity form, via transport of vortex elements – and provides us with the instantaneous velocity field $\boldsymbol{\mathsf{v}}$ and strength of the bound vortex sheet $\boldsymbol{\mathsf{f}}$ on any bodies. This approach, which satisfies the incompressibility constraint by expressing velocity in the null space of the divergence operator (i.e. as curl of a streamfunction), obviates the need for computing pressure, the Lagrange multiplier for this constraint. However, now that we seek pressure, we use the Euler equations in their velocity form to provide it.

It is important to note that the immersed boundary treatment enriches the Euler equations’ standard form with surface terms (Eldredge Reference Eldredge2022). Written in their spatially discrete form, these immersed boundary Euler equations are

(2.78)\begin{align} \rho \frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}t} + \rho (\boldsymbol{\mathsf{w}} + {\mathcal{N}}{}{\boldsymbol{\mathsf{R}}_\mathcal{S}} \boldsymbol{\mathsf{f}} ) \times \boldsymbol{v} &={-} \boldsymbol{\mathsf{G}} \left(\boldsymbol{\mathsf{p}} + \frac{1}{2}\rho|\boldsymbol{v}|^{2} \right) \nonumber\\ &\quad + {\mathcal{F}}{}{\boldsymbol{\mathsf{R}}_\mathcal{V}} [\boldsymbol{\mathsf{n}}_{} \circ \Delta\boldsymbol{\mathsf{p}} + \rho \Delta \boldsymbol{v} ({\mathcal{V}}{}{\boldsymbol{\mathsf{E}}_\mathcal{F}}\boldsymbol{v} - \boldsymbol{v}_{b} )\boldsymbol{\cdot}\boldsymbol{\mathsf{n}}_{}], \end{align}

where $\rho$ is the fluid density and $ {\mathcal {V}}{}{\boldsymbol{\mathsf{E}}_\mathcal {F}}$ interpolates grid data from $\mathcal {F}$ to the space of vector-valued surface data $\mathcal {V}^{N_s}$. Note the appearance of a few terms that are not typically seen in the Euler equations. First, we note the final bracketed pair of surface terms on the right-hand side, containing the jumps in pressure and velocity ($\Delta \boldsymbol{\mathsf{v}}$) across the surface, as well as the difference between the normal components of the fluid velocity (interpolated onto the surface) and the surface velocity itself, $\boldsymbol{\mathsf{v}}_{b}\in \mathcal {V}^{N_s}$. (This difference is zero by virtue of the no-flow-through condition, but we keep it here since it combines with other terms in later manipulations.) These surface terms are immersed into the grid by the operator $ {\mathcal {F}}{}{\boldsymbol{\mathsf{R}}_\mathcal {V}}$. In fact, if we had chosen to solve the Euler equations in velocity – instead of by streamfunction–vorticity – form, then we would have used the pressure jump in this term as a Lagrange multiplier for enforcing the no-flow-through condition on $\boldsymbol{\mathsf{v}}$. The second new term is the bound vortex sheet strength, appearing alongside the fluid vorticity on the left-hand side. This term emerges because the curl of the velocity field in the immersed boundary method generates both of these: $\boldsymbol{\mathsf{C}}^{T} \boldsymbol{\mathsf{v}} = \boldsymbol{\mathsf{w}} + {\mathcal {N}}{}{\boldsymbol{\mathsf{R}}_\mathcal {S}} \boldsymbol{\mathsf{f}}$, where $\boldsymbol{\mathsf{C}}^{T}: \mathcal {F} \mapsto \mathcal {N}$. (Our notation for this operator is consistent with that of previous authors, such as Colonius & Taira Reference Colonius and Taira2008.)

These quantities, $\Delta \boldsymbol{\mathsf{p}}$ and $\boldsymbol{\mathsf{p}}$, can be solved for simultaneously from the Schur complement system that arises from solving the Euler equations (2.78) and the associated constraints of divergence-free velocity and no-flow-through condition. Indeed, the approach we outline here is the natural outcome of that system. However, rather than present a detailed derivation, we present the equations with an intuitive explanation.

First, we develop an equation for $\Delta \boldsymbol{\mathsf{p}}$ by taking the discrete curl $\boldsymbol{\mathsf{C}}^{T}$ of (2.78) to eliminate the gradient term – since $\boldsymbol{\mathsf{C}}^{T} \boldsymbol{\mathsf{G}}$ is identically zero – and obtain a vorticity form of the immersed boundary Euler equations. With some manipulation to account for the motion of the surface (embodied in the time variation of $ {\mathcal {N}}{}{\boldsymbol{\mathsf{R}}_\mathcal {S}}$), these equations can be written as

(2.79)\begin{align} \rho \left( \frac{\mathrm{d}\boldsymbol{\mathsf{w}}}{\mathrm{d}t} - \boldsymbol{\mathsf{C}}^{T} (\boldsymbol{v}\times \boldsymbol{\mathsf{w}})\right) + \rho \,{\mathcal{N}}{}{\boldsymbol{\mathsf{R}}_\mathcal{S}} \frac{\mathrm{d}\boldsymbol{\mathsf{f}}}{\mathrm{d}t} = \boldsymbol{\mathsf{C}}^{T}\, {\mathcal{F}}{}{\boldsymbol{\mathsf{R}}_\mathcal{V}} \boldsymbol{\mathsf{n}}_{} \circ\left[ \Delta\boldsymbol{\mathsf{p}} + \rho ({\mathcal{V}}{}{\boldsymbol{\mathsf{E}}_\mathcal{F}}\boldsymbol{v}-\boldsymbol{v}_{b})\boldsymbol{\cdot}(\,\boldsymbol{\mathsf{f}}\times \boldsymbol{\mathsf{n}}_{}) \right], \end{align}

where $\mathrm {d}\boldsymbol{\mathsf{f}}/\mathrm {d}t$ represents the time derivative of each element of $\boldsymbol{\mathsf{f}}$ while following a point moving with velocity $\boldsymbol{\mathsf{v}}_{b}$. We assume that each immersed point moves with this local surface velocity, so this is simply the time derivative of the vector $\boldsymbol{\mathsf{f}}$. We have written the equation in this intermediate form on purpose in order to make a few key points. First, it is important to note that we have already satisfied the vorticity equation in the fluid – the first terms in parentheses on the left-hand side – by advecting the point vortices. Thus, we can set these terms to zero, leaving only those terms associated with the surface. Aside from $\Delta \boldsymbol{\mathsf{p}}$, these remaining terms involve only known quantities, and we could solve them in the current form for $\Delta \boldsymbol{\mathsf{p}}$ using similar techniques to the ones we will describe below. However, we will first write the equation in a more familiar form, and define some helpful operators and quantities to enable this.

In (2.79) we see a composite of the curl of the regularization operator; let us write this more compactly as $\boldsymbol{\mathsf{C}}_\mathcal {S}^{T}: \mathcal {S}^{N_s} \mapsto \mathcal {N}$,

(2.80)\begin{equation} \boldsymbol{\mathsf{C}}_\mathcal{S}^{T} \boldsymbol{\sigma} = \boldsymbol{\mathsf{C}}^{T} \,{\mathcal{F}}{}{\boldsymbol{\mathsf{R}}_\mathcal{V}} \left(\boldsymbol{\mathsf{n}}_{} \circ \boldsymbol{\sigma} \right), \end{equation}

for some surface scalar data, $\boldsymbol{\sigma } \in \mathcal {S}^{N_s}$. We refer to this as a surface curl operator. Its transpose, $\boldsymbol{\mathsf{C}}_\mathcal {S}: \mathcal {N} \mapsto \mathcal {S}^{N_s}$, also arises in what follows, and can also be described as a surface curl operator. It is defined as

(2.81)\begin{equation} \boldsymbol{\mathsf{C}}_\mathcal{S} \boldsymbol{\mathsf{s}} = \boldsymbol{\mathsf{n}}_{}\boldsymbol{\cdot} {\mathcal{V}}{}{\boldsymbol{\mathsf{E}}_\mathcal{F}} \boldsymbol{\mathsf{C}} \boldsymbol{\mathsf{s}}, \end{equation}

for $\boldsymbol{\mathsf{s}} \in \mathcal {N}$. The operator $\boldsymbol{\mathsf{C}}_\mathcal {S}$ obtains the normal component of velocity on an immersed surface for a given streamfunction $\boldsymbol{\mathsf{s}}$. Before we explain the role of its transpose, it is useful to remember that any potential flow generated by (or about) an impenetrable surface can be equivalently described by either a distribution of vortices (a vortex sheet, with strength $\boldsymbol{\mathsf{f}}$) or a distribution of dipoles (a double layer) on the surface. In the latter case, the strength of the double layer is given by the negative of the jump in scalar potential across the surface, $-\boldsymbol{\Delta} {\boldsymbol{\phi} } \in \mathcal {S}^{N_s}$. In fact, the two distributions’ strengths can be related to each other, either by using Stokes’ theorem or by the properties of the generalized functions that underpin the immersed boundary method (Eldredge Reference Eldredge2022), leading to

(2.82)\begin{equation} {\mathcal{N}}{}{\boldsymbol{\mathsf{R}}_\mathcal{S}} \boldsymbol{\mathsf{f}} ={-} \boldsymbol{\mathsf{C}}_\mathcal{S}^{T} \boldsymbol{\Delta {\phi}}. \end{equation}

Thus, $\boldsymbol{\mathsf{C}}_\mathcal {S}^{T}$ produces the equivalent bound vorticity distribution (immersed into the grid) associated with a given jump in scalar potential on the surface. To calculate $\boldsymbol{\Delta {\phi}}$ in terms of $\boldsymbol{\mathsf{f}}$, we apply $\boldsymbol{\mathsf{C}}_\mathcal {S} \boldsymbol{\mathsf{L}}^{-1}$ to both sides of (2.82) to equate the normal velocity induced on the surface by each distribution. The composite operator $-\boldsymbol{\mathsf{C}}_\mathcal {S}\boldsymbol{\mathsf{L}}^{-1} \boldsymbol{\mathsf{C}}_\mathcal {S}^{T}$ is positive semi-definite, and its null space can be shown to consist only of uniform values on the surface (i.e. the null space has an equivalent bound vorticity equal to zero). Thus, the jump in scalar potential associated with a vortex sheet of strength $\boldsymbol{\mathsf{f}}$ is

(2.83)\begin{equation} \boldsymbol{\Delta {\phi}} ={-}(\boldsymbol{\mathsf{C}}_\mathcal{S}\boldsymbol{\mathsf{L}}^{{-}1} \boldsymbol{\mathsf{C}}_\mathcal{S}^{T} )^{{-}1} \boldsymbol{\mathsf{C}}_\mathcal{S}\boldsymbol{\mathsf{L}}^{{-}1} {\mathcal{N}}{}{\boldsymbol{\mathsf{R}}_\mathcal{S}} \boldsymbol{\mathsf{f}}, \end{equation}

to which we can add any constant value without affecting the result.

Now, armed with this insight, we can return to (2.79), and rewrite the surface terms on the left-hand side in terms of $\boldsymbol{\Delta {\phi}}$

(2.84)\begin{equation} \boldsymbol{\mathsf{C}}_\mathcal{S}^{T} \left[ \Delta\boldsymbol{\mathsf{p}} + \rho ({\mathcal{V}}{}{\boldsymbol{\mathsf{E}}_\mathcal{F}}\boldsymbol{v}-\boldsymbol{v}_{b})\boldsymbol{\cdot}(\,\boldsymbol{\mathsf{f}}\times \boldsymbol{\mathsf{n}}_{}) + \rho \frac{\mathrm{d} \boldsymbol{\Delta {\phi}}}{\mathrm{d}t}\right] = 0, \end{equation}

where $\mathrm {d} \boldsymbol{\Delta {\phi}} /\mathrm {d}t$ denotes the time derivative of $\boldsymbol{\Delta {\phi}}$ associated with a particular immersed point, assumed to be moving with velocity $\boldsymbol{\mathsf{v}}_{b}$. This equation implies that the expression in brackets must be equal to a uniform value, which we can take to be zero without loss of generality. We can immediately write an immersed boundary Bernoulli equation

(2.85)\begin{equation} \Delta\boldsymbol{\mathsf{p}} + \rho ({\mathcal{V}}{}{\boldsymbol{\mathsf{E}}_\mathcal{F}}\boldsymbol{v}-\boldsymbol{v}_{b})\boldsymbol{\cdot}(\,\boldsymbol{\mathsf{f}}\times \boldsymbol{\mathsf{n}}_{}) + \rho \frac{\mathrm{d} \boldsymbol{\Delta {\phi}}}{\mathrm{d}t} = 0. \end{equation}

This equation is the discrete equivalent of a continuous version that appears in previous works, e.g. Jones (Reference Jones2003) and Eldredge (Reference Eldredge2019). At each time step, we use (2.84) to compute the instantaneous jump in scalar potential associated with the vortex sheet $\boldsymbol{\mathsf{f}}$, and then use (2.85) to find $\Delta \boldsymbol{\mathsf{p}}$.

Now we can substitute $\Delta \boldsymbol{\mathsf{p}}$ from (2.85) into the Euler equations (2.78) and operate on these equations with the discrete divergence operator, so that $\boldsymbol{\mathsf{L}} = \boldsymbol{\mathsf{D}}\boldsymbol{\mathsf{G}}$ acts on $\boldsymbol{\mathsf{p}}$. We solve the resulting equation, obtaining

(2.86)\begin{equation} \boldsymbol{\mathsf{p}} ={-}\frac{1}{2} \rho |\boldsymbol{v}|^{2} - \rho \boldsymbol{\mathsf{L}}^{{-}1} \left[ \boldsymbol{\mathsf{D}} (\boldsymbol{\mathsf{w}} \times \boldsymbol{v}) + {\mathcal{F}}{}{\boldsymbol{\mathsf{R}}_\mathcal{V}}\left( \boldsymbol{\mathsf{n}}_{} \circ \frac{\mathrm{d}\boldsymbol{\Delta {\phi}}}{\mathrm{d}t} - \boldsymbol{v}_{b}\times\boldsymbol{\mathsf{f}} \right) \right]. \end{equation}

It should be noted that any uniform value can be arbitrarily added to this expression. Also, we note in passing that the final set of terms (with the inverse Laplacian acting on the term in brackets) is equivalent to the time derivative of the scalar potential field, rendering the overall equation equivalent to a Bernoulli equation on the grid. However, there is no particular advantage in writing the equation in that form. In the current form, the first term in brackets (the divergence of the Lamb vector) represents the direct force exerted on the fluid by the fluid vorticity. The remaining terms in brackets collectively constitute the effects of surface motion and of the surface's modification of the flow induced by fluid elements (e.g. vorticity, free stream).

We now demonstrate the pressure calculations with two examples that have an analytical solution for the pressure distribution on the surface. Figures 12 and 13 depict the pressure field and surface distribution for a vortex near a cylinder and a plate, respectively. The pressure inside the cylinder is close to the exact solution of a uniform value, except for some noise at the side near the vortex, which is visible on the interior surface distribution of the pressure. The exterior surface pressure distribution on the cylinder shows good agreement with the analytical solution. In the example of the flat plate, shown in figure 13, the bottom and top surface pressure distributions again show good agreement with the analytical solution except for at the edges, which is expected. Note that, in these examples, we used a smaller value for $\Delta S/\Delta x$ than before, which was necessary to prevent the low pressure from the vortex from leaking through the surface. As a result, high-frequency components of the surface pressure are incorrectly amplified through the regularization and interpolation operators (Goza et al. Reference Goza, Liska, Morley and Colonius2016), especially at the edges.

Figure 12.

(a) Contours (——) of the discrete pressure for a point vortex ($\bullet$) with strength $\varGamma _{v,1}$ at $(R_v,0)$ near a cylinder consisting of $N_s$ points with radius $R_c$ and a bound circulation $-\varGamma _{v,1}$. (b) The scaled discrete pressure at the exterior (——) and interior ($\cdot \cdot \cdot \cdot \cdot \,\cdot$) of the cylinder. Overlaid is the exact continuous solution (- - - -, red dashed) for the exterior pressure. The simulation is performed with $R_v/R_c = 3/2$, $\Delta x/R_c=0.05$, $\Delta S/\Delta x=1.4$ and $\Delta t \varGamma _{v,1}/R_c^{2}=0.01$.

Figure 13.

(a) Contours (——) of the discrete pressure for a point vortex ($\bullet$) with strength $\varGamma _{v,1}$ at $(0,R_v)$ near a flat plate consisting of $N_s$ points with chord length $c$ and a bound circulation $-\varGamma _{v,1}$. The scaled discrete pressure (——) at (b) the top side and (c) the bottom side of the plate. Overlaid is the exact continuous solution (- - - -, red dashed). The simulation is performed with $R_v/c = 0.25$, $\Delta x/c=0.01$, $\Delta S/\Delta x=1.4$ and $\Delta t \varGamma _{v,1}/c^{2}=0.005$.

2.6.2. Impulse-based calculations of force and moment

The integral of the pressure distribution over the surface (plus any edge-suction parameters in the case of sharp edges) will be equal to the force on the surface. However, in this section we provide an alternative means of calculating the force and moment on the body through the negative rate of change of impulse in the fluid. The continuous expressions for linear and angular impulse (about the origin) are, in two dimensions (Saffman Reference Saffman1993; Eldredge Reference Eldredge2019)

(2.87)\begin{gather} \boldsymbol{P} = \int_{\varOmega} \boldsymbol{x}\times \boldsymbol{\omega}\,\mathrm{d}V + \int_{\mathcal{S}} \boldsymbol{x}\times \left( \boldsymbol{n}\times\boldsymbol{u}\right)\mathrm{d}S , \end{gather}

(2.88)\begin{gather}\boldsymbol{\varPi}_{{O}} = \tfrac{1}{2}\int_{\varOmega} \boldsymbol{x}\times \left(\boldsymbol{x}\times \boldsymbol{\omega}\right)\mathrm{d}V + \tfrac{1}{2} \int_{\mathcal{S}} \boldsymbol{x}\times \left[ \boldsymbol{x}\times\left( \boldsymbol{n}\times\boldsymbol{u}\right)\right]\mathrm{d}S . \end{gather}

If there is only a single body, then the force and moment (about the origin) exerted by the fluid on that body are given by

(2.89a,b)\begin{equation} \boldsymbol{F} ={-}\rho \frac{\mathrm{d}\boldsymbol{P}}{\mathrm{d}t}, \quad \boldsymbol{M}_{{O}} ={-}\rho\frac{\mathrm{d}\boldsymbol{\varPi_O}}{\mathrm{d}t}, \end{equation}

where $\rho$ is the fluid density. In the two-dimensional applications of this paper, the angular impulse and the moment have only a single component, e.g. $\boldsymbol {\varPi_{\mathrm {O}} =\varPi_{O}}\boldsymbol {e}_z$, where $\boldsymbol {e}_z$ is the unit vector out of the plane.

It should be observed that, by definition, the bound vortex sheet strength $\gamma$ is equal to the jump in tangential velocity between the fluid and the surface, $\boldsymbol {n}\times \boldsymbol {u}= \gamma \boldsymbol {e}_z + \boldsymbol {n}\times \boldsymbol {u}_{b}$, where $\boldsymbol {n}$ is the unit surface normal vector directed into the fluid, $\boldsymbol {u}$ is the fluid velocity, and $\boldsymbol {u}_{b}$ is the velocity of the surface. Thus, the surface integrals in (2.87) and (2.88) can be re-written in terms of the vortex sheet strength and the body motion.

We can easily develop discrete forms of the integrals (2.87) and (2.88) with the solutions and notation described in this paper. For the volume integrals, let us denote diagonal matrices containing the coordinates of the grid nodes by $\boldsymbol{\mathsf{D}}_{\boldsymbol{x}_{\mathcal{N}}}$ and $\boldsymbol{\mathsf{D}}_{\boldsymbol{y}_{\mathcal{N}}}$. Thus, the expressions in (2.87) and (2.88) can be written in discrete form as

(2.90)\begin{gather} P_{x} = {\langle \boldsymbol{y}_{\mathcal{N}},\boldsymbol{\mathsf{w}} \rangle}_{\mathcal{N}} + {\langle \boldsymbol{r}_{y},\boldsymbol{\mathsf{f}} + \boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{n}}_{x}} \boldsymbol{v}_{b,y} - \boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{n}}_{y}} \boldsymbol{v}_{b,x} \rangle}_{\mathcal{S}^{N_s}}, \end{gather}

(2.91)\begin{gather}P_{y} ={-}{\langle \boldsymbol{x}_{\mathcal{N}},\boldsymbol{\mathsf{w}} \rangle}_{\mathcal{N}} - {\langle \boldsymbol{r}_{x},\boldsymbol{\mathsf{f}} + \boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{n}}_{x}} \boldsymbol{v}_{b,y} - \boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{n}}_{y}} \boldsymbol{v}_{b,x} \rangle}_{\mathcal{S}^{N_s}}, \end{gather}

and

(2.92)\begin{align} \varPi_{{O}} ={-}\tfrac{1}{2} {\langle \boldsymbol{\mathsf{D}}_{\boldsymbol{x}_{\mathcal{N}}}\boldsymbol{x}_{\mathcal{N}} + \boldsymbol{\mathsf{D}}_{\boldsymbol{y}_{\mathcal{N}}}\boldsymbol{y}_{\mathcal{N}},\boldsymbol{\mathsf{w}} \rangle}_{\mathcal{N}} -\tfrac{1}{2} {\langle \boldsymbol{\mathsf{D}}_{\boldsymbol{r}_{x}}\boldsymbol{r}_{x} + \boldsymbol{\mathsf{D}}_{\boldsymbol{r}_{y}}\boldsymbol{r}_{y},\boldsymbol{\mathsf{f}} + \boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{n}}_{x}} \boldsymbol{v}_{b,y} - \boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{n}}_{y}} \boldsymbol{v}_{b,x} \rangle}_{\mathcal{S}^{N_s}}. \end{align}

The overall force and moment exerted on the body are obtained from calculating these impulses and computing their rates of change in (2.89a,b). Part of this force and moment is attributable to the dynamics of vorticity in the fluid. The remaining part is due to surface motion relative to the fluid, and we will discuss this in the next section.

To illustrate the accuracy of the impulse-based calculation of force, we apply the method to two examples. In the first example, we simulate the trajectories of two point vortices of opposite strength, in which case each vortex is convected past a cylinder due to the presence of the other vortex. The time stepping in this example is carried out using a fourth-order Runge–Kutta scheme. Figure 14 shows the trajectories and the $x$ component of the impulse together with their exact solutions, which show good agreement with the exact solution. In the second example, we compare our simulation of the first instants of the unsteady, fully separated flow around a flat plate after impulsively starting a uniform flow with the Biot–Savart method from Darakananda & Eldredge (Reference Darakananda and Eldredge2019), using the same positioning rules to insert point vortices and the same time step. The vortex positions and the corresponding impulse and lift are compared in figure 15 and show good agreement as well.

Figure 14.

(a) Numerically simulated trajectories (——) of two point vortices ($\bullet$) of opposite strengths $\varGamma _{v,1}$ and $\varGamma _{v,1}=-\varGamma_{v,}[2]$ being convected past a circular cylinder with radius $R$, and (b) the $x$ component of the associated, numerically simulated impulse (——) in the fluid. Overlaid are the exact continuous trajectories and impulse (- - - -, red dashed). The simulation is performed with $\Delta x/R=0.04$, $\Delta S/\Delta x=2$ and $\Delta t \varGamma _{v,1}/R^{2}=0.1$.

Figure 15.

Comparison of the simulated vortex shedding behind a flat plate of length $c$ at $60^{\circ }$ in a uniform flow using the method in this paper ($\bullet$ and ——) and using the Biot–Savart method of Darakananda & Eldredge (Reference Darakananda and Eldredge2019) ($\circ$, red open circle and - - - -, red dashed). (a) The positions of the shedded point vortices, (b) the lift coefficient and (c) the moment coefficient, one convective time after impulsively starting the uniform flow. At both edges, the Kutta condition is enforced. The simulation is performed with $\Delta x/c=0.01$, $\Delta S/\Delta x=2$ and $\Delta t U_{\infty }/c=0.05$.

2.6.3. Added mass

The added mass tensor provides a measure of the inertial influence of the fluid on the body in response to changes in the body's translational or rotational motion. The coefficients of the added mass tensor of a body are obtained by computing the impulse components associated with a unit-valued component of motion. The motion's influence is both direct, via the surface velocity, and indirect, in the bound vortex sheet that develops on the surface.

For example, suppose that we consider translation at unit velocity in the $x$ direction, for which the motion is described by $\boldsymbol{\mathsf{v}}_{b,x} = \boldsymbol{\mathsf{1}}_{\mathcal{S}}$, $\boldsymbol{\mathsf{v}}_{b,y} = 0$ and $\boldsymbol{\mathsf{s}}'_{b} = \boldsymbol{r}_{y}$, and the associated bound vortex sheet – obtained without the Kutta condition by solving the basic problem (2.25) – is $\boldsymbol{\mathsf{f}} = \boldsymbol{\mathsf{S}}^{-1}\boldsymbol{\mathsf{P}}^{\varGamma }\boldsymbol{\mathsf{s}}'_{b} = \boldsymbol{\mathsf{S}}^{-1}\boldsymbol{\mathsf{P}}^{\varGamma }\boldsymbol{r}_{y}$. The added mass coefficients corresponding to this motion are derived by substituting these into the impulse formulas (2.90)–(2.92)

(2.93)\begin{gather} P_{x}^{(x)} = {\langle \boldsymbol{r}_{y},\boldsymbol{\mathsf{S}}^{{-}1}\boldsymbol{\mathsf{P}}^{\varGamma}\boldsymbol{r}_{y} - \boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{n}}_{y}}\boldsymbol{\mathsf{1}}_{\mathcal{S}} \rangle}_{\mathcal{S}^{N_s}}, \end{gather}

(2.94)\begin{gather}P_{y}^{(x)} ={-}{\langle \boldsymbol{r}_{x},\boldsymbol{\mathsf{S}}^{{-}1}\boldsymbol{\mathsf{P}}^{\varGamma}\boldsymbol{r}_{y} - \boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{n}}_{y}}\boldsymbol{\mathsf{1}}_{\mathcal{S}} \rangle}_{\mathcal{S}^{N_s}}, \end{gather}

(2.95)\begin{gather}\varPi_{{O}}^{(x)} ={-}\tfrac{1}{2} {\langle \boldsymbol{\mathsf{D}}_{\boldsymbol{r}_{x}}\boldsymbol{r}_{x} + \boldsymbol{\mathsf{D}}_{\boldsymbol{r}_{y}}\boldsymbol{r}_{y},\boldsymbol{\mathsf{S}}^{{-}1}\boldsymbol{\mathsf{P}}^{\varGamma}\boldsymbol{r}_{y} - \boldsymbol{\mathsf{D}}_{\boldsymbol{\mathsf{n}}_{y}} \boldsymbol{\mathsf{1}}_{\mathcal{S}} \rangle}_{\mathcal{S}^{N_s}}. \end{gather}

Thus, the components of the added mass coefficients tensor associated with translation in the $x$ direction are

(2.96)\begin{gather} m^{F}_{xx} = \rho P_{x}^{(x)}, \end{gather}

(2.97)\begin{gather}m^{F}_{xy} = \rho P_{y}^{(x)}, \end{gather}

(2.98)\begin{gather}m^{M}_{x} = \rho \left(\varPi_{{O}}^{(x)}-\boldsymbol{X}_c \times \boldsymbol{P}^{(x)} \right), \end{gather}

where $\boldsymbol {X}_c$ is the centroid of the body, which can be calculated using (B8), and the superscripts $F$ and $M$ are used to denote the coefficient for the force and moment, respectively.

A similar approach can be used to obtain the added mass coefficients due to unit translation in the $y$ direction, for which $\boldsymbol{\mathsf{v}}_{b,x} = 0$, $\boldsymbol{\mathsf{v}}_{b,y} = \boldsymbol{\mathsf{1}}_{\mathcal{S}}$ and $\boldsymbol{\mathsf{s}}'_{b} = -\boldsymbol{r}_{x}$. The coefficients due to unit rotation follow from taking $\boldsymbol{\mathsf{v}}_{b,x} = -\boldsymbol{r}_{y}$, $\boldsymbol{\mathsf{v}}_{b,y} = \boldsymbol{r}_{x}$ and $\boldsymbol{\mathsf{s}}'_{b} = -\tfrac {1}{2}(\boldsymbol{\mathsf{D}}_{\boldsymbol{r}_{x}}\boldsymbol{r}_{x}+\boldsymbol{\mathsf{D}}_{\boldsymbol{r}_{y}}\boldsymbol{r}_{y})$.