Vortex cusps

TL;DR

This paper models vortex sheet cusps using a reduced Birkhoff-Rott equation, revealing conditions for their existence, and highlights the importance of accurate numerical schemes to avoid spurious results.

Contribution

It derives a reduced model for vortex cusp formation, providing formulas for cusp exponents and identifying conditions for their existence based on similarity exponent and strain.

Findings

Piecewise quadratic vortex sheet approximations agree with the model.

Linear schemes produce spurious results and violate mass conservation.

Vortex cusps only exist under certain conditions related to circulation and strain.

Abstract

We consider pairs of self-similar 2d vortex sheets forming cusps, equivalently single sheets merging into slip condition walls, as in classical Mach reflection at wedges. We derive from the Birkhoff-Rott equation a reduced model yielding formulas for cusp exponents and other quantities as functions of similarity exponent and strain coefficient. Comparison to numerics shows that piecewise quadratic and higher approximation of vortex sheets agree with each other and with the model. In contrast piecewise linear schemes produce spurious results and violate conservation of mass, a problem that may have been undetected in prior work for other vortical flows where even point vortices were sufficient. We find that vortex cusps only exist if the similarity exponent is sufficiently large and if the circulation on the sheet is counterclockwise (for a sheet above the wall with cusp opening to the right), unless a sufficiently positive strain coefficient compensates. Whenever a cusp cannot exist a spiral-ends jet forms instead; we find many jets are so narrow that they appear as false cusps.

Equations

$$\omega (t,\mathbf{x})=t^{-1}\omega (1,t^{-\mu }\mathbf{x})$$

$$\omega (0,\mathbf{x})=r^{-\kappa }\mathring{\omega }(\theta )$$

$$\mathbf{v}(t,\mathbf{x})=\mathbf{v}(t^{-1}\mathbf{x}).$$

$$v^x-i{v}^y=\frac{1}{2\pi i}\frac{d{\Gamma }^{\prime }}{z_0-z(t,{\Gamma }^{\prime })}.$$

$$w(z_0)=\text{p.v.}\int \frac{d{\Gamma }^{\prime }}{2\pi i(z_0-z(t,{\Gamma }^{\prime }))},$$

$$\frac{1}{2\pi i}\int \frac{-d{\Gamma }^{\prime }}{z_0-z^{\ast }(t,{\Gamma }^{\prime })},$$

$$z_t(t,\Gamma )=w^{\ast },w=\text{p.v.}\int \frac{1}{z(t,\Gamma )-z(t,{\Gamma }^{\prime })}-\frac{1}{z(t,\Gamma )-z^{\ast }(t,{\Gamma }^{\prime })}\frac{d{\Gamma }^{\prime }}{2\pi i}$$

$$z(t,\Gamma )=t^{\mu }z(1,t^{1-2\mu }\Gamma ),$$

$$0=w^{\ast }-\mu z+(2\mu -1)\Gamma z_{\Gamma },w=\text{p.v.}\int \frac{1}{z(\Gamma )-z({\Gamma }^{\prime })}-\frac{1}{z(\Gamma )-z^{\ast }({\Gamma }^{\prime })}\frac{d{\Gamma }^{\prime }}{2\pi i}.$$

$$0=(2\mu -1)\Gamma z_{\Gamma }-\mu z$$

$$\Gamma =\text{const}\cdot |z{|}^{2-1/\mu },\frac{z}{|z|}=\text{const}.$$

$$\frac{d{\Gamma }^{\prime }}{2\pi i(z_0-x^{\prime }-i{y}^{\prime })}$$

$$\frac{-d{\Gamma }^{\prime }}{2\pi i(z_0-x^{\prime }+i{y}^{\prime })};$$

$$dw$$

$$0$$

$$\approx \frac{1}{2}{\Gamma }_x-\mu x+(2\mu -1)\Gamma x_{\Gamma }.$$

$$e(z)=e_0+e_1(z-x_{\ast })+o(|z-x_{\ast }|)$$

$$\mathrm{Re}e=e_0+e_1\xi +...$$

$$\mathrm{Im}e(z)=e_{1}y+...$$

$$0=\frac{1}{2}{\Gamma }_x+\underset{\approx \mathrm{Re}e}{\underbrace{e_0+e_1\xi }}-\mu x+(2\mu -1)\Gamma x_{\Gamma }$$

$$\nabla \cdot \mathbf{q}=\nabla \cdot (\mathbf{v}-\mu \mathbf{x})=-2\mu$$

$${\int }_{-y(\xi )}^{y(\xi )}{q}^{x}ds=-2\mu \cdot \text{cusp triangle area}$$

$$0>{\underline{q}}^x\ge -\mu \xi ,$$

$$0$$

$${\underline{q}}^x$$

$$0=\frac{1}{2}({\underline{q}}^x{)}^2-\frac{1}{2}(e_0+e_1\xi -\mu x{)}^2+(2\mu -1)\Gamma .$$

$$(2\mu -1){\Gamma }_{\ast }=\frac{1}{2}(e_0-\mu x_{\ast }{)}^2\ge 0$$

$$e_0=\mu x_{\ast },$$

$${\underline{q}}^x$$

$$0$$

Full text

1. Introduction

This article studies planar self-similar incompressible vortex sheets meeting in cusps (fig. 1 right, fig. 9). Self-similar flows, also called “pseudo-steady” or “quasi-steady”, have vorticity

[TABLE]

for some similarity exponent $\mu$. The flow is “similar” at all times $t>0$, with any spatial distance $L$ dilated by a factor $t^{\mu}$ and any velocity $\mathbf{v}$ scaled accordingly. Such flows arise from initial data of type

[TABLE]

where $r,\theta$ are polar coordinates. Since vorticity has dimensions of inverse time, the exponent $-\kappa$ of the initial singularity fixes $t^{-1}\sim L^{-\kappa}$ and thus determines the similarity exponent $\mu=1/\kappa$. To ensure velocity is locally integrable at the origin we only consider $\mu>\frac{1}{2}$. The case $\mu=1$ is commonly studied for compressible flow as well since it corresponds to

[TABLE]

If $\omega$ is single-signed, or if the other sign is present but small enough to be “dominated”, then strong rotation at the origin winds $\omega$ (in particular its $\omega=0$ level sets) into algebraic spirals of type $r\sim\theta^{-\mu}$ (fig. 2 center). For vortex sheets algebraic spiral flows have been modelled in [Kad31], [Rot56], [Ste56], [Bir62], [MW67], [Moo75]; recently [Ell16] gave a mathematically rigorous existence proof for a class of mixed-sign smooth $\mathring{\omega}$. (For logarithmic spirals with entirely different scaling see e.g. [Pra24], [Ale71], [Kan89], [EG19].)

It is natural to ask what happens in the borderline case where neither sign of vorticity dominates. In this article a pair of vortex sheets of equal strength but opposite sign is considered. It is shown that in some cases they may form a spiral-free cusp of type $y\sim\xi^{\alpha}$ (with $y$ distance to symmetry axis, $\xi$ distance to cusp along axis); otherwise a spiral-ends jet appears.

While “equal strength” may at first glance appear to be a less important borderline case, it is precisely the case of a self-similar vortex sheet merging into a slip-condition wall (fig. 3 center), since eliminating the wall by a reflection yields a mirror sheet of opposite vorticity. Vortex-wall merging is commonly observed, for example in Mach reflection (fig. 3; see [MW75], [Neu43], [Hor86], [BD92]; experimental observations e.g. [vD82] fig. 236, 237). [HVBDE03] and [VBDEH04] have begun a theoretical and numerical study on when the sheet merges smoothly rather than flipping backwards and forming a spiral-end jet near the wall.

