Introduction

The so-called liquid bridge (sometimes referred to as LB through this manuscript), namely, a small amount of liquid held by surface tension between two co-axial solid supports displaying rotational symmetry (i.e. simple disks or circular rods) was originally introduced at the end of the 70s as a vehicle to investigate the dynamics of surface-tension driven flows in well-defined conditions (Schwabe and Scharmann 1979; Chun and Wuest 1979). Notably, with such configuration, convection can simply be produced through the application of differential heating along its axis and insights into the emerging flow can be obtained accordingly via direct visualisation of tracer particles dispersed into the considered transparent (work) liquid. This problem has always been widely open and has attracted the attention of distinct research groups with various backgrounds and motivations (Shevtsova 2011). This general interest has led over the years to establishing a common elegant theoretical framework that is currently more or less universally known as the liquid-bridge problem (Lappa 2009).

Although several hundreds of papers have been devoted to this subject (the liquid bridge is indeed the key element of a range of technological processes, Hurle (1994)), we will content ourselves with giving only some flavour of the most important findings and provide a limited set of representative samples for the existing literature (the reader interested in a more thorough discussion may consult the book by (Lappa 2009), where a large cross section of past findings have been described in detail).

Perhaps the simplest way to elaborate a historical perspective on this subject is to start from the remark that, with the progressive rise in its theoretical relevance and suitability as a proper workhorse for the investigation of the fundamental properties of surface-tension-driven (Marangoni) flows, the initially largely “terrestrial” utilisation of this configuration (mostly based on the use of small scale systems to emphasise Marangoni flow in comparison to buoyancy effects) received a continuation with new studies conducted in the weightless conditions provided by sounding rockets (relevant examples along these lines being the works by Monti 2000; Monti et al. 1998; Savino 2001; Schwabe 2002; Schwabe 2005).

Space experiments have gained additional momentum in the recent years owing to the advent of multiple orbiting platforms able to provide a microgravity environment of better quality (and for an ‘extended’ time) with respect to that attainable with Sounding rockets. Such new studies have been spurred, on the one-hand, by the need to move a step forward with respect to the overwhelming earlier research conducted using millimetre-sized liquid bridges and, on the other hand, to gain a better picture of the possible dynamics of Marangoni flow through extensive parametric investigation (not allowed in the frame of the “one-shot” Sounding-rocket based approach). It is known that this type of flow is initially very regular, but it becomes oscillatory and three-dimensional if the applied temperature difference is increased beyond a given threshold. The related critical conditions depend on a variety of influential factors, including (but not limited to) the geometrical aspect ratio AR (ratio of the height and diameter of the liquid bridge), the so-called volume ratio S (ratio of the effective amount of liquid held between the supporting disks and a cylindrical volume having the same height of the liquid bridge and diameter corresponding to the supporting disks) and the temperature of the gas surrounding the free liquid surface.

Though a fairly coherent and exhaustive picture of the richness of possible scenarios for the supercritical states of Marangoni convection in this configuration has been elaborated over the years (see again Lappa 2009), aspects exist that still elude us.

An example of new findings is represented by the MEIS (Marangoni Experiment in Space) series of experiments conducted using the FPEF, i.e. the Fluid Physics Experiment Facility, a Japan Aerospace Exploration Agency (JAXA) subrack located in the Ryutai (fluid) rack inside the Kibo module of the International Space Station (ISS). A synthesis of the main outcomes of such recent studies can be found in Nishino et al. (2015), where the critical Marangoni number for liquid bridges made of different liquids (5 cSt silicone oil and 20 cSt silicone oil) has been reported for a relatively wide interval of aspect ratios (up to \(AR=2.0\)) and an almost cylindrical liquid-gas interface, i.e. \(S\cong 0.95\). These experiments have shown that for a fixed value of S the critical Marangoni number can behave as a non-monotonic function of the aspect ratio. A similar concept also applies to the frequency of the emerging oscillatory flow, which can even display a significant jump when the aspect ratio is varied in the range between \(AR=1.25\) and \(AR=1.5\).

Paralleling these efforts are equivalent studies recently carried out onboard the Chinese space laboratory Tiangong-2. Using a dedicated facility, Kang (2019) have conducted a rich set of space experiments for liquid bridges of 5 cSt silicone oil with different aspect and volume ratios. It has been found that if the system dynamics are examined in the reduced space of parameters accounting for geometrical effects only, i.e. the plane (AR, S), two dominant behaviors can be distinguished on the basis of the magnitude (low or high, respectively) of the frequency of oscillation.

