4.1. Alternative A1: forced responses

There is a relationship between the solutions of (4.1) at squared forcing frequency $\omega ^2=\sigma _r^2=1/(1+r^2)$, associated with waveform function $f_r$, and those at $1-\sigma _r^2=r^2/(1+r^2)=1/(1+(1/r)^2)$, associated with $f_{1/r}$. The details are provided in Appendix A. It is therefore sufficient to consider $r\in (0,1)$, i.e. $\sigma _r^2\in (1/2,1)$.

First, we identify waveform functions $f_r$ satisfying (4.12) for selected rational values $r=n/m$ with $m>n$. Evaluating (4.12) at $\boldsymbol {x}=[x_j]_{1\le\, j\le m}$ with $x_j=j/(2m)$, $j=1,\ldots ,m$, and using the symmetries approximately $0$ ($\,f_r$ odd) and $1/2$ (4.10), yields an $m\times m$ linear system in $\boldsymbol {f}_r(\boldsymbol {x})=[\,f_r(x_j)]_{1\le\, j\le m}$,

(4.14)\begin{equation} \boldsymbol{\mathsf{A}}\,\boldsymbol{f}_r(\boldsymbol{x})=2\boldsymbol{x}. \end{equation}

Harlander & Maas (Reference Harlander and Maas2007) derived a similar system using a boundary collocation approach for the solution of the Poincaré equation in two-dimensional (2-D) containers.

When $m$ and $n$ have opposite parities, this system is non-singular and has a unique solution $\boldsymbol {f}_r(\boldsymbol {x})$. The construction of the system and its solution is illustrated in Appendix B for the case with $m=11$ and $n=4$. The linear spline obtained by connecting the points $[x_j,f_r(x_j)]$ satisfies (4.10) and (4.12) for all $x\in [-(1+r)/2,(1+r)/2]$. This solution $f_r$ is plotted in red in the first row of figure 6 for the cases $n=1$ and $m=2,4,6,8$ and 10. A unique 2-periodic extension, obtained by enforcing (4.10) and $f_r(-\zeta )=-f_r(\zeta )$ for all real $\zeta$, is shown in blue. This extension satisfies $f_r(\pm 1)=0$. Note that additional periodic extensions exist for $0 < r < 1$, $[-(1+r)/2,(1+r)/2]\subset [-1,1]$.

Figure 6. Waveform function $f_r$ in $[-1.5,1.5]\times [-1.25,1.25]$ and the associated scaled forced responses $\eta _r$ at phase 0 and $\theta _r$ at phase ${\rm \pi} /2$, from (4.15) at $r$ as indicated. The blue curve is the 2-periodic extension of the red curve $f_r$ defined in the interval $[-(1+r)/2,(1+r)/2]$.

The value of the waveform function $f_r$ in $[-(1+r)/2,(1+r)/2]$ completely determines the response flow of the forced linear system (4.1) everywhere in the square cavity. This response flow is

(4.15)\begin{equation} \left. \begin{gathered} \begin{bmatrix} u_r \\ w_r \end{bmatrix} =c_r\sigma_r\begin{bmatrix}r[\,f_r(x+rz)-f_r(x-rz)]\\ 2x-f_r(x+rz)-f_r(x-rz)\end{bmatrix}\cos\sigma_rt,\\ \begin{bmatrix} p_r \\ \theta_r \end{bmatrix} ={-}c_r\begin{bmatrix} \lambda_r[2r xz-F_r(x+rz)+F_r(x-rz)]\\ 2x-f_r(x+rz)-f_r(x-rz) \end{bmatrix} \sin\sigma_rt. \end{gathered} \right\} \end{equation}

Its magnitude is linearly proportional to the forcing amplitude $\alpha$.

The forced response (4.15) is now considered for three different sequences of rationals $r=n/m$ with $n$ and $m$ of opposite parities, which converge to $r=0$, $r=1/\sqrt {2}$ (an irrational number) and $r=1$. These sequences reveal qualitatively different behaviours.

The vorticity $\eta _r$ and temperature deviation $\theta _r$ for the sequence $r=1/m$, $m=2,4,6,\ldots$ converging to $r=0$, associated with $\sigma _r^2=1/(1+r^2)=1$, are illustrated in figure 6 for the first five cases $m=2,4,6,8$ and $10$. For $r=1/2$, $\eta _r$ is piecewise constant with the pieces delimited by the characteristic lines emanating from the corners forming a pattern known as a harlequin print, and $\theta _r$ is piecewise linear in the form of two pyramids, one positive and the other negative, either side of $x=0$. As the even integer $m$ increases, the 2-periodic extension of $f_r$ looks more and more like a triangular wave with $f_r(x)=x$ for $|x|\le 0.5$ and the response flows consist of $m/2$ copies of the $r=1/2$ response with $x\to (m/2)x$.

We now consider a sequence of rationals $r=n/m$ with $m$ and $n$ of opposite parities, converging to the irrational value $r=1/\sqrt {2}$, corresponding to $\sigma _r^2=1/(1+r^2)=2/3$. Such a sequence can be obtained from successive partial sums of the continued fraction representation

(4.16)\begin{equation} \dfrac{1}{\sqrt{2}}=\dfrac{\sqrt{2}}{2}= \dfrac{1}{2}\left[1+\dfrac{1}{2+\dfrac{1}{2+\cdots}}\right], \end{equation}

according to the iteration $r\to (1+r)/(1+2r)$, starting with the case $r=1/2$. Figure 7 shows the waveform function $f_r$ and the associated response $\eta _r$ and $\theta _r$ for five consecutive iterates of $r$. The value $f_r(1/2)$ is observed to rapidly approximate $1$ as $r\to 1/\sqrt {2}$, with $f_r$ exhibiting increasingly higher frequency spatial variations, and finer and finer structure in the corresponding $\eta _r$ and $\theta _r$. The curve $f_{1/\sqrt {2}}$, obtained at the limit $r=1/\sqrt {2}$, appears to be fractal.

Figure 7. Waveform function $f_r$ in $[-1.5,1.5]\times [-1.25,1.25]$ and associated scaled forced responses $\eta _r$ at phase 0 and $\theta _r$ at phase ${\rm \pi} /2$ for a series of rational $r$ values converging to the irrational $1/\sqrt {2}$, corresponding to $\sigma _r^2=2/3$. The blue curve is the 2-periodic extension of the red curve $f_r$ defined in the interval $[-(1+r)/2,(1+r)/2]$.

Figure 8 shows $f_r$ and $\eta _r$ at a rational value of $r$ approximating $1/\sqrt {2}$, obtained from the fourteenth iterate of the sequence. The figures appear qualitatively similar to lower iterates shown in figure 7; the higher iterates allow us to zoom-in and explore the fractal nature of the response in the limit $r\to 1/\sqrt {2}$, from which the similarity relation

(4.17)\begin{equation} f_r(0.5s+[1-s]x) = s\,f_r(0.5)+[1-s]f_r(x), \end{equation}

with $s=2r/(1+r)=2\sqrt {2}-2\approx 0.828$, can be inferred. The relation (4.17), equivalent to

