On Determining the Critical Velocity in the Shot Sleeve of a High-Pressure Die Casting Machine Using Open Source CFD

: This paper investigates the critical plunger velocity in high-pressure die casting during the slow phase of the piston motion and how it can be determined with computational ﬂuid dynamics (CFD) in open source software. The melt-air system is modelled via an Eulerian volume-of-ﬂuid approach, treating the air as a compressible perfect gas. The turbulence is treated via a Reynolds-averaged Navier Stokes (RANS) approach that uses the Menter SST k - ω model. Two different strategies for mesh motion are presented and compared against each other. The solver is validated via analytical models and empirical data. A method is then presented to determine the optimal velocity using a two-dimensional (2D) mesh. As a second step, it is then discussed how the results are in line with those obtained for an actual, industrially relevant, three-dimensional (3D) geometry that also includes the ingate system of the die.


Introduction
High-pressure die casting (HPDC) is an important process for the manufacture of high-volume and low-cost automotive components, such as automatic transmission housings, crank cases and gear box components [1][2][3]. Liquid metal, generally aluminium or magnesium, is poured into a shot sleeve chamber and then injected through a complex system of gates and runners and into the die at speeds of between 50 and 100 ms −1 at the ingate, and under pressures as high as 100 MPa [4]. The normal high-pressure die casting process typically consists of three stages: pre-filling, die-filling (the shot), dwell-pressure. The content of this paper will revolve around the first stage only, a schematic for which is shown on Figure 1. This shows a chamber, typically cylindrical, which is also often referred to as a shot sleeve, that is initially partially filled with melt via a pouring hole; subsequently, a plunger, or piston, pushes the melt to the right, forcing it to exit via an outlet, and thence through an ingate to a die cavity (not shown in the figure).
An important aspect of this stage is the nature of the flow that occurs in the chamber, as this affects the soundness of the casting process. It is, by now, well-established that different flow patterns can occur, depending on the combination of plunger velocity profile as a function of time, plunger diameter, initial melt height in the chamber and chamber length; the principal possibilities are summarized in Figure 2. A plunger moving too slowly, as in Figure 2a, results in the right-moving wave reaching the chamber wall on the right and reflecting back; thus, there are waves travelling along the length of the chamber in both directions, preventing proper expulsion of air into the die cavity. On the other hand, if the plunger moves too quickly, as in Figure 2c, it creates large overturning waves on the surface of the metal that facilitate turbulent mixing of the air into the bulk metal. In either case, the outcome is excessive porosity in the final casting. In an ideal situation, as suggested in Figure 2b, the slope of the metal free surface should be directed away from the plunger everywhere along the length of the cylinder and, at the same time, should not be too steep. Moreover, to minimize air entrapment, the molten metal injection process usually consists of slow-and fast-shot phases; in the slow phase, the plunger accelerates slowly until the molten metal fills the upper section of the shot sleeve, and then moves at constant speed until the sleeve is completely filled with molten metal. There turns out to be exactly one critical constant speed that can be measured experimentally [5][6][7][8], or, under certain assumptions, determined analytically [9][10][11][12]; in particular, this requires: assuming a two-dimensional (2D) geometry, an isothermal flow and that the horizontal length scale is much greater than the vertical length scale; and neglecting the air flow, the effects of melt viscosity and turbulence. Under these conditions, it is possible to derive the shallow-water equations and analytical, or at worst quasi-analytical, solutions can be obtained using the method of characteristics. However, through numerical simulation employing computational fluid dynamics (CFD), it has more recently become more possible to relax some of these assumptions, albeit while still using the early analytical results for benchmarking purposes; a fairly up-to-date review is given in [13]. Even so, there appears to be no experimentally validated model which simultaneously considers turbulent, non-isothermal, two-phase flow in three dimensions (3D), which is what we shall provide here. Furthermore, all of the work that has taken into account some of these effects has tended to use commercial software, principally ANSYS Fluent [13][14][15] and Flow-3D [16,17], although MAGMA 5 is also often used for the modelling of high-pressure die casting [18]. On the other hand, there is some work that has used non-commercial software to produce a 2D model [19], but our purpose here is to use open source software, in particular OpenFOAM [20][21][22], to develop a numerical 3D model for the shot-sleeve dynamics; moreover, this work extends the suite of models that we have been developing in OpenFOAM for investigating different aspects of the high-pressure die casting process [18,[23][24][25].
(a) (b) (c) Figure 2. Impact of increasing plunger propagation speed according to [1]. The plunger velocity increases on going from (a-c).
The structure of this paper is as follows. In Section 2, we present the model equations for the 3D, non-isothermal, turbulent, two-phase flow in a shot sleeve, whereas in Section 3 we detail how the numerical solver was developed and tested. In Section 4, we present some indicative results that we have obtained using the new solver, and conclusions are drawn in Section 5.

Theory
As in our earlier CFD modelling of high-pressure die casting [18,23,24], we model the two-phase flow of molten metal and air in the shot sleeve of a high-pressure die casting machine, as shown in Figure 1, by using the volume-of-fluid (VOF) method [26], wherein a transport equation for the VOF function, γ, of each phase is solved simultaneously with a single set of continuity and Navier-Stokes equations for the whole flow field; note also that γ, which is advected by the fluids, can thus be interpreted as the liquid fraction. Considering the molten melt and the air as Newtonian, compressible and immiscible fluids, the governing equations can be written as [18,27,28] where t is the time, U the mean fluid velocity, p the pressure, g the gravity vector, F s the volumetric representation of the surface tension force and T denotes the transpose. In particular, F s is modelled as a volumetric force by the Continuum Surface Force (CSF) method [29]. It is only active in the interfacial region and formulated as F s = σκ∇γ, where σ is the interfacial tension and κ = ∇ · (∇γ/|∇γ|) is the curvature of the interface. The term U r is a supplementary velocity field for compressing the phase-interface introduced by the solving scheme for the γ-field (MULES) [27,30]. The material properties ρ and µ are the density and the dynamic viscosity, respectively, and are given by where the subscripts g and l denote the gas and liquid phases, respectively. We take ρ l , µ g and µ l to be constant, but assume the air to be a compressible perfect gas, i.e., its density changes with pressure and temperature. Hence, the equation of state for our model reads where M is the molecular weight for air, R s is the specific gas constant and T is the temperature, which is governed by the energy equation, written as where K = 1 2 U · U is the kinetic energy, c v g and c v l denote the specific heat capacities at constant volume for the gas and liquid phases, respectively, and α e f f is given by where λ g and λ l denote the thermal conductivities of the gas and liquid phases, respectively, and σ tur is the turbulent Prandtl number, whose value is set to 0.9 [31]. Furthermore, µ tur in Equations (2) and (8) denotes the turbulent eddy viscosity, which will be calculated via the Menter k-ω-SST model [32]; this turbulence model has been shown to be robust for a wide range of applications, and also to provide results that are in excellent agreement with experimental data [33]. A compact summary of the equations for this turbulence model is given in [23]. At this point, it is also worth mentioning why we have included the effects of turbulence, even though earlier models have not and even though we will later on benchmark against quasi-analytical solutions of the laminar and inviscid shallow-water equations. Although the latter may serve as good approximations at early times, this will not be the case for the situation in Figure 2c, or even in Figure 2a, once the reflecting wave hits the plunger surface. Furthermore, the air cannot simply be considered as a passive entity in these situations and must itself be modelled if one is to evaluate the extent of entrapment. We also point out that the model being developed here can considered as the front-end to our earlier models for the flow in the die [18,23,24], which is most certainly turbulent. Lastly, we should note that, in view of the properties of the melt, the size of the geometry and the plunger speed -all of which are quantitified in Section 3-the typical value of the Reynolds number would be at least in the order of 10 5 ; it is therefore natural to expect turbulence to be present. Equations (1)-(7) require boundary conditions. In the situation under consideration, there are essentially two types of boundaries: an outlet and a wall. At the former, which is labelled by "Gate to the die cavity" in Figure 1, we set zero normal gradient conditions for the velocity, the temperature and the liquid fraction; thus, respectively, where n O denotes the outward unit normal vector at the outlet. We can comment that Equations (9)-(11) would be quite natural conditions for the case of an outlet to a long channel [31] (pp. 39, 76), but are perhaps harder to motivate in the present case; however, as we are not here modelling the die cavity itself, this seems to be the simplest option. In addition, for the turbulence-model quantities, we set where k is the turbulence kinetic energy and ω is the specific rate of dissipation of the turbulence kinetic energy. As regards the wall-type boundary conditions, we have no-slip and no normal-flow conditions at the fixed chamber walls, so that in addition, the handling of the boundary conditions for k and ω is documented in [32,34]. At the plunger wall, we set where U p (t) denotes time-dependent plunger velocity and e x denotes the unit vector in the x-direction; note that the plunger wall is located at it is necessary to set conditions for the turbulence quantities relevant to the turbulence model, i.e., k and ω, at this boundary; this is done via the turbulence intensity and length scale for the largest eddies, which were set at 0.01% and 2 mm, respectively, with the rationale for these values being as follows. Although we would not expect much turbulence at the plunger surface, we cannot simply set the turbulence intensity to be zero there: for example, doing so for the k-ε turbulence model leads to µ tur being indeterminate. On the other hand, CFD codes normally set values of between 1% and 6% [31] for the value of the turbulent intensity at an inlet; thus, we have set a compromise value which is significantly lower that 1%, yet greater than zero. As for the length scale, we have merely used the same value as in our earlier models for other aspects of the high-pressure die casting process [18,23]. Lastly, for all walls, we set zero normal gradient conditions for γ and T; thus, where n W denotes the inward unit normal vector at a wall. For the second of these, we assume for simplicity that all walls are adiabatic, which is in line with previous modelling studies that have made this assumption [15,35]. Lastly, we also require initial conditions. We assume that the chamber initially contains melt to a certain height, h 0 , and above that only air. Both phases are at rest and the air is at room temperature, T air . The metal is at the melting temperature, T melt , of liquid molten aluminium alloy 4600, i.e., 823 K. Thus, we set Moreover, it is useful to relate h 0 to the initial fill fraction, ϕ. Referring to Figure 1, we have where H = 2h 0 /H, noting that, for a 2D problem, H is interpreted as the height of the chamber, whereas for a 3D problem, H is interpreted as the diameter of the chamber.
The non-geometrical model parameters are given in Table 1; note that it is usually c p g and c p l , the specific heat capacities at constant pressure, that are tabulated, rather than c v g and c v l . However, for air, the value of the isentropic expansion factor or heat capacity ratio, c p g /c v g , can be assumed to be constant and equal to 1.4 for air in the given temperature range [36,37] with acceptable accuracy. On the other hand, the aluminium melt alloy was treated as incompressible in this model; thus, c v l = c p l . Table 1. Non-geometrical model parameters. The parameters for gas are those for air; those for metal are for the alloy AlSi 9 Cu 3 .

Symbol
Value Units

Development
Given the complexity of the physics involved and also the niche nature of the application, the C++ toolbox OpenFOAM [20][21][22] was used to implement the model, as it is freely available, amenable to extension by the user and scales up very well for industrial application, as extra cores in parallelisation do not require additional licenses.
As regards the solver, it is necessary to distinguish between the physics solver, that solves the partial differential equations (PDEs) given in Section 2, and how to handle the ever-shrinking computational domain. For the former, we may use the solver inter-DyMFoam or compressibleInterDyMFoam, both of which are available inside the standard foundation release of OpenFOAM. Note that the equations that were documented in Section 2 are the more complex ones for the solver compressibleInterDyMFoam. The solver interDyMFoam assumes incompressibility for both phases and is therefore able to simplify the model equations significantly; in this case, all time derivatives of ρ are set to zero, no equation of state is needed, and therefore also no additional partial differential equation for temperature. The solver that was used for the mesh motion in this paper is available in the OpenFOAM foundation release and in a slightly altered form also inside the foam-extend package; details of the mesh motion strategies that we attempted are given Appendix A. The code's overall solution procedure is shown in Figure 3; basically, it follows the same procedure as the classic PIMPLE algorithm, which is a combination of the SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) [38] and PISO (Pressure Implicit with Splitting of Operator) [39] solution procedures, with the addition that, at the beginning of each loop, the PDE for the phase fraction γ is being solved. As outlined in Section 2, already the melt-air system is a rather complex fluid flow problem. However, there is a solver for solving such a problem inside the foundation release of OpenFOAM as well as in the release of the foam-extend package, although the implementation of the two is different. While the compressible flow solver compressibleInter-DyMFoam inside the foundation release uses the equations as stated in Section 2, the solver inside the foam-extend release uses the simplification of modelling the gas as a barotropic fluid, i.e., the equations of state (6) changes to by omitting the influence of temperature and modelling the density of air dependent on pressure via a compressibility factor ψ. Due to the simplified equation of state, it is not necessary for the foam-extend solver to solve the PDE for temperature, Equation (7).

Testing
Upon successful compilation, the solver was benchmarked against results from existing literature. For this, we considered the problem introduced in [9] for a 2D chamber for which the plunger velocity profile solution is given by where the values for α, β, h 0 and H, as well as the chamber length L, are given in Table 2. We have given two sets of values for h 0 , H and L; the first set is for the geometry considered in [35,40], whereas the second is more reflective of the geometry dimensions and initial fill fraction that occur in industrial practice and which we adopt in Section 4 for an actual 3D geometry.
In [9], the solution to this problem was obtained in quasi-analytical form by formulating in terms of the shallow-water theory equations and using the method of characteristics until the wave hits and reflects from the wall at x = L, after which the problem must be solved numerically [12], and via CFD using the PHOENICS code [40]; subsequently, these results were used for CFD code validation in [13][14][15]35], and we will do likewise. First of all, we note that although this problem has mathematical similarities with that of a onedimensional flow of gas in a semi-infinite cylindrical pipe terminating in a plunger [41], as pointed out in [9], it is actually physically more similar to that of bore formation described in [42]; however, in bore formation, there is no limit on the height that the wave can take, and the geometry is of semi-infinite extent in the horizontal direction. In addition, we should observe that the profile (21) was chosen so that the maximum height of the melt is no greater than the chamber height, provided that there is no interaction with the wave that reflects from x = L. Another feature is that t h is much smaller than the value of t, call it t s , when a shock would occur as a consequence of the non-constant part of U p (t), i.e., αβ e αt − 1 . To see this, recall that, as shown in [9], a shock would form at where c 0 = gh 0 , with g as the acceleration due to gravity, i.e., 9.81 ms −2 . Then, since c 0 > 2αβ from Table 2, we arrive at t s = 1.695 s, whereas t h = 1.473 s; thus, the possibility of a shock occurring does not constrain the quasi-analytical solution obtained via the method of characteristics. However, a constraint on its validity does come from the fact that the wave hits and reflects from the wall at x = L at t = t L := L/c 0 (≈ 1.13 s). Thereafter, the quasi-analytical solution is no longer valid for all of the spatial region that remains, which is X p (t) ≤ x ≤ L. As pointed out in [11], the solution must be determined numerically for t > t L , although the quasi-analytical form is still valid for X p (t) ≤ x ≤ 2L − c 0 t, but only until i.e., when t = t * , where this gives t * = 1.5236 s. Moreover, it is until this time that the analytical expression for the height of the melt at the plunger surface, h p (t), remains valid, where as quoted in [43]; Equation (26) can also be obtained via the method of characteristics. For this benchmarking study, we use the first set of geometry and fill fraction data; for these computations, a uniform 220 × 51 mesh was used, and a comparison is made against the results in [35]. The results are given in Figure 4 and show the position of the melt-air interface inside the shot sleeve chamber at various times. From this, it is clear that, at all four documented time steps, the symbols are more or less on top of each other, although one slight deviation is observable: the OpenFOAM result seems to indicate a slightly greater melt height at the plunger interface, and also the melt seems to propagate slightly more quickly through the chamber. However, we have reason to believe that our result is actually more correct. This is because, of the two curves for t = 1.47 s, ours is closer to the chamber ceiling at the plunger surface, recalling that shallow-water theory indicates that this happens at 1.473 s. We therefore conclude that here our solution is slightly more accurate. Note that we also carried out the same computation with a 60 × 30 mesh, which is the same size as the one used originally in [35], but there was no discernible difference in the results; hence, they can be considered to be mesh-independent. Also, and with reference to Figure 4, we can summarize the changing dynamics as follows: • at t = 1.13 s, the wave first hits the wall at x = L and reflects; • at t = 1.473 s, the wave at the plunger surface hits the ceiling of the chamber-from this time onwards, the plunger no longer continues to accelerate, but moves at constant speed; • at t = 1.5236 s, the reflected wave reaches the plunger surface; • at t ≈ 1.57 s (or slightly less), the melt first exits the outlet.  As we were interested in benchmarking the solver on a geometry and initial fill fraction closer to those of ultimate industrial interest, a 2D geometry was created with the second set of values for α, β, h 0 , H and L shown in Table 2. The mesh spacing for this study was 2 mm in the horizontal direction and 1 mm in the vertical direction, corresponding to a 507 × 130 mesh. With these new parameters, we have t L = 1.115 s, t h = 0.865 s, t * = 1.681 s t s = 1.619 s; for this case also, there will be no shock, since t h < t s . The result of this study can be seen in Figure 5, which shows h p as a function of time for 0 ≤ t ≤ t h , after which h p (t) = H; it is clear that the result from numerical simulation of the OpenFOAM model is in excellent agreement with the analytical formula in Equation (26) analytical solution (eq. (26)) OpenFOAM solution Figure 5. The melt height at the plunger boundary; analytical model according to [43] vs. result of the numerical simulation.