In section 2 we discuss several heuristic ODE models derived from the Birkhoff-Rott equation. Section 2.2 approximates axis-parallel on-sheet velocity $v^{{x}}$ as a function of $\xi$ alone, deriving a formula $v^{{x}}\sim\xi^{\beta-1}$ (see (24)) with exponent $\beta$ depending only on $\mu$ and strain parameter $e_{1}$, a single real number representing the limit of $\partial v^{{x}}/\partial\xi$ outside the cusp. Section 2.3 derives a model for the cusp shape as ${y}\sim\xi^{\alpha}$, predicting cusp exponents $\alpha$ which are also functions of $\mu,e_{1}$ alone (see (31)).

In section 3 we test our model against numerical data. Section 3.1 describes the numerics used for calculating cusps and the new obstacles compared to previous sheet calculations. Section 3.2 compares piecewise linear reconstruction to quadratic and higher degree; while the latter produce matching results that agree with our preferred model, we find that linear reconstruction produce a cusp exponent matching oversimplified conservation-violating models. (This is in contrast to the case of vortex spirals for which even point vortex methods were successful [Pul78, Pul89].) Section 3.2.3 and 3.2.4 confirm that our $\beta,\alpha$ formulas agree with numerical observations, with correct dependence on $\mu$ and $e_{1}$. Section 3.2.5 shows that clockwise circulation (on the upper sheet for a cusp opening to the right) is sometimes possible if $e_{1}$ is sufficiently large, although generally counterclockwise circulation is observed, which matters because that is the circulation produced by Mach reflections (see [Neu43, HM98, Ser07, Ell18]). In cases where cusps do not appear we generally observe spiral-ends jets instead, as in section 3.2.6 where we find that even for parameters far from the jet-cusp transition the jets can be rather small, mimicking cusps, suggesting that jetting is frequently obscured by viscous or kinetic effects in numerics or experiments, or simply small enough to be overlooked.

We conclude that our model is correct, that vortex cusps exist and can be computed accurately using vortex segments of at least quadratic degree.

2. Modelling

2.1. Birkhoff-Rott equation

We are interested in a vortex sheet in the upper halfplane meeting the wall in a “cusp”. We may assume the sheet approaches the cusp point from the right since approach from the left reduces to this case by mirror reflection. Equivalently we may think of a pair of vortex sheets in the whole plane (fig. 4), mirror-symmetric across the horizontal axis, one above one below, with circulation of opposite sign. $\Gamma$ is the circulation on the upper sheet, $z(t,\Gamma)$ with $z=x+iy$ the location of the corresponding sheet point at time $t$. The infinitesimal segment $\Gamma^{\prime}$ to $\Gamma^{\prime}+d\Gamma^{\prime}$ of the upper sheet, regarded as an infinitesimal point vortex, induces in another point $z_{0}$ the complex velocity

[TABLE]

The entire upper sheet induces a complex velocity

[TABLE]

where “p.v.” indicates principal value, the arithmetic average of the velocity limits on each side, i.e. of $w(z_{\pm})$ as $z_{\pm}$ approach $z_{0}$ from each side of the sheet. The lower sheet induces

[TABLE]

where $-$ before $d\Gamma^{\prime}$ is due to lower-sheet circulation of opposite sign, whereas $z^{*}$ is the complex conjugate, the mirror image $x-iy$ of $z=x+iy$. The upper sheet evolves according to the Birkhoff-Rott equation

[TABLE]

where subscripts indicate partial derivatives. The integration domain is chosen based on the case at hand. The Birkhoff-Rott equation is equivalent to the incompressible Euler equations under mild assumptions (see e.g [LFLS07]).

We are interested in self-similar vortex sheets: since circulation $\Gamma$ has dimensions of length squared over time, it is scaled by a factor $t^{2\mu-1}$. We pass to dimensionless $z,{\Gamma}$ with the ansatz

[TABLE]

resulting in the self-similar Birkhoff-Rott equation

[TABLE]

$w^{*}-\mu z$ is the complex form of pseudo-velocity $\mathbf{q}=\mathbf{v}-\mu{\mathbf{x}}$. ${\Gamma}\mapsto z({\Gamma})$ parametrizes the upper sheet, so $z_{\Gamma}={x}_{\Gamma}+i{y}_{\Gamma}$ is a tangent. Hence by (1) the pseudo-velocity $w^{*}-\mu z=q^{{x}}+iq^{{y}}$ is everywhere tangent. This represents the pseudo-velocity on the sheet; the limits on each side of the sheet differ by adding or subtracting $\frac{1}{2}{\Gamma}_{s}\mathbf{s}$ where $\mathbf{s}$ is a sheet unit tangent, subscript $s$ indicates derivative with respect to arc length; hence the one-sided limits of pseudo-velocity are also tangential. $\mathbf{q}\cdot\mathbf{n}=0$ corresponds to $\mathbf{v}\cdot\mathbf{n}=\sigma$ for general unsteady vortex sheets, $\sigma$ being normal speed of the sheet.

When considering flows in the entire 2d plane we assume the velocity $w$ is bounded or at least $o(|z|)$ near infinity (as $|z|\rightarrow\infty$), so that velocity is uniquely determined by vorticity. If so, then the term $-\mu z$ in (1) dominates $w$ near infinity, leading to an approximation

[TABLE]

which has a simple solution

[TABLE]

Hence the vortex sheets become asymptotic to straight lines near infinity. Those asymptotes correspond precisely to the initial data (fig. 1 left) in original $t,x$ coordinates: $t\searrow 0$ there corresponds in the similarity plane to zooming out, $t\nearrow\infty$ to zooming into the origin.

2.2. Horizontal velocity approximation

2.2.1. Cusp contributions

We consider cusp solutions of (1). A key idea (see fig. 4 left) is that the velocity induced in a point $z_{0}$ by a not too close point $z^{\prime}$ is almost completely cancelled by the mirror image ${z^{\prime}}^{*}$. More precisely, an infinitesimal segment $d{\Gamma}^{\prime}$ at $z^{\prime}$ on the upper sheet induces

[TABLE]

whereas its mirror image induces

[TABLE]

the sum is

[TABLE]

${y}^{\prime}$ is very small near the cusp, so the sum is much smaller than the summands unless $z_{0}$ is very close to $z^{\prime}$.

We emphasize the importance of this cancellation: the velocity integral, an operator that is generally non-local and thus complicated, is in the cusp case dominated by its local part. Of course “local” nurtures hopes of “differential”, and indeed this convenient reduction will soon become apparent.

In addition, for $z_{0}$ on the upper sheet there is cancellation between the “near” parts to the left and to the right as well (see fig. 4 right): for positive ${\Gamma}_{x}$ the left part induces almost upward velocities and the right one almost downward ones. So the dominant contribution to velocity is from the nearby mirror image.

To calculate its effect the following argument is convenient: in near-cusp points ${x}+i{y}$ the sheet tangent is almost horizontal (fig. 5 left); hence we approximate the sheet and its mirror image by two straight lines with constant velocity jump ${\Gamma}_{x}({x})$ (fig. 5 right). Such an approximation has velocity $\overline{\mathbf{v}}=(0,0)$ above the upper sheet (and below the lower one), $\underline{\mathbf{v}}=({\Gamma}_{x},0)$ between the two sheets; on the upper sheet the principal value is the arithmetic average $\mathbf{v}=(\frac{1}{2}{\Gamma}_{x},0)$. For the real (horizontal) part of (1):

[TABLE]