(4.18)\begin{equation} \begin{bmatrix}\zeta\\\, f_r(\zeta)\end{bmatrix} =\begin{bmatrix}0.5\\ \,f_r(0.5)\end{bmatrix} +(1-s)\begin{bmatrix}x-0.5\\ \,f_r(x)-f_r(0.5) \end{bmatrix}, \end{equation}

shows that the graph of $f_{1/\sqrt {2}}$ is invariant under the homothety with similarity factor $1-s=3-2\sqrt {2}\approx 0.172$ and centre $(0.5,f_{1/\sqrt {2}}(0.5))=(0.5,1)$. Figure 8 illustrates how the curve $f_{1/\sqrt {2}}$ in the range $[s/2, 1-s/2]$, of size $1-s$ centred around $0.5$, is a scaled-down copy of itself from the range $[0,1]$. The supplementary movie 3 shows a continuous zoom into the region in the neighbourhood of the point $(x,z)=(0.5,0)$ of the response $\eta _r$ shown in figure 8. The relation (4.17) is not trivial and can be verified to hold for other quadratic irrational values of $r$, with exactly the same definition of $s(r)$ when $r^2=1-1/m$, $m=3,4,\ldots$, and when $r^2=1+1/m$, $m=1,2,\ldots$, or with a modified rational relationship for other quadratic irrational values of $r$.

Figure 8. Zoom-in in the neighbourhood of $(0.5,0)$ for $\eta _r$ at phase 0; (a–c) waveform function $f_{1/\sqrt {2}}$ in windows $[-1.5,1.5]\times [-1.25,1.25]$ (a,d), $[0,1]\times [0,1]$ (b,e) and $[s/2, 1-s/2]\times [s, 1]$ (c, f) with $s=2\sqrt {2}-2$; (d–f) $\eta _r$ in regions $[-0.5, 0.5]\times [-0.5, 0.5]$ (square cavity; a,d), $[0, 1/2]\times [-1/4, 1/4]$ (b,e) and $[s/2, 1/2]\times [-(1-s)/4,(1-s)/4]$ (c,f). The supplementary movie 3 shows a continuous zoom-in of the response $\eta _r$ about the point $(x,z)=(0.5,0)$.

The third sequence of rationals considered is $r=n/m=(m-1)/m=1-1/m$, $m=2,3,4,\ldots$ converging to $r=1$, which is associated with $\sigma _r^2=1/(1+r^2)=1/2$. Figure 9 illustrates the behaviour of the solution of the functional equations, $f_r$, and the corresponding forced response $\eta _r$ and $\theta _r$. As $m$ increases, the maximum value of $f_{1-1/m}$, which is attained at $0.5$, increases unboundedly as ${O}(m)={O}(1/(1-r))$. As a consequence, the forced response also grows in magnitude as $m\to \infty$ and $r\to 1$, with the piecewise constant vorticity $\eta _r$ showing a finer and finer harlequin pattern converging towards a piecewise linear vorticity profile in the shape of an inverted pyramid, and the temperature deviation $\theta _r$ converges towards a piecewise quadratic field with isocontours in the form of perfect (arcs of) circles centred at $(\pm 0.5,0)$ and hyperbolas centred at $(0,\pm 0.5)$ meeting at the limiting characteristics $x\pm z=0$; see figure 9. The unbounded response at the limit $r=1$ corresponds to resonance, which is discussed below under the umbrella of the Fredholm alternative A2. The unbounded growth of the magnitude of the response being inversely proportional to the distance to the limiting critical frequency also conforms to the standard result from forced undamped mechanical linear oscillators.

Figure 9. Scaled waveform function $f_r$ in $[-1.5,1.5]\times [-1.25,1.25]$ and the associated scaled forced responses $\eta _r$ at phase 0 and $\theta _r$ at phase ${\rm \pi} /2$ for a series of rational $r$ values converging to $r=1$, corresponding to $\sigma _r^2=1/2$. The blue curve is the 2-periodic extension of the red curve $f_r$ defined in the interval $[-(1+r)/2,(1+r)/2]$.

4.2. Alternative A2: intrinsic modes

At exactly $r=1$, the Fredholm alternative A1 fails: the conditions $f(1)=f(0)$ ((4.10) with $\zeta =0.5$) and $f(1)+f(0)=1$ (alternative A1 with $x=0.5$) lead to a contradiction with $f(0)=0$ ($\,f$ odd). Equivalently, the system (4.14) reduces to the inconsistent single equation $[0]f(0.5)=2[0.5]$. The inconsistency of the system (4.14) is also illustrated for $r=3/5$ in Appendix B. In these cases, the Fredholm alternative A2 must hold. For $r=1$, alternative A2 reduces to (4.10); for every odd waveform function $f=f(\zeta )$ which is even around $\zeta =0.5$, there is a corresponding solution of the unforced homogeneous problem (4.1) constructed from (4.6). Examples of $f$ that naturally satisfy these symmetry conditions include the family of shifted Euler splines $f(\zeta )=S_q(\zeta -0.5)$, of regularity class $C^{q-1}$(Olver et al. (Reference Olver, Lozier, Boisvert and Clark2010), (24.17.3)). These splines are defined from the Euler polynomials $E_q$ (Olver et al. (Reference Olver, Lozier, Boisvert and Clark2010), § 24). They are $1$-antiperiodic (Olver et al. (Reference Olver, Lozier, Boisvert and Clark2010), (24.2.12)) and, hence, $2$-periodic. They also belong to the more general family of exponential Euler splines (Schoenberg (Reference Schoenberg1973), Lecture 3, § 5). This family includes the $C^0$ linear triangular wave $S_1$, the $C^1$ quadratic spline $S_2$ defined by $S_2(\zeta -0.5)=4\zeta (1-|\zeta |)$ for $|\zeta |\le 1$, and the $C^\infty$ waveform function $f(\zeta )=\sin ({\rm \pi} \zeta )=\lim _{q\to \infty }S_q(\zeta -0.5)$. Figure 10 shows the structure of $S_1$, $S_2$, $S_3$ and $S_\infty$; note the remarkable similarity between $S_q$ for $q\ge 2$.

Figure 10. Shifted Euler splines $S_q$, of regularity class $C^{q-1}$, shown in a window $[-1.5,1.5]\times [-1.25,1.25]$: $(a)$ $S_1(\zeta -0.5)=E_1(\zeta +0.5)/E_1(0)=2\zeta$ for $-0.5<\zeta <0.5$, $(b)$ $S_2(\zeta -0.5)=E_2(\zeta )/E_2(0.5)= 4\zeta (1-\zeta )$ for $0<\zeta <1$, $(c)$ $S_3(\zeta -0.5)= -E_3(\zeta +0.5)/E_3(0)=\zeta (3-4\zeta ^2)$ for $-0.5<\zeta <0.5$, and $(d)$ $S_\infty (\zeta -0.5)=\sin ({\rm \pi} \zeta )$, where $E_q$ is the Euler polynomial of degree $q$.

The choice $f_1(\zeta ) =S_\infty (\zeta -0.5)=\sin ({\rm \pi} \zeta )$ leads to the $C^\infty$ intrinsic mode