Another example of the renewed interest in these dynamics is the so-called JEREMI experiment. The acronym JEREMI stands for “Japanese European Research Experiments on Marangoni Instabilities”. This ambitious project is part of an agreement between the European Space Agency (ESA) and the Japanese Space Agency (JAXA). Started in 2001, it is now entering its final stage of preparation (Shevtsova 2011; Kuhlmann et al. 2014) and is currently endorsed to fly on the International Space Station (ISS). It will consider both cylindrical and non-cylindrical liquid bridges. In this regard, it may be seen as the successor of MEIS. Like its predecessor, JEREMI will be based on the FPEF facility (Shevtsova 2008). Unlike previous efforts, it however, will target the exploration of “new phenomena” (indeed, this field continues to burgeon and bring surprises to this day).

In particular, an important objective of JEREMI relates to the identification of the cause-and-effect relationships at the root of the so-called “Partice Accumulation Structures” (PAS). It stems from the relatively recent discovery (Tanaka et al. 2006) that, under the effect of Marangoni flow, solid particles initially distributed uniformly in a liquid bridge can demix spontaneously from the surrounding fluid and form highly ordered three-dimensional aggregates. The resulting cluster or pattern formed by particles looks like a spatially extended “closed wire” or “circuit” (having the shape of a flower with several petals). This fascinating structure has been observed to float inside the liquid and rotate in space with constant angular velocity, thereby giving the illusion of a freely-floating “rotating solid body” (Gotoda 2016; Melnikov 2011; Melnikov et al. 2013; Melnikov and Shevtsova 2017; Pushkin et al. 2011; Ueno 2021). As an alternative, a completely different (stationary) pattern can be formed, where the aforementioned one-dimensional particle circuit is replaced by the two-dimensional accumulation planes, which partition the liquid bridge into different sub-regions “like the pieces of a cake cut by a knife” (see, e.g., Schwabe and Mizev 2011; Lappa 2014 for related terrestrial experiments and numerical results, respectively).

Due to such (and other similar) studies, eventually, it has been recognised that these two different attractors, i.e. loci of particle self-organisation simply reflect a well-known dichotomy affecting the possible supercritical spatio-temporal states of Marangoni flow in liquid bridges, i.e. the ability of this flow to produce “travelling waves” or “standing waves”.

Given these premises, the main objective of the present work is the provision of meaningful data for the optimisation of this new series of space experiments. Though the liquid (5 cSt silicone oil) and particles to be used for JEREMI have already been defined, the effective conditions to be considered for the flight (in terms of aspect ratio, volume ratio and applied temperature difference) have not been rigidly fixed yet.

Stripped to its basics, our inquiry scheme envisions a parametric numerical investigation of the response of the liquid bridge to the application of an imposed temperature difference for different values of AR via solution of the balance equations for mass, momentum and energy in their complete time-dependent form. Dispersed particles are tracked separately in the framework of a Lagrangian approach. Conditions for the onset of oscillatory flow (equivalent to the neutral curves, which would be obtained in the context of a linear stability analysis) are determined through extrapolation of the disturbance growth rates obtained via solution of the non-linear governing equations.

Problem Characterisation

As anticipated in the introduction, we consider the instability of thermocapillary flows and ensuing PAS formation in liquid bridges made of a high Prandtl liquid under microgravity conditions.

The system under examination is constituted by a liquid fluid column held between two co-axial solid disks of diameter, D (set equal to 3 cm throughout the present study to mimic an actual experimental setup planned for the JEREMI project), placed at mutual distance H apart and maintained at different temperatures, \(T_c\) and \(T_h\), as indicated in Fig. 1. The interface is represented by a rigid straight adiabatic boundary in which thermocapillary stresses are modeled by imposing a suitable boundary conditions (refer to Flow Dynamics for the details about its formal definition). In these conditions, the flow field is fully characterised by three dimensionless parameters, namely the Marangoni number, \(Ma= -\sigma _T H\Delta T/\rho \nu \alpha\) (here, \(\sigma _T=\partial \sigma /\partial T\), where \(\sigma\) is the interfacial tension, \(\Delta T=T_h-T_c\), while \(\rho\) and \(\nu\) are the fluid density and kinematic viscosity, respectively), the Prandtl number, \(Pr=\nu /\alpha\), where \(\alpha\) is the thermal diffusivity, and the aspect ratio, \(AR=H/D\), of the bridge. In line with the set of experiments planned for execution in space, several values of AR are considered in the present work. When the Marangoni number reaches a critical value, \(Ma_c=f\left( AR, Pr \right)\), the flow becomes unstable and if the conditions are favorable (see e.g., Schwabe et al. (2007)), small particles dispersed in the liquid column can accumulate giving rise to well-defined patterns and structures. Since a microgravity environment is considered, the Stokes number, \(St=2/9d_p^2/H^2\), and the density ratio, \(\xi =\rho _p/\rho\), (here \(\rho _p\) represents the density of the tracers), are the only parameters required to characterise the dynamics of particles. Put simply, the overall problem relating to Marangoni flow and PAS can be completely characterised in terms of five independent dimensionless parameters only, i.e.: Ma, AR, Pr, St and \(\xi\).