In one case our approximations are appropriate not only near the cusp point but everywhere: let $2\phi_{\infty}$ be the angle the straight-line asymptotes of vortex sheet and mirror image enclose near infinity (fig. 1 left). If the infinity angle $\phi_{\infty}$ is small, then we may expect that the sheet tangents are nearly horizontal everywhere.

2.2.2. Non-cusp contributions; symmetry

Of course we have neglected that there are velocity contributions from the outer (non-cusp) part of the sheets, as well as those of other sheets, vortex patches or any other form of vorticity present in the 2d similarity plane. Since they are uniformly distant from the cusp we may model the velocity induced by them as a complex-analytic function $e=e(z)$. We use its Taylor expansion

[TABLE]

where ${x}_{*}$ is the cusp location which is real (horizontal axis); the coefficients $e_{0},e_{1}$ are likewise real since velocity is horizontal at the wall by the slip condition. In section 4.1 we argue that expansions beyond $e_{1}$ cannot always be justified.

The constant $e_{0}$ does not have a direct influence on the cusp shape; it merely corresponds to a constant shift in the similarity plane. The strain coefficient $e_{1}$ corresponds to a saddle flow $\mathbf{v}=(e_{1}({x}-{x}_{*}),-e_{1}{y})$; for $e_{1}>0$ this flow is expanding along the symmetry axis of the cusp (horizontal axis), compressing in the perpendicular direction.

The horizontal part of $e$ is

[TABLE]

where $\xi={x}-{x}_{*}$ is horizontal distance to the cusp; we neglect (at this point) the vertical part

[TABLE]

since ${y}$ decays rapidly in a cusp. Now (7) is replaced by the more accurate model

[TABLE]

In a few special cases we may make additional assumptions. E.g. for small infinity angle $\phi_{\infty}$, and without non-cusp flow features present, we may assume that $e$ becomes small (possibly in contrast to $\frac{1}{2}{\Gamma}_{x}$).

Another special case is $N$-fold symmetry (see fig. 1): instead of a single pair of equal-strength opposite-circulation sheets consider $N$ pairs, with equal angles $2\pi/N$ between the symmetry axes of adjacent pairs, and cusps in the origin (as suggested by numerical experiments). Then by symmetry the Taylor expansion of $e$ vanishes up to and including $z^{N-1}$. For $N=2$ that means $e_{0}=0$; for $N\geq 3$ we additionally have $e_{1}=0$. Symmetry is a convenient device for suppressing $e_{0},e_{1}$ so that other influences shaping the cusp, such as dependence on $\mu$ or $\phi_{\infty}$, can be studied in isolation.

Symmetry also helps avoiding problems at infinity: the integral in the self-similar Birkhoff-Rott equation (1) behaves like $\int|z|^{-1/\mu}d|z|$ near infinity which generally diverges for any $\mu\geq 1$. In special cases it does converge; for instance with $N\geq 2$ symmetry the antipodal vortex sheets cancel out sufficiently to improve the integrand to $\sim|z|^{-1-1/\mu}$, convergent for all $\mu<\infty$. Even so infinity becomes dominant like $O(\mu)$ as $\mu\nearrow\infty$, but for $N\geq 3$ a further cancellation avoids even that effect.

2.2.3. Conservation arguments

(9) is not only a differential equation but an ordinary one; we have modelled $v^{{x}}$ as a function of $\xi$ alone. The ODE is not difficult to solve, but we minimize effort on irrelevant solutions by first considering conservation of mass. $\nabla\cdot\mathbf{v}=0$ has a self-similar form

[TABLE]

where $2$ comes from the number of dimensions. Integrate this over a small “cusp triangle” (formed in fig. 6 by solid and dashed vortex sheet and dashed-dotted “mouth”) from the cusp ${x}_{*}$ to some ${x}>{x}_{*}$: since the normal pseudovelocity $\mathbf{q}\cdot\mathbf{n}$ is zero at each side of a vortex sheet, only the “mouth” boundary term remains:

[TABLE]

On the left we use $q^{{x}}\approx\underline{q}^{x}(\xi)$. The right-hand side is always negative, but for a “cusp” the area is reasonably estimated as less than the area of the straightened-side triangle (dotted in fig. 6) which is $\xi y(\xi)$. Hence

[TABLE]

and as $\xi$ is increased $\underline{q}^{x}$ is strictly decreasing. Velocity $v^{{x}}$ must be a constant plus $O(\xi)$, and same for ${\Gamma}_{x}$ which differs from $v^{{x}}$ by a smooth function $\operatorname{Re}e$.

2.2.4. Constraints

Multiply both sides of (9) by ${\Gamma}_{x}$ to obtain

[TABLE]

completing the square for ${\Gamma}_{x}$ yields

[TABLE]

The first parenthesis is recognized to be $\underline{q}^{x}$ since (recall fig. 5 right) the self-induced velocity is ${\Gamma}_{x}$ between the sheets, so

[TABLE]

So we obtain

[TABLE]

Mass conservation was shown to require $\underline{q}^{x}\nearrow 0$ as we approach the cusp ($\xi\searrow 0$), which implies (for $*$ indicating values in the cusp limit)

[TABLE]

Hence ${\Gamma}_{*}=0$ or ${\Gamma}_{*}>0$; the latter appears to happen only in borderline cases, so we defer it to section 4.2.

2.2.5. ${x}$ ODE solutions in case ${\Gamma}_{*}=0$

By (13) ${\Gamma}_{*}=0$ requires

[TABLE]

so (11) and (12) simplify to

[TABLE]

and

[TABLE]

After $\partial_{\xi}$ of the latter, producing

[TABLE]

the former substitutes the final ${\Gamma}_{\xi}$ to yield

[TABLE]

This ODE turns autonomous if we eliminate the $\xi$ dependence by the substitution

[TABLE]

(which is also natural in light of $\underline{q}^{x}=O(\xi)$) to obtain

[TABLE]

Divide by $\xi>0$ and factor using $(\mu-1+e_{1})+(\mu-e_{1})=2\mu-1$:

[TABLE]

Our conservation-of-mass arguments showed that $Q=\underline{q}^{x}/\xi$ is negative but bounded in the cusp limit $-\log\xi\rightarrow+\infty$ (see (10)). This is possible only if $Q$ converges to a nonpositive root of the numerator. The two roots are

[TABLE]

The two roots are equal if and only if $e_{1}=\frac{1}{2}$, a borderline case we ignore to keep the discussion concise and avoid $\log\xi$ terms.

First consider solutions that converge to one of the roots from one side. Then that root must be asymptotically stable on that side. The right-hand side of (19) has positive denominator for $Q<0$ whereas the leading $-Q^{2}$ part of the quadratic numerator makes it positive between the roots, negative outside. Hence if the roots are distinct, the smaller root is either unstable on both sides or not negative and hence irrelevant, so we may ignore it. On the other hand if the larger root is negative, then it is asymptotically stable on both sides. The corresponding solutions will indeed turn out to lead to meaningful vortex cusps.

If $e_{1}<\frac{1}{2}$, then $Q=1-e_{1}-\mu$ root is the larger root. It is negative if and only if $e_{1}>1-\mu$. We only consider $\mu>\frac{1}{2}$ (see introduction), so $1-\mu<\frac{1}{2}$, meaning the range ${]1-\mu,\frac{1}{2}[}$ for $e_{1}$ is never empty. However, the constraint $1-\mu<e_{1}$ is significant especially for smaller $\mu$; for $\mu<1$ it does not permit cusps with $e_{1}$ arbitrarily small. The root corresponds to $G=1-2e_{1}+...$, so we obtain solutions

[TABLE]