Parameter Studies with a 2D Shot Sleeve
Having tested the solver against results from numerical simulations in the literature, it was decided to use it to try to interpret experimental data from the literature. For this, we turn to the study of the optimal velocity for the slow stage in the die casting machine conducted by Garber [5]. Although the original experiments were carried out for an actual 3D geometry, here we will attempt to account for the experimental results with a 2D model. For this purpose, we consider the second of the two geometries already discussed in Section 3.2.
With this shot sleeve, a parameter study was conducted and the results evaluated. An array of 100 different slow phase velocities was tested and the residual air in the sleeve at the point where the melt front reaches the ingate was evaluated. The result of these simulations is shown in Figure 6. The abscissa features the different velocities in ms −1 and the ordinate shows the corresponding amount of trapped air in the sleeve. The time for the evaluation differs as the melt front reaches the outlet at different times for different velocities.
The tested velocity in the slow stage was increased from 0 to 1 ms −1 in increments of 0.01 ms −1 . One sees that, in the beginning, the differences in the results between several of the tested velocities can be quite large. The reason for this pattern is believed to be the random nature of the wave proceeding inside the chamber. If the peak of the wave is somewhere near the plunger when the first melt front leaves the chamber, the residual air inside the shot sleeve will be quite neatly pushed out in front of the metal. This is a desirable state also in the serial process, as the air can then be vented out of the die via the venting system. However, as Figure 6 indicates, already small changes in the velocity can turn the result towards the other side where the wave peak is near the outlet when the melt front propagates through the inlet of the casting. By doing so, it trapped a rather large bubble of air between the plunger and the melt front. This process state is not desirable, as this melt is basically stuck there as it cannot be vented via the de-airing system of the die either. The melt will flow into this venting system first and seal it for the residual air in the process. This process state is therefore to be avoided.  Table 2.
At about 0.13 ms −1 , one observes a stable local minimum in the amount of trapped air. This is due to a state of the system where plunger speed is still much slower than the wave propagation speed in the chamber. But given the system constraints, such as the length of the shot sleeve, the wave can now more or less travel to and fro in the chamber, thus resulting in only a minimal amount of air being trapped behind the melt front, with the majority being ejected in front of the melt-air interface. This stable regime lasts until a plunger velocity of about 0.17 ms −1 is reached.
After this value, the trapped air increases again very steeply until again a stable plateau is reached, starting from the velocity value of about 0.23 ms −1 . There, one observes the opposite effect. The plunger advances now too quickly to allow the wave to travel back and forth inside the shot sleeve. On the contrary, it is interrupted on this journey while the peak is still in a forward location. Thus, a large amount of melt is sealed between the wave peak that is now near the outlet and the plunger surface. This state is more or less stable until a velocity value of 0.32 ms −1 is reached.
After that, the plunger enters into the regime where the optimal velocity is found. It starts at about 0.38 ms −1 , as seen in Figure 6. There, one sees a constantly low level of trapped air that is close to zero. The flow pattern that corresponds to this result is a wave peak that builds up in front of the plunger until the ceiling of the shot sleeve is reached. By applying this pattern, the casting operator can ensure that all of the air that is inside the shot sleeve gets pushed out in front of the melt and hence no melt is sealed behind the outlet. The optimal range continues up to a plunger speed of 0.46 ms −1 . These values match existing experimental data very well, since Garber [5] and Brunnhuber [2] report 0.46 ms −1 as the critical velocity for the given configuration.
After this value, the processes that are occurring change completely. While the amount of air that was trapped behind the wave peak previously determined how much air was left in the chamber, there is no longer a wave travelling through the shot sleeve from this point on. This process goes on until a critical limit is reached and the wave breaks, thus entraining air in the process that is later recorded as residual air in the CFD model. It is an interesting observation that the amount of air that is entrained by this process increases only very gradually. The recommendation for process engineers would therefore be that it is less harmful to let the plunger move a bit faster than the exact critical velocity, as the entrainment process is much less critical compared to the other process of locking an air pocket behind the melt peak. An additional observation further underscores this thought. The contaminated melt is injected into the casting in front of the clean melt. It is therefore a part of the melt that will end up in the de-airing channels of the die. On the other hand, an air bubble that is behind the melt front will most likely end up somewhere inside the casting.