Fig. 1
figure 1

Schematic representation of the liquid bridge showing its geometrical constrains and the two heated disks

Full size image

The properties of the liquid and the particles that will be used for the JEREMI experiments are reported in Table 1.

Table 1 Properties of the KF-96L 5 cSt silicone oil at the working temperature of 25 C and characteristics of the Soken Chemical & Engineering Co., Ltd. tracers that will be used for the flight experiments
Full size table

Flow and Particle Dynamics Modeling

Flow Dynamics

The flow motion is governed by the incompressible Navier-Stokes and the energy equations, which in dimensional form read:

$$\begin{aligned} \mathbf {\nabla }\cdot \mathbf{u }= & {} 0,\end{aligned}$$
(1)
$$\begin{aligned} \frac{\partial \mathbf{u} }{\partial t} +\mathbf{u} \cdot \mathbf {\nabla }\mathbf{u }= & {} -\frac{1}{\rho }\mathbf {\nabla }{p} + \nu {\nabla }^2\mathbf{u },\end{aligned}$$
(2)
$$\begin{aligned} \frac{\partial T}{\partial t} +\mathbf{u} \cdot \mathbf {\nabla }{T}= & {} \alpha {\nabla }^2{T}. \end{aligned}$$
(3)

Here, \(\mathbf{u}\) and t are the velocity vector and time, respectively, p is the pressure.

The kinematic viscosity is modeled as a decreasing function of the temperature by adopting the following law, which was specifically determined for the Shin-Etsu KF-96 silicone oil by (Nishino et al. 2015)

$$\begin{aligned} \frac{\nu }{\nu _{25}}=\mathrm {exp}\left( 5.892\frac{25-T}{273.15+T} \right) , \end{aligned}$$
(4)

being \(\nu _{25}\) the kinematic viscosity measured at 25 C, and the temperature T is also given in degree Celsius.

Closure of the mathematical problem also requires the introduction of the tangential stress boundary condition needed to account for the thermocapillary force arising as a consequence of thermally induced surface tension imbalance. Assuming an inviscid surrounding gaseous phase this condition simply reads (see, e.g., (Capobianchi and Lappa 2020))

$$\begin{aligned} \mu \left[ \mathbf {\mathbf {\mathbf {\nabla }}}\mathbf{u } + \mathbf {\nabla }\mathbf{u ^{\mathrm {T}}} \right] \cdot \mathbf{n} = \sigma _T \mathbf {\nabla }_{s}{T}, \end{aligned}$$
(5)

where \(\mu\) is the dynamic viscosity, \(\mathbf{n}\) is the unit normal vector at the liquid-gas interface, and the operator \(\mathbf {\nabla }_{s}\) represents the surface gradient, i.e., the projection of the vector gradient on the interface of the liquid bridge.

Particle Dynamics

Once the flow field has been computed, particles are tracked individually in a Lagrangian framework assuming negligible back influence of the particle on the flow field as well as disregarding their mutual interactions. This is equivalent to the so-called one-way coupling approach. Existing studies in the literature have shown that a judicious utilization of this approximation has to be based on “a priori” evaluation of the degree of dilution of the considered dispersion (Bodnir 2017) and that the threshold value beyond which it ceases to be valid (or it might lead to unphysical results) strongly depends on the considered flow. As an example, assuming the volume fraction of the dispersed phase (\(\phi\)) as a proper measure of the aforementioned level of dilution (defined as the volume globally occupied by the solid particles in a certain sub-region of the considered fluid domain divided by the volume of the considered region itself, namely \(\phi ={N_{part}\pi d_p^3}/{6\Omega }\), where \(N_{part}\) is the number of solid particles and \(\Omega\) is the volume of the considered portion of fluid), Elghobashi (1994) found the one-way coupling to be suitable for isotropic turbulence only for\(\phi <10^{-6}\). For laminar flow, experimental evidence exists that the limiting value is much higher. In particular, for the specific case of time-periodic Marangoni flow in liquid bridges the threshold can be increased to \(\phi \approx 10^{-4}\) (\(\phi\) ratio of the total volume of particles and the total volume of carrier liquid), as demonstrated by the good agreement between numerical simulations carried out in the framework of the one-way coupling simplified model and the experiments by Schwabe et al. (2007). In the light of these arguments, the applicability of the present numerical results to effective experiments in space should be limited to cases for which the amount of particles does not exceed the aforementioned limiting conditions (based on the global volume fraction of the solid phase). Although a possible objection to this way of thinking might be based on the argument that when the particles accumulate in the PAS, the volume fraction does locally (in a restricted neighborhood of the PAS) exceed the \(\phi \approx 10^{-4}\) condition, experimental evidence has demonstrated that interparticle forces or other “back influence” effects are not enough intense to force particles to leave the attractor (thereby making high-fidelity simulation accounting for two- and four-way effects unnecessary, Capobianchi and Lappa 2021).