(4.19) \begin{equation} \left. \begin{gathered} \begin{bmatrix} u_h \\ w_h \end{bmatrix} \propto \sigma_{1} \begin{bmatrix} \cos{\rm \pi} x \: \sin{\rm \pi} z \\ -\sin{\rm \pi} x \: \cos{\rm \pi} z \end{bmatrix} \cos\sigma_1t, \\ \begin{bmatrix} p_h \\ \theta_h \end{bmatrix} \propto \begin{bmatrix} (\lambda_{1}/{\rm \pi}) \sin{\rm \pi} x \: \sin{\rm \pi} z \\ \sin{\rm \pi} x \: \cos{\rm \pi} z \end{bmatrix} \sin\sigma_1t, \end{gathered} \right\} \end{equation}

with $\sigma _1=1/\sqrt {2}$ and $\lambda _1=1/2$, which was shown to be excited by purely vertical forcing oscillations of the container (Yalim et al. Reference Yalim, Lopez and Welfert2018). In contrast, for the horizontal forcing considered here, the limit of the forced responses associated with the sequence $r=1-1/m$ with $m\ge 2$ corresponds to the quadratic spline $S_2$. Specifically,

(4.20)\begin{equation} \lim_{m\to\infty} \dfrac{f_{1-1/m}(\zeta)}{f_{1-1/m}(0.5)} = 4\zeta(1-|\zeta|)=S_2(\zeta-0.5), \quad |\zeta|\le1; \end{equation}

see Appendix C for a proof. The eigenmode components $\eta _1$ and $\theta _1$ are visually indistinguishable from the bottom row of figure 9, obtained for $r=1-1/m$ with $m$ large.

The resonance at $\omega ^2=1/2$ (i.e. $r=1$) corresponds to the largest peak in the Bode response diagram in figure 2$(a)$. Similar resonances occur for all rational values $r=n/m$ with both $m$ and $n$ odd integers, leading to other peaks in the Bode diagram at the corresponding values of $\omega ^2=\sigma _r^2=1/(1+r^2)$. These are all the cases for which Fredholm alternative A2 holds. The scaling relation

(4.21)\begin{equation} f_{n/m}(\zeta)=f_1(m\zeta), \end{equation}

shown in Appendix D, implies that the resulting mode, denoted $M_{m:n}$, can be simply constructed from $M_{1:1}$, with vorticity and temperature deviation components $\eta _r$ and $\theta _r$ arranged in a regular $m\times n$ tiling of the scaled components of $M_{1:1}$.

Figure 11 shows the vorticity $\eta _r$ and temperature deviation $\theta _r$ for the eigenmode $M_{5:3}$ resonating at squared frequency $\sigma _r^2=25/34\approx 0.7353$ (corresponding to $r=3/5$). The response consists of an exact $5\times 3$ tiling of copies of the components of the mode $M_{1:1}$, with piecewise linear $\eta _1$ in a pyramid pattern and piecewise quadratic $\theta _1$. The responses at slightly detuned frequencies on either side (slightly smaller and larger values of $r$) of the resonant modal case are also shown. These forced responses, being associated with rational values of $r$ that are approximations of the resonant case $r=3/5$, have piecewise constant $\eta _r$ and piecewise linear $\theta _r$. At the macroscopic scale of the figure, they appear to be smooth with one order of regularity higher than they actually possess. This is due to the natural visual (Riemann) integration of the small scales associated with the underlying network of characteristic lines emanating from the corners, resulting from the large values of $m$ and $n$ in the irreducible representation of $r=n/m$. Whereas the eigenmode has velocity and temperature fields of regularity class $C^1$ and satisfies the homogeneous linear inviscid system (4.1) in a pointwise or strong sense, the forced response at a detuned frequency has velocity and temperature fields only of class $C^0$ globally, as a result of the nature of the forcing, and are weak solutions.

Figure 11. Scaled waveform function $f_r$ in $[-1.5,1.5]\times [-1.25,1.25]$ and associated scaled forced response ($\eta _r$ at phase 0 and $\theta _r$ at phase ${\rm \pi} /2$) for two values $r=(300\pm 1)/500$ close to $3/5$ (a,e,i,d,h,l) and eigenmode at $r=3/5$ (b,f,j,c,g,k, associated with opposite scaling). The blue curve is the 2-periodic extension of the red curve $f_r$ defined in the interval $[-(1+r)/2,(1+r)/2]$.

Note that the waveform functions $f_r$ shown in figure 11 are normalized, with extrema at ${\pm }1$. In fact, the extrema of $f_r$ are, as in the case at $r=1$, unbounded at the resonant frequency. As the forcing frequency is varied across the resonant frequency, there is an inversion $f_r\to -f_r$. At resonance, the intrinsic mode of the unforced system is given by

(4.22)\begin{equation} \left. \begin{gathered} \begin{bmatrix} u_r \\ w_r \end{bmatrix} =\sigma_r\begin{bmatrix}r[\,f_r(x+rz)-f_r(x-rz)]\\ -f_r(x+rz)-f_r(x-rz)\end{bmatrix} \cos\sigma_rt, \\ \begin{bmatrix} p_r \\ \theta_r \end{bmatrix} =\begin{bmatrix} \lambda_r[F_r(x+rz)-F_r(x-rz)]\\ f_r(x+rz)+f_r(x-rz) \end{bmatrix} \sin\sigma_rt, \end{gathered} \right\} \end{equation}

obtained by letting $\alpha =0$ in (4.7). So, the change in sign in $f_r$ is equivalent to a half-period phase shift in time. This is evident from comparing figures 11(a–d) and 11(i–l). Interestingly, this shift occurs in the high spatial frequency component of $f_r$, not the low spatial frequency component. This shift is emphasized in figures 11(b,f,j) and 11(c,g,k), representing the eigenmode, scaled to $\pm 1$ according to which side of the resonant frequency and the associated parameter $r=3/5$ they correspond. Also note how the cellular grid modal pattern of the piecewise quadratic temperature deviation $\theta _r$ becomes distorted along the characteristic directions and the cells merge upon detuning away from the resonance.

4.3. Spectral representations of modal and forced responses

The visual similarity between the quadratic spline $S_2$ and the sine wave $S_\infty$ in figure 10 can be quantified via the sine Fourier series expansion (Olver et al. Reference Olver, Lozier, Boisvert and Clark2010, (24.8.4))

(4.23)\begin{align} S_2(\zeta-0.5)=\frac{32}{{\rm \pi}^3}\,\sum_{k\ge0} \frac{\sin(2k+1){\rm \pi}\zeta}{(2k+1)^3}\, =\frac{32}{{\rm \pi}^3}\,\left[S_\infty(\zeta-0.5)+ \frac{1}{3^3}S_\infty(3\zeta-0.5)+\ldots\right]. \end{align}

By approximating the factor $32/{\rm \pi} ^3\approx 1.032\approx 1$ in (4.23), the resonating mode $M_{m:n}$ can be expanded as a superposition of an infinite number of $C^\infty$ intrinsic modes $I_{(2k+1)m:(2k+1)n}$. These are the odd spatial harmonics of $I_{m:n}$, all of which have the same intrinsic frequency $\sigma _{(2k+1)n/(2k+1)m}=\sigma _{n/m}$,