Application to a 3D Shot Sleeve
The solver was then also applied to an actual industrially relevant geometry, i.e., the shot-sleeve of the EA211 crank case at Volkswagen. The diameter of the shot sleeve in this case is 130 mm, whilst the sleeve is 1.015 m long and the fill fraction is 65%; note that this data is identical to that in the original experiments in [5], which then means that we can compare the results of our 3D simulations against actual experimental results, as far as the optimal plunger velocity is concerned. For these computations, we used an initial grid consisting of approximately 15 million 2.5 mm polyhedral cells; as explained in Appendix A, the cell count decreases with time when using the layer-removal meshing strategy. The velocity was here increased by increments of 0.05 ms −1 , starting from 0.15 ms −1 . The plunger velocity as a function of distance is shown in Figure 7a; in fact, this is a profile that is often referred to in the literature [6,17,44,45]. The resulting amount of trapped air is plotted in Figure 7b; in view of the data given in the caption, the figure indicates that a poor choice of plunger velocity can lead to as much as 13% of the initial air volume in the chamber being trapped in the final casting. We may also note that, in the velocity interval 0.2-0.45 ms −1 , the profiles in Figures 6 and 7b are qualitatively very similar. This may be linked to the observations made in [12,43] that the simplication from a 3D cylinder to a 2D rectangle is justified for initial fill fractions of 40-60%; conversely, this does all indicate the need for 3D simulations for initial fill fractions that are outside of this range.
The data shown in Figure 7b was taken from a 3D simulation that also included the particular ingate geometry of the casting; its shape can be seen in Figure 8. The figure shows the profiles for γ at the instant when the melt first enters the ingate of the die for four different first-stage terminal plunger speeds; from these, we clearly see that the amount of air trapped behind the melt front decreases with increasing plunger speed. It should also be noted that this is normally critical for the casting's quality, as this amount of air cannot be extracted from the die; neither a vacuum nor any other venting system can carry this out, as they are typically automatically closed as soon as the melt front approaches.
Finally, Figure 9 shows what happens when a terminal speed is used that is far from the critical one for minimal air entrainment. While Figure 9a shows a process setup where the plunger speed is far lower than the critical speed, Figure 9b shows the other extreme; here, the plunger propagates too fast. In Figure 9a, the initial right-moving wave has had to reflect back from the right wall, and the profile seen is a juxtaposition of the right-moving and left-moving waves. In Figure 9b, the melt accumulates in front of the plunger much faster than the propagating wave can transport the metal away from the melt-plunger interface; the result is a breaking of the wave as soon as the air-melt interface hits the ceiling of the round chamber.