Our decision to disregard the Saffman lift force, the Faxen corrections and the Basset force in the particle tracking equation also requires a proper rationale. Relevant information about the possibility to neglect the Basset force (also known as history term) can be found in the earlier studies by Hjelmfelt and Mockros (1966) and Lappa (2013) where a relevant criterion was formulated as \(St, St_{\omega }<2/9\xi ^{1/2}\), where St is the classical particle Stokes number and the additional parameter \(St_{\omega }\) (see, e.g., Coimbra and Rangel 2001; Capobianchi and Lappa 2020) reads \(St_{\omega }=\omega _{dim}d_p^2/18\nu =\omega St\) where \(\omega _{dim}\) is the dimensional angular frequency of oscillation of the velocity field (\(\omega _{dim}=2\pi f_c^{dim}\)). With regard to the conditions considered in the present work, \(\xi\) and the particle Stokes number can be evaluated considering the data reported in Table 1. However, while for the particles MX-3000N and ME-200N, \(\xi _{MX-3000N} \cong 0.754\) and \(\xi _{ME-200N}\cong 0.721\), respectively, the quantitative evaluation of the particle Stokes number also requires the liquid bridge height, which (for a fixed diameter of the supporting disks) depends on the aspect ratio AR. Assuming two representative values of this parameter (\(AR=0.25\) and \(AR=0.9\) in line with the results presented in Results) one gets, \(St^{AR=0.25}_{MX-3000N}\cong 3.55\times 10^{-6}\), \(St^{AR=0.25}_{ME-200N}\cong 1.58\times 10^{-4}\) and \(St^{AR=0.9}_{MX-3000N}\cong 2.74\times 10^{-7}\), \(St^{AR=0.9}_{ME-200N}\cong 1.2\times 10^{-5}\), respectively. Substitution of these data into the previously defined inequality immediately shows that it is largely satisfied in terms of particle Stokes number. Verification of the second inequality (in terms of \(St_{\omega }\)) requires the value of the frequency of oscillation related to the two considered representative aspect ratios. This can be gathered together with the corresponding critical value of the Marangoni number from Fig. 2. Given the specific definition of this non-dimensional frequency, i.e. \(f_c=\left( D/2\right) ^2/\left( \alpha \sqrt{Ma_c}\right) f_c^{dim}\) (Fig. 2b), where \(f_c^{dim}\) is the frequency in dimensional form, with the present reference quantities the non-dimensional \(\omega\) can be expressed as \(2\pi f_c \sqrt{Ma_c} AR^2\) \(\rightarrow\) \(St_{\omega }=St 2\pi f_c \sqrt{Ma_c} AR^2\) where \(Ma_c\) can be extracted from Fig. 2a. Accordingly, for \(AR=0.25\rightarrow f_c\cong 1\rightarrow \omega \cong 100\) and for \(AR=0.9\rightarrow f_c\cong 0.4\rightarrow \omega \cong 350\). Substitution of these data into the second inequality shows that the condition for which the Basset force can be neglected it is still largely met, thereby lending reliability to the present approach. Following Kuhlmann et al. (2014), we also neglect the Saffman lift force and the Faxén corrections (known to be very small for values of the Stokes number such as those considered in the present work). Moreover, evidence exists that numerical simulations carried out without these terms can adequately capture PAS for conditions corresponding to the experiments.

Accordingly, the particle motion is described by the following reduced form of the Maxey-Riley equation (refer, e.g., to the work of Babiano et al. (2000))

$$\begin{aligned} \rho _p\frac{d\mathbf{v} }{dt} = \rho \frac{D\mathbf{u} }{Dt} -\frac{9\mu }{d_p^2}\left( \mathbf{v} -\mathbf{u} \right) - \frac{\rho }{2} \left( \frac{d\mathbf{v} }{dt} - \frac{D\mathbf{u} }{Dt}\ \right) \end{aligned}$$
(6)

where \(\rho _p\), \(\mathbf{v}\) and \(d_p\) are, in the order, the density, velocity and diameter of the particles, while D/Dt is the material derivative operator.

Numerical Methodology

All the simulations reported in the present work have been carried out by means of the transient incompressible “buoyantBoussinesqPimpleFoam” solver available in the OpenFOAM package, while the thermocapillary stress condition was implemented anew and thoroughly tested by our group (see e.g., Capobianchi and Lappa (2020)). Moreover, the boundary/interface-particle interaction schemes available in the native OpenFOAM releases have been modified to take into account the actual dimension of the particles (see further below for a more detailed discussion of this model).