(4.24)\begin{equation} M_{m:n} \propto \sum_{\ell\ge0} (2\ell+1)^{{-}3}\, I_{(2\ell+1)m:(2\ell+1)n}. \end{equation}

The components $(u_h, w_h, p_h, \theta _h)$ of the modes $I_{(2\ell +1)m:(2\ell +1)n}$ have the form

(4.25)\begin{equation} \left. \begin{gathered} \begin{bmatrix}u_h\\ w_h\end{bmatrix} \propto \sigma_{n/m}\begin{bmatrix}\dfrac{n}{m} \cos(2\ell+1)m{\rm \pi} x \: \sin(2\ell+1)n{\rm \pi} z\\ -\sin(2\ell+1)m{\rm \pi} x \: \cos(2\ell+1)n{\rm \pi} z\end{bmatrix} \cos\sigma_{n/m}t, \\ \begin{bmatrix}p_h\\ \theta_h\end{bmatrix} \propto \begin{bmatrix} \dfrac{\lambda_{n/m}}{(2\ell+1)m{\rm \pi}} \sin(2\ell+1)m{\rm \pi} x \: \sin(2\ell+1)n{\rm \pi} z\\ \sin(2\ell+1)m{\rm \pi} x \: \cos(2\ell+1)n{\rm \pi} z\end{bmatrix} \sin\sigma_{n/m}t. \end{gathered} \right\} \end{equation}

Figure 12 shows the vorticity and temperature deviation associated with the first few partial sums in the series expansion of $M_{5:3}$ in (4.24), as well as of $M_{5:3}$, illustrating the rapid convergence of the series (4.23) and, thus, of (4.24). It is difficult to visually distinguish the primary fields, such as the temperature deviation, of $M_{m:n}$ from those of $I_{m:n}$. The contribution of terms beyond the first term in the expansion (4.24) is more visible for the vorticity due to the larger factor, $(2k+1)^{-2}$ vs. $(2k+1)^{-3}$, resulting from the term-wise differentiation of the velocity field.

Figure 12. Plots of the vorticity $\eta _r$ and temperature deviation $\theta _r$ associated with the partial sums in the series expansion of $M_{5:3}$ in (4.24), for the number of indicated terms.

The expansion of forced responses in terms of odd harmonics of the sine wave can be obtained from the sine Fourier expansion

(4.26)\begin{equation} S_1(\zeta-0.5) ={-}\frac{1}{4}S'_2(\zeta) = \frac{8}{{\rm \pi}^2}\sum_{k\ge0}\frac{({-}1)^k}{(2k+1)^{2}}\sin(2k+1){\rm \pi}\zeta. \end{equation}

Expressing the linear spline as a superposition of $m$ shifted triangular waves $S_1(\cdot -\zeta _j)$,

(4.27)\begin{equation} f_{n/m}(\zeta)=\sum_{j=1}^m a_j S_1(\zeta-\zeta_j), \end{equation}

where $\zeta _j=(2j-1)/(2m)$, $j=1,\ldots ,m$ and $a_j={\pm }0.5$ are coefficients whose signs depend on $j$, $m$ and $n$ in a non-trivial way (see Appendix E), and substituting (4.26) into (4.27) leads to

(4.28)\begin{equation} f_{n/m}(\zeta) = \frac{4}{{\rm \pi}^2} \sum_{k\ge0}\frac{s_k}{(2k+1)^2}\sin(2k+1){\rm \pi}\zeta, \end{equation}

with

(4.29)\begin{equation} s_k = 2\sum_{j=1}^m a_j\sin(2k+1){\rm \pi}\zeta_j, \quad k\ge0. \end{equation}

The expression (4.29) amounts to a standard discrete sine transform of the coefficients $a_j$ (the discrete sine transform denoted DST-IV in Britanak, Yip & Rao (Reference Britanak, Yip and Rao2007)). The closed form expression

(4.30)\begin{equation} s_k = ({-}1)^k \sec\left(\frac{2k+1}{2}\right)r{\rm \pi}, \end{equation}

with $r=n/m$, is obtained and illustrated for the case $m=11$, $n=4$, $k=0$ in Appendix F. Note that $s_k$ remains finite as long as $(2k+1)r$ is not an odd integer, i.e. $r$ is not rational with $m$ and $n$ both odd. Using a density argument, the expression (4.30) extends to irrational values of $r$ (in which case the sequence $\{s_k\}_{k\ge 0}$ never repeats), so that

(4.31)\begin{equation} f_r(\zeta) = \frac{4}{{\rm \pi}^2}\sum_{k\ge0} \frac{({-}1)^k}{(2k+1)^2} \sec\left(\frac{2k+1}{2}\right)r{\rm \pi} \sin(2k+1)\zeta{\rm \pi}. \end{equation}

The waveform function is verified to be odd, $f_r(-\zeta )=-f_r(\zeta )$. The relation $(-1)^k\sin (2k+1)\zeta {\rm \pi}=\cos (2k+1)(\zeta -0.5){\rm \pi}$ emphasizes the parity of $f_r$ around $\zeta =0.5$. Moreover, the direct verification that $f_r$ indeed satisfies (4.12) in the Fredholm alternative A1 reduces to a trivial trigonometric exercise, once the right-hand side of (4.12), $2x=S_1(x-0.5)$, is expanded for $|x|\le 0.5$ using (4.26). While the spectral approach employed in Maas & Lam (Reference Maas and Lam1995) and Wotherspoon (Reference Wotherspoon1995), valid for any $r$ different from resonant values $r=n/m$ with $m$ and $n$ odd integers, can be used to derive the infinite series expansion (4.31) directly from the functional relation (4.12), evaluating the series requires truncation to a finite number of terms and it is not clear how the truncation depends on $r$. On the other hand, the spatial approach leads to a linear system that provides an exact solution $f_r$, but only for rational $r=n/m$ with $m$ and $n$ of opposite parities. For irrational $r$, this system becomes unbounded and again a truncation, now to a finite system, must be made.

The value of (4.31) at $\zeta =0.5$ is related to a secant zeta function introduced in Lalín, Rodrigue & Rogers (Reference Lalín, Rodrigue and Rogers2014) that motivated the recent study in Welfert (Reference Welfert2020), where it was shown that $f_{1/\sqrt {2}}(0.5)=1$, as observed in figure 8.

Replacing the quantity $2x$ appearing in the components $w_r$, $p_r$ and $\theta _r$ of the forced response (4.15), hereby denoted $R_f$, by its expression in terms of the function $f_r$ from (4.12), and then substituting (4.31), shows that $R_f$ can be decomposed as the superposition of responses $R_{f_k}$,

(4.32)\begin{equation} R_f = \frac{4}{\rm \pi}\sum_{k\ge0} \frac{({-}1)^k}{2k+1}\, R_{f_k}, \end{equation}

of the sequence of linear forced systems