If on the other hand $e_{1}>\frac{1}{2}$, then $Q=e_{1}-\mu$ is the larger root. It is negative if and only if $e_{1}<\mu$; again the range ${]\frac{1}{2},\mu[}$ is nonempty but especially for smaller $\mu$ a significant constraint. The root corresponds to $G=0+...$, so the leading term of ${\Gamma}$ has to be determined by linearizing (19) there: passing to

[TABLE]

the linearization at $G=0$ is

[TABLE]

with solutions $G\sim\xi^{(2e_{1}-1)/(\mu-e_{1})}$.

In summary:

[TABLE]

Finally consider solutions that equal one of the roots, so that asymptotic stability is not necessary. $Q=e_{1}-\mu$ corresponds to ${\Gamma}_{x}=0$, the trivial case of no sheets at all which we ignore. $Q=1-e_{1}-\mu$ for $e_{1}<\frac{1}{2}$ is a special case of the discussion above; for $e_{1}>\frac{1}{2}$ it does not produce cusps, as will become apparent at the end of section 2.3.3.

2.3. Vertical velocity modelling

2.3.1. Simple model

Given the relationship between $\xi$ and ${\Gamma}$, it remains to derive a model for the upper sheet height ${y}={y}(\xi)$, using the imaginary (vertical) part of the self-similar Birkhoff-Rott equation (1):

[TABLE]

Returning to our idea for approximating the lower and upper sheet as locally straight and horizontal (fig. 5), the simplest model is to adopt the resulting cusp-induced velocity $v^{{y}}=0$. Together with the non-cusp part $-\operatorname{Im}e$ from (8) we obtain

[TABLE]

(23) yields ${\Gamma}/{\Gamma}_{\xi}=\xi/\beta$, so near the cusp the approximation

[TABLE]

follows, with solution

[TABLE]

for cusp exponent

[TABLE]

Substitute (24) to obtain

[TABLE]

2.3.2. Improved model

However, we will see later that (27) is not the correct exponent if $e_{1}<\frac{1}{2}$. While the horizontal-straight-sheet approximation is appropriate for the ${x}$ ODE, ${y}$ decays faster so that relative errors are larger.

A better approximation is obtained by the following simple argument (fig. 7): the jump of velocity from below to above the upper sheet is tangential, hence must be

[TABLE]

($-$ due to counterclockwise circulation for positive ${\Gamma}_{s}$). We allow nonzero but still small slopes ${y}_{x}$, so to leading order the square root is $1$ and also $-{\Gamma}_{s}\approx-{\Gamma}_{x}$. So we still obtain the same approximation $-{\Gamma}_{x}$ for the $v^{{x}}$ jump, but now $v^{{y}}$ jumps by $-{\Gamma}_{x}{y}_{x}$. To respect the slip condition we take $\underline{v^{{y}}}=0$ on the lower side of the upper sheet, yielding $\overline{v}^{{y}}=-{\Gamma}_{\xi}{y}_{\xi}$ on the upper side and an arithmetic average

[TABLE]

on the sheet. Now (26) yields the improved model

[TABLE]

Compare the new ${y}_{\xi}$ term on the right-hand side to the left-hand side one which has coefficient $\sim\xi$. For $e_{1}<\frac{1}{2}$ the new coefficient ${\Gamma}_{\xi}$ is $\sim\xi$, hence its term non-negligible, and (25) yields

[TABLE]

with solution

[TABLE]

That this “improved” exponent is not the “simple” exponent in (27) already confirms that our previous approximation omitted a significant term.

2.3.3. Correct model

Both numerical experiments and heuristic arguments will show later that the improved exponent $\alpha$ is still incorrect. The simplest way to see the inadequacy is to invoke again the conservation of mass arguments from section 2.2.3. Again integration of $-2\mu=\nabla_{\mathbf{x}}\cdot\mathbf{q}$ over the cusp triangle (fig. 6) yields

[TABLE]

omitting all but the leading terms and using ${y}(\xi)=C\xi^{\alpha}+...$ and the $\underline{q}^{x}$ asymptotics from (15) and (23) we see

[TABLE]

We solve for

[TABLE]

Hence, given ${\Gamma}(\xi)$, conservation of mass already determines ${y}(\xi)$!

We emphasize that $e_{1}$ is not a fudge factor that can be used to realize any observed cusp exponent $\alpha$. For example in case of $N\geq 3$ symmetry $e_{1}$ is zero so that the prediction (31) depends only on $\mu$. Even without symmetry, $e_{1}$ can be measured independently, for example in numerics by calculating the velocity field derivative in an evaluation point near the cusp but outside and not too close so that the cusp part of the sheet pair almost cancels. Besides, the ${\Gamma},{x}$ relationship can also be measured; the formula (24) for $\beta$ is likewise dependent on $e_{1}$ and $\mu$, providing a second constraint for the two. $e_{1}$ depends on the vorticity in every part of the similarity plane; it cannot be modelled from local considerations.

Our restrictions

[TABLE]

have physical interpretations. The lower bound is due to $\underline{q}^{x}<0$ inside the cusp. Upper bound: consider a flow with no vortex sheets at all, but a velocity $e=e_{1}z$. Then the pseudo-velocity $\mathbf{q}=\mathbf{v}-\mu{\mathbf{x}}=((e_{1}-\mu){x},(-e_{1}-\mu){y})$ has a zero ${\mathbf{x}}=0$ that is a sink if $e_{1}<\mu$, but a saddle if $e_{1}>\mu$. Naturally a background pseudo-velocity that is expanding along the horizontal axis does not favor sheets approaching the origin along that axis, let alone forming cusps, unless the sheets are sufficiently strong to overcome the background with their self-induced velocity.

At this point we may collect the loose end left in section 2.2.5 where we ignored the constant solution $Q=1-e_{1}-\mu$, i.e. ${\Gamma}_{x}=(1-2e_{1}){x}$, in case $e_{1}>\frac{1}{2}$. Repeating our conservation of mass argument we obtain

[TABLE]

But a proper ${y}\sim{x}^{\alpha}$ cusp, in particular our small-slopes assumption $|{y}_{x}|\ll 1$, requires $\alpha>1$. Although this $Q$ represents an appropriate solution of the ${x}$ ODE, it leads to an unusable solution of the ${y}$ ODE. Hence we were justified in rejecting it.

A full “correct” model ODE yielding $\alpha$ as in (31) can be derived by considering again conservation in the cusp triangle (fig. 6):

[TABLE]

Take $\partial_{\xi}$, corresponding to moving the “mouth”:

[TABLE]

Solve for

[TABLE]

Substitute (15):

[TABLE]

For $e_{1}<\frac{1}{2}$ (23) shows

[TABLE]

which generates the desired exponent. The ODE can also be used to obtain terms beyond the leading ones, but as the asymptotic expansion deepens physical effects we neglected during derivation of the model may become significant.

3. Numerical validation

3.1. Methods

To test our models we compare their predictions to numerical data obtained from discretizing the Birkhoff-Rott equation itself. Besides general methods for unsteady vortex sheet rollup or other free-surface evolutions (e.g. [Kra91] or [MBO82]), numerics adapted to self-similar flow using point vortices were successfully used to compute algebraic vortex spirals in [Pul78, Pul89]. The vortex sheet was approximated by replacing the continuous sheet with discrete point vortices, choosing evaluation points for the Birkhoff-Rott equation midway between adjacent vortices to mitigate the strong singularity of the point vortex velocity field.