In particular, the Navier-Stokes and energy Equations (1)-(3) have been discretised in space on the basis of the Finite Volume Method (FVM). The PISO algorithm of Issa (1986) with a collocated grid arrangement of the variables has been used to integrate the momentum equations and enforce mass conservation, while the energy Equation (3) has been solved subsequently in a segregated manner. Integration in time of the governing equations (for both the flow and the particles) has been performed adopting the backward Euler scheme while convective terms have been discretised using the second-order accurate, linear central-differences scheme.

The PISO algorithm can be seen as a specific realization of a more general class of numerical techniques known as projection methods. These methods rely on a common set of principles and ensuing modus operandi. More specifically, they are articulated in three main steps. In the first step a provisional velocity field (\(\varvec{V}^*\)) is determined through solution of a simplified version of the momentum equation where the gradient of pressure is neglected. In a second step, this velocity is formally corrected reintroducing the gradient of pressure (the effective velocity is expressed as the difference between the provisional one and the gradient of pressure, i.e. \(\varvec{V}=\varvec{V}^*-\Delta t\varvec{\nabla }p\) where \(\Delta t\) is the time step used to integrate the simplified momentum equation). The corrected velocity is then forced into the continuity equation. In this way a Poisson equation \({\nabla }^2p=\varvec{\nabla }\cdot \varvec{V}^*/\Delta t\) is obtained that can be used to evaluate numerically the (otherwise unknown) pressure. Such a splitting relies mathematically on the Hodge decomposition theorem, that states that any vector field (\(\varvec{V}^*\)in this case) can be decomposed into a solenoidal part and the gradient of a scalar quantity. It also takes advantage of the theorem of inverse calculus, by which a vector field is fixed when its divergence and curl are assigned. The introduction of the Poisson equation reveals the intrinsic “elliptic” nature of the problem from a spatial point of view. With OpenFoam, this equation is solved iteratively in the framework of a Generalized Geometric-Algebraic Multi-Grid (GAMG) method with homogeneous Neumann boundary conditions.

In such a context, we also wish to highlight that Eq. (6) intrinsically relies on the idea that tracers can be considered geometric points in the sense that the particles centre of mass can occupy any position (the same position can also simultaneously be occupied by different particles) within the domain. In this sense, the model lacks the ability to account for the effective particle dimension, which however is known to influence the dynamics of PAS formation when they reach the boundaries of the domain (Hofmann and Kuhlmann 2011). In order to circumvent such limitation, a particle-wall interaction model has been implemented by modifying an existing OpenFOAM library on the basis of a paradigm already implemented in similar numerical experiments (see, e.g., (Hofmann and Kuhlmann 2011)). This is based on the idea that the particle velocity component normal to the wall must be annihilated once the particle periphery (i.e., its surface) reaches the boundaries of the domain. The interested reader is referred to Capobianchi and Lappa (2020) for more detailed information about the rationale underlying this way of thinking and its practical implementation in OpenFOAM.

Non uniform grids have been used for the simulations. Starting from a structured mesh with uniform spatial step along the axial coordinate (consisting of a number of points \(40\times 40\times 40\) in the radial, azimuthal and axial directions, respectively, hereafter simply referred to as \(M_1\) grid), the mesh spacing has been progressively reduced (by a factor 1.5 for each iteration of the grid assessment study) in order to obtain a local refinement in proximity to the heated rods sufficient to capture the thin thermal boundary layers expected to emerge there. The outcomes of this refinement study are summarized in Table 2 for two representative values of the liquid bridge aspect ratio, i.e. \(AR=0.9\) and \(AR=0.33\). As quantitatively substantiated by this table, the case \(AR=0.9\) has required three levels of local refinement in order to obtain a sufficiently small relative percentage error. For the case \(AR=0.33\), a good rate of convergence has been obtained after only two increments in the resolution.

Table 2 Mesh convergence study for \(AR=0.9\) and \(AR=0.33\). The number of grid points for each case has been kept constant while adjusting the spacing \(\Delta y_{min}\) in proximity to the two disks
Full size table

Results

Determination of the Critical Conditions

Here, we report a preliminary study dedicated to the determination of the critical conditions for the onset of instability for liquid bridges with different aspect ratios. The knowledge of the critical conditions before starting to deal with the investigation of PAS formation is extremely useful since, as indicated by available experimental studies, accumulation structures usually appear in ranges of Marangoni number which are about twice the critical one, \(Ma_c\) (Schwabe et al. 2007).