(4.33)\begin{equation} \left. \begin{gathered} \partial_t u_k ={-}\partial_x p_k + z\alpha\cos(2k+1){\rm \pi} x \sin\omega t,\\ \partial_t w_k ={-}\partial_z p_k + \theta_k,\\ \partial_t \theta_k ={-}w_k, \\ \partial_x u_k + \partial_z w_k = 0. \end{gathered} \right\} \end{equation}

An explicit expression of the components $(u_k,w_k,p_k,\theta _k)$ of $R_{f_k}$ is given in Appendix G. The pointwise convergence of (4.32) away from characteristics inside the cavity also follows directly from the Fourier cosine expansion

(4.34)\begin{equation} S_0(x)=\frac{4}{\rm \pi}\sum_{k\ge0} \frac{({-}1)^k}{2k+1}\,\cos(2k+1){\rm \pi} x, \end{equation}

of the 2-periodic square wave $S_0$, whose value is $(-1)^\ell$ for $x\in (\ell -0.5,\ell +0.5)$, $\ell \in \mathbb {Z}$, and, in particular, $1$ for $x\in (-0.5,0.5)$. The magnitude of the velocity components and temperature deviation of $R_{f_k}$ itself decreases inversely proportionally to $k$ (quadratically for the pressure); see (G5). In other words, the response to (4.33) vanishes for spatially highly oscillatory forcings. As a result, (4.32) converges absolutely, whereas (4.34) does not.

Figures 13 and 14 show the partial sums of (4.32) for the vorticity $\eta$ and temperature deviation $\theta$, respectively, with an increasing number of terms, for the same $\sigma _r^2=\omega ^2$ values as those from figures 3 and 4. Note that in order to avoid $\sec (2k+1){\rm \pi} r=\infty$, the partial sums at $\sigma _r^2=0.1^+, 0.5^+, 0.9^-$ were computed using $\sigma _r^2$ values slightly detuned by a relative amount equivalent to ten times the machine epsilon ($10\times 2^{-52}$). Adding terms to the partial sums improves the approximation. The convergence depends on whether $\sec (2k+1){\rm \pi} r/2$ is large or not. For example, there is a big change for $\sigma _r^2=0.4$ when adding the fifth term and for $\sigma _r^2=0.6$ when adding the sixth term. This is particularly noticeable when $\sigma _r$ is very close to resonances. These near resonances occur when $(2k+1)r\approx 2\ell +1$, i.e. for every $k=m$ terms if $r\approx n/m$ with $m$ and $n$ both odd. For $\sigma _r^2=0.9^-$, ($r\approx 1/3=n/m$), the changes occur at terms $2$ ($k=1$, $\ell =0$), $5$ ($k=4$, $\ell =1$), etc.

Figure 13. Plots of the vorticity $\eta _r$ associated with the partial sums in the series expansion of $R_f$ in (4.32), at $\sigma _r^2$ and a number of terms as indicated.

Figure 14. Plots of the temperature deviation $\theta _r$ associated with the partial sums in the series expansion of $R_f$ in (4.32), at $\sigma _r^2$ and a number of terms as indicated.

For the most part, the comparison between the linear inviscid theory using 100 terms in the partial sums shown in figures 13 and 14 shows excellent agreement with the numerical results at the largest $R_N=10^7$ shown in figures 3 and 4.

5.1. Inviscid versus viscous global enstrophy and a reduced-order viscous model

The global mean enstrophy $\langle \mathcal {E}\rangle _\tau$, defined by (3.1), of the linear inviscid responses is now compared with that of the viscous responses obtained at increasing values of $R_N$ from figure 2$(a)$. The inviscid enstrophy is given by (for details, see Appendix H)

(5.1)\begin{equation} \langle\mathcal{E}\rangle_\tau = \frac{\alpha^2}{r^2\sigma_r^2} \sum_{k\ge0} \lambda_k^{{-}2} \left\{1+\frac{1+r^2}{2}\left[(1+r^2)\sec^2\lambda_k -(3-r^2) \frac{\tan\lambda_k}{\lambda_k}\right]\right\}. \end{equation}

The Dirichlet series (5.1) diverges at resonances which, for the square cavity, occur for every rational $r=n/m$ with $m$ and $n$ both odd integers. Near a specific resonance, $\cos \lambda _k\approx 0$ for some $k$, (5.1) shows that the scaled global inviscid enstrophy, $(r/\alpha ^2)\langle \mathcal {E}\rangle _\tau ={O}(r^{-3}\sigma _r^{-2}[1+r^2]^2) ={O}([r+1/r]^3)$, is invariant under the transformation $r\leftrightarrow 1/r$, i.e. $\sigma _r^2\leftrightarrow 1-\sigma _r^2$. Figure 15$(a)$ illustrates this balance in the response magnitude between the $\sigma _r^2<0.5$ and $\sigma _r^2>0.5$ regimes. The graph is obtained by truncating the series (5.1) after 10 terms ($k=9$). This truncation only retains a finite subset of rational values $r=n/m$ at which the inviscid response becomes unbounded. However, truncating the infinite series after a fixed number of terms independent of $r$ allows for a larger set of resonances to be identified in the lower forcing frequency regime, $\sigma _r^2<0.5$, i.e. $r>1$, compared with the higher frequency regime, $\sigma _r^2>0.5$, because of the higher frequency in the variations of $\cos \lambda _k$. In contrast, truncating after $1+{O}(1/\sqrt {r})$ terms results in a comparable representation of the spatial variations of $\eta$ in the $x$ direction at $\sigma _r^2\approx 1$ and in the $z$ direction at $\sigma _r^2\approx 0$, making the graph and the various resonances visually symmetric approximately $\sigma _r^2=0.5$; see figure 15$(b)$. This modification in the number of terms increases the number of resonances detected by the truncated series in the $\sigma _r^2>0.5$ regime. For example, it decreases the frequency gap between consecutive resonances at $\sigma _r^2\approx 0.8$, with an inverse effect in the $\sigma _r^2<0.5$ regime, as can be observed from the increased gap around $\sigma _r^2\approx 0.2$. Increasing the proportionality constant of the ${O}(1/\sqrt {r})$ term allows for the representation of more oscillatory components throughout the $\sigma _r^2$ parameter inertial range, leading to many more resonances at $r=n/m$ with larger $n$ and $m$, including those closer to $\sigma _r^2=0.5$ ($r=1$), and defining an inviscid response ‘envelope’, shown in figure 15$(c)$. This figure also includes the synchronous ($\sigma _r^2=\omega ^2$) response curves from DNS at various $R_N$ from figure 2$(a)$, adjusted for viscous effects, $\delta (r/\alpha ^2)\langle \mathcal {E}\rangle _\tau$, where

(5.2)\begin{equation} \left. \begin{gathered} \delta = 0.8{\rm \pi}^2R_N^{{-}1/2} (0.25\lambda^{1.5}+0.75\lambda^6)^{{-}1}, \\ \lambda = 4\omega^2(1-\omega^2) =4\sigma_r^2(1-\sigma_r^2)=\left(\frac{2r}{1+r^2}\right)^2. \end{gathered} \right\} \end{equation}

