This paper investigates the applicability of smoothed particle hydrodynamics in simulations of superfluid helium-4. We devise a new approach based on Hamiltonian mechanics suitable for simulating thermally driven and weakly compressible flows with free surfaces. The method is then tested in three cases, including a simulation of the fountain effect. We obtain remarkable agreement with referential and theoretical results. The simulations provide new physical insight, such as the pressure and temperature fields in a vessel experiencing the fountain effect.

Quantum fluids are governed by laws of quantum mechanics rather than classical ones and represent a fairly recent addition to the family of various types of fluids with relevance both in fundamental science and in technological applications. Prime examples can be found in superfluid phases of helium isotopes or in Bose–Einstein condensates of rarefied gases. These fluids possess remarkable properties, such as the possibility of frictionless flow, the existence of numerous sound modes, the formation of superfluid films, quantized circulation around discrete vortex lines, or the ability to transport heat convectively with zero local mass flow. In addition to the fields of low-temperature physics and fluid dynamics, this behavior makes quantum fluids appealing, via certain analogies,62 for researchers across multiple disciplines, including astrophysicists13,44 or cosmologists,22,65 for which they represent an accessible model system.

In this work, we deal specifically with the isotope helium-4. The liquid phase, historically called He I, exists below 4.2 K at standard pressure. When cooled further by pumping on its vapors, helium-4 undergoes a second-order phase transition at 2.17 K and enters the superfluid phase, called He II. Most physical properties of helium are well-known and tabulated.6 