Having this in mind, the critical conditions (the above mentioned \(Ma_c\) and the associated critical frequency, \(f_c\), of the hydrothermal wave) have been computed through an extrapolation technique relying on the observation of the rate of growth of the maximum value of the azimuthal velocity, i.e., the velocity component associated with the the instability of Marangoni flow (which typically breaks the initial axisymmetry of this type of convection in liquid bridges). These results are summarised in Fig. 2, where they are also compared with the findings of previous experiments of (Nishino et al. 2015) and a linear stability analysis performed by the same authors assuming non-adiabatic interfaces. As evident in this figure, our results are in qualitative agreement with available data, although the critical Marangoni number is overestimated for all the aspect ratios considered. We ascribe such discrepancy to the fact that contrarily to the analysis of (Nishino et al. 2015), in the present model convective heat exchange at the interface is totally neglected. Previous works (see e.g., Xun et al. 2009; Melnikov and Shevtsova 2014; Melnikov et al. 2015; Yano et al. 2017) have shown, in fact, that moderate heat transfer can have a destabilizing effect. The corresponding critical dimensionless frequency,\(f_c=\left( D/2\right) ^2/\left( \alpha \sqrt{Ma_c}\right) f_c^*\), where \(f_c^*\) is the frequency in dimensional form, on the contrary, appears in better agreement with the available data, suggesting that this quantity is less sensitive to the thermal behavior of the interface.

Fig. 2
figure 2

Critical conditions as a function of the aspect ratio, AR, for \(Pr=68\) liquid bridges in microgravity conditions: critical Marangoni number, \(Ma_c\)(left), and dimensionless critical frequency (right). The present results (red triangles) are compared with previous investigations (linear stability analysis and experiments by Nishino et al. (2015))

Full size image

Particle Accumulation in Low-aspect-ratio Liquid Bridges

Having completed the determination of the critical conditions (as a pre-requisite for the exploration of the PAS formation phenomena in the supercritical region of the space of parameters), in the present section we turn to analysing particle aggregation phenomena.

In particular, first we concentrate on shallow liquid bridges. Then taller liquid columns are considered.

In this regard, we start from the observation that, in line with earlier studies on the subject (Lappa et al. 2001), we have found pulsating flow (i.e. standing waves) to be the dominant mode of convection for relatively small values of the aspect ratio.

Despite the stability of this specific waveform is still a matter of debate in the literature, the present simulations (though based on a completely different computational platform) confirm the earlier findings by Lappa et al. (2001), namely, that pulsating flows are the preferred mode of convection in shallow liquid bridges (these waveforms have also been found in experiments where they have been observed as long lasting modes of convection (Matsugase et al. 2015; Schwabe and Mizev 2011)).

As shown by the present computations, in such conditions particles self-organise in centrosymmetric patterns, which can ideally be subdivided in a number of circular sectors equal to \(n=2m\) (“cake-like accumulations”), where m is the azimuthal wavenumber. Similar phenomena were previously observed experimentally by (Schwabe and Mizev 2011).

The present results for \(AR=0.25\) (\(m=3\)) and \(Ma/Ma_c=2.25\) are summarised in Fig. 3 More precisely, while Fig. 3a illustrates the outcomes of the numerical simulations for the 30 \(\upmu\)m particles (MX-3000N), Fig. 3b shows the equivalent pattern obtained in the presence of larger particles (MR-200N). It can be seen that the two particle distributions display similar features. In particular, \(n=6\) circular sectors are evident along the circumferential extension of the liquid bridge. They repeat at regular azimuthal intervals of \(\pi /3\) (such subdivisions appear more evident for the case with larger particles). Moreover, a depleted region (clear fluid) is created in a well-defined area surrounding the axis of the liquid bridge. It has the shape of a hexagon (dashed lines have been added to highlight this geometrical structure) whose extension depends on the tracer type. In particular, the depletion area observed for the case with MX-3000N particles (\(\rho _p=690\, \text {kg} \text {m}^{-3}\)) appears more extended than the one formed for the case with the MR-200N particles (\(\rho _p=660\, \text {kg} \ m^{-3}\)). As the percentage difference in terms of particle density is almost negligible (\(\cong 4\%\), see Table 1), we argue that the cause for the different behavior must be ascribed to the particle size (\(\cong 600\%\) difference in terms of diameter). Another distinguishing mark relates to the variation of the radial extension of this internal clear (depleted) area with the axial coordinate (see Fig. 4). While for the smaller particles changes in size are almost negligible (Fig. 4a, b), for the larger particles the radial extension undergoes significant modifications along the axis of the liquid bridge (compare Fig. 4c and Fig. 4d).

Fig. 3
figure 3

Accumulation patterns for \(AR=0.25\) and \(Ma/Ma_c=2.25\) (as seen through one of the supporting disks): (a) MX-3000N particles, (b) MR-200N particles

Full size image
Fig. 4
figure 4

Particle locations for the case \(AR=0.25\) and \(Ma/Ma_c=2.25\) and contours of the temperature field in the liquid bridge midplane: (a, b) combined particle-temperature distribution as seen from above and from below (through the supporting disks), respectively for the MX-3000N particle type, (c, d) for the MR-200N particle type

Full size image