The ${O}(R_N^{-1/2})$ correction leads to the asymptotic collapse (away from lower-order resonant conditions) of the response curves from figure 2$(a)$ onto a single curve as $R_N$ increases, with the proportionality constant calibrated by the response at $\omega ^2=0.5$ ($\lambda =1$). The factor in $\delta$ inside the parentheses is a correction away from the $r=1$ case, in the spirit of the modal viscous reduced-order model used by Yalim et al. (Reference Yalim, Welfert and Lopez2019a), to obtain an approximate critical instability curve for the linearly stratified base state subjected to purely vertical oscillatory forcing of the container. The $\lambda ^{1.5}$ term controls the shape of the curves at $\omega ^2\approx 0$ or $1$, while the $\lambda ^6$ term controls their relative ‘flatness’ in the intermediate ranges around the values $0.25$ and $0.75$. Since $\lambda \le 1$, this mode-dependent correction factor is larger at $\omega ^2\approx 0$ and $1$, where the differences between the horizontal and vertical spatial scales of the responses are large, associated to the nearly vertical/horizontal characteristic directions, and for which the viscous boundary layers have a larger impact.

Figure 15. Graph of the enstrophy $\langle \mathcal {E}\rangle _\tau$ from (5.1) scaled by $\alpha ^2/r$, evaluated at $10^6$ values $\omega ^2\in (0,1)$ using $(a)$ $1+9=10$ terms, $(b)$ $\lfloor 1+9/\sqrt {r}\rfloor$ term(s) or $(c)$ $\lfloor 1+99/\sqrt {r}\rfloor$ term(s) in the series, where $\lfloor x\rfloor$ represents the integer part of $x$. Also shown in $(c)$ are the scaled curves obtained from DNS for the different values of $R_N$ from figure 2 using the same colours, and a viscous correction factor $\gamma$ from (5.2). The supplementary movie 4 shows a frequency sweep of 9999 inviscid solutions uniformly distributed between $\omega ^2\in (0,1)$ at phase 0 for $\eta _r$ and phase ${\rm \pi} /2$ for $\theta _r$.

The 9999 frames in the supplementary movie 4 correspond to snapshots of $\eta _r$ and $\theta _r$ at phases $0$ and ${\rm \pi} /2$, respectively, of the linear inviscid response flows at square forcing frequencies $\omega ^2\in (0,1)$. The animation provides a visual catalog of the Fredholm alternative A1 responses obtained from the truncated series (5.1) with $\lfloor 1+9/\sqrt {r}\rfloor$ term(s), where $\lfloor x\rfloor$ represents the integer part of $x$.

5.2. Flow field comparisons

The nonlinear viscous DNS responses from figures 3 and 4 at $R_N=10^7$ are now compared with the inviscid responses predicted by the analysis from § 4 at similar forcing frequencies $\omega$. These comparisons are shown in figures 16 and 17. All the inviscid responses shown were obtained using the discrete approach in physical space, described in Appendix B, rather than using the truncated spectral representations discussed in § 4.3. Responses associated with quadratic irrationals $r$ (e.g. for $\sigma _r^2=0.3,0.4,0.6$ and $0.7$) were obtained via approximations by a high-order convergent of appropriate continued fractions, similar to (4.16). Inviscid responses associated with rationals $r$ obtained from a lower-order convergent, which correspond to slightly detuned responses, are also included for comparison with viscous results obtained at smaller $R_N$. For example, the case $r=22/27$ in the vicinity of $\sigma _r^2=0.6$ is obtained from $r=4/5$ using the iteration

(5.3)\begin{equation} \begin{bmatrix} n\\m \end{bmatrix} \leftarrow \boldsymbol{\mathsf{R}}\begin{bmatrix} n\\m \end{bmatrix}, \quad \boldsymbol{\mathsf{R}}=\begin{bmatrix} 3 & 2 \\ 3 & 3 \end{bmatrix}, \end{equation}

while the case $r=27/22$ in the vicinity of $\sigma _r^2=0.4$ is obtained from $r=5/4$ using a similar iteration based on $\boldsymbol{\mathsf{R}}^T$ instead. Note that the ratio of coefficients of the right eigenvector of $\boldsymbol{\mathsf{R}}$ is $r=\sqrt {2/3}$, while the same ratio for the left eigenvector is $r=\sqrt {3/2}$, values that are associated with the exact forcing frequencies.

Figure 16. Vorticity $\eta$ from DNS at $R_N=10^7$ $(a{,}c{,}e{,}g{,}i{,}k{,}m{,}o{,}q)$ and $\eta _r$ from the inviscid theory $(b{,}d{,}f{,}h{,}j{,}l{,}n{,}p{,}r)$, all at phase $0$, at the indicated squared forcing frequencies $\omega ^2$ (with the corresponding $r$ values indicated). The inviscid responses at $\sigma _r^2=0.3$, $0.4$, $0.6$ and $0.7$ are fractal and are obtained as limits of sequences with rational $r$.

Figure 17. Temperature deviation $\theta$ from DNS at $R_N=10^7$ $(a{,}c{,}e{,}g{,}i{,}k{,}m{,}o{,}q)$ and $\theta _r$ from the inviscid theory $(b{,}d{,}f{,}h{,}j{,}l{,}n{,}p{,}r)$, all at phase ${\rm \pi} /2$, at the indicated squared forcing frequencies $\omega ^2$ (with the corresponding $r$ values indicated). The inviscid responses at $\sigma _r^2=0.3$, $0.4$, $0.6$ and $0.7$ are fractal and are obtained as limits of sequences with rational $r$.

For ease of comparison, the viscous (DNS) results at large $R_N=10^7$ from figures 3 and 4 have been reproduced in figures 16$(a{,}c{,}e{,}g{,}i{,}k{,}m{,}o{,}q)$ and 17$(a{,}c{,}e{,}g{,}i{,}k{,}m{,}o{,}q)$, albeit at the fixed phases of the forcing, $0$ for the vorticity $\eta$ and ${\rm \pi} /2$ for the temperature deviation $\theta$. For rational $r=n/m$ with sufficiently small $m$ and $n$, the DNS are very similar to the results predicted from the inviscid analysis shown in figures 16$(b{,}d{,}f{,}h{,}j{,}l{,}n{,}p{,}r)$ and 17$(b{,}d{,}f{,}h{,}j{,}l{,}n{,}p{,}r)$. For quadratic irrational $r$, such as $r=\sqrt {7/3}$, $\sqrt {3/2}$, $\sqrt {2/3}$ and $\sqrt {3/7}$, the DNS results at $R_N=10^7$ are unable to fully capture the fractal nature of the inviscid response. For finite $R_N$, the viscous effects smear out the small-scale details. The DNS responses appear more like the inviscid responses obtained at nearby rational values of $r=n/m$ with moderate $m$ and $n$, especially for $\eta$.