In our experience this approach is quite effective for vortex spiral calculations, well beyond the cases discussed by [Pul78, Pul89]. But in the case of vortex cusps or similar cases it is clearly ineffective, due to the cancellation between positive and negative vorticity (see fig. 8 left): if the upper sheet is resolved with a moderate number of point vortices they will near the cusp inevitably be at wide horizontal spacing relative to the vertical distance to the mirror vortex. Then the velocity is highly singular in their vicinity, but nearly zero velocity away from each pair. Clearly such a field has no resemblance to the real one. To obtain a reasonable approximation the horizontal spacing must be small compared to the sheet-image gap. But this is rather expensive because in a cusp the sheet and mirror image approach very rapidly.

Hence we must use low horizontal resolution and increase the degree instead. A natural concern is that physically relevant wavelengths may be suppressed, leading to spurious solutions. But oscillations are not expected in self-similar (as opposed to unsteady) vortex sheets; e.g. it has been shown that vortex sheet with some minimum regularity that exist for a positive time must a fortiori be analytic ([Leb02], [Wu06]). It is also easy, especially after a solution has already converged at coarse resolution, to test by selecting a few segments and refining them far more than the rest, to the point where horizontal resolution is much finer than the vertical gap; none of our experiments suggested that additional resolution was required.

While point vortices have obvious problems, it is more surprising that piecewise linear vortex segments are still inadequate. The strongest hint is in our “correct” ${y}$ ODE (33): it contains a term ${\Gamma}_{\xi\xi}$ that cannot be resolved by a piecewise linear relationship between $\xi$ and ${\Gamma}$. Indeed our numerical results will show clearly that a continuous piecewise linear representation produces cusps with the “improved” cusp exponent (28), not the “correct” one (31). In fact the linear approximation violates conservation of mass like the “improved” model did, as can be observed in numerical solutions by testing $-2\mu=\nabla\cdot\mathbf{q}$ via numerically evaluated mass flux boundary integrals. The defect is not cured by moderate refinement. It turns out that piecewise quadratic or higher approximation is needed to obtain accurate results.

This is rather inconvenient: not only are higher-degree representations segment-by-segment more expensive to evaluate than linear ones, which are already more costly than point vortices, but in addition higher order numerics are prone to oscillations and stability problems. Fortunately higher order is only needed for segments that are close to the evaluation point.

In detail our approximation works as follows: we parametrize the upper sheet by a function $z$ of ${\Gamma}$. The numerical scheme solves for variable $z_{i}$ approximating $z({\Gamma}_{i})$ in finitely many ${\Gamma}_{i}$. (For point vortices it is important that the ${\Gamma}_{i}$ are evenly spaced, but as the degree of approximation increases the influence of spacing decreases.) The sheet segment from $z_{i}$ to $z_{i+1}$ is approximated by a polynomial interpolating these two points and a few adjacent one, for example $z_{i+2}$ or $z_{i-1}$ for quadratic and both for cubic, etc.

Although the approximated sheet is continuous, it may have corners at each $z_{i}$ even if the degree is high. It is of course possible to require tangents and possible higher derivatives to match, but this turned out to have unpleasant side effects in our experiments. Asking for continuity of as many derivatives as possible leads to natural cubic splines and similar interpolation schemes, which are non-local, slowing the calculation, although locality can be obtained by compromising on top-order continuity. But more importantly the higher the spline smoothness the higher the tendency for oscillations at nonsmooth parts of the sheet. Consider for example the fork point in fig. 9 second from top, or for a more extreme example the ODE solution for $\mu=1+\epsilon$ which is almost ${y}(\xi)=0$ for $x$ inside the cusp, ${y}(\xi)=\text{const}\cdot\xi$ outside, i.e. a kink). Trying to approximate such shapes by natural cubic splines would cause strong oscillations near the kink, with upper-lower sheet intersection very hard to avoid, among other problems.

On the other hand if the sheet approximation has corners, then the velocity integral produces logarithmic singularities. In our experiments this did not appear to be important. As a sufficiently smooth part of the sheet is refined, the coefficient of each logarithmic singularity decreases. Moreover the evaluation points are chosen in the midpoints between the corners, far from the singularities and possibly with an extra order of cancellation, similar to [Pul78, Pul89] where point vortices had far more severe $1/r$ singularities in place of the mild $\log r$.

For a piecewise polynomial function $z$ of ${\Gamma}$ the Birkhoff-Rott integral $\int(Z-z)^{-1}d\Gamma$ can be done by partial fractions, yielding formulas with several complex logarithms. For higher degree this requires polynomial root finding. For the purpose of this paper, more concerned with correctness of the modelling rather than optimizing speed of numerical schemes, we simply used the usual quadratic or Cardano formulas.

The Birkhoff-Rott equation (1) is discretized by evaluation in the arithmetic averages ${\Gamma}_{i+\frac{1}{2}}$ of ${\Gamma}_{i},{\Gamma}_{i+1}$. On the left-hand sides $z$ and its derivatives are replaced by the segment $(i,i+1)$ polynomial approximation. This produces finitely many nonlinear algebraic equations with complex values depending on complex variables $z_{i}$. Newton iteration is applied to the system, solving the linearization by direct methods using the exact derivative matrix. Although the cubic cost growth limited the number of vortices to a few thousand, this was more than sufficient for our purposes so that it was not necessary to explore more efficient alternatives.

While for the exact sheet ${\Gamma}$ ranges from [math] (cusp) to $\infty$ (outer limit), the approximation uses inner and outer cutoffs. At the outer cutoff several fixed $z$ are calculated by the straight-sheet asymptotes (3). The outer cutoff must be taken large enough to avoid large errors, especially for $N\leq 2$ and large $\mu$. At the inner cutoff, for the innermost segment $z_{0},z_{1}$ the polynomial degree was limited to quadratic, using $z_{2}$ as third interpolation point.

The innermost point and mirror image have a small nonzero gap. There the sheet and image having a tendency to split apart if the calculation becomes unstable (e.g. if overly large Newton iteration steps are taken or if the initial approximation is poor). We experimented with extending the sheet by fitting $y\sim\xi^{\alpha}$ cusp shapes and generating additional $z$ points on them, but stability seemed to worsen rather than improve.

Instead of calculating “global” cusps, especially for $N=1$ (no symmetry) it is interesting to study “local” ones, with a small outer cutoff; this is reasonable since in many applications the vortex cusp is merely part of a larger incompressible flow, or a quasi-incompressible region of a compressible flow. Changing the outer cutoff adds vorticity at a uniform distance from the cusp, hence induces a holomorphic velocity contribution, so the effect of cutoff changes can be modelled by adjusting our $e$. In addition, to stabilize the solution against outer-sheet changes, it is possible to subtract from all calculated $w$ the value computed in a fixed point, say the origin. A constant added to $w$ can in (1) be absorbed by a constant shift added to $z$, since every term except $\mu z$ is either a difference of two $z$ or a derivative of $z$. Hence only the location of the cusp is changed, but not its shape or circulation distribution.

We do not specify the strength of the vortex sheets since it can be scaled arbitrarily by scaling of time $t$ and space $x$, corresponding to scaling the similarity coordinate $t^{-\mu}x$. The sheet strength determines the scale at which the sheet pair transitions from near-straight lines to a cusp. This scale could be defined as (to give only one example) the radius at which the cusp point angle between sheet points $z$ and their image $z^{*}$ has decreased from $2\phi_{\infty}$ to $\phi_{\infty}$ .

For calculating spiral-end jets we combine our new methods with the techniques of Pullin: in spiral parts of sheets we approximate ${\Gamma}$ and spiral center distance $r$ as functions of traversed angle $\beta$, choosing $\beta$ step width to be $2\pi$ divided by an integer, so that interpolation points lie on several common rays, with evaluation points at half-angles between these rays. The transition between coordinates $\beta$ inside the spiral and ${\Gamma}$ outside is troublesome, but less unstable than in pure point vortex calculations; adaptive refinement, which is easily implemented in 1d, is helpful in avoiding coordinate singularities.