As a concluding remark for this aspect ratio, we wish to highlight that the abovementioned hexagonal area reflects the pattern of the temperature contour. This correspondence is not limited to purely topological (static) aspects (i.e. ”the six-node shape of the polygon line delimiting the internal depleted area”). It also concerns the spatio-temporal behavior of the overall velocity and temperature fields.

As the pulsating temperature field evolves in time, the sides of the hexagon (this polygon is made visible in the temperature field by cold fluid i.e. by the blue color) are periodically replaced by corners (where sides meet) and vice versa. The particle accumulation pattern keeps constantly reorganising accordingly (as we will illustrate in more detail for the case with aspect ratio \(AR=0.33\)).

The results for \(AR=0.33\) (mode \(m=2\)) and \(Ma/Ma_c\approx 1.2\) are shown in Fig. 5 and Fig. 6. The flow still behaves as a standing wave. In particular, the top figures provide the accumulation structure obtained considering the MX-3000N particles at three different times equally spaced within the oscillation period (taken approximately at \(\mathrm {T}/4\) apart, where \(\mathrm {T}\) is the period of the hydrothermal wave). Similarly, the bottom snapshots show the time evolution when the larger particles are considered (MR-200N particle type). As the reader will realize by inspecting these figures, for both conditions, the particles show the tendency to gather in four different sectors (note that in this case \(n=4\)) repeated at regular intervals of \(\pi /2\). Correspondingly, a rectangular depletion area is formed in the central region. Interestingly, the major side of this 4-side polygon does not keep a fixed orientation; rather it alternatively aligns with one or the other two fundamental directions which correspond to the sides of the rectangle (in other words, the quadrilateral stretches continuously but it does not rotate, in agreement with the main features of the standing wave).

Notably, unlike what has been observed for \(m=3\), the central depletion region exhibits four well defined branches that extend in the radial direction up to the periphery of the bridge (as a result, the four sectors where particles accumulate behave as well separated regions).

If considered with the previous results shown for \(AR=0.25\) and those obtained by Lappa (2014), these findings are extremely interesting as they lead to the conclusion that the boundaries separating the circular sectors can behave as loci of attraction or vice versa as loci of repulsion for particles. Moreover, this change is accompanied by an inversion in the relationship between the size of the internal depleted region and the diameter of particles (now larger particles correspond to the situation with a more extended internal area of clear fluid).

The identification of the precise cause-and-effect relationships underlying these modifications will be the main subject of investigation of a forthcoming dedicated study (it being beyond the scope of the present analysis). Here we limit ourselves to pointing out that for \(AR=0.33\) the accumulations structures have not been observed to vary appreciably over the axis of the liquid bridge.

Fig. 5
figure 5

Patterning behavior for \(AR=0.33\) and \(Ma/Ma_c\approx 1.2\): From (a) to (c) three different snapshots are shown (taken approximately at \(\mathrm {T}/4\) apart). The top panels (black tracers) illustrate the evolution in time for the case with the MX-3000N particles, the bottom panels (red tracers) refer to the situation with the MR-200N particles

Full size image
Fig. 6
figure 6

Combined view showing particle clustering effects and fluid velocity field at two times equally spaced in the oscillation period (\(AR=0.33\) and \(Ma=58800\))

Full size image

Particle Accumulation in a High-aspect-ratio Liquid Bridges

In this section we finally concentrate on the accumulation structure observed for \(AR=0.9\) (\(m=1\)) (obtained for \(Ma/Ma_c\approx 1.75\)). In this case hydrothermal disturbances have been found to behave as a traveling wave.

The related particle accumulation pattern for the case with the MR-200N particles type (the same structure emerged also when we seeded the 30 micron particles) is reported in Fig. 7 (where the ”typical” thread-like structure can be recognised, the red continuous line being representative of the KAM tori, which behaves as a “template” for the accumulation of particles in such cases (Romanò and Kuhlmann 2018)).

Despite other structures of such a kind have already been reported in the literature for azimuthal wavenumber \(m=1\) (e.g. Schwabe et al. (2007)), the PAS presented here shows some peculiar features that should be regarded as an interesting distinguishing factor with respect to other known numerical simulations or experiments. The two loops which contribute to the spiraling nature of the resulting structure are essentially oriented along planes perpendicular to the axis of the liquid bridge (unlike the classical PAS observed by Schwabe et al. (2007), which is constituted by a central region that coils up and then extends from the base of the liquid bridge toward the opposite side following a path that for an external observer would appear inclined with respect to the axis of the liquid bridge).

In this regard, the present phenomenon is much more similar to the \(m=1\) particle circuits observed in previous on-the-ground experiments by Sasaki et al. (2004), later reproduced numerically by Capobianchi and Lappa (2021). In line with the observations of Sasaki et al. (2004), this PAS is not centrosymmetric. It resembles a helicoid (the similarity with their experimental results becomes very evident by comparing the axial view of their PAS with the outcomes of the present numerical simulations shown in Fig.7e).

Fig. 7
figure 7

PAS for \(AR=0.9\), \(Ma=47500\) and MR-200N particles: (a)-(d) Lateral views obtained for different rigid rotations of the entire domain (by constant angular increments of \(\pi /2\)), (e) Top view corresponding to the lateral view shown in panel (a), and (f) “perspective” view

Full size image

As a heretofore unseen feature (see Fig. 7a, d), however, the present helicoid is confined in the upper part of the bridge (in correspondence of the heated rod); surprisingly, its axial extension is smaller than half of the liquid bridge height (i.e. the distance between the two supporting disks), which may indicate that this kind of structure tends to be axially compressed when the Prandtl number is increased (compare, e.g., with Capobianchi and Lappa (2021)).

Finally, we wish to highlight that numerical simulations about particle accumulation for \(AR=0.5\) and \(AR=2\) (for both cases we observed the \(m=1\) mode), which have been included in the analysis of the critical conditions, did not lead to the observation of any accumulation structures in the range of values of the Marangoni number explored in the present work (\(Ma/Ma_c\approx 2.1\) for both the \(AR=0.5\) and \(AR=2\) cases).

Conclusions

In the present work we have investigated numerically the particle structure formation process in a high-Pr liquid liquid bridge considering microgravity conditions. Various bridge aspect ratios have been examined and different behaviors in terms of emerging (supercritical) flow have been observed accordingly.

In particular, for relatively low aspect ratios (\(AR=0.25\) and \(AR=0.33\)) the dominant waveform is represented by a standing wave (stable over the entire observation time (more than 100 periods of the HTW). In these conditions, “cake-like” accumulation patterns are the typical outcome of the particle self-organization mechanisms at work. These patterns are characterised by recognizable regions repeated circumferentially at regular intervals, bounded internally by a central area of clear fluid.

For \(AR=0.25\), for which the supercritical flow is represented by a standing wave with \(m=3\), particles gather in six different, essentially seamless, circular sectors. The central area of depletion (which displays the shape of a hexagon), however exhibits a different radial extension depending on the size of the considered particles. The particle diameter has also a non-negligible impact on the axial morphology of the overall internal depletion volume. For relatively small particles, this volume is delimited laterally by a straight (parallel to the liquid bridge axis) surface, whereas this property is lost when larger particles are considered.

For \(AR=0.33\), the accumulation pattern displays four sectors repeated at regular intervals (as expected for a standing wave with \(m=2\)). The peculiar signature of this case, however, is the emergence of four well-defined depletion areas of small (but finite) azimuthal thickness protruding from the central depletion area towards the surface of the liquid bridge. If considered with earlier numerical results on the subject where sharp radial particle-dense strips were found in place of depletion elongated areas, this finding is particularly interesting as it indicates the possible existence of a “dichotomy” in the mechanism leading to the emergence of the cake-like particle patterns, which would deserve additional investigation using stability techniques or other relevant numerical strategies.

Finally, we were able to observe particle accumulation also for the liquid bridge having aspect ratio \(AR=0.9\), for which we came across a rotating wave having mode \(m=1\). In this case, we have observed a kind of PAS, which to the best of our knowledge, was never encountered before. This particular type of PAS displays features somehow similar to those seen in earlier on-the-ground experiments (Sasaki et al. 2004). More precisely, the two types of structures, the present one and the one found by Sasaki et al. (2004), display notable elements of similarity only if observed from an axial perspective (i.e. through one of the supporting disks, by which the observer would get the illusion of a double-loop circuit). Putting aside this analogy, the peculiarity of the present PAS resides in its extremely limited axial extension (it is entirely confined in the half upper part of the bridge, whereas the structure of Sasaki et al. (2004) develops over the entire extension of the bridge).

The present set of results confirm that particle self-organization driven by laminar (time-periodic) flows may reflect different underlying mechanisms. As a concluding remark, we wish to recall that these time consuming simulations (taking up to one week for some specific cases) have been conducted with the specific intent to provide a basis to support the preparation and interpretation of the forthcoming JEREMI experiments. In this regard, the following clear indications for the JEREMI Science Team can be inferred from the present study: 1) Shallow liquid bridges made of a high-Pr fluid tend to favor standing waves; 2) accordingly, towards the end to clarify the mechanisms underlying one-dimensional PAS formation, special attention should be paid to liquid bridges with relatively high aspect ratios (for which traveling disturbances are the preferred supercritical mode of convection); 3) regardless of whether standing or travelling waves are produced, it is recommended to visualize the particle behavior through one of the supporting disks (while for standing waves this is the only way to recognize the peculiar arrangement of particles in circular sectors, for traveling waves, this may lead to the successful identification of heretofore unseen structures like that shown in Fig. 7, characterized by well-distinguishable loops in the axial view).