Over long periods of time, even small detunings in the forcing frequency can lead to sizeable phase differences in the responses and peaks in an inviscid response at one forcing frequency that may correspond to a minimal response at a nearby forcing frequency. In the viscous DNS responses, viscous detuning at small forcing amplitudes may lead to responses that are not perfect standing waves. Figure 18 illustrates such behaviour at $\omega ^2=0.6$ and $R_N=10^7$. At phases $0$ and ${\rm \pi}$ (corresponding to fractions 0 and 0.50 of a forcing period for $\eta$) or ${\rm \pi} /2$ and $3{\rm \pi} /2$ (i.e. fractions $0.25$ and $0.75$ of a forcing period for $\theta$) the DNS response, already shown in figures 16 and 17, respectively, matches that expected from the inviscid prediction at that forcing frequency. However, at other phases, structures associated with inviscid responses at nearby $\sigma _r^2$ corresponding to $r=n/m$ with relatively small $m$ and $n$ appear to dominate. In particular, the structure associated with the resonant response at $r=9/11$ becomes clearly visible at phases where a standing wave response would be expected to vanish. Which inviscid response is identified in the DNS at a particular phase depends on the relative energies of the various responses at nearby (i.e. viscously detuned) forcing frequencies in relation to viscous effects as measured by $R_N$.

Figure 18. (g–v) Snapshots of $\eta$ and $\theta$ from the DNS limit cycle at $R_N=10^7$ and $\omega ^2=0.6$ over one forcing period (see supplementary movies 1 and 2), while (a–f,w–z,aa,ab) show inviscid standing wave responses at maximal amplitude at $r$ and $\sigma _r^2$ values as indicated.

5.3. Horizontal profile comparisons

This subsection compares the horizontal profiles at $z=0$ of the vorticity $\eta$ and temperature deviation $\theta$ from the DNS at $R_N=10^7$ with the linear inviscid theory. The comparison is performed at the maximal phase for the respective fields. In order to do this, the corresponding inviscid $\eta _r$ and $\theta _r$ need to be evaluated from the waveform function $f_r$. How these are evaluated depends on whether the response is forced or resonant.

For the forced responses, (4.15) results in

(5.4)\begin{equation} \alpha^{{-}1}\begin{bmatrix} \eta_r \\ \theta_r \end{bmatrix} = (1-\sigma_r^2)^{{-}1} \begin{bmatrix} \sigma_r^{{-}1}(\,f_r^\prime(x)-\sigma_r^2) \\ f_r(x)-x \end{bmatrix}, \end{equation}

with $\eta _r$ affinely related to the derivative of $\theta _r$. For the resonant responses, (4.22) leads to

(5.5)\begin{equation} \begin{bmatrix} \eta_r \\ \theta_r \end{bmatrix} \propto 2 \begin{bmatrix} \sigma_r^{{-}1}f_r^\prime(x) \\ f_r(x) \end{bmatrix}, \end{equation}

with $\eta _r$ now directly proportional to the derivative of $\theta _r$.

The profiles from the DNS and linear inviscid theory are shown in figure 19. The resonant responses are in figures 19(a,b) and 19(i,j) and the forced responses in figures 19(c,d) through 19(g,h). For $\omega ^2=\sigma _r^2=0.5$ with $r=1/1$, the response is resonant and $f_r(x)=4x(1-|x|)$ for $|x|\le 0.5$ from (4.20). This waveform function $f_r$ is piecewise quadratic. Likewise, for $\omega ^2=\sigma _r^2=0.9$ with $r=1/3$, the response is resonant and (4.21) yields the piecewise quadratic $f_r(x)=4(3x)(1-3|x|)$ for $|x|\le 1/6$, and is extended by symmetry and periodicity for $|x|>1/6$. For the forced response at $\omega ^2=\sigma _r^2=0.8$ with rational $r=1/2$, $f_r(x)=|x+0.25|-|x-0.25|$ for $|x|\le 0.5$. This piecewise linear $f_r$ is obtained either by evaluating the spectral series (4.31) for $f_{1/2}$ in closed form, or by solving the system (4.14) and collecting the linear pieces.

Figure 19. Horizontal profiles of the (scaled) vorticity $\eta$ and temperature deviation $\theta$ at $z=0$ and $\omega ^2$ as indicated. The black profiles correspond to the inviscid prediction based on the expression (4.22) (for a,b,i,j) or (4.15) (for c–h). The waveform function $f_r$ used in these expressions is obtained either explicitly for the cases $\omega ^2=0.5$, $0.8$ and $0.9$ associated with rational values of $r$, or via the series (4.31) truncated to the indicated number of terms for the cases $\omega ^2=0.6$ and $0.7$ associated with irrational values of $r$. An appropriate factor is used to scale the resonant cases (eigenmodes; a,b,i,j) in order to facilitate the comparison with the DNS response obtained at $R_N=10^7$ from figure 5 (here shown in cyan).

For the forced responses at $\omega ^2=\sigma _r^2=0.6$ and $0.7$, for which $r$ is irrational, the evaluation of $f_r$ is no longer straightforward. The inviscid profiles that are illustrated in figures 19(c,d) and 19(e,f) were obtained by truncating the series representation of $f_r$, equation (4.31), as indicated in the figure. The number of terms retained is determined such that the resulting fields capture the oscillatory behaviour of the DNS responses without introducing high spatial frequency content due to large coefficients $\sec [(2k+1)/2]{\rm \pi} r$ for arbitrarily large $k$ when $r$ is irrational. While the truncation captures the oscillations in the DNS, it overestimates their amplitude. This is partly due to the inviscid theory not having viscous damping.

In order to improve the comparison, filtering techniques, such as the use of mollifiers (Tadmor Reference Tadmor2007), may be applied to the series representation of $f_r$. The use of a mollifier softens the effects of the truncation and provides a much better quantitative account of the viscous results in the DNS. This is illustrated in figure 20, where the mollifier $\exp (-[(2k+1){\rm \pi} /40]^2)$ is used for $\omega ^2=\sigma _r^2=0.6$ and the mollifier $\exp (-[(2k+1){\rm \pi} /70]^4)$ is used for $\omega ^2=\sigma _r^2=0.7$. The factors 40 and 70 are fitting parameters, on the order of $R_N^{1/4}\approx 56$. These factors are consistent with the $R_N^{-1/2}$ factor appearing in the viscous dissipation correction (5.2) at the global enstrophy level. The power $2$ in the first mollifier corresponds to the Gaussian filtering associated with Laplacian dissipation. The power $4$ in the second mollifier corresponds to a biharmonic filtering, known for its scale selectivity. It is used extensively in ocean models (Griffies Reference Griffies2004), and its tendency to overshoot large changes in the data (Delhez & Deleersnijder Reference Delhez and Deleersnijder2007) is a feature exploited here to better capture the DNS results.

Figure 20. Horizontal profiles of the (scaled) vorticity $\eta$ and temperature deviation $\theta$ at $z=0$ similar to figures 19(c,d) and 19(e,f), but using a mollifier rather than hard truncation of the spectral series (4.31) of the waveform function $f_r$ to account for viscous effects.

6.1. Summary