Conclusions
In this paper, we have used freely available open-source software, in the form of OpenFOAM, to model the computationally rather complex flow problem of shot-sleeve dynamics in the slow stage of plunger motion in high-pressure die casting. The documented results are well in line with previous CFD simulations on the matter and also analytical models derived from the shallow-water equations. Two different mesh motion strategies were presented, with a preference for the layer-removal approach. The solver was also ultimately capable of solving the melt flow in 3D shot sleeves with regular industrial ingate systems attached. It was in general possible to reproduce also the results expected for 3D geometries in practical tests. The economical benefits of increasing the velocity in the slow stage of the plunger motion, not only for the sake of less air entrapment, were also pointed out.

Data Availability Statement:
The data presented in this study are available on request from the first author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Mesh Motion Strategies
Two different mesh motion strategies have been applied in order to include the motion of the plunger and the therefore shrinking fluid domain into the model. Both simulations start with the base mesh at t = 0, with Figure A1 showing a simplified schematic for this. The xand y-axes point in the horizontal and vertical directions, respectively, and the boundary adjacent to the plunger is on the left-hand side of the rectangle. The two strategies that shall be explained here, briefly, are mesh compression and layer removal. After a certain interval of time ∆t has elapsed, the mesh looks as in Figure A2 if the compressing mesh motion strategy is applied. The entire domain has shrunk to some fraction of its original length and volume; for ease of presentation, we have taken the fraction to be a half here. One can see that all cells themselves have also shrunk to half of their original horizontal lengths. On the other hand, the alternative strategy-layer removal-compresses only the layer of the mesh closest to the plunger; the basic principle is illustrated in Figure A3. The width of this layer is reduced until a user-defined minimum layer thickness is reached; for illustrative purposes, we have taken it here to be δx/2, where δx is the initial mesh size in the x-direction. Figure A3a shows the situation just before this layer is removed. If the width of this last layer is about to decrease below the user-specified limit, the mesh motion solver removes that layer and adds its residual width to the next layer by adjusting the points of that layer, as shown in Figure A3b. In this paper, the strategy of layer removal was preferred over mesh compression, as it proved to be more stable. Also, a computational speed-up was observed towards the end of the simulation using the layer removal strategy, as the solver had fewer and fewer cells to handle during the later time steps.