He II exhibits all of the exotic phenomena mentioned above; for instance, it flows easily through narrow capillaries without friction and displays two-fluid behavior (in the sense of Landau's model, see further below). The latter leads to the mechano-caloric effect, or the famous fountain effect,1 where a heater causes the motion of helium in the form of a fountain-like jet.2 Additionally, various forms of turbulent motion may exist in He II, including a state of two-fluid “double turbulence;” for a phenomenological treatment of quantum turbulence, see Ref. 47.

Historically, experimental investigations of dynamics in He II may be divided into studies of thermally generated flows, such as the fountain effect, counterflow in a channel driven by a heater placed at its closed end,52 or mechanically driven flows, mostly by submerged oscillating structures.14,20,43 A superfluid wind tunnel experiment has been constructed,41 and a large-scale Kármán flow experiment is now in operation.40 Methods for flow visualization using tracer particles25,26 or He excimers42 have been developed recently and applied to both thermally and mechanically driven flows.

The motion of superfluid helium-4 can be simulated by solving the Gross–Pitaevskii equation,46 the Hall–Vinen–Bekarevich–Khalatnikov (HVBK) models,5,15 the one-component models,35 or by the vortex-filament method.10,19 However, many experiments are notoriously difficult to simulate numerically using mesh-based techniques such as finite element or finite volume methods, as they contain either a free surface (fountain effect) or a moving boundary (oscillating bodies). Therefore, it is advantageous to use a method based on Smoothed Particle Hydrodynamics (SPH),31,33,58 which is mesh-free and thus more suitable for these problems.

This paper presents a new numerical mesh-free method for simulations of superfluid helium-4. An SPH method for helium II has already been developed by Tsuzuki,54–56 who investigated vortices emerging in a rotating cylinder. Our approach is different because we discretize helium II using only one type of fluid particle with two densities as opposed to the work of Tsuzuki, where normal and super particles are distinguished. Our method comes at the cost of including additional convective terms, which consider that entropy, for instance, is advected by normal velocity and not by the collective or coflow velocity. However, it liberates us from assuming that superfluid behaves like a mixture and facilitates the implementation of heating bodies (a crucial feature in thermally driven flows, like the fountain effect). Additionally, we use a slightly different form of equations, a Hamiltonian variant of the two-fluid model that goes beyond the Landau–Tisza approach. We explain this in detail in  Appendix A. The Hamiltonian characteristic of the equations is advantageous also in the SPH numerical scheme.23 

The first macroscopic models of superfluid helium-4 were proposed by Tisza et al.16,27,49,51 Their final model21 consists of four evolution equations, namely, the continuity equation, balance of momentum, evolution of superfluid velocity, and entropy balance. The Landau–Tisza model can also be derived from the Gross–Pitaevskii equation.53 However, the model has several limitations. First, it does not allow for non-zero superfluid vorticity (quantum vortices).

Second, the Landau–Tisza model is formulated in terms of five quantities (superfluid density ρs, normal density ρn, superfluid velocity v s, normal velocity v n, and entropy density s), despite having only four evolution equations. This inconsistency is overcome by setting a dependence of the ratio ρ n / ρ on temperature ( ρ = ρ s + ρ n being the total mass density). Although this setting closes the evolution equations, it goes against the nature of superfluid helium-4, which is not a mixture of two fluids, but rather a single fluid with two motions, as expressed by Landau:27,29 “It must be particularly stressed that we have here no real division of the particles of the liquid into ‘superfluid’ and ‘normal’ ones….” The dependence ρ n / ρ ( T ) actually defines how free energy depends on temperature, which has to be taken into account because otherwise compatibility with Hamiltonian mechanics would be violated.50 The HVBK model resolves the problem of four equations for five quantities by requiring both v n and v s to be divergence-free (incompressible).4,38 Another solution was proposed by Zilsel,64 who introduced a continuity equation for both ρs and ρn, despite that there are no real two components in the superfluid helium (as noted by Landau). Yet another solution is the so-called one-fluid model,35 where entropy exhibits an extra motion instead of density.

Here, we follow a route that builds upon Hamiltonian structures of the evolution equations for superfluid helium-4. The Hamiltonian structure of Landau–Tisza equations was derived from quantum commutators between the mass density and phase of the wave function.7,27 Moreover, Volovik and Dotsenko extended the Hamiltonian structure to a model related to HVBK dynamics by also considering Hamiltonian motion of quantum vortices.59–61 Alternatively, the Hamiltonian form of equations can be derived from the one-fluid model.50 This is the structure we use in the current manuscript. Another Hamiltonian formulation of HVBK dynamics was obtained by Holm and Kuperschmidt,11,17,18 which contains additional state variables (vector and scalar potentials of the superfluid velocity).

This manuscript is organized as follows. Section II contains both the continuous model of superfluid helium-4 and the discrete SPH counterpart. Section III describes numerical simulation of the fountain effect and its comparison with experimental data.

This section contains a two-fluid model of superfluid helium-4 that extends the Landau–Tisza model and its discretization within SPH.

The reversible part of our model consists of the following four evolution equations:
D ρ D t = ρ · v , D s D t = 1 ρ · ( ρ s χ s v n s ) , D v D t = 1 ρ · ( ρ χ n χ s v n s v n s + p I ) , D v s D t = χ n v n T v n s p ρ + s T ,
(1)
where “unknowns” ρ , s , v , and v s are the density, specific entropy, coflow velocity, and superflow velocity, respectively. We denote the temperature as T and pressure as p. Dimensionless variables χ n , χ s ( 0 , 1 ) are mass fractions of normal and super components and satisfy χ n + χ s = 1. The coflow velocity v and the counterflow velocity v n s satisfy the following relations:
v = χ n v n + χ s v s ,
(2)
v n s = v n v s .
(3)
We consider T , χ n , χ s, and p to be smooth functions of the unknowns. Equations (1) contain convective derivatives with respect to the overall coflow velocity v
D φ D t : = φ t + v · φ = 1 ρ ( ( ρ φ ) t + · ( ρ φ v ) ) , φ .
(4)
This is the difference with respect to the standard Landau–Tisza model,49,53 where the superfluid velocity is convected only by itself ( v s) and not by the whole coflow velocity. Equations (1) should, however, be superior to the Landau–Tisza model as they contain terms that at least partly take into account effects of quantum vorticity, see  Appendix A.
Equations (1) represent only the reversible part of superfluid dynamics. It can be beneficial to add some irreversibility to the system for increased realism, enhanced numerical stability, and noise suppression. We will consider two types of dissipative processes: viscosity of normal flow and parabolic thermal conduction. In future, the latter could be replaced with a complete model of quantum turbulence, which includes mutual friction. We include related entropy production terms to conserve the total energy. Thus, we write
D ρ D t = ρ · v , D s D t = 1 ρ · ( ρ s χ s v n s ) + β ρ Δ T + ζ ρ T , D v D t = 1 ρ · ( ρ χ n χ s v n s v n s + p I ) + 2 μ ρ · D n , D v s D t = χ n v n T v n s p ρ + s T ,
(5)
where μ > 0 is the dynamic viscosity, and β > 0 is a diffusion parameter,
D n = 1 2 ( v n + v n T )
(6)
is a normal velocity deformation tensor, and
ζ = β | T | 2 + 2 μ | D n | 2
(7)
is the dissipative power.

Although SPH can compute reversible flows, the stabilization of entropy helps to prevent the self-emergence of disorder. Numerical noise is a recognized problem in explicit SPH, which makes velocity converge to Boltzmann distribution and eventually leads to the unphysical gas-like behavior of simulated particles.23 We stabilize by adding a “small Laplacian” to the right-hand side of the entropic balance. Particle shifting32 would be an alternative remedy.

Closing the system (5) requires the knowledge of functions p, T, and χn. To the best of our knowledge, there is no theory that would express the energy of helium-4 for wide range of s , ρ , v n s 2 and be in agreement with experiments. However, we point out that due to extremely good heat conducting properties of helium-4,30 large temperature gradients are unlikely to exist in practice. For instance, the superfluid fountain, which is the main interest of this paper, can be powered by temperature difference between a heater and a reservoir as small as 10 3 K.2 Thus, it should be sufficient to use a linearized model, which is valid in a vicinity of certain referential temperature T0. For χn, χs, we use
χ n = χ n 0 + χ ( s s 0 ) , χ s = χ s 0 χ ( s s 0 ) ,
(8)
where χ , χ n 0 , χ s 0 , and s 0 are constant values at T0.
Yet, how does the energy depend on the state variables s, ρ, v n s, and v? By means of Galilean invariance, it can be derived that differential of the specific total energy e must be
d e = v n · d v χ s v n s · d v s + p ρ 2 d ρ + T d s
(9)
and that energy has to be in the form of
e = 1 2 v 2 + χ s 2 χ n ( v v s ) 2 + e 0 ( ρ , s ) ;
(10)
see  Appendix A for details.
In the absence of any flow ( v = v n s = v s = 0), the last undetermined term in the formula for energy (10) can be found by comparison with experimental data. Its differential reads
d e 0 = T d s + p ρ 2 d ρ .
(11)
Due to the low compressibility and thermal expansion coefficient,57 we can estimate the pressure as an affine function of density
p = u 1 2 ( ρ ρ 0 ) ,
(12)
where ρ 0 and u 1 are a referential density and the first speed of sound.
Similarly, in the absence of any flow, we can use
T = ( 1 + s s 0 C ) T 0 ( for v s = v n = 0 ) ,
(13)
where C is a referential heat capacity. Below T λ, heat capacity is connected to the second speed of sound u2 by27,
u 2 2 = χ s 0 T 0 s 0 2 χ n 0 C .
(14)
Finally, we integrate (12) and (13) and substitute them in Eq. (10) to obtain the following approximate total-energy formula:
e = 1 2 v 2 + 1 2 χ s χ n ( v v s ) 2 + u 1 2 ( ln ( ρ ρ 0 ) + ρ 0 ρ 1 ) + T 0 ( s s 0 ) + T 0 2 C ( s s 0 ) 2 .
(15)
Table I contains values of C , χ n 0 , χ s 0 , s 0 , χ , ρ 0 , and u 1 for a given T0 obtained from tabulated experimental data. For a comprehensive collection of measured values, we refer to Donnelly and Barenghi.6 
TABLE I.

Parameters used in simulations. We list only values relevant to the simulation output.

Cavity Second sound Fountain effect
T0    1.9 K  1.65 K 
ρ0  145.5 kg / m 3  145.2 kg / m 3 
s0    725.5 J / ( kg K )  335.0 J / ( kg K ) 
χ n 0  0.4195  0.1934 
χ   5.697 × 10 4 kg K / J  5.851 × 10 4 kg K / J 
u1  20  40 m / s  40 m / s 
u2    18.83 m / s  20.37 m / s 
C    3902 J / ( kg K )  1861 J / ( kg K ) 
μ  1 Re  10 4 kg / ( m s ) 
β  d ρ 0 s 0 2 200 u 2 
δr  1 N  L N  d 40 
h  3 δ r  3 δ r  2.8 δ r 
δt  h 10 u 1  h M u 1  h 5 u 1 
t end  40  2 L u 1  0.3 s 
H      2 × 10 3 m 
Cavity Second sound Fountain effect
T0    1.9 K  1.65 K 
ρ0  145.5 kg / m 3  145.2 kg / m 3 
s0    725.5 J / ( kg K )  335.0 J / ( kg K ) 
χ n 0  0.4195  0.1934 
χ   5.697 × 10 4 kg K / J  5.851 × 10 4 kg K / J 
u1  20  40 m / s  40 m / s 
u2    18.83 m / s  20.37 m / s 
C    3902 J / ( kg K )  1861 J / ( kg K ) 
μ  1 Re  10 4 kg / ( m s ) 
β  d ρ 0 s 0 2 200 u 2 
δr  1 N  L N  d 40 
h  3 δ r  3 δ r  2.8 δ r 
δt  h 10 u 1  h M u 1  h 5 u 1 
t end  40  2 L u 1  0.3 s 
H      2 × 10 3 m 
In SPH, derivatives are approximated by transferring them onto a convenient smoothing kernel function w = w ( r ), which satisfies
d w ( r ) d r = 1 ,
(16)
where d is the geometric dimension of a simulation, usually 1 or 2 (for translationally symmetrical problems) or 3. In this paper, we will use the Wendland's kernel function:
w ( r ) = { q d h d ( 1 r 2 h ) 4 ( 1 + 2 r h ) , r 2 h , 0 , r 2 h ,
(17)
where
{ q 2 = 7 4 π , q 3 = 21 16 π .
(18)
Here, constant h is called the smoothing length. Wendland's kernel is often advantageous to the Gaussian kernel because it has compact support and is defined via a piece-wise polynomial functions, which are easy to evaluate.
In our numerical method, superfluid is represented by N particles with positions r a for a = 1 , 2 , , N. Each particle has a certain velocity v a, superflow velocity v s , a density ρa, entropy sa, and temperature Ta and moves as a material point. At initial time, they can be arranged in a grid, such that particle a is centered in a cell of mass ma (often, ma can be chosen the same for every a). Replacing continuous derivatives using SPH techniques yields a system of ordinary differential equations, which are given as follows:
ρ ̇ a = b m b w a b r a b r a b · v a b , s ̇ a = b m b w a b r a b ( 1 ρ a 2 j a · r a b + 1 ρ b 2 j b · r a b 2 β T a b ρ a ρ b ) + ζ a T a ρ a , v ̇ a = b m b w a b r a b ( 1 ρ a 2 Π a + 1 ρ b 2 Π b ) r a b + b m b w a b r a b 2 ( d + 2 ) μ ρ a ρ b v n , a b · r a b r a b 2 + η 2 r a b , v ̇ s , a = b m b w a b r a b χ n , a ρ a ( v n , a b · v n s , a ) r a b b m b w a b r a b ( p a ρ a 2 + p b ρ b 2 + s a ρ a T a b ) r a b , r ̇ a = v a ,
(19)
where for every particle indices a, b
r a b = r a r b , r a b = | r a b | , w a b = w ( r a b ) , w a b = d w d r ( r a b ) , v a b = v a v b , v n , a b = v n , a v n , b , T a b = T a T b ,
(20)
the entropy flux and the momentum flux are
j a = ρ a s a χ s , a v n s , a , Π a = p a I + ρ a χ n , a χ s , a v n s , a v n s , a ,
(21)
the normal and counterflow velocity can be expressed via
v n , a = v a χ s , a v s , a χ n , a , v n s , a = v a v s , a χ n , a ,
(22)
and the dissipative power on particle a is
ζ a = β b m b ρ b w a b r a b T a b 2 2 μ ( d + 2 ) b m b ρ b ( v n , a b · r a b ) 2 r a b 2 + η 2 w a b r a b .
(23)
The entropy production term in this form is also present in the SDPD (Smoothed Dissipative Particle Dynamics) method.8 A numerical parameter η = h 10 prevents division by zero when two particles overlap. For detailed explanation of how Eqs. (19) are obtained, we refer to  Appendix B, where we also prove that they conserve energy and momentum and satisfy the entropic inequality. An important aspect of our method is that the do-nothing condition represents a free adiabatic surface. This means that, conveniently, free surfaces do not require any implementation and particles on this boundary do not need to be identified. Helium vapors are not modeled in this approach. Note that by this, we neglect any aerodynamic effects induced by vapors on liquid helium.
Instead of updating density iteratively, we prefer to use an equivalent closed formula
ρ a = b m b w a b + C a ,
(24)
where Ca is an integration constant specified by the initial condition. This avoids accumulation of time discretization errors.

Simple explicit integrators, like leap-frog, are commonly used in SPH codes. In this paper, we employ the following scheme, which is very similar to leapfrog, except we update entropy and density using mid-time positions:

  1. (initial step only) find rate of v and v s

  2. update v and v s by δ t 2 step

  3. update x by δ t 2 step

  4. recompute the cell list and find ρ

  5. find the rate of s

  6. update s by δ t 2 step

  7. update x by δ t 2 step

  8. recompute the cell list and find ρ

  9. find the rate of v and v s

  10. update v and v s by δ t 2 step

In each step, rates are evaluated using the most recently computed values of s , ρ , v , v s , and x.

Before we jump to examples involving actual superfluidity, we assess the validity of our model for classical fluid (Navier–Stokes equations), which can be understood as a special case of system (1) for χ s = 0. Clearly, our discrete equations (19) must work for classical fluids; otherwise, there would be no hope to use them on more complex problems.

Lid-driven cavity is usually formulated as a dimensionless problem: Incompressible, viscous flow is confined in a box ( 0 , 1 ) × ( 0 , 1 ). No slip boundary is prescribed at the left, bottom, and right walls. At the top, velocity is v = ( 1 , 0 ) T, density is ρ 0 = 1, and viscosity μ = 1 Re. Reynolds number Re has various values.

Walls are implemented using a layer of dummy particles, which are treated normally with the exception that their velocities are constantly zero. The lid is also implemented this way, but with velocity v a = ( 1 , 0 ) T appearing in the evaluation of viscous forces. The list of simulation parameters is displayed in Table I. We measure transverse velocities along x and y centerlines and compare them to a reference solution by Ghia et al.12 The results are shown in Figs. 1–3 and convergence curve in Fig. 4. Error grows approximately linearly with δ r = 1 N.

FIG. 1.

SPH result compared to the referential solution for Re = 100 and N = 336. Left: y-component of velocity for y = 0.5. Right: x-component of velocity for x = 0.5.

FIG. 1.

SPH result compared to the referential solution for Re = 100 and N = 336. Left: y-component of velocity for y = 0.5. Right: x-component of velocity for x = 0.5.

Close modal
FIG. 2.

SPH result compared to the referential solution for Re = 400 and N = 336. Left: y-component of velocity for y = 0.5. Right: x-component of velocity for x = 0.5.

FIG. 2.

SPH result compared to the referential solution for Re = 400 and N = 336. Left: y-component of velocity for y = 0.5. Right: x-component of velocity for x = 0.5.

Close modal
FIG. 3.

The streamline plot for Re = 400. Left: Our SPH result. Image created using line integral convolution in Paraview.48 Right: Reference solution.9 

FIG. 3.

The streamline plot for Re = 400. Left: Our SPH result. Image created using line integral convolution in Paraview.48 Right: Reference solution.9 

Close modal
FIG. 4.

Convergence curves for Re = 100 (left) and Re = 400 (right). Error is measured as a maximum distance from reference data. Regression slopes were 0.998 and 1.189, respectively, suggesting approximately linear convergence.

FIG. 4.

Convergence curves for Re = 100 (left) and Re = 400 (right). Error is measured as a maximum distance from reference data. Regression slopes were 0.998 and 1.189, respectively, suggesting approximately linear convergence.

Close modal
Contrary to the previous case, we will now investigate a situation where coflow velocity v is almost zero. Let us consider a standing wave of second sound, written in terms of entropy as
s ( x , t ) = s 0 + A s 0 sin ( π x L ) sin ( π y L ) cos ( 2 π u 2 t L )
(25)
for x Ω = ( 0 , L ) × ( 0 , L ) and L = 1 cm. This is an approximate solution of continuous two-fluid equation (1) linearized around v = v s = 0 , ρ = ρ 0, and s = s0. As a test of consistency, we would like to verify that our numerical model approaches (25) for fine resolution and small data ( A 0). We impose initial condition
ρ ( x , 0 ) = ρ 0 , s ( x , 0 ) = s 0 + A s 0 sin ( π x L ) sin ( π y L ) , x Ω , v ( x , 0 ) = v s ( x , 0 ) = 0
(26)
and free adiabatic boundaries. Particles are initially arranged on a Vogel spiral with spatial step δ r = L N. Time step is δ t = L M u 1, and we end the simulation at t = 2 L u 2. Let us define two types of errors. First, a normalized l l 2 error
ϵ 1 = max n a | s a ( t n ) s exact ( x a , t n ) A s 0 | 2 .
(27)
Second, a dimensionless energy error
ϵ 2 = max n | a e a ( t n ) a e ( 0 ) 1 | .
(28)
Figure 5 shows the results for various values of A, N, and M. We observe that energy error is very small and we can also achieve satisfactory ϵ1 error below 1% (see Fig. 6). It is likely that the convergence rate could be improved with renormalized operators.37 
FIG. 5.

Left: ϵ1 error for M = 50 and different values of A and N. Right: ϵ2 error for N = 200 and different values of M and A.

FIG. 5.

Left: ϵ1 error for M = 50 and different values of A and N. Right: ϵ2 error for N = 200 and different values of M and A.

Close modal
FIG. 6.

Evolution of entropy at x = y = L 4 in a simulation compared to the formula (25). Result for A = 0.01, M = 50, and N = 300.

FIG. 6.

Evolution of entropy at x = y = L 4 in a simulation compared to the formula (25). Result for A = 0.01, M = 50, and N = 300.

Close modal

Unlike the two previous examples, in superfluid fountain effect, coflow and counterflow are equally significant. Let us begin by describing our idealized setting. A cell (with shape indicated in Fig. 7) is submerged in liquid helium-4 at temperature T 0 = 1.65 K. In its center, there is a point heat source of constant power W ̇. The cell is insulated by adiabatic walls and a thin superleak at the bottom.

FIG. 7.

Scheme of the simulated apparatus. All lengths are expressed in millimeters. Shape is composed of line segments and circular sections with indicated radii.

FIG. 7.

Scheme of the simulated apparatus. All lengths are expressed in millimeters. Shape is composed of line segments and circular sections with indicated radii.

Close modal
In (19), we add gravitational acceleration g, superleak friction force f sup , a, and dissipative power from heater ζ heater , a,
s ̇ a = + 1 T a ρ a ζ heat , a , v ̇ a = + g + f sup , a , v ̇ s , a = + g ,
(29)
where
ζ heat , a = W ̇ ( w H * δ { x = x heat , y = y heat } ) ( r a ) , f sup , a = χ n , a v n , a τ ( w H * δ { y = y sup , | x | r sup } ) ( r a ) .
(30)
Function wH is the kernel (17) with smoothing length H and * denotes the convolution operator. Mollification by wH is necessary because it is impossible to work with Dirac distributions in SPH. To ensure that the superleak is practically impenetrable for normal component, relaxation time τ must be as small as 1 μ s. Therefore, it is better to split superleak friction from other forces and compute it implicitly.

Adiabatic walls of the cell are implemented using static dummy particles, which take part only in density and pressure computation. At the walls of the cryostat, let us prescribe no-slip for v, free condition for v s, and Dirichlet condition S = S0 for entropy (which is almost the same as T = T0). Consequently, the wall acts as a cooler, preventing the temperature from growing indefinitely. This condition is implemented using dummy particles, which we force to have constant entropy and zero coflow velocity, but we allow them to conduct heat.

In the initial state, v = v s = 0 and T = T0 everywhere. A fountain is quickly generated. The speed of the jet is measured by counting particles that escape from the cell in a given time window. We then compare these numbers to a theoretical formula from the book by Landau and Lifshitz,29 which we adapted to a two-dimensional flow
v jet , ideal = W ̇ T s ρ d ,
(31)
where T , s , and ρ are values of temperature, entropy density, and mass density inside the cell (which are nearly constant). We should emphasize that (31) is merely an approximation as it does not take into account heat loss through capillary and dissipative effects. It is not very clear whether the corrected value of v jet should be above or below the theoretical estimate. Assuming stationary state, the balance of entropy yields
0 = d d t cell s d m = W ̇ heater T + 1 T ρ cell ζ d m s ρ d v n , jet ,
(32)
where the dissipative power ζ is positive. This would, paradoxically, predict that v jet , ideal is a lower estimate when dissipation is considered. However, a counterflow can occur through capillary, which means that the normal speed v n , jet can be different from v jet. Additionally, heat can partially escape through walls and superleak in practice. Measurements by Amigó et al.2 report values both slightly lower and higher than (31).
The Reynolds number can be given as follows:
Re = ρ 0 d v jet , ideal μ = W ̇ T s μ .
(33)
With μ = 1.296 × 10 6 Pa s, this yields Re in the range of hundreds of thousands for reasonably high values of power. Therefore, classical turbulence appears in the simulation, which is difficult to resolve without a dedicated turbulent model. We sidestep this problem by choosing an artificial μ, keeping the Reynolds number around 1500. A parameter analysis reveals that the relative jet speed does not vary significantly with Re for values in range 100–1500, which are accessible to our model—see Fig. 11.

The simulation never reaches a stationary state, but the jet speed is relatively stable after a brief transitional period. Figure 8 shows that the computed jet speed is very close to the value (31) but slightly lower on average. Figure 9 shows the shape of the fountain, velocity, and pressure. The simulation remains stable even during the violent impact of jetted helium. Figure 10 shows the temperature in the first millisecond of simulated time. After this short period, second sound waves fill the cell, and temperature gradients vanish, except for a steep jump at the superleak. The subsequent evolution is shown in Figs. 11 and 12. When the fountain connects with colder helium in the cryostat, a temporary drop in temperature occurs. Nonetheless, the effect of the second sound traveling upstream and entering the cell through the capillary was not observed. Table I lists all simulation parameters for reproduction purposes.

FIG. 8.

Evolution of the jet speed in comparison to theoretical prediction (31). The left image is for W ̇ = 60 W / m and the right for W ̇ = 80 W / m. The discrepancy for t ( 0 , 0.1 s ) is transitional effect caused by the evolution from the initial state with constant temperature and zero velocity field. Calculating the average for t > 0.2, we find that the SPH prediction is smaller than (31) by 1.95% (left image) and 2.38% (right image).

FIG. 8.

Evolution of the jet speed in comparison to theoretical prediction (31). The left image is for W ̇ = 60 W / m and the right for W ̇ = 80 W / m. The discrepancy for t ( 0 , 0.1 s ) is transitional effect caused by the evolution from the initial state with constant temperature and zero velocity field. Calculating the average for t > 0.2, we find that the SPH prediction is smaller than (31) by 1.95% (left image) and 2.38% (right image).

Close modal
FIG. 9.

Coflow velocity plot (up) and pressure plot (down) at t = 0.3 s for different values of heating power. From left to right: W ̇ = 60 , 80 , and 100 W / m. Each simulation consists of approximately 200 000 particles. The pressure is very noisy, especially in the capillary and the plunge pool.

FIG. 9.

Coflow velocity plot (up) and pressure plot (down) at t = 0.3 s for different values of heating power. From left to right: W ̇ = 60 , 80 , and 100 W / m. Each simulation consists of approximately 200 000 particles. The pressure is very noisy, especially in the capillary and the plunge pool.

Close modal
FIG. 10.

Plot of temperature inside the cell at t = 0.25, 0.5, 0.75, and 1 ms (from left to right).

FIG. 10.

Plot of temperature inside the cell at t = 0.25, 0.5, 0.75, and 1 ms (from left to right).

Close modal
FIG. 11.

Left: Evolution of temperature in time for W ̇ = 80 W / m. The dashed orange line indicates when the fountain connects with the cryostat. Right: Relative jet speed v ̂ = v jet / v jet , ideal (averaged over time) for various values of Re.

FIG. 11.

Left: Evolution of temperature in time for W ̇ = 80 W / m. The dashed orange line indicates when the fountain connects with the cryostat. Right: Relative jet speed v ̂ = v jet / v jet , ideal (averaged over time) for various values of Re.

Close modal
FIG. 12.

Plot of temperature distribution and fountain shape for W ̇ = 80 W / m at equally spaced time frames from 37.5 to 300 ms.

FIG. 12.

Plot of temperature distribution and fountain shape for W ̇ = 80 W / m at equally spaced time frames from 37.5 to 300 ms.

Close modal

We described a new energy-conserving SPH-based numerical method for superfluid helium. We verified its validity in three essential test cases, demonstrating that the technique can accurately capture coflow and thermally driven counterflow. Our fountain effect simulation verifies the applicability of the formula for fountain jet speed by Landau and Lishitz.

The numerical scheme discretizes a Hamiltonian one-fluid formulation generalizing the Landau–Tisza two-fluid model. We also provide a closed form of the energy functional in terms of density, coflow velocity, entropy density, and superfluid velocity, which has been missing in the literature.

In the future, we would like to enrich our discrete system by implementing the dissipative effects of quantum vorticity. Particle simulation of Rollin films is another promising research direction, as well as simulation of helium droplets,24 or simulation of inertial particles in superfluid helium.45 Moreover, the mesh-free nature of SPH is advantageous in multiphase flow simulations, including phase change.34,63 Consequently, we expect that this method could be extended to simulate solidification and evaporation in superfluid helium.

O.K. and D.S. were supported by Project No. START/SCI/053 of Charles University Research program. M.P. was supported by Czech Science Foundation, Project No. 23‐05736S.

The authors have no conflicts to disclose.

Ondřej Kincl: Investigation (equal); Methodology (equal); Software (equal); Visualization (lead); Writing – original draft (lead). David Schmoranzer: Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). Michal Pavelka: Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal).

The data that support the findings of this study are available within the article.

This appendix contains a derivation of model (1) from a Hamiltonian formulation of superfluid-helium dynamics.50,59 The field of superfluid density interacts with itself partly due to the convective term and partly due to the presence of quantum vortices. When the reversible (non-dissipative) part of interactions mediated by quantum vortices is taken into account, we obtain a system of four Hamiltonian equations for four state variables (mass density ρ, momentum density m, volumetric entropy density s ¯, and superfluid velocity v s) as follows:

(A1)
ρ t = k ( ρ E m k + E v s k ) ,
(A1a)
m i t = j ( m i E m j ) j ( v s i E v s j ) ρ i E ρ m j i E m j s ¯ i E s ¯ v s k i E v s k + i ( E v s k v s k ) ,
(A1b)
s ¯ t = k ( s ¯ E m k ) ,
(A1c)
v s k t = k E ρ k ( v s j E m j ) + ( k v s j j v s k ) ( E m j + 1 ρ E v s j ) ,
(A1d)
where the total energy E with subscripts stands for functional derivatives of the energy with respect to the fields in the subscript. In particular, Es = T is the temperature, E m = v n is the normal velocity, and E ρ is the generalized chemical potential μ. The derivative E v s has no particular name, but it is investigated below. Equations (A1) represent a Hamiltonian system generated by a Poisson bracket that is derived from quantum commutators and dynamics of vortices7,60 or by the requirement of unconditional validity of Jacobi identity.50 Note that these equations differ from those of Landau and Tisza in the terms coupling the evolution of v s with the derivative of energy with respect to v s, which are a consequence of the vortex dynamics or of the Jacobi identity.

In order to close evolution equations (A1), we have to provide an energy functional E = e ¯ d r. We assume that the volumetric energy density e ¯ ( ρ , m , s , v s ) = ρ e ( ρ , v , s , v s ) is a smooth function of the state variables (not a function of gradients of the state variables). This assumption excludes, for instance, any explicit dependence of energy on superfluid vorticity ω = × v s.7,17 Moreover, from the requirement of Galilean invariance,28,39 it follows that the energy must be in the form
e ¯ = ϵ ¯ ( ρ , s , m ρ v s ) + ( m ρ v s ) · v s + 1 2 ρ v s 2 ,
(A2)
where ϵ ¯ ( ρ , s , m 0 ), denoting m 0 = m ρ v s, is the volumetric energy density in the frame of reference co-moving with the superfluid velocity (to be further specified below).
To obtain model (1), however, we have to transform the momentum density and volumetric entropy density to the overall velocity v = m / ρ and entropy density per mass, s = s ¯ / ρ. Derivatives of the energy density with respect to momentum density and superfluid velocity are
v n = e ¯ m = ϵ ¯ m 0 + v s .
(A3)
Because the energy ϵ ¯ depends only on the norm of m 0, not on its direction,7 we can introduce an auxiliary field called normal density ρn,
ρ n = def m 0 v n v s ,
(A4)
because the nominator and denominator are collinear due to Eq. (A3). This, in turn, allows to define superfluid density as ρ s = def ρ ρ n, which gives the usual formula m = ρ s v s + ρ n v n. It should be kept in mind, however, that ρn and ρs are just auxiliary variables determined from the actual form of energy, not state variables. This can be seen for instance from that there are two scalar equations in Eqs. (A1) (for the overall density and entropy density), and so if auxiliary field ρn were taken as a state variable, one scalar evolution equation would be missing. Nevertheless, it is customary in the two-fluid models of superfluid helium-4 to discuss superfluid and normal densities,29 so we formulate our model in terms of the auxiliary fields for better readability.
When we have experimental data on the dependence ρ n ρ ( T ), which is often the case,3 Eq. (A3) actually represents a differential equation for e ¯ 0
ρ n ( ρ , s ) ρ = m 0 ρ ( ϵ ¯ m 0 ) ρ , s .
(A5)
A general solution to this equation is
ϵ ¯ = m 0 2 2 ρ n ( ρ , s ) + e ¯ 0 ( ρ , s ) ,
(A6)
where e ¯ 0 ( ρ , s ) is an integration constant (volumetric energy density in the absence of any flow). In terms of the auxiliary densities (ρn and χn), the energy density becomes
ϵ ¯ = 1 2 χ n ( ρ , s ) v n 2 + 1 2 χ s ( ρ , s ) v s 2 + e ¯ 0 ( ρ , s ) = 1 2 ρ χ n ( ρ , s ) ( v n v s ) 2 + e ¯ 0 ( ρ , s ) ,
(A7)
and the overall energy density per mass is then
e ( ρ , v , s , v s ) = 1 2 v 2 + 1 2 χ s χ n ( v v s ) 2 + e 0 ( ρ , s ) ,
(A8)
where e 0 = e ¯ 0 ( ρ , ρ s ) / ρ (energy density per mass in the absence of any flow) is to be determined from thermodynamic data.
The derivative of energy with respect to the superfluid density reads
ρ e v s = e ¯ v s = ρ ϵ ¯ m 0 + m ρ v s = ρ s v n s ,
(A9)
where v n s = v n v s is the counterflow velocity. The remaining two derivatives (with respect to ρ and s) transform as
e ¯ ρ = e + ρ ( e ρ e s s ρ e v · v ρ ) ,
(A10)
e ¯ s ¯ = T = e s .
(A11)
Formulas (A3), (A9), (A10), and (A11) permit us to write the differential of energy density [Eq. (9)],
d e = v n · d v χ s v n s · d v s + p ρ 2 d ρ + T d s .
(A12)
Hamiltonian equations (A1) then transform to model (1).

In summary, the reversible part of our model [Eqs. (1)] together with the identification of derivatives of energy (9) and a form of the energy itself (determined up to the internal energy) are implied by Hamiltonian mechanics of superfluid helium-4 [Eqs. (A1)] and by the requirement of Galilean invariance. The model is more precise than the standard Landau–Tisza model because it takes into account at least the reversible interaction between the field of superfluid velocity with itself due to the quantum vortices. This results, for instance, in the presence of convective derivative t + v · in evolution equations for all the state variables, instead of having the superfluid velocity v s convected only by itself.53 

In this appendix, we describe a heuristic derivation of our discrete model from the continuous equations and prove the conservation of energy and momentum and entropic inequality. To simplify, we split the ordinary differential system (B19) into reversible and irreversible parts and establish conservation laws for each separately. This procedure works since a conservation law for a generic ordinary differential equation
y ̇ = f ( y )
(B1)
can be formulated as a geometric condition that f belongs to the tangent space of a particular manifold. If two different right-hand sides f 1 and f 2 satisfy the condition, then so does their sum. A similar argument can be used for entropic inequality when we replace tangent space with a half-space.

1. Reversible part

We shall begin with the reversible part of dynamics
D ρ D t = ρ · v , D s D t = 1 ρ · ( ρ s χ s v n s ) , D v D t = 1 ρ · ( ρ χ n χ s v n s v n s + p I ) , D v s D t = χ n v n T v n s p ρ + s T .
(B2)

Let us imagine a classical solution to these equations, which lives in a time-dependent domain Ωt and such that the boundary δ Ω t is an adiabatic free surface, where we assume v n s · n = 0 and p = 0. (Pressure is only relevant up to an additive constant. We fix this constant so that vapor pressure is zero.) The adiabatic free surface is our choice of do-nothing boundary condition. (Implementation of other boundary conditions, such as at the walls of a cryostat, necessitates specific treatment.)

Multiplying equations in (B2) by arbitrary smooth test functions ρ ̃ , s ̃ , v ̃ , and v ̃ s and integrating by parts, we infer
Ω t D ρ D t ρ ̃ d m = Ω t ρ ρ ̃ · v d m , Ω t D s D t s ̃ d m = Ω t s χ s v n s · s ̃ d m , Ω t D v D t · v ̃ d m = Ω t χ n χ s ( v n s v n s ) : v ̃ d m + Ω t p ρ · v ̃ d m , Ω t D v s D t · v ̃ s d m = Ω t χ n ( v n s v ̃ s ) : v n d m + ( p ρ · v ̃ s + s T · v ̃ s ) d m ,
(B3)
where d m = ρ d V. In the next step, we approximate all integrals using particles as quadrature nodes
Ω t φ d m a m a φ a , φ
(B4)
and replace all derivatives by a discrete SPH operator
i φ | r = r a { i φ } a ,
(B5)
which approximates i using particle positions and field values and which we specify later. Hence, we find the following “discrete weak formulation:”
a m a ρ ̇ a ρ ̃ a = a m a ρ a ρ ̃ a { · v } a , a m a s ̇ a s ̃ a = a m a s a χ s , a v n s , a · { s ̃ } a , a m a v ̇ a · v ̃ a = a m a χ n , a χ s , a ( v n s , a v n s , a ) : { v ̃ } a + a m a p a ρ a { · v ̃ } a , a m a v ̇ s , a · v ̃ s , a = a m a χ n , a ( v n s , a v ̃ s , a ) : { v n } a + a m a ( p a ρ a { · v ̃ s } a + s a v ̃ s , a · { T } a )
(B6)
for all values ρ ̃ a , s ̃ a , v ̃ a , and v ̃ s , a. The next logical step is to perform a discrete analogy of integration by parts and localization and express ρ ̃ a , s ̃ a , v ̃ a , and v ̃ s , a explicitly. Before we do that, however, we would like to use (B6) to prove conservation properties.
Definition 1 (order of consistency). A discrete operator { i · } a is k-order consistent if
{ i φ } a = i φ ( r a )
(B7)
for every polynomial of degree at most k.
Proposition 1 (conservation laws). If { i · } a is 0-order consistent for all a, then the solution of (B6) satisfies conservation of energy
d d t ( a m a e a ) = 0 ,
(B8)
entropy
d d t ( a m a s a ) = 0 ,
(B9)
and momentum
d d t ( a m a v a ) = 0 .
(B10)
If { i · } a is 1-order consistent for all a, then the solution also conserves angular momentum:
d d t ( a m a r a × v a ) = 0 .
(B11)
Proof. For conservation of entropy, we set s ̃ a = 1 , a in (B6). The right-hand side of entropy balance vanishes, since { 1 } a = 0. For momentum, we proceed similarly, but instead, we choose
v ̃ = ( 1 0 0 ) a
(B12)
to show conservation of the first component of momentum, and so on. For energy, we remind that
d e = v n · d v χ s v n s · d v s + p ρ 2 d ρ + T d s .
(B13)
Therefore, when we substitute in v ̃ = v n , v ̃ s = χ s v n s , ρ ̃ = p ρ 2 , and s ̃ = T and sum the four equations, we obtain the desired conservation law.
If our discrete derivative is also first order consistent, then we set
v ̃ = ( 0 z a y a ) a .
(B14)
The right-hand side of momentum equation vanishes due to the anti-symmetry of
{ v ̃ } a = ( 0 0 0 0 0 1 0 1 0 ) a
(B15)
and the fact that the double dot product of symmetric and anti-symmetric matrix equals zero. This demonstrates conservation of x-component of angular momentum. Conservation of y and z components follows analogically. □
In this paper, in sake of simplicity, we use a simple and standard SPH operator
{ i φ } a : = 1 ρ a b m b w a b r a b x a b i
(B16)
(xi denotes ith component of vector r), which is merely 0-order consistent. Thus, we conserve all quantities above except for angular momentum. However, operators of order 1 and higher exist and can be found in the literature.37 
Proposition 2 (dual operator). For any pair of discrete variables φ , ψ and any coordinate index i,
a m a ρ a { i φ } a ψ a = a m a ρ a φ a { i ψ } a * ,
(B17)
where { i φ } a is given by (B15,) and
{ i ψ } a * : = ρ a b m b ( ψ a ρ a 2 + ψ b ρ b 2 ) w a b r a b x a b i .
(B18)
Proof. Can be found in Ref. 58. □
Using the dual operator, we can easily solve the algebraic system (B6) and we immediately obtain
ρ ̇ a = ρ a { · v } a , s ̇ a = 1 ρ a { · ( ρ s χ s v n s ) } a * , v ̇ a = 1 ρ a { · ( ρ χ n χ s v n s v n s + p I ) } a * , v ̇ s , a = χ n , a { v n T } a v n s , a 1 ρ a { p } a * + s a { T } a ,
(B19)
and this is precisely the reversible component of (B19).

2. Irreversible part

Next, we turn our attention to the irreversible part of dynamics, which is
ρ t = 0 , s t = β ρ Δ T + ζ ρ T , v t = 2 μ ρ · D n , v s t = s T .
(B20)

We will need two additional operators, both of which are well known and established.

Definition 2 (discrete Laplace operator36).
{ Δ φ } a : = 2 b m b ρ b w a b r a b φ a b .
(B21)
Definition 3 (Monaghan's viscosity operator33).
{ · D } a M : = ( d + 2 ) b m b ρ b v n , a b · r a b r a b 2 + η 2 .
(B22)

In analogy to (B18), we define

Definition 4 (dissipation duals).
{ | φ | 2 } a : = b m b ρ b w a b r a b ( φ a b ) 2 ,
(B23)
{ | D | 2 } a M : = d + 2 2 b m b ρ b ( v n , a b · r a b ) 2 r a b 2 + η 2 w a b r a b
(B24)
such that
Proposition 3.
a m a ρ a { Δ φ } a φ a = a m a ρ a { | φ | 2 } a ,
(B25)
a m a ρ a v n , a · { · D } a M = a m a ρ a { | D | 2 } a M .
(B26)
Now, we are able to write the discrete irreversible part of dynamics as follows:
ρ ̇ a = 0 , s ̇ a = β ρ a { Δ T } a + ζ a T a ρ a , v ̇ a = 2 μ ρ a { · D } a M , v s ̇ a = 0 ,
(B27)
where
ζ a = β { | T | 2 } a + 2 μ { | D | 2 } a M .
(B28)
This gives the remaining irreversible component of (B19).
Proposition 4 (conservation laws—irreversible part). Solution of irreversible equations (B27) and (B28) satisfies conservation of energy
d d t ( a m a e a ) = 0 ,
(B29)
momentum
d d t ( a m a v a ) = 0 ,
(B30)
and angular momentum
d d t ( a m a r a × v a ) = 0.
(B31)
Also, it satisfies entropic inequality:
d d t ( a m a s a ) 0.
(B32)
Proof. Energy conservation law follows immediately from
e ̇ a = v n , a · v a ̇ χ s , a v n s , a · d v s , a + p a ρ a 2 ρ ̇ a + T a s ̇ a ,
(B33)
Eqs. (B27), and the duality relationships (B25) and (B26). As in classical SPH method for Navier–Stokes equations, conservation of momentum follows from
a m a ρ a { · D } a M = 2 ( d + 2 ) a , b m a m b ρ a ρ b v n , a b · r a b r a b 2 + η 2 a b w = 0
(B34)
(the summand is anti-symmetric in a, b). Similarly, conservation of angular momentum is obtained from
a m a ρ a r a × { · D } a M = 2 ( d + 2 ) a , b m a m b ρ a ρ b v n , a b · r a b r a b 2 + η 2 ( r a b × a b w ) = 0.
(B35)
Same technique can be used to show that
a m a β ρ a { Δ T } a = 0 .
(B36)
Since w a b is always non-positive, ζ a 0 for every a. Therefore, we obtain the entropic inequality. □
1.
J. F.
Allen
and
H.
Jones
, “
New phenomena connected with heat flow in helium II
,”
Nature
141
(
3562
),
243
244
(
1938
).
2.
M. L.
Amigó
,
T.
Herrera
,
L.
Neñer
,
L.
Peralta Gavensky
,
F.
Turco
, and
J.
Luzuriaga
, “
A quantitative experiment on the fountain effect in superfluid helium
,”
Eur. J. Phys.
38
(
5
),
055103
(
2017
).
3.
C. F.
Barenghi
,
L.
Skrbek
, and
K. R.
Sreenivasan
, “
Introduction to quantum turbulence
,”
Proc. Natl. Acad. Sci. U. S. A.
111
,
4647
4652
(
2014
).
4.
C. F.
Barenghi
,
Introduction to Superfluid Vortices and Turbulence
(
Springer
,
Berlin/Heidelberg
,
2001
).
5.
I. L.
Bekarevich
and
I. M.
Khalatnikov
, “
Phenomenological derivation of the equations of vortex motion in He II
,”
Sov. Phys. JETP
13
(
3
),
67
(
1961
).
6.
R. J.
Donnelly
and
C. F.
Barenghi
, “
The observed properties of liquid helium at the saturated vapor pressure
,”
J. Phys. Chem. Ref. Data
27
(
6
),
1217
1274
(
1998
).
7.
I. E.
Dzyaloshinskii
and
G. E.
Volovick
, “
Poisson brackets in condensed matter physics
,”
Ann. Phys.
125
(
1
),
67
97
(
1980
).
8.
P.
Español
and
M.
Revenga
, “
Smoothed dissipative particle dynamics
,”
Phys. Rev. E
67
,
026705
(
2003
).
9.
J.
Evans
and
T.
Hughes
, “
Isogeometric divergence-conforming B-splines for the steady Navier-Stokes equations
,”
J. Comput. Phys.
241
,
141
167
(
2013
).
10.
L.
Galantucci
,
A. W.
Baggaley
,
C. F.
Barenghi
, and
G.
Krstulovic
, “
A new self-consistent approach of quantum turbulence in superfluid helium
,”
Eur. Phys. J. Plus
135
(
7
),
547
(
2020
).
11.
F.
Gay-Balmaz
and
T. S.
Ratiu
, “
The geometric structure of complex fluids
,”
Adv. Appl. Math.
42
(
2
),
176
275
(
2009
).
12.
U.
Ghia
,
K. N.
Ghia
, and
C. T.
Shin
, “
High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method
,”
J. Comput. Phys.
48
(
3
),
387
411
(
1982
).
13.
G.
Greenstein
, “
Superfluid turbulence in neutron stars
,”
Nature
227
(
5260
),
791
(
1970
).
14.
A.
Guthrie
,
S.
Kafanov
,
M. T.
Noble
,
Y. A.
Pashkin
,
G. R.
Pickett
,
V.
Tsepelin
,
A. A.
Dorofeev
,
V. A.
Krupenin
, and
D. E.
Presnov
, “
Nanoscale real-time detection of quantum vortices at millikelvin temperatures
,”
Nat. Commun.
12
(
1
),
2645
(
2021
).
15.
H. E.
Hall
,
W. F.
Vinen
, and
D.
Shoenberg
, “
The rotation of liquid helium II II. The theory of mutual friction in uniformly rotating helium II
,”
Proc. R. Soc. London, Ser. A
238
(
1213
),
215
234
(
1956
).
16.
H.
London
, “
Thermodynamics of the thermomechanical effect of liquid He II
,”
Proc. R. Soc. London, Ser. A
171
,
484
496
(
1939
).
17.
D. D.
Holm
,
Introduction to HVBK Dynamics
(
Springer
,
Berlin/Heidelberg
,
2001
), pp.
114
130
.
18.
D. D.
Holm
and
B. A.
Kupershmidt
, “
Superfluid plasmas: Multivelocity nonlinear hydrodynamics of superfluid solutions with charged condensates coupled electromagnetically
,”
Phys. Rev. A
36
,
3947
3956
(
1987
).
19.
R.
Hänninen
and
A. W.
Baggaley
, “
Vortex filament method as a tool for computational visualization of quantum turbulence
,”
Proc. Natl. Acad. Sci.
111
,
4667
4674
(
2014
).
20.
T.
Kamppinen
and
V. B.
Eltsov
, “
Nanomechanical resonators for cryogenic research
,”
J. Low Temp. Phys.
196
(
1–2
),
283
292
(
2019
).
21.
I. M.
Khalatnikov
,
An Introduction to the Theory of Superfluidity
(
Avalon Publishing
,
1989
).
22.
T. W. B.
Kibble
, “
Topology of cosmic domains and strings
,”
J. Phys. A
9
(
8
),
1387
1398
(
1976
).
23.
O.
Kincl
and
M.
Pavelka
, “
Globally time-reversible fluid simulations with smoothed particle hydrodynamics
,”
Comput. Phys. Commun.
284
,
108593
(
2023
).
24.
K.
Kolatzki
,
M. L.
Schubert
,
A.
Ulmer
,
T.
Müller
,
D.
Rupp
, and
R. M. P.
Tanyag
, “
Micrometer-sized droplets from liquid helium jets at low stagnation pressures
,”
Phys. Fluids
34
(
1
),
012002
(
2022
).
25.
M.
La Mantia
, “
Particle dynamics in wall-bounded thermal counterflow of superfluid helium
,”
Phys. Fluids
29
(
6
),
065102
(
2017
).
26.
M.
La Mantia
,
P.
Svancara
,
D.
Duda
, and
L.
Skrbek
, “
Small-scale universality of particle dynamics in quantum turbulence
,”
Phys. Rev. B
94
(
18
),
184512
(
2016
).
27.
L.
Landau
, “
Theory of the superfluidity of helium II
,”
Phys. Rev.
60
,
356
358
(
1941
).
28.
L. D.
Landau
and
E. M.
Lifschitz
,
Statistical Physics, Part 1 in Course of Theoretical Physics
(
Pergamon Press
,
1969
).
29.
L. D.
Landau
and
E. M.
Lifshitz
,
Fluid Mechanics
(
Elsevier Science
,
2013
).
30.
P.
Lebrun
and
L.
Tavian
, “
Cooling with superfluid helium
,” arXiv:1501.07156 (
2015
).
31.
S. J.
Lind
,
B. D.
Rogers
, and
P. K.
Stansby
, “
Review of smoothed particle hydrodynamics: Towards converged Lagrangian flow modelling
,”
Proc. R. Soc. A
476
(
2241
),
20190801
(
2020
).
32.
S. J.
Lind
,
R.
Xu
,
P. K.
Stansby
, and
B. D.
Rogers
, “
Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves
,”
J. Comput. Phys.
231
(
4
),
1499
1523
(
2012
).
33.
J. J.
Monaghan
, “
Smoothed particle hydrodynamics
,”
Annu. Rev. Astron. Astrophys.
30
(
1
),
543
574
(
1992
).
34.
J. J.
Monaghan
,
H. E.
Huppert
, and
M. G.
Worster
, “
Solidification using smoothed particle hydrodynamics
,”
J. Comput. Phys.
206
(
2
),
684
705
(
2005
).
35.
M. S.
Mongiovì
,
D.
Jou
, and
M.
Sciacca
, “
Non-equilibrium thermodynamics, heat transport and thermal waves in laminar and turbulent superfluid helium
,”
Phys. Rep.
726
,
1
71
(
2018
).
36.
J. P.
Morris
,
P. J.
Fox
, and
Y.
Zhu
, “
Modeling low Reynolds number incompressible flows using SPH
,”
J. Comput. Phys.
136
(
1
),
214
226
(
1997
).
37.
G.
Oger
,
M.
Doring
,
B.
Alessandrini
, and
P.
Ferrant
, “
An improved SPH method: Towards higher order convergence
,”
J. Comput. Phys.
225
(
2
),
1472
1492
(
2007
).
38.
O.
Outrata
,
M.
Pavelka
,
J.
Hron
,
M.
La Mantia
,
J. I.
Polanco
, and
G.
Krstulovic
, “
On the determination of vortex ring vorticity using Lagrangian particles
,”
J. Fluid Mech.
924
,
A44
(
2021
).
39.
S. J.
Putterman
,
Superfluid Hydrodynamics
(
North Holland, Amsterdam
,
1974
).
40.
B.
Rousset
,
P.
Bonnay
,
P.
Diribarne
,
A.
Girard
,
J. M.
Poncet
,
E.
Herbert
,
J.
Salort
,
C.
Baudet
,
B.
Castaing
,
L.
Chevillard
,
F.
Daviaud
,
B.
Dubrulle
,
Y.
Gagne
,
M.
Gibert
,
B.
Hébral
,
T.
Lehner
,
P. E.
Roche
,
B.
Saint-Michel
, and
M.
Bon Mardion
, “
Superfluid high Reynolds von Karman experiment
,”
Rev. Sci. Instrum.
85
(
10
),
103908
(
2014
).
41.
J.
Salort
,
C.
Baudet
,
B.
Castaing
,
B.
Chabaud
,
F.
Daviaud
,
T.
Didelot
,
P.
Diribarne
,
B.
Dubrulle
,
Y.
Gagne
,
F.
Gauthier
,
A.
Girard
,
B.
Hebral
,
B.
Rousset
,
P.
Thibault
, and
P. E.
Roche
, “
Turbulent velocity spectra in superfluid flows
,”
Phys. Fluids
22
(
12
),
125102
(
2010
).
42.
H.
Sanavandi
,
S.
Bao
,
Y.
Zhang
,
R.
Keijzer
,
W.
Guo
, and
L. N.
Cattafesta
, III
, “
A cryogenic-helium pipe flow facility with unique double-line molecular tagging velocimetry capability
,”
Rev. Sci. Instrum.
91
,
053901
(
2020
).
43.
D.
Schmoranzer
,
M. J.
Jackson
,
S.
Midlik
,
M.
Skyba
,
J.
Bahyl
,
T.
Skokankova
,
V.
Tsepelin
, and
L.
Skrbek
, “
Dynamical similarity and instabilities in high-Stokes-number oscillatory flows of superfluid helium
,”
Phys. Rev. B
99
(
5
),
054511
(
2019
).
44.
J.
Schumacher
and
K. R.
Sreenivasan
, “
Colloquium: Unusual dynamics of convection in the Sun
,”
Rev. Mod. Phys.
92
(
4
),
041001
(
2020
).
45.
S.
Shukla
,
A. K.
Verma
,
V.
Shukla
,
A.
Bhatnagar
, and
R.
Pandit
, “
Inertial particles in superfluid turbulence: Coflow and counterflow
,”
Phys. Fluids
35
(
1
),
015153
(
2023
).
46.
V.
Shukla
,
R.
Pandit
, and
M.
Brachet
, “
Particles and fields in superfluids: Insights from the two-dimensional Gross-Pitaevskii equation
,”
Phys. Rev. A
97
,
013627
(
2018
).
47.
L.
Skrbek
,
D.
Schmoranzer
,
S.
Midlik
, and
K. R.
Sreenivasan
, “
Phenomenology of quantum turbulence in superfluid helium
,”
Proc. Natl. Acad. Sci. U. S. A.
118
(
16
),
e2018406118
(
2021
).
48.
A.
Henderson Squillacote
,
J.
Ahrens
,
C.
Law
,
B.
Geveci
,
K.
Moreland
, and
B.
King
,
The Paraview Guide
(
Kitware
,
Clifton Park, NY
,
2007
), Vol.
366
.
49.
C. F.
Squire
,
Macroscopic Theory of Superfluid Helium, Superfluids
(
Wiley/Chapman & Hall
,
New York/London
,
1954
), Vol.
2
;
C. F.
Squire
Science
121
(
3146
),
554
555
(
1955
).
50.
M.
Sỳkora
,
M.
Pavelka
,
M.
La Mantia
,
D.
Jou
, and
M.
Grmela
, “
On the relations between large-scale models of superfluid helium-4
,”
Phys. Fluids
33
(
12
),
127124
(
2021
).
51.
L.
Tisza
, “
Transport phenomena in helium II
,”
Nature
141
(
3577
),
913
913
(
1938
).
52.
J. T.
Tough
,
Progress in Low Temperature Physics
(
North-Holland Publ. Co
.,
1982
), Vol.
8
.
53.
M.
Tsubota
,
K.
Fujimoto
, and
S.
Yui
, “
Numerical studies of quantum turbulence
,”
J. Low Temp. Phys.
188
(
5
),
119
189
(
2017
).
54.
S.
Tsuzuki
, “
Particle approximation of the two-fluid model for superfluid 4He using smoothed particle hydrodynamics
,”
J. Phys. Commun.
5
(
3
),
035001
(
2021
).
55.
S.
Tsuzuki
, “
Reproduction of vortex lattices in the simulations of rotating liquid helium-4 by numerically solving the two-fluid model using smoothed-particle hydrodynamics incorporating vortex dynamics
,”
Phys. Fluids
33
(
8
),
087117
(
2021
).
56.
S.
Tsuzuki
, “
Theoretical framework bridging classical and quantum mechanics for the dynamics of cryogenic liquid helium-4 using smoothed-particle hydrodynamics
,”
Phys. Fluids
34
(
12
),
127116
(
2022
).
57.
S. W.
Van Sciver
,
Helium Cryogenics
, International Cryogenics Monograph Series (
Springer
,
New York
,
2012
).
58.
D.
Violeau
,
Fluid Mechanics and the SPH Method: Theory and Applications
(
Oxford University Press
,
Oxford, UK
,
2012
).
59.
G. E.
Volovik
, “
Poisson bracket scheme for vortex dynamics in superfluids and superconductors and the effect of the band structure of the crystal
,”
J. Exp. Theor. Phys. Lett.
64
(
11
),
845
852
(
1996
).
60.
G. E.
Volovik
and
V. S.
Dotsenko
, Jr.
, “
Poisson brackets and continuous dynamics of the vortex lattice in rotating He II
,”
JETP Lett.
29
(
10
),
5
(
1979
).
61.
G. E.
Volovik
and
V. S.
Dotsenko
, Jr.
, “
Hydrodynamics of defects in condensed media, using as examples vortices in rotating He II and disclinations in a planar magnet
,”
Sov. Phys.-JETP
51
(
1
),
65
(
1980
).
62.
G. E.
Volovik
,
The Universe in a Helium Droplet
, International Series of Monographs on Physics (
OUP
,
Oxford
,
2009
).
63.
X.
Yang
and
S.-C.
Kong
, “
Smoothed particle hydrodynamics method for evaporating multiphase flows
,”
Phys. Rev. E
96
(
3
),
033309
(
2017
).
64.
P. R.
Zilsel
, “
Liquid helium II: The hydrodynamics of the two-fluid model
,”
Phys. Rev.
79
,
309
313
(
1950
).
65.
W. H.
Zurek
, “
Cosmological experiments in superfluid-helium
,”
Nature
317
(
6037
),
505
508
(
1985
).
Published open access through an agreement with Univerzita Karlova Matematicky ustav