A static vertically stratified fluid, with temperature $T=z$, in a rectangular container is viscously stable. In order to have a sustained non-trivial response flow, the system needs to be continuously forced. The nature of the response flow depends on the nature of the forcing. Here, we consider the effects of subjecting the container to small horizontal oscillations that are harmonic in time, leading to a horizontally modulated effective gravity. The static stratified state is not an equilibrium of the forced system, and the response flow is non-trivial for all non-zero forcing frequencies $\omega$ and amplitudes $\alpha$. The flows obtained via DNS of the Navier–Stokes–Boussinesq equations are synchronous with the forcing and retain the spatio-temporal symmetries of the forced system. As viscous effects are reduced (by increasing $R_N$, the buoyancy frequency relative to the viscous time scale), these flows tend to be standing waves with spatial structure that is increasingly less regular, with nearly piecewise constant or piecewise linear vorticity distributions, depending on the forcing frequency.

The high $R_N$ flows are well described in the inviscid limit by a perturbation analysis of the unforced equilibrium using the forcing amplitude as the small perturbation parameter. The forcing term in the full nonlinear problem, $\alpha T \sin \omega t$, reduces to $\alpha z \sin \omega t$ in the first-order perturbed system, which is linear and non-homogeneous with Dirichlet boundary conditions. This non-homogenous linear boundary value problem is reduced to a Poincaré equation for the temperature deviation, which is independent of the forcing.

The Poincaré equation is solved by a general standing wave ansatz for the temperature deviation, motivated by the DNS results. Substituting this into the first-order perturbation system leads to corresponding velocity components and pressure. Enforcing symmetries and boundary conditions on these then leads to a system of functional equations for the waveform function $f$, one of which depends on the forcing. The solutions depend on whether the forcing resonates or not with the intrinsic modes of the static vertically stratified fluid in the container, as categorized by a Fredholm alternative for $f$. Forced responses are obtained at forcing frequency $\omega$ such that $\omega ^2=1/(1+r^2)$, with $r$ being irrational, or $r=n/m$ being rational with integers $m$ and $n$ of opposite parities. The linear inviscid system has then a unique response flow (alternative A1) which scales with the forcing amplitude $\alpha$. For rational $r=n/m$ with $m$ and $n$ of opposite parities, a linear spline $f$ can be determined by solving a non-homogeneous linear system. The resulting response has piecewise constant vorticity which is discontinuous across characteristic lines and forms a regular harlequin pattern. For irrational $r$, the waveform function $f$ is better expressed in the form of a spectral Dirichlet series. At least for quadratic irrationals $r$, the sequence of rational cumulants of the continued fraction representation of $r$ leads to self-similar properties of $f$ and the harlequin pattern of the vorticity becomes fractal.

Resonant responses are obtained at forcing frequencies $\omega ^2=1/(1+r^2)$ for $r=n/m$ with integers $m$ and $n$ both odd (alternative A2). These responses correspond to solutions of the unforced, homogeneous inviscid linear perturbation equations and are intrinsic eigenmodes of the unforced system. The spatial structure of the specific modes excited by the particular forcing is obtained by considering the limit of forced responses as the forcing frequency approaches a resonant frequency. The response flow becomes unbounded, with piecewise quadratic velocity and temperature deviation, and, thus, a piecewise linear vorticity field, and can be represented as a superposition of an infinite set of odd spatial harmonics of smooth intrinsic modes of the unforced system.

Models of viscous dissipation both at the global enstrophy level and at the local level in a horizontal cut via mollifiers were introduced in the inviscid model to account for viscous effects at finite $R_N$. The use of filtering strategies in post-processing is in general necessary when differentiating low regularity solutions.

6.3. Extensions and further remarks

In the present study we have considered the square cavity. However, the results are straightforwardly extended to rectangular cavities of width-to-depth aspect ratio $\varGamma$. The dispersion relation (4.5) then becomes $\sigma _r^2=\varGamma ^2/(\varGamma ^2+r^2)$, with responses according to Fredholm alternatives A1 and A2 depending on $r$ in a similar fashion as when $\varGamma =1$. If the container is not rectangular, for example, due to walls at oblique angles to gravity, or is tilted, the responses may become non-trivial as a result of the peculiar reflections of internal waves, which may lead to internal wave attractors or focusing into edges or corners. While some aspects of this has been previously studied (e.g. Maas & Lam Reference Maas and Lam1995; Wotherspoon Reference Wotherspoon1995), there is still much to learn (e.g. Bajars et al. Reference Bajars, Frank and Maas2013; Pillet et al. Reference Pillet, Ermanyuk, Maas, Sibgatullin and Dauxois2018; Pillet, Maas & Dauxois Reference Pillet, Maas and Dauxois2019), and approaches from the study of how inertial waves reflect in the obliquely rotating cube (Wu et al. Reference Wu, Welfert and Lopez2020) can be adapted for the internal wave problem.

In order to determine whether the low regularity patterns in the 2-D solutions persist in a 3-D container, two cases, corresponding to $\omega ^2=0.2$ and $\omega ^2=0.5$, were simulated at $R_N=10^6$ and $\alpha =1.75\times 10^{-6}$ in a box $[-0.5,0.5]\times [-0.2,0.2]\times [-0.5,0.5]$ with no-slip, insulated walls at $y=\pm 0.2$. The spanwise aspect ratio $0.4$ is the same as that used in the experiments of Benielli & Sommeria (Reference Benielli and Sommeria1998) and the simulations of Yalim et al. (Reference Yalim, Lopez and Welfert2020) for the vertically forced system. The $\omega ^2=0.2$ case is a prototypical example of a forced response in the 2-D setting (Fredholm alternative A1), and the $\omega ^2=0.5$ case corresponds to a resonant response associated with the main peak in the Bode diagram (figure 2) in the 2-D setting (Fredholm alternative A2). Figure 21(a–d) shows the 2-D DNS solutions from figures 3 and 4 for comparison with the figure 21(e–h) showing the vorticity and temperature deviation at the spanwise midplane $y=0$ of the 3-D DNS solutions. For $\omega ^2=0.2$, the solutions are qualitatively the same and quantitatively the 3-D vorticity level is approximately $77\,\%$ that of the 2-D vorticity, while the 3-D temperature deviation is approximately $2\,\%$ larger. For $\omega ^2=0.5$, the solutions are also qualitatively the same, but the 3-D vorticity and temperature deviation levels are both approximately $46\,\%$ those of the corresponding 2-D solution. From the 2-D study, it appears that resonant responses tend to be more affected by viscous effects than the forced responses are. The viscous dissipation in the boundary layers along the spanwise walls of the 3-D container are also responsible for dampening the response. This warrants a more detailed separate investigation. As noted by Benielli & Sommeria (Reference Benielli and Sommeria1998), the effects of friction in parametrically forced stratified confined flows are not well understood, and, hence, neither are those associated with variations of the Prandtl number.

Figure 21. Panels (a–d) show the 2-D $\eta$ and $\theta$ at $R_N=10^6$ and $\omega^2$ as indicated (repeated from figures 3 and 4 for ease of comparison) and panels (e–g) show the 3-D $\eta$ and $\theta$ in the spanwise midplane using the same colour levels as in the 2-D cases (a–d). Panels (i–l) show the corresponding isosurfaces of the 3-D solutions (using a different colour map), with the boundary layer regions clipped.