We emphasize that our observations leading to the correct ODE model apply to all flows with sheet and mirror image that are close and nearly parallel, even if they do not terminate in a joint cusp, so for calculations of narrow jets quadratic or higher degree is still essential, e.g. for the outside-spiral parts of fig. 9 bottom.

Like Pullin’s point vortex calculations ours are very sensitive due to the inherent instability of vortex sheets. An added concern is the proximity of sheet and image; intersection must be avoided carefully, although for special cases it is possible to modify branch cuts in logarithms so that the approximate BR integral remains smooth even as upper-sheet vortices pass below lower-sheet ones. (Of course such intersections are permissible only during iteration, not in a converged solution.)

It is important to start with a good initial approximation, such as those obtained from our ODE models, and to limit iteration step size until a good basin of attraction is reached. In particular at the inner cutoff the sheets are prone to separate, causing uncontrollable oscillations and self-intersection. For large infinity angles $\phi_{\infty}$ we found it helpful to start with numerical solutions for small ones, then gradually increasing the angle. This numerical homotopy was also useful for other parameters like $\mu$ or $e_{1}$.

3.2. Results

Vector fields in all diagrams use a truncated logarithmic scaling to avoid very large or very small arrows. Numerical spirals shown in this paper may use an inner cutoff to a central point vortex, leaving an empty core in the figures as in fig. 9 bottom, physical spirals without cutoff look like fig. 2 right, with filled-in core.

Our numerical results (see fig. 9 top) show that vortex cusps exist under some conditions: if $\mu$ is sufficiently large, and if the data is not too “large”, which means in the special case of global cusps that $\phi_{\infty}$ is sufficiently small, below a limit angle that depends on $\mu$ and other parameters. In cases where cusps do not appear spiral-end jets can be observed instead.

3.2.1. Linear and higher degree

We find that the linear reconstruction is insufficient, qualitatively incorrect and also numerically unstable: consider the case of global cusps with $\phi_{\infty}=1^{\circ}$ using $N=3$ symmetry, hence $e_{1}=0$ so that its effects can be suppressed temporarily, so the “correct” cusp exponent $\alpha$ is $(\mu+1)/(\mu-1)$ (see (31)), the “improved” one $\mu/(\mu-1)$ (see (28)). In figure 11 we vary $\mu$ in a wide range from $1$ to $5$ and plot the numerical cusp exponent $\alpha$.

Cusp exponents were estimated by fixing two sufficiently distant points well inside the cusp region of the sheet pair but not too close to the inner cutoff to suppress spurious influences. Since for $N=3$ we know the cusp location is the origin, fitting ${y}=C_{1}{x}^{\alpha}$ (or ${\Gamma}=C_{2}{x}^{\beta}$) to the data is after logarithms a simple linear equation for $\log C_{i},\alpha,\beta$.

Numerics with linear reconstruction consistently produces numerical solutions with “improved” exponent cusps (dotted curve in fig. 11) when $\mu$ is large. This is reasonable since the linear approximation is unable to resolve the 2nd derivative ${\Gamma}_{{x}{x}}$ appearing as the extra term in (33) that makes the difference between “improved” and “correct” exponent. In addition it was checked numerically that the sheets produced by the piecewise linear approximation violate conservation of mass $\nabla\cdot\mathbf{q}=-2\mu$; which is not surprising since conservation of mass arguments also forced us to accept the “correct” exponent in (31).

For $\mu$ below $2$ calculation with linear segments becomes increasingly unstable. Numerical instabilities for decreasing $\mu$ have several causes. First, cusp exponents $\alpha$ increase to infinity; then ${y}$ decreases much faster than ${x}$ in the inner cusp so that calculations involving both are prone to large roundoff errors. This problem can be reduced by more careful evaluation or higher precision arithmetic. Secondly, the limit angle for $\phi_{\infty}$ at which cusps change to jets decreases to zero as $\mu$ approaches $1$. These two issues are observed for higher degree as well, but thirdly linear approximation also seems inherently more unstable, causing more breakdowns in the Newton iteration and requiring smaller time step, contrary to the standard expectation that lower-order numerics are more robust.

Quadratic reconstruction, as well as higher-degree ones, are more robust and yield consistent results, in particular the same “correct” cusp exponent (solid curve in fig. 11, closely matching the “x” for quadratic numerics). The match is good except near $\mu=1$ where large exponent and cusp-jet transition start to take effect.

Although we are not concerned with fast numerics here, it is interesting to point out that accurate exponents can be obtained by using quadratic or higher degree only for the sheet segments near the evaluation point, but linear segments everywhere else. Again this matches our modelling since the cusp exponents and shape are largely governed by an ordinary differential equation from the local part of the BR integral, with cancellations largely eliminating non-local contributions. This allows considerable numerical acceleration, although for complex flows with additional sheets or other forms of vorticity it may be necessary to detect proximity to them, changing degree as needed, which then causes differentiability problems affecting Newton iteration and other stability issues.

In any case we use either quadratic or quadratic-nearby-linear-elsewhere reconstruction for the rest of this paper; experiments did not show significant differences between them or cubic or higher degrees.

3.2.2. Limit angle

The limit angle for $N=3$ global cusps is shown in figure 13. Clearly below $\mu=1.3$ the limit angles become rather restrictive, converging to [math] as $\mu\searrow 1$. For $\mu\leq 1.05$, near the left side, large cusp exponents cause unreliable calculations. The limit angle diagram for $N=2$ and other cases is similar; in those cases $e_{1}\rightarrow 0$ along with $\mu\searrow 1$ so that it is not possible to obtain cusps at or below $\mu=1$.

However, it is possible to see local cusps for $\mu\leq 1$ if a large strain coefficient $e_{1}$ is generated by other vorticity in the similarity plane or by externally imposed harmonic velocity fields. In figure 13 we add to the BR integral velocity an additional harmonic velocity field to obtain $e_{1}\approx 0.3$. The two straight segments on the right of fig. 13 are fixed, only the curved parts solve the BR equations. (The pronounced corner between straight and curved segments was inconsequential.) It is possible to achieve cusps at $\mu=0.9$ and lower, at fairly “large” angles between the straight segments.

The problems of linear reconstruction are also important for theoretical reasons. It is natural to look for initial data that permits more than one self-similar solution, for example two jets or jets and cusps coexisting in a “hysteresis” region of the parameter space. Sufficiently inaccurate numerical methods do suggest that jets may exist below the limit angle at which cusps disappear.

3.2.3. ${y}$,${x}$ exponents with strain

We also wish to demonstrate the effect of the “saddle coefficient” $e_{1}$ on cusp exponents. Calculations using $N=1$, i.e. no symmetry, cause the velocity integral to diverge if $\mu\geq 1$, which is a dilemma since we do not generally expect cusps for most $\mu<1$. The divergence could be controlled by numerically subtracting the velocity induced at $z=0$ from every other velocity, i.e. adjusting the velocity integral by a large constant that becomes infinity as the outer cutoff is increased to infinity. But the results did not differ much from the alternative, using $N=2$, i.e. eliminating the divergence by symmetry, which is the only case we present here (fig. 14).

We used a cusp with $\phi_{\infty}=10^{\circ}$; since $e_{1}$ is mostly generated by the outer non-cusp part of the sheet, a smaller angle would cause a relatively small and harder-to-estimate $e_{1}$, while a larger angle causes cusps to cease to exist for $\mu$ far above $1$ (which is one of the reasons for the inaccurate results seen near $\mu=1$ in the diagram, the other reason being that the larger the cusp exponent the more their numerical estimation is affected by numerical roundoff errors from tiny ${y}$). The coefficient $e_{1}$ was estimated numerically by evaluating the velocity integral derivative in some point near the cusp but outside the sheet pairs where their cusp parts undergo strong cancellation and most of the velocity is induced by far-away sheet parts.

We evaluate the “correct” cusp exponent $\alpha=(\mu+1-e_{1})/(\mu-1+e_{1})$ from (31) once with $e_{1}=0$ (dotted curve in fig. 14) and once with the numerically estimated $e_{1}$ (solid curve); none of these $e_{1}$ exceeded $\frac{1}{2}$. The numerical measurements fit the solid curve well, so we conclude that our modelling likely yielded the correct $e_{1}$ dependence of the cusp exponent. For calculating the diagram we used an unncessarily coarse discretization and smallish outer cutoff to display the errors more clearly. The numerical cusp exponents fit the solid curve closely, except near the $\mu\searrow 1$ end where large cusp exponents and transition to jetting cause now-familiar inaccuracy; the fit at large $\mu$ can be improved further by taking larger outer cutoffs. This corresponds to $e_{1}$ being more sensitive to outer parts of the sheets when $\mu$ is large. To improve the fit from $\mu=1.5$ to $\mu=3$ it was sufficient to refine the discretization alone.

3.2.4. ${\Gamma}$,${x}$ exponents with strain

Finally, in figure 15 we consider the exponent $\beta$ in the ${\Gamma}=C\xi^{\beta}+..$ relationship in the cusp limit. Formula (24) distinguishes the cases $e_{1}<\frac{1}{2}$ and $e_{1}>\frac{1}{2}$. The strain coefficient $e_{1}$ induced by the non-cusp parts of the sheet increases monotonically as the infinity angle $2\phi_{\infty}$ between the outer parts increases; $e_{1}$ does exceed $\frac{1}{2}$ if $\mu$ is chosen sufficiently large (with $\mu=2.5$ chosen for computing fig. 15), allowing not only a larger angle but also using that for $\mu\nearrow\infty$ the effects of the outer sheet in the velocity integral become ever stronger.

The results, again using a numerically estimated $e_{1}$, show a very close match between formula (24) (solid curve) and the numerically determined ${\Gamma},{x}$ exponent $\beta$ (dashed and dotted curves). Near $e_{1}=\frac{1}{2}$ the results are almost $0.02$ apart. Errors decrease rapidly at a distance from $e_{1}=\frac{1}{2}$ as the inner cutoff is taken smaller and simultaneously the initial $-d\xi/\xi$ stepsize is reduced (the dotted curve has inner cutoff improved by a factor $1/10$ and resolution by $1/2$). But convergence is very slow near $e_{1}\approx\frac{1}{2}$, which is natural: not only is it difficult to distinguish a power law exponent $\beta=0.5$ from $\beta=0.49$ numerically, but as $e_{1}\nearrow\frac{1}{2}$ the coefficient $C$ in ${\Gamma}=C{x}^{0.5}+...$ also vanishes (see (25)), meaning the leading term is dominated by the next one for all but astronomically small inner cutoffs.

We conclude that our formulas (24) and (31) for ${\Gamma},{x}$ and ${y},{x}$ exponents are probably the physically correct exponents, with close match between numerical data and heuristic derivations. We emphasize again that the ${y},{x}$ exponent is completely determined by the conservation of mass argument, so that the point of attack for criticism is confined to the ${\Gamma},{x}$ exponent and the ODE model and simplifying assumptions that led to it.

3.2.5. Clockwise upper circulation

For $e_{1}<\frac{1}{2}$ our ${x}$ ODE discussion in section 2.2.5 yielded solutions ${\Gamma}_{\xi}=(\frac{1}{2}-e_{1})\xi+...$; note $\frac{1}{2}-e_{1}>0$ so that ${\Gamma}_{\xi}>0$ near the cusp, corresponding to counterclockwise circulation on the upper sheet near cusps opening to the right. For clockwise circulation jets rather than cusps will be observed, as in the local flow in figure 16 ($\mu=1.3$, $N=1$).

However, if $e_{1}>\frac{1}{2}$, then section 2.2.5 permitted solutions ${\Gamma}_{\xi}=C\xi^{\beta-1}+...$ for $\beta-1>1$ (see (24)), with arbitrary $C$. In particular we may choose $C<0$, allowing clockwise circulation on the upper sheet at the cusp. Numerically it can be confirmed that such solutions are possible: figure 17 shows a local $N=2$ symmetric cusp with $\mu=1.5$ and $e_{1}$ adjusted to $0.7$ by adding an external harmonic velocity field $w=cz$ with suitable constant $c$. Again the straight outer segments are fixed chosen data, whereas the curved inner parts are numerical Birkhoff-Rott solutions forming a cusp in the origin.

Generally the case of global flows with clockwise upper sheet circulation and infinity angle $\phi_{\infty}$ is equivalent to counterclockwise circulation and infinity angle $360^{\circ}/N-\phi_{\infty}$. For example the jets in fig. 10 with $\phi_{\infty}=48^{\circ}$ and $\phi_{\infty}=77^{\circ}$ can be rotated $90^{\circ}$ to yield solutions for the clockwise upper sheet circulation case for $\phi_{\infty}=42^{\circ}$ and $\phi_{\infty}=13^{\circ}$.

3.2.6. Jetting and near-limit-angles

As $\phi_{\infty}$ approaches the limit angle from below the cusp attains the $Y$ shape in fig. 9 second from top. Between the sheets near the meeting point of the Y the pseudo-velocity $\mathbf{q}$ approaches zero. We were not able to discerm by numerical experiments whether it reaches zero or whether transition to jets occurs earlier. This would be of interest for obtaining formulas or criteria predicting the limit angle.

A natural idea is to consider where the cusp fails to exist because the denominator $\mu-1+e_{1}$ of the cusp exponent $\alpha$ in (31) crosses zero. But this “local” criterion is easily refuted once we consider $N=3$ where $e_{1}=0$ due to cancellation so that cusp-jet transition would always occur at $\mu=1$, clearly at odds with our numerical observations (see fig. 13).

Hence cusp-jet transition is a more non-local effect, influenced strongly by the ambient flow especially near the meeting point rather than the cusp. Better criteria may be based on conservation of mass arguments, possibly using $\underline{q}^{x}$ or similar variables approaching zero. It is not clear whether we should expect a single criterion to be suitable for large classes of flows.

Above the limit angle the cusp ceases to exist; sheet and mirror image separate with ends flipping rightward in between them and forming a spiral-ends jet; see fig. 9 bottom two as well as fig. 10. As $\phi_{\infty}$ is increased this flow pattern persists, with sheets separating further. On the other hand as $\phi_{\infty}$ converges to the limit angle from above the distance of the two sheets converges to zero from above.

For $\mu=1.3$ the limit angle is near $\phi_{\infty}=25.5^{\circ}$, but already at $\phi_{\infty}=48^{\circ}$ (fig. 10 top), and even more so at $\phi_{\infty}=35^{\circ}$ (fig. 9 bottom two) which is still almost $10^{\circ}$ above the limit angle, the jets are very small. Indeed numerical calculations show that distance between the spiral centers decreases quadratically in the distance of $\phi_{\infty}$ from limit angle. Due to the instability of vortex sheets in general and the smallness of the jet in particular, computing jet solutions near the limit angle is even more challenging than computing vortex cusps.

3.2.7. Relation to Mach reflection

[HVBDE03] and [VBDEH04] have begun a theoretical and numerical study on jetting in Mach reflection of shock waves for compressible flow, whose quasi-incompressible regions correspond to $\mu=1$ in our context. It is well-established (see [HM98], [Ser07] and [Ell18]) that the vortex sheets generated by single Mach reflection (fig. 3 right) have counterclockwise circulation. This is very important because it allows, by our analysis above, the presence of cusps even if the ambient flow generates only weak or negative-sign strain coefficients $e_{1}$.

By regarding sheet-wall interaction as an essentially incompressible phenomenon and adapting vortex sheet numerics to obtain better resolution, we found that when jets occur they may be rather small, possibly small enough to be obscured by numerical or physical dissipation. Besides, for $\mu=1$ we find that cusps are possible only if the ambient flow generates a positive (but not too large) strain coefficient $e_{1}$. Even then it is necessary that the outer (non-cusp) part of the sheet forms a rather small angle $\phi_{\infty}$ to the wall, increasingly smaller as $e_{1}$ approaches zero from above. This suggests that many flows that appear to be cusps are really jets in disguise. For instance $N$-symmetric Mach reflections with $N\geq 3$ must have $e_{1}=0$ for cusps in the origin which is not possible if $\mu=1$. The problem is compounded by the rather large cusp exponents that result from near-zeros of the denominator $\mu-1+e_{1}$ in (31) when $\mu=1$ and $e_{1}>0$ is small; such exponents cause most of the cusp part to be very close to the wall, making it appear at large scale that the sheet meets the wall at a positive angle.

Shock capturing methods on non-adapted grids are much less accurate than the vortex methods employed here. In calculations small jets are easily obscured by numerical or physical boundary layers and by general inaccuracy. Of course in an exactly self-similar flow the inviscid features of the flow often outgrow the viscous ones; whether the jet-cusp distinction is important would depend on the intrinsic length and time scales of the application at hand. Small-scale cusp-jet differences may also make themselves felt in attachment vs. separation of boundary layers that cause differences at a much larger scale.

A natural and relatively cheap prediction for cusp vs. jet is to measure in the output of shock capturing calculations the sheet-wall angle $\phi_{\infty}$ and the velocity $e$ outside the hypothetical cusp; this requires sufficiently high resolution for a significant region relatively free of boundary effects and with Mach number low enough that flow can be treated as incompressible. A negative $e_{1}$ immediately rules out cusps; given a positive value either the model ODE or (for larger $\phi_{\infty}$ where our approximations are inaccurate outside the inner cusp region) numerical solution of the full Birkhoff-Rott equation can determine whether jetting occurs. This approach avoids modification of existing numerical software which is more time-consuming even when the source code is available.

4. Additional considerations

4.1. $e$ asymptotics

During modelling starting with (4) we have neglected the velocities induced in points $z$ by “non-near” other points $z^{\prime}$ when both are located in the cusp region of the sheets. Of course despite strong cancellation these contributions are also nonzero; here we perform a crude check of the decay rate.

Assume for simplicity the cusp ${x}_{*}$ has been shifted to [math]. We regard $z^{\prime}$ as “not near” $z$ if $|{x}^{\prime}-{x}|\geq\delta{x}$ for some constant $\delta>0$, with $0<{x},{x}^{\prime}<\overline{x}$ where $\overline{x}$ marks the “boundary” of the cusp region, and consider the ${x},{x}^{\prime}\rightarrow{x}_{*}=0$ limit. First consider the $0<{x}^{\prime}<(1-\delta){x}$ part of the Birkhoff-Rott integral:

[TABLE]

We obtained solutions of type ${\Gamma}\sim{x}^{\beta}$ and ${y}\sim{x}^{\alpha}$ with $\alpha>1$, so ${y},{y}^{\prime}$ in the denominator are dominated by ${x},{x}^{\prime}$ (if non-near):

[TABLE]

Non-near means $({x}-{x}^{\prime})^{2}\sim({x}^{\prime})^{2}$, so

[TABLE]

assuming the last exponent of ${x}$ is positive. Our modelling results yielded $2\leq\beta<\infty$; proper cusps have $1<\alpha<\infty$, so the exponent is in fact $>1$.

However, if $\beta=2$ and $\alpha=1+\epsilon$, then the exponent is $1+\epsilon$, just barely above $1$. Such $\beta,\alpha$ do occur; consider for example $N\geq 3$ symmetry so that $e_{1}=0$ and thus $\beta=2$, then take $\mu\rightarrow\infty$ so that exponent $\alpha+\beta-2=(\mu+1-e_{1})/(\mu-1+e_{1})\searrow 1$; already for $\mu>3$ the exponent is below $2$.

For the integral over $(1+\delta){x}<{x}^{\prime}<\overline{x}$ the only difference is a term from the $\overline{x}$ boundary that is constant in ${x}$, hence can be regarded as part of $e_{0}$.

We conclude that our approximation $e(z)=e_{0}+e_{1}cz+...$ was reasonable, but that expanding to

[TABLE]

or deeper is not always reasonable since the quadratic term may be dominated by effects we already neglected.

4.2. ${\Gamma}_{*}>0$ cusps

In section 2.2.4 we deferred the case ${\Gamma}_{*}>0$; we revisit it here, showing that such cusps may occur but are “non-generic”. In $t,\mathbf{x}$ coordinates $\Gamma$ is merely a parameter that can be changed by an additive constant without consequences. But the $t$-scaled ${\Gamma}$ in in similarity coordinates is not as arbitrary since it corresponds the rate at which circulation passes through a point through self-similar movement alone. Recall (13) that ${\Gamma}_{*}>0$ means $\mu{x}_{*}\neq e_{0}$. Take $\partial_{\xi}$ of (12): again we use

[TABLE]

to calculate

[TABLE]

The middle terms decay as $\xi\searrow 0$, but the last term is nonzero unless the scalar constraint

[TABLE]

is satisfied. Assume it is not. Then (36) yields

[TABLE]

for a nonzero constant $C$ (necessarily positive as evident below). Integrate:

[TABLE]

where $\underline{q}^{x}\rightarrow 0$ as $\xi\searrow 0$ eliminated the integration constant; $-$ is chosen since $\underline{q}^{x}<0$ (derived from conservation of mass). Then (38) shows

[TABLE]

Now consider our “correct” model (32) for the upper sheet ${y}={y}({x})$:

[TABLE]

By (40) the numerator is dominated by $\underline{q}^{x}_{\xi}$; using also (39) we find

[TABLE]

Hence a cusp cannot form; ${y}$ necessarily blows up as the cusp is approached.

We emphasize that this blowup is relatively mild and easily overlooked in low-accuracy numerical calculations where it may appear falsely that a cusp might form upon further refinement. Besides, blowup is entirely due to the ${\Gamma}_{\xi\xi}$ term (corresponding to $\underline{q}^{x}_{\xi}$ above) which is only present in the “correct” model, but not in the other two; numerics with less than quadratic approximation would not recognize blowup, producing false cusps!

On the other hand if the constraint (37) is satisfied, there may be a chance for cusps to form. A candidate for such cusps are the flows arising as the maximal-$\phi_{\infty}$ limit, see fig. 9 right: in that limit the “standard” ${\Gamma}_{*}=0$ cusp collapses to two coinciding horizontal lines; the region where they separate appears to be a cusp as well, but necessarily with ${\Gamma}_{*}>0$ as observed in numerics. However, further modelling requires deeper Taylor expansions of $e$ which are not always justified, as explained in the previous section. Besides, the theoretical condition (37) imposes a scalar constraint on the parameters; similarly $\phi_{\infty}$ being the limit angle is a “borderline” case, not robust under parameter changes. Since ${\Gamma}_{*}>0$ cusps do not obviously appear in other circumstances our investigation appears to have reached the point of diminishing returns.