Oscillatory large-scale circulation in liquid-metal thermal convection and its structural unit

Abstract In Rayleigh–Bénard convection, the size of a flow domain and its aspect ratio $\varGamma$ (a ratio between the spatial length and height of the domain) affect the shape of the large-scale circulation. For some aspect ratios, the flow dynamics includes a three-dimensional oscillatory mode known as a jump rope vortex (JRV); however, the effects of varying aspect ratios on this mode are not well investigated. In this paper, we study these aspect ratio effects in liquid metals, for a low Prandtl number ${{Pr}}=0.03$. Direct numerical simulations and experiments are carried out for a Rayleigh number range $2.9 \times 10^4 \leq {{Ra}} \leq 1.6 \times 10^6$ and square cuboid domains with $\varGamma =2$, $2.5$, $3$ and $5$. Our study demonstrates that a repeating pattern of a JRV encountered at aspect ratio $\varGamma \approx 2.5$ is the basic structural unit that builds up to a lattice of interlaced JRVs at the largest aspect ratio. The size of the domain determines how many structural units are self-organised within the domain; the number of the realised units is expected to scale as $\varGamma ^2$ with sufficiently large and growing $\varGamma$. We find the oscillatory modes for all investigated $\varGamma$; however, they are more pronounced for $\varGamma =2.5$ and $5$. Future studies for large-aspect-ratio domains of different shapes would enhance our understanding of how the JRVs adjust and reorganise at such scaled-up geometries, and answer the question of whether they are indeed the smallest superstructure units.


Introduction
Thermal convection manifests not only in various geo-and astrophysical systems, but also in smaller-scale phenomena ranging from industrial processes to our daily lives such as household heating.Given its importance, natural thermal convection has been the subject of intensive research for over a century (Bénard 1900;Lord Rayleigh 1916).Investigation of thermal convection in low-Prandtl-number fluids (Prandtl numbers Pr 1) is of particular importance for a better understanding of convection on surfaces of stars, where Pr can be as low as 10 −8 (Spiegel 1962;Hanasoge et al. 2016), and, in case of liquid metals, for numerous technical applications, e.g. the advancement of cooling technology (see, e.g., Scheel et al. 2013;Frick et al. 2015;Schumacher et al. 2015;Scheel & Schumacher 2016;Heinzel et al. 2017;Teimurazov et al. 2017;Zürner et al. 2019;Pandey et al. 2022;Zwirner et al. 2022).
Natural thermal convection occurs in a fluid layer due to a temperature difference imposed at its surfaces.Here, the orientation of the fluid layer surfaces with respect to the gravity vector plays an important role (see, e.g., Shishkina & Horn 2016;Teimurazov & Frick 2017;Zürner et al. 2019;Zwirner et al. 2020a;Teimurazov et al. 2021).One of the classical and probably the most intensively investigated configurations of natural thermal convection is Rayleigh-Bénard convection (RBC) (Bénard 1900;Lord Rayleigh 1916;Bodenschatz et al. 2000;Ahlers et al. 2009;Chillà & Schumacher 2012).In RBC, the heated and cooled surfaces are placed orthogonal to the gravity vector, and the fluid layer is heated from below and cooled from above.Thermal expansion causes warm fluid to rise and cool fluid to sink.At sufficiently large Rayleigh number, Ra ≡ Δ 3 /(), the resulting turbulent convective flow self-organises through an inverse energy cascade from small-scale thermal turbulence to large flow structures.Here,  is the isobaric thermal expansion coefficient,  is the kinematic viscosity,  is the thermal diffusivity, Δ is the temperature difference between the heated and cooled surfaces,  is the distance between these surfaces (i.e., the height of the container) and  denotes the acceleration due to gravity.The energy of small scales is directly transferred to large scales via three-dimensional modes, and is different from the classical two-dimensional inverse energy cascade (Ecke & Shishkina 2023;Boffetta & Ecke 2012).
At sufficiently large Ra, the flow is self-organised into a large-scale circulation (LSC), or a turbulent thermal wind, the concept of which is an important ingredient in the heat and momentum transport theory (Grossmann & Lohse 2000, 2001, 2011), and boundary-layer theory for natural thermal convection (Shishkina et al. 2015;Ching et al. 2019;Tai et al. 2021).The resulting flow structures strongly depend on the Rayleigh number Ra, which is a measure of the thermal forcing that drives convection in the system, and on the Prandtl number Pr ≡ /, which describes the diffusive properties of the considered fluid (Ahlers et al. 2009).In addition, the geometric characteristics of the container, especially the shape of the container and, in particular, the aspect ratio  of its spatial length  and height ,  ≡ /, influence the global flow structure and the mean characteristics of the flow (Shishkina 2021;Ahlers et al. 2022).
Turbulent RBC in a cylindrical container with equal height and diameter (aspect ratio  = 1) is the most extensively studied.For containers with  ≈ 1, the principle structure of the LSC can be delineated as follows.There exists a vertical plane (called the LSC-plane), in which the LSC is observed as a big domain-filling roll with two smaller secondary rolls in the corners (Sun et al. 2005a;Ahlers et al. 2009;Chillà & Schumacher 2012), while in the orthogonal vertical plane, the LSC for this geometry of the container is seen as a four-roll structure, with an inflow at mid-height (Shishkina et al. 2014).Not only the LSC is generally unsteady, but also the LSC path can exhibit dynamic behaviour.Thus in containers with  ≈ 1, the LSC can display various modes of periodic or chaotic oscillations which can take the form of sloshing, precession, and torsion (Cioni et al. 1997;Xi et al. 2004;Funfschilling & Ahlers 2004;Sun et al. 2005b;Xi et al. 2006;Brown & Ahlers 2006, 2007;Xi & Xia 2007, 2008;Funfschilling et al. 2008;Zhou et al. 2009;Brown & Ahlers 2009;Sugiyama et al. 2010;Assaf et al. 2011;Stevens et al. 2011;Wagner et al. 2012;Sakievich et al. 2016Sakievich et al. , 2020;;Cheng et al. 2022).The sloshing mode is associated with the motion of the LSC-plane in the radial direction, while the precession and torsion modes are related to the azimuthal motion of the LSC-plane (Cheng et al. 2022;Horn et al. 2022).In the precession mode, the entire LSC-plane drifts in the azimuthal direction, while in the torsion mode, the azimuthal motion of the LSC-plane in the upper half of the container is generally in the opposite direction to the motion of the LSC-plane in the lower half of the container.
In slender containers with the aspect ratio  < 1, a single big-roll structure of the LSC is By contrast, for wide containers with  > 1, the more rolls of the LSC mean the more efficient heat transport (van der Poel et al. 2011(van der Poel et al. , 2012;;Wang et al. 2020), also in highly turbulent cases.For  > 1, the rolls or roll-like structures are attached to each other and aligned in horizontal directions (Hartlep et al. 2003;von Hardenberg et al. 2008;Emran & Schumacher 2015).In the two-dimensional case, the range of possible aspect ratios of particular convective rolls and, hence, the total number of the rolls in any confined domain, are restricted, and there exist quite accurate theoretical estimates for the lower and upper bounds of possible aspect ratios of the rolls (Wang et al. 2020;Shishkina 2021).
For three-dimensional domains, the typical length scales of the self-organised coherent turbulent flow structures are not yet well-studied and their accurate prediction remains an unsolved problem so far.These flow structures can be identified as turbulent superstructures (Stevens et al. 2018;Pandey et al. 2018;Green et al. 2020;Krug et al. 2020;Berghout et al. 2021), since their lifetime is much larger than the free-fall time, and their length scales are generally larger than the typical length scale in RBC, which is the height of the container .Several studies suggest that the characteristic length scale of these coherent turbulent large-scale flow structures increases with growing Ra, see, e.g., Fitzjarrald (1976);Hartlep et al. (2003); Pandey et al. (2018);Akashi et al. (2019); Krug et al. (2020).Depending on the considered parametric space of the Ra-, Prand  in different studies, different preferable length scales of the turbulent superstructures are reported, which are always larger that the container height .Thus the values of order 10 (Busse 1994), or between 6𝐻 and7𝐻 (Hartlep et al. 2003;Pandey et al. 2018;Stevens et al. 2018) were proposed.Although the typical horizontal wavelengths of the turbulent superstructures generally grow with Ra, they tend to decrease with decreasing Prandtl number (Pandey et al. 2018).This fact is pretty remarkable, since decreasing Pr is usually associated with even stronger turbulence and therefore one might expect a certain similarity to the situation when Ra is increased.
Recent laboratory and numerical experiments show that in an intermediate range of moderate aspect ratios,  1.4, the LSC displays a low-frequency oscillatory dynamics (Vogt et al. 2018;Horn et al. 2022;Akashi et al. 2022;Cheng et al. 2022).The precession, torsional and sloshing dynamics of the LSC, which dominates at  = 1, is replaced by a mode which can be described as a jump rope vortex (JRV).In this flow pattern, a curved vortex is formed, which swirls around the cell centre in the direction opposite to the LSC direction, resembling the motion of a swirling jump rope, see figure 1e.This phenomenon was first demonstrated for liquid metal convection in a cylinder with aspect ratio  = 2 (Vogt et al. 2018).Numerical simulations showed that the JRV exists also for a cylindrical container of the aspect ratio  = √ 2 and that the JRV structure is present not only in low-Pr liquid metal convection, but also in water at Pr = 4.8.This has been confirmed in several other experiments and simulations of comparable aspect ratios for both water and liquid metal (Horn et al. 2022;Cheng et al. 2022;Li et al. 2022).
Flow measurements in containers of different shapes such as a cuboid domain with  = 5 (Akashi et al. 2022) showed that the strongly oscillating velocity and temperature fields could also be attributed to the presence of the JRV-like structures.However, instead of only one vortex, four JRVs interlaced in that case.The ends of the JRVs cross perpendicularly at a certain point in space (see figure 1a).Here, the detached (opposite) JRVs oscillate  out of phase, whereas adjacent JRVs do so with a lag of /2.Akashi et al. (2022) demonstrated that the JRVs can form a lattice structure of different vortices, which determines a fundamental flow mode that for the considered combinations of Ra and Pr can dominate the dynamics at moderate, and possibly also at very large aspect ratios.
As such, liquid metals are well suited to investigate the JRV-like flow dynamics.
The objective of the present work is to investigate in more detail the aspect ratio and geometry dependence of the three-dimensional oscillatory JRV-like large-scale circulation in liquid-metal thermal convection.In particular, how increasing aspect ratios result in a lattice of oscillatory flow pattern via the formation of JRVs, starting from the smallest structural building block to that of the more interlaced JRVs at a higher aspect ratio.To this end, we study the LSC dynamics in RBC of liquid metal with Pr = 0.03 in square cuboids with different aspect ratios, which vary from 2 to 5, using both experimental and numerical approaches.

Methods
2.1.Direct numerical simulations Thermal convection under the assumption of the Oberbeck-Boussinesq approximation is described by the following Navier-Stokes, energy, and continuity equations: (2.1) Here,   denotes the substantial derivative,  = (  ,   ,   ) is the velocity vector field,  is the reduced kinematic pressure,  the temperature,  0 = ( + +  − )/2 is the arithmetic mean of the top ( − ) and bottom ( + ) temperatures,   is the unit vector that points upward.The considered domain is a square cuboid with the height  and equal width  and length ,  = , so that the domain aspect ratio equals  ≡ /.The system (2.1)-( 2.3) is closed by the following boundary conditions: no-slip for the velocity at all boundaries,  = 0, constant temperatures at the end-face of the box, i.e.,  =  + at the bottom plate at  = 0 and  =  − at the top plate at  = , and adiabatic boundary condition at the side walls, / = 0, where  is the vector orthogonal to the surface.Equations (2.1)-( 2.3) are non-dimensionalised by using the height , the free-fall velocity    , the free-fall time    , and the temperature difference between the heated plate and the cooled plate, Δ, as units of length, velocity, time and temperature, respectively.The resulting dimensionless equations are solved numerically using the latest version (Reiter et al. 2022(Reiter et al. , 2021) ) of the direct numerical solver (Shishkina et al. 2015;Kooĳ et al. 2018), which applies a fourth-order finite-volume discretisation on staggered grids.Three-dimensional direct numerical simulations (DNS) were performed for square cuboid domains with the aspect ratios  = 2, 2.5, 3, and 5.The utilised staggered computational grids, which are clustered near all rigid walls, are sufficiently fine to resolve the Kolmogorov microscales (Shishkina et al. 2010), see Tables 1 and 2.

Experimental set-up
A schematic of the experimental set-up is presented in figure 2 along with the measuring positions of ultrasound probes.The set-up consists of a cuboid vessel with a base area of  ×  = 200 × 200 mm 2 and a height  = 66 mm, resulting in an aspect ratio  ≈ 3. The top and bottom plates of this vessel are made of copper plates, whereas the side walls are made of polyvinyl chloride (PVC) of 30 mm thickness.This vessel is filled with an eutectic liquid metal alloy GaInSn of Gallium, Indium, and Tin, that serves as the working fluid in the experiment.Thermophysical properties of GaInSn are reported in Plevachuk et al. (2014).In particular, the melting point of GaInSn is 10.5 • C and the Prandtl number equals Pr ≈ 0.03.
The liquid layer enclosed within the vessel is heated from the bottom and cooled from the top by adjusting the temperature of water flowing through channels in the copper plates.The temperature of water in these channels is held constant at set temperatures via two external thermostats.To minimise heat losses, the tubes transporting the hot and cold water and the entire vessel are wrapped in about 30 mm thick insulating foam tubes and additional envelope.Two platinum resistance thermometers (Pt-100) (accuracy of ±0.005) have been utilised to accurately monitor the temperatures of water entering (  ) and leaving (  ) the hot and cold plates, respectively.These temperature readings are essential in measuring the nondimensional convective heat transport, the Nusselt number Nu, expressed as Nu = Φ/ Φ  .Here, Φ  =  2 Δ/ is the conductive heat flux, with  being the thermal conductivity of the liquid metal.Φ =    (  −   ) is the total heat flux exchanged in the set-up, whereas   is the isobaric heat capacity of water, and  is the flow rate of the circulating water determined via an axial turbine flow sensor at the cooling outlet of the set-up.
Prior to measurements, calibrations are performed.To account for the measurement uncertainty and heat losses of the set-up; hose split valves are used to split cold and hot water outlets respectively.One set of a cold and hot pair is used to feed the top plate and the other set to that of the bottom plate while ensuring that the temperature of both the plates remained at a set temperature of 20 • C using the external thermostats.Once the temperature in the plates reaches an equilibrium; an hour-long time series of temperature readings from both sets of thermocouples are recorded.Using the least square method, offsets from each of these thermocouples are extracted which are then used to correct the temperature measurements.This procedure gives a lower threshold of temperature difference attainable for the setup, measurements below Δ 0.22   Each ultrasound transducer is marked with a letter that indicates the distance to the bottom ("T" -close to the top, "M" -matching the middle plane, "B" -close to the bottom) followed by a number.All dimensional distances in (a-c) are given in mm.Blue and red colours indicate the cooled and heated plates respectively.are recorded after the temperature difference between the hot and the cold plates reached a constant value, when the system attains thermal equilibrium.
Principles of Ultrasound Doppler Velocimetry (UDV), a technique widely used for opaque flow diagnostics, are implemented to determine the fluid velocity (Tsuji et al. 2005;Eckert et al. 2007).Nine UDV transducers (TR0805SS, Signal Processing SA) are installed in a direct contact with the fluid.Each of these transducers acquire an instantaneous velocity profiles sequentially along the measuring lines as shown in figure 2 using multiplexing.The velocity measurements are performed with a resolution of about 0.5 mm/s and a sampling frequency of 1 Hz.
For the numerical results, statistical equilibrium or convergence is reached after several hundreds of free-fall time units.Throughout this paper, the length, velocity, and time are made non-dimensional using the cell height , the free-fall velocity    , and the free-fall time unit    ≡ /   , respectively, see equations (2.4).

Phase averaging procedure
To analyse the 3D flow dynamics from the experimental data, the whole field mapping of the velocity flow field is required, which is currently not possible using the UDV techniques.However, this sort of flow-field measurements can be assessed via the numerical techniques.The flow pattern consists of oscillatory coherent structures over several range of scales.To visualise the coherent structures, it is advisable to remove the background turbulent fluctuations using statistical means.Pandey et al. (2018) implemented an averaging method, which was later adopted by Akashi et al. (2022) in the form of a phase averaging algorithm.In this algorithm, one complete oscillation period,   = 1/   , is equally divided into certain (e.g.16) intervals or phases.Averaging of the temperature and velocity field data is carried out within each of these phases.This method reveals the underlying coherent structures in a flow field with high oscillations, such as that encountered in the three-dimensional cellular regime by Akashi et al. (2022).
Vogt et al. (2018) used conditional averaging to showcase the 3D structures of the JRVs in a cylinder.The method of conditional averaging is similar to that of the phase-average process, with the only difference in the choice of the conditioning intervals.In the conditional averaging approach, the intervals for one complete cycle are divided into seven intervals bounded by multiples of standard deviations of the average temperature of the fluid.
In the present study, the phase averaging method is applied to the simulation data, which cover 16 oscillation periods for the cases  = 2.5,  = 3,  = 5, and 8 oscillation periods for the case  = 2. Every oscillation period is divided into 16 phases and each phase is represented by 20 snapshots of all flow fields.Then the corresponding snapshots are averaged within each phase.Finally, the conditional averaging is applied to the flow fields within each phase and for all oscillation periods, which gives a phase-averaged temporal evolution of all flow fields during the period.

Results
The results of all conducted DNS and experiments are summarised in Table 2.For all experimental and numerical data for sufficiently large Ra, an oscillatory behaviour of the LSC was identified.Like in the case of  = √ 2 and  = 2 of a cylindrical container (Vogt et al. 2018), the JRV-like oscillatory structures leave imprints on almost all flow characteristics, for the considered ranges of Ra and  of a cuboid domain.The oscillatory behaviour of the LSC is reflected in temporal evolution of the temperature and particular components of the velocity fields and is also seen in the vertical heat flux temporal evolution.
Once the dominating frequency  0 is evaluated (we will discuss this in more detail later), one can analyse the mean flow dynamics within the time period that lasts   = 1/  0 .For that, the temporal evolution of the flow fields, which are obtained in the DNS, are split into separate periods, according to the dominating frequency  0 , and then a phase-averaged temporal evolution of all flow fields during the period is calculated.
Our DNS for Ra = 10 6 and Pr = 0.03 and two different aspect ratios,  = 5 and  = 2.5, show a very remarkable similarity of the global flow structure and its dynamics.In figure 3, phase-averaged instantaneous temperature distributions in horizontal cross-sections are presented, which are considered at distances  = 0.5 (figure 3  counterclockwise in the time interval [0, 0.5  ] (see supplementary movies), suggesting the presence of oscillatory flow dynamics which periodically changes the flow topology.For fixed values of Ra and the cell height , the spatial length of the convection cell in the case  = 5 is twice larger than in the case  = 2.5 for the same Ra and .Therefore, for any fixed , one can expect a similarity of the flow pattern in the horizontal cross-section at the height  in the case  = 2.5 with the flow pattern in the 1/4 of the area of the horizontal cross-section at the same height  in the case  = 5.Indeed, figure 3 shows that the temperature distribution in the region marked with black dashed lines for  = 5 (figure 3 a-d) is very similar to the temperature distribution in the corresponding cross-sections for  = 2.5 (figure 3 e-h) if considered at the same phase.
To gain more evidence for this similarity, we evaluate the horizontal components of the velocity,   and   , along the lines marked T1 and T2 in figure 3 (c, d) ( = 5) and compare them with the corresponding horizontal components of the velocity along the lines marked T1 and T2 in figure 3 (g, h) ( = 2.5).The temporal evolutions of these velocity components for  = 5 and  = 2.5 are compared in figure 4 for Ra = 10 6 and  = 0.85.One can see that the lower halves of the spatio-temporal velocity maps in figure 4 (a, b), which correspond to the measurements along the lines T1 and T2 within the 1/4 area that is marked in figure 3 (c, d) with the black dashed lines, mimic the spatio-temporal velocity maps in figure 4 (c, d), which correspond to the measurements along the lines T1 and T2 in figure 3 (g, h).
Qualitatively, the signals for  = 5 and  = 2.5 are similar, however, the frequency of the oscillations in the latter case is slightly lower than that in the former case, with six versus five oscillations during the same time interval.Also at  = 5 the signal seems to be less stable than in the case  = 2.5. Figure 5 shows a comparison of the experimental and simulation results for the same  = 3 and Ra = 10 6 .Here a comparison of the temporal evolution of the horizontal components of the velocity is made exactly at the same locations in the DNS and in the experiment.One can see a good qualitative agreement between the experimental and simulation data.However, the dominant frequency obtained in the experiment is slightly higher compared to the frequency evaluated from the simulation data.More precisely, we obtained on average 11 oscillations in the experiment versus 9 oscillations in the DNS for the same time interval.

Three-dimensional cellular flow dynamics
Figure 6 shows phase-averaged streamlines in Rayleigh-Bénard convection for Pr = 0.03, Ra = 10 6 , as obtained in the direct numerical simulations for cuboid domains for all considered aspect ratios  at the beginning ( = 0) and at the middle ( = 0.5  ) of the oscillation period.There are four interlacing JRVs in the case  = 5 (figure 6 a, b).The flow structure resembles a cellular structure which previously was observed in Akashi et al.
(2022) for  = 5 and Ra ≈ 1.2 × 10 5 .There are only two JRVs for the aspect ratios  = 3 (figure 6 c, d), and  = 2.5 (figure 6 e, f ), which is in contrast to a lattice of four JRVs in the  = 5 case.What is more striking is that the JRVs in the  = 3 and 2.5 cells represent a quadrant of the JRV lattice of the  = 5 cell (see also figure 1).Only one vortex is observed in the convection cell with  = 2.This dynamic interplay between changing aspect ratio and organisation of the JRVs highlights the influence of shape and size of the geometry of the container.It also rises an important question as to whether there is a certain hierarchy in these systems when it comes to the reorganisation of JRVs within a container.
A closer look at the 3D flow structure for  = 2.5 shows how the two vortices connect to each other in the central part of the box (see supplementary movies).In this case, the JRVs are connected with two vortices in the upper part and two vortices in the lower part of the  domain.Figure 7 shows the phase-averaged streamlines for the same  = 0.25  , but the connecting vortices in the upper central part of the domain are highlighted with red colour in figure 7 a, while the connecting vortices in the lower part of the domain are highlighted with blue colour in figure 7 (b).Thus colours here reflect the distance  from the bottom plate, i.e., the vertical coordinate of the structure.For clarity, all other streamlines are shown transparently.It is also worth noting that for the same Ra, the JRVs are more stable and better pronounced for  = 2.5 compared to  = 3.
Figure 8 shows in detail the specifics of the spatial-temporal velocity and temperature maps for  = 2.In this case (see figure 6 e, f ) there is only one vortex that rotates in the direction opposite to the LSC direction.Note that for the considered Ra = 10 6 , the LSC is oriented along a vertical wall of the container rather than diagonally.The flow pattern is similar to that obtained for the  = 2 cylinder (Vogt et al. 2018).To examine this similarity, we evaluate at the mid-height, at  = 0.5 (figure 8a), the horizontal component of the velocity   and the temperature along the straight line marked in the figure "M" and along the circle marked "C", which correspond to the lines along the diameter and along the midplane circumference of an inscribed cylinder (as it was considered in Vogt et al. 2018), respectively.
The dominant frequency  0 is visible in the spatio-temporal maps of both velocity (figure 8b) and temperature (figure 8c). Figure 8b  The spatio-temporal maps of the temperature along the circle "C" (figure 8d) also look similar to those in figure 4(a) in Vogt et al. (2018) and indicate the presence of a dominant frequency.It is worth noting that the oscillations of the LSC orientation, which are characterised by the azimuthal angle   , and computed here using the single-sinusoidal fitting method by Cioni et al. (1997) (indicated with the green line in figure 8d) are strong.Although, in cuboid domain, the LSC direction is expected to be more stable compared to that of a cylindrical domain.This is clearly demonstrated by longer time series, which are presented in figure 9.In contrast to the relatively short time interval with only 9 oscillation periods, when the LSC orientation was mostly stable (figure 8), the longer time series reveal quite strong temporal oscillations of the azimuthal angle   (figure 9c).Spatio-temporal maps for both velocity and temperature (figure 9 a, b) show that intervals with relatively stable regular oscillations alternate with intervals with less stable signals.This leads to difficulties for conditional averaging in this case.Therefore, as mentioned in section 2.3, we used averaging for only 8 oscillation cycles for  = 2, whereas for other values of  the averaging was performed for 16 oscillation cycles.

Oscillation frequency
Any periodic oscillation of the JRV has a certain dominant frequency  0 .For each studied case, these frequencies were extracted both from the velocity and temperature time series.The frequencies were non-dimensionalised using the dissipative time scale.Two length scales were used for the dissipative time scale: the height of the domain  and the overall path length of the LSC, following Cheng et al. (2022), is coarsely approximated as 2,  ≡  +  , where  equals to 2, 2.5, 3 and 5/2 for  = 2, 2.5, 3 and 5 respectively.Note that in the latter case  = 5/2, since for the  = 5 container, the flow pattern consists of two JRV building blocks, repeated in both horizontal directions, as it was shown above.The diffusion times are then   =  2 / and    = ( +  ) 2 / with the corresponding diffusion frequencies   = 1/  = / 2 and    = 1/   = /( +  ) 2 respectively.Figure 10 shows the values of the dominant frequencies  0 versus Ra for all considered values of .The data from the present study are shown in red and blue colours, all other data are shown in grey colour.The values that the frequencies take in different flow configurations are presented in table 2. Variation of  0 , normalised with the dissipative time scale (+ ) 2 /, is shown in figure 10a, as a function of Ra.One can see that our new experimental data for  2022) for  = 5 (grey circles) collapse onto one master scaling line.Note that for  = 3 the oscillatory JRV mode occurs at higher Ra compared to the case of  = 5.The so-called roll regime at lower Ra does not have clear dominant frequencies, therefore only data for the cases where an oscillatory mode exists are presented in the figure .A comparison between the experimental and simulation results for  = 3 shows that the oscillation frequency values obtained in the experiment (as it was already shown in figure 5) is slightly higher than that evaluated from the simulation data for  = 3.
The numerically obtained normalized frequencies for  = 3, 2.5, and 2 are very close to each other.For  = 5, which is the only case with two JRV building blocks, the normalized frequency is generally higher than the frequencies for other  (cf.crosses and asterisk at Ra = 10 6 ).The numerically obtained dimensionless frequency for  = 5 at Ra = 10 6 is in very good agreement with the scaling line for  = 5 reported in Akashi et al. (2022) and with the new experimental data for  = 3.For a lower Rayleigh number, Ra = 1.2 × 10 5 , the numerically obtained dimensionless frequency for  = 5 is also in very good agreement with the numerical and experimental data from Akashi et al. (2022).Experimental data for a  = 2 cylinder from Vogt et al. (2018) andCheng et al. (2022) give scaling relations with slightly lower exponent values and for the considered Ra range; they locate slightly below the fitting lines for  = 3 and 5.The frequency value from our  = 2 box simulations is close to that for  = 2 cylinder experimental data from Vogt et al. (2018).The numerical data for all considered  are located between the fitting lines obtained in the experiments for  = 3, 5 box and cylinder  = 2.
To sum up, all the experimental and numerical data show a very similar frequency dependence, as one can see in a  0 versus Ra plot, across all aspect ratios with the frequency normalisation based on the path length .In figure 10b we normalise the frequency  0 with the thermal diffusion time   =  2 /.In that case, without taking into account the spatial length of the vortex path, the deviation between the data points for  = 5 and 2.5 remains the same (as the length scale is the same for these two cases), while the data points for  = 3 move down and the points for  = 2 move significantly up.We conclude that the spatial length of the domain is an important control parameter, which together with the height of the fluid layer determines the relevant length and the scaling relations for the oscillation frequency.

Heat transport
In this section, the effect of the flow dynamics on the heat transport is discussed.The volumeaveraged Nusselt number Nu  can be evaluated from the simulation data as follows: (3.1) where Ω  is a component of the full heat flux vector  ≡ ( − ∇)/(Δ/) along the vertical axis and •  , denotes the the time-volume average.In the experiments, the Nusselt numbers Nu are computed as discussed in section 2.2.The global heat transport scaling across various Ra is shown in figure 11.The flow dynamics does not seem to have any dramatic effect on the heat transport.This is true for all studied aspect ratios.The cases without oscillations are shown in the figure with open symbols.The fitted curve gives a scaling relation Nu = 0.22 × Ra 0.23 , which differs slightly from that reported by Vogt et al. (2021): Nu = 0.166 × Ra 0.25 .However, this difference can be attributed to the difference in the geometry of the cell.
An interesting feature of the JRV regimes is that the Nusselt numbers, which are computed using the phase-average method as discussed in section 2.3, demonstrate an oscillatory behaviour during the JRV cycle.This is demonstrated for one period of oscillation in figure 12. Qualitatively, the oscillatory behaviour of the local vertical heat flux Nu() during the JRV cycle is the same for all considered .It is clear that the Nu shows oscillatory behaviour with distinct peaks of maxima and minima.This sort of behaviour was also reported in previous study of Akashi et al. (2022) for  = 5.However, the amplitude of the oscillations for a certain given Ra is different: it decreases with increasing .
In addition to figure 12 that shows the volume-averaged Nusselt number Nu  during one  oscillation period, we present in figure 13  Figure 14 shows phase-averaged isosurfaces of the full heat transport vector  as obtained in the direct numerical simulations for cuboid domains at the beginning ( = 0) and at the middle ( = 0.5  ) of the oscillation period.The isosurfaces of the full heat transport vector  follow the JRV flow structure at all considered  values (cf.figure 6).
Figure 15 demonstrates -isosurfaces together with the distribution of the magnitude, ||, in the horizontal cross-section at  = 0.5 for  = 5 and  = 2.5.The heat flow is mainly realized in the gaps between the isosurfaces that envelop the JRVs.Thus, the JRVs are not efficient in transporting the heat and are located in the areas of minimum heat flux.
Figure 16 shows the vertical component of the local heat flux, Ω  , at  = 0.5 for  = 5 and  = 2.5.Analogously to figure 3 with the temperature distributions, here one can see a similarity of the Ω  distribution pattern in the case  = 2.5 with the pattern in the 1/4 of the area in the case  = 5 (see supplementary movies).But the resemblance is not complete, probably because of the influence of the sidewalls.How sidewalls affect the movements of vortices is a subject for future study.

Discussion
We have presented a combined numerical and experimental investigation of a liquid metal convection flow in different geometries.The Prandtl number in these investigations is Pr ≈ 0.03 and the Rayleigh numbers are between 2.9 × 10 4 Ra 1.6 × 10 6 .The investigations focus on the influence of the size of the flow domain (via its aspect ratio) on the dominant oscillation modes of the large-scale circulation.Results for four different cuboid domains with varying spatial length-to-height aspect ratios  = 5,  = 3,  = 2.5 and  = 2 were compared with the results of a cylindrical  = 2 cell.
The results show that the oscillations in all aspect ratios investigated are due to the presence of jump rope vortices.A jump rope vortex forms at the centre of the large-scale circulation, and moves analogous to a swirling jump rope.However, the direction of motion of the JRV is opposite to the direction of flow of the LSC.The JRV, which was first discovered in a cylindrical  = 2 convection cell (Vogt et al. 2018), also forms in a square cuboid domain of aspect ratio,  = 2, as demonstrated in this work.The appearance of the JRV is almost identical in both the cylindrical and cuboid domain of same aspect ratio.If a cylinder is numerically cut out from the rectangular cell, the similarity becomes more pronounced, also with respect to the JRV-induced sidewall temperature distribution.
In domains with larger spatial length, the appearance of the JRV vortices changes.For domains with aspect ratios of  = 2.5 and  = 3, the vortices form an orthogonal cross that periodically rotates alternately clockwise and counterclockwise.In a  = 5 cell, a lattice of four JRVs interlace each other, which oscillate in a synchronised manner.Therefore, a key finding of this work is that the JRV is an extremely robust flow feature that adapts and reorganises depending on different aspect ratio of a domain, with ability to form intricate lattice of repetitive flow structures in large aspect-ratio containers.Moreover, our findings further reinforce that the shape of the domain does matter: we encounter the presence of a JRV in a square cuboid of  = 3, whereas, (Cheng et al. 2022) did not find any evidence of a JRV in a cylinder of the same aspect ratio.The frequency of the oscillations show a consistent scaling for the different aspect ratios with a good agreement between numerics and experiment.Slight deviations between the different aspect ratios are likely due to the non-uniform path length of the LSC for the different aspect ratios.
The heat transport scaling relations show only minor (if any) deviations between different flow pattern regimes.The data from the regime close to the onset, convection roll dominated regime and the turbulent JRV regime collapse on a master curve.However, the oscillations of the JRV are clearly visible in the time evolution of the Nusselt number.The frequency of the Nu oscillations is thereby twice as high as that of the JRVs.The maxima of the Nusselt numbers occur when the horizontal velocity components reach a minima during the JRV cycle (see Akashi et al. 2022).
Questions, which are difficult to answer based on both experimental and numerical approaches, are whether the JRV structures have an upper Ra number limit and whether  they are subsequently displaced by other structures as soon as Ra exceeds a certain critical value.In the previous experiments, JRVs detected were stable over two orders of magnitude in Ra (see Vogt et al. 2018).Since the flows in these measurements and simulations are already in a turbulent state, one might expect that the JRV-like oscillatory structures can occur for even larger Ra.It is worth noting that the JRVs occur not only for low Prandtl numbers like that we studied here, but have also been detected in a  = 2 cylinder for water, which has approximately two orders of magnitude higher Pr (Horn et al. 2022).
Our study poses few more questions for future studies that could potentially be investigated, such as: how do the JRVs behave in even larger containers with even higher spatial length domains, and what role do they play in formation of convective turbulent superstructures?The present study suggests that in the case of large , the global structure of the oscillatory mode can be thought of as a lattice of interlaced JRV-like building blocks found for the aspect ratio  ≈ 2.5, repeated spatially.However, such investigations come with their own challenges.Numerical cost increases with the square of the domain aspect ratio, whereas the stabilizing influence of the sidewalls decreases with increasing aspect ratios, giving the flow more degrees of freedom, which results in JRVs that are less stable.This makes intractable the detection of the JRVs by known experimental techniques or by numerical techniques such as conditional averaging.Current ongoing research efforts at the HZDR aim to tackle this problem head-on by experimentally investigating the dynamics of oscillatory liquid metal thermal convection in a square cuboid with a large aspect ratio of  = 25, which is under construction at the time of writing this paper.Supplementary data.Supplementary material and movies are available at ...

Figure 1 :
Figure 1: Phase-averaged streamlines in Rayleigh-Bénard convection for Pr = 0.03, Ra = 10 6 , as obtained in direct numerical simulations for different aspect ratios (a)  = 5, (b)  = 3, (c)  = 2.5 and (d)  = 2 for square cuboid domains (new simulations) and (e, adapted from Vogt et al. (2018), available under a CC BY-NC-ND 4.0) a cylindrical domain with  = 2.These streamlines envelope the oscillating vortex, and the colour scale is according to the vertical velocity component   .Blue (red) colour corresponds to a negative (positive) value of   , indicating a downward (upward) direction of the flow.The structures in the lower aspect ratio cases ( = 2, 2.5 and 3) are building units of the structure formed within the largest aspect ratio case ( = 5).

Figure 2 :
Figure 2: Schematics of the experimental set-up: (a-c) three projections showing (a) the top view and (b,c) two side views and (d) a three-dimensional sketch, illustrating the positions of all ultrasound transducers.Each ultrasound transducer is marked with a letter that indicates the distance to the bottom ("T" -close to the top, "M" -matching the middle plane, "B" -close to the bottom) followed by a number.All dimensional distances in (a-c) are given in mm.Blue and red colours indicate the cooled and heated plates respectively.

Figure 3 :
Figure 3: Phase-averaged snapshots of the temperature at a distance  from the bottom: (a, b, e, f )  = 0.5 and (c, d, g, h)  = 0.85, for different container aspect ratios (a-d)  = 5 and (e-h)  = 2.5, as obtained in the simulations for Ra = 10 6 and Pr = 0.03 (see supplementary movies).The virtual probe lines T1 and T2 (see also figure 4) are indicated with dashed white lines.The black squares indicate the areas that correspond to the areas of the container with  = 2.5.

Figure 4 :
Figure 4: Spatio-temporal velocity maps for Ra = 10 6 at  = 0.85, as obtained in the direct numerical simulations at the virtual probe lines T1 (a, c) and T2 (b, d), for the aspect ratios  = 5 (a, b) and  = 2.5 (c, d).The black dashed lines in (a, b) correspond to the measurements in the cuboid with  = 2.5 (c, d), respectively.

Figure 5 :
Figure 5: Spatio-temporal velocity maps for Ra = 10 6 , as obtained in the direct numerical simulations (left column) and in the experimental measurements (right column), for the aspect ratio  = 3.The numerical data are probed exactly at the same locations where the UDV sensors are located in the experiment: T3 (a, b), B3 (c, d), T2 (e, f ), B2 (g, h), T1 (i, j), B4 (k, l) and M3 (m, n).

Figure 7 :
Figure 7: Phase-averaged streamlines in Rayleigh-Bénard convection for Pr = 0.03, Ra = 10 6 , as obtained in the direct numerical simulations for  = 2.5,  = 0.25  .Vortices in the upper part of the domain are highlighted in (a), vortices in the lower part of the domain are highlighted in (b).For convenience, all other streamlines outside the center of the domain are shown transparent.Colours correspond to the vertical coordinate of the structure .

Figure 8 :
Figure 8: Data for Ra = 10 6 for  = 2, as obtained in the direct numerical simulations.Measurement positions in the midplane at  = 0.5 are shown in (a): the central straight line marked "M" and the circle marked "C" are shown in black and blue colours, respectively.Spatio-temporal (b) velocity and (c) temperature maps along the M-line.(d) Spatio-temporal temperature map along the C-circle.The instantaneous position angle of the LSC is marked with the green line (cf.figure 4 in Vogt et al. (2018), for a cylinder with the same  and Ra).

Figure 9 :
Figure 9: Data for Ra = 10 6 and  = 2, as obtained in the direct numerical simulations.A similar figure to figure 8 (b-d), but the time series is longer.Spatio-temporal (a) velocity and (b) temperature maps along the M-line.(c) Spatio-temporal temperature map along the C-circle.The instantaneous position angle of the LSC is marked with the green line.

Figure 12 :
Figure12: Phase-averaged Nusselt number Nu() during one oscillation period, as it is evaluated from the direct numerical simulations for Ra = 10 6 and different aspect ratios  of a parallelepiped container.

Figure 13 :
Figure 13: Evolution of the phase-averaged Nusselt number Nu() during the one oscillation period, as it is evaluated from the direct numerical simulations for Ra = 10 6 and  = 2 (a),  = 2.5 (b),  = 3 (c) and  = 5 (d).Nu is calculated over the top and bottom plates, over the horizontal cross-section in the middle plane at  = 0.5 and over the entire volume.

Figure 14 :
Figure 14: Isosurfaces of the magnitude of the full heat transport vector  ≡ ( − ∇)/(Δ/), for  = 5 (a, b) and  = 3 (c, d),  = 2.5 (e, f ),  = 2 (g, h), as obtained in the simulations for Ra = 10 6 .The surfaces are coloured by the temperature: blue (red) colour corresponds to the temperature below (above) the arithmetic mean of the top and bottom temperatures.

Figure 15 :
Figure 15: Isosurfaces of the magnitude of the full heat transport vector  are shown together with the || distribution at  = 0.5 for  = 5 (a) and  = 2.5 (b) for Ra = 10 6 .The JRV-like vortex structures are associated with the minimal heat flux.

Figure 16 :
Figure 16: Phase-averaged vertical component of the local heat flux Ω  at  = 0.5 for  = 5 (a, b) and  = 2.5 (c, d) for Ra = 10 6 .The black squares indicate the areas that correspond to the areas of the container with  = 2.5 (see supplementary movies).

Table 1 :
Details on the conducted DNS, including the number of nodes   ,   ,   in the directions ,  and , respectively; the number of the nodes within the thermal boundary layer, N  , and within the viscous boundary layer, N v , the relative thickness of the thermal boundary layer,   /, and the viscous boundary layer  v /; the Kolmogorov microscale, ℎ K , and the relative mean grid stepping, ℎ DNS /ℎ K .
the values of Nu, which are computed over the surfaces: Nu  at the bottom plate, Nu   at the top plate and Nu  over the horizontal cross-section in the middle plane at  = 0.5.Compared to Nu  , for all studied values of , there is a shift between Nu evaluated at the plates (Nu  , Nu   ) and the volume-averaged Nusselt number Nu  .The maxima and minima of Nu, calculated at the horizontal walls, occur always later than they appear in the Nu  evolution.Nu  and Nu   are synchronised with each other.Nu  seems to be less smooth and gives a larger difference between the maximum and minimum values compared to Nu  and Nu   .