Numerical Investigation on Instability Flow Behaviors of Liquid Oxygen in a Feeding Pipeline with a Five-Way Spherical Cavity

The hydrodynamic information of liquid oxygen in the conveying pipeline of cryogenic launch vehicles directly determines the reliability of the operation of the turbopump. A 0.09 MPa anomalous pressure fall phenomenon in the feeding system has been observed during the flight and run test of a cryogenic rocket with four parallel engines. In previous work, we set up a full-scale experimental system with liquid oxygen as media. The anomalous pressure fall was successfully reproduced. Experimental studies of this phenomenon suggest that the problem might be associated with vortices into the five-way spherical cavity structure. The objective of this study was to determine the three-dimensional instability flow by computational methods to identify and better understand the anomalous pressure fall phenomenon. A numerical model developed by the turbulent conservation equations was validated by experimental data. The generation and evolution of vortices into the five-way spherical cavity of feeding pipelines was captured. It was found that the root cause of the instability flow causing the unusual pressure fall is the formation of a spindle-like vortex into the five-way spherical cavity due to disturbance of the inlet liquid oxygen. The results showed that there is a mirror-symmetrical four-vortices structure in the absence of disturbance, in which the liquid oxygen pressure fall with the rise of the Reynolds number is in good agreement with the predicting values calculated by a set of traditional empirical correlations. In the case of the specific operating conditions, it is also consistent with the experimental results. When the disturbance occurs at the inlet of the spherical cavity, the mirror-symmetrical four-vortices structure gradually evolves into the mirror-symmetrical two-vortices structure. When the disturbance is further enhanced, the mirror-symmetrical two-vortices structure merge with each other to form a spindle-like vortex, which is similar to the Rankine vortex structure. The pressure fall on the corresponding side of the spindle-like vortex core reduces abnormally, and is about 0.07 MPa, which is consistent with the experimental data under certain disturbance conditions. Moreover, it was found that the spindle-like vortex is a stable eddy structure, and would continue to exist once it is formed, which could also not disappear with the removal of the disturbance.


Introduction
The distribution behavior of flow in a pipeline system is vital for fluid transportation, especially in cryogenic rockets, which require liquid oxygen oxidizer to be fed from a propellant tank into multiple engines in parallel during flight with a constant mass flow rate. Flow pressure loss must be small systems, which was combined with the pressure fall model for individual parallel channels with a pump curve in a system flow network. The results showed that below some critical inlet, the uniform flow distribution is always stable, and maldistribution does not occur, regardless of the heat flux and flow rate. Chen [20] derived and compared vortex identification methods for detecting vortices in planar velocity fields. It was found that while all methods ( , Q, λ ci , λ 2 criteria) are able to extract strong vortices, their efficiencies in identifying weaker vortices are not necessarily the same. Zhou [21,22] studied the viscous evolution of a single hairpin vortex-like structure with strength and a core diameter in a mean turbulent field of a low Reynolds number channel flow using direct numerical simulation. It was found that the initial vortical structure always develops into a classical hairpin-shaped vortex with a Ω-shaped head and a pair of counter-rotating long quasi-streamwise vortex legs. Roth [23] presented a novel method to extract vortical structures from 3D CFD vector fields automatically. He distinguished locating the core line from calculating attributes of the strength and quality. Jeong [24] proposed a definition of a vortex in an incompressible flow in terms of the eigenvalues of the symmetric tensor S 2 + Ω 2 ; here S and Ω are respectively the symmetric and antisymmetric parts the velocity gradient tensor, u. The definition could capture the pressure minimum in a plane perpendicular to the vortex axis at high Reynolds numbers.
For the propellant feeding system of cryogenic launch vehicles, much work has been performed on the optimal design and unsteady flow characteristics [26][27][28][29][30][31][32][33][34][35][36][37][38]. For H-IIB launch vehicles in Japan, both liquid hydrogen (LH 2 ) feed lines and liquid oxygen (LOX) feed lines each consist of two independent feed pipes from a tank [26]. This configuration was selected to avoid interference between two feed lines. Ahuja [28] reported that conveying components could experience cavitation-based instability flow, thus leading to large-scale shedding of vapor clouds and pressure oscillations. They presented simulations of the diverse unsteady phenomenon associated with the valve and feed systems that included a valve stall, valve timing studies, as well as two different forms of cavitation instabilities. For the propulsion feeding system of Korea Sounding Rocket-III, a feeding pipeline was designed to feed a certain amount of propellant from the propellant tank to the turbopump manifold inlet during combustion [29]. They found that the pressure in the combustion chamber changed abruptly from ambient to 200 psi during the ignition period. The propellant flow rate that flowed into the combustor could be changed greatly in this period according to the pressure difference, which might cause feeding system instability. For X-34 managed by Marshall Space Flight Center, which is a single-stage LOX/rocket propellant 1(RP-1)-fueled launch vehicle, the LOX feed subsystem design is an Inconel tube with two dual axis gimbals, a z-axis pinned gimbal, and a pressure-compensating elbow [30]. Vu [31] illustrated the X-34 main propulsion system configuration, which included liquid oxygen and rocket propellant #1 feedlines. They conducted a CFD analysis to determine the flow field in the new feedlines. The results showed that the new duct design does not have a dramatic effect on the flow at the outlet of the feedline. The introduction of the neck-down and the tighter bend did not grossly affect the flow and since the flow is fully developed at the engine interface, there is no need to use the flow straighteners at the engine inlet. Dill [35] discussed the LOX feed system of the National Launch System-2 vehicle, which has dual LOX primary feedlines. Each line feeds three space transportation main engines. They developed six preliminary designs for a single LOX primary feedline configuration, and discovered that the spider concept is promising. Chen [36] observed a 4000 Hz vibration phenomenon during the test firings of several space shuttle main engines. They suggested that the typical phenomenon might relate to vortex shedding from the vanes within the LOX tee manifold.
Although many investigators have studied the optimal design of the feeding pipeline in cryogenic rockets, instability flow features in the feeding pipeline with the five-way spherical cavity structure are barely reported in public references. In order to better understand the instability flow and reveal the root reason causing the anomalous pressure to fall in the feedline, a 3D numerical model for the 1:1 full-scale feeding pipeline with the five-way spherical cavity structure was established. The validity of the numerical calculation model was verified by the experimental results of our previous study. By much trial calculation, the generation and evolution of the vortex into the five-way spherical cavity were accurately captured. The pressure decrease along the flow direction in the feedline and physical field were analyzed. These results will help to the fault location of the anomalous pressure decrease in the feeding systems of cryogenic multiple engines in parallel.

Description of Instability Flow Phenomenon in LOX Feedlines
According to the telemetry data, it was found that during the operation of cryogenic rockets, the outlet pressure of symmetrical two-branch pipelines suddenly falls abnormally at a certain time, as shown in Figure 1. Before falling down in Figure 1, the outlet pressure of four-branch pipelines remained stable basically. At the moment corresponding to the blue circle, an anomalous pressure decrease of about 0.09 MPa occurred at the inlet of engine II/IV, which was arranged symmetrically. The pressure differences would continue until the four-branch pipelines are shut down, then the outlet pressure of each branch pipeline is restored to the same level.
Energies 2020, 13, x FOR PEER REVIEW 4 of 21 in the feedline and physical field were analyzed. These results will help to the fault location of the anomalous pressure decrease in the feeding systems of cryogenic multiple engines in parallel.

Description of Instability Flow Phenomenon in LOX Feedlines
According to the telemetry data, it was found that during the operation of cryogenic rockets, the outlet pressure of symmetrical two-branch pipelines suddenly falls abnormally at a certain time, as shown in Figure 1. Before falling down in Figure 1, the outlet pressure of four-branch pipelines remained stable basically. At the moment corresponding to the blue circle, an anomalous pressure decrease of about 0.09 MPa occurred at the inlet of engine II/IV, which was arranged symmetrically. The pressure differences would continue until the four-branch pipelines are shut down, then the outlet pressure of each branch pipeline is restored to the same level. To address the problem and achieve fast fault location, in a previous study [1], we reproduced the anomalous pressure fall phenomenon using full-scale liquid oxygen experiments. The experimental results indicated that the pressure fall under certain conditions exists objectively, and is regular. To further understand the instability flow behaviors causing the anomalous pressure fall, the numerical simulation presented in this paper was conducted.

Physical Model
The actual feeding system consisted of a main pipeline, a flowmeter, a five-way spherical cavity, four-branch pipelines, valves, pressure sensors, temperature sensors, and the corrugated structure of each branch pipeline, which was to evenly feed liquid oxygen from a tank to four engines in parallel. Considering the complexity of the system, the geometry structure needed to be simplified to better study the instability flow of liquid oxygen, and capture the evolution of vortex in the flow field, so the physical object was simplified, including a main pipe, a five-way spherical cavity, and four branch pipes, as shown in Figure 2. The main geometric dimensions remained unchanged, such as the diameter, pipe length, five-way spherical cavity size, elbow size, etc., as shown in Table 1. Moreover, the corrugated pipes were simplified as straight pipes equivalently. To address the problem and achieve fast fault location, in a previous study [1], we reproduced the anomalous pressure fall phenomenon using full-scale liquid oxygen experiments. The experimental results indicated that the pressure fall under certain conditions exists objectively, and is regular. To further understand the instability flow behaviors causing the anomalous pressure fall, the numerical simulation presented in this paper was conducted.

Physical Model
The actual feeding system consisted of a main pipeline, a flowmeter, a five-way spherical cavity, four-branch pipelines, valves, pressure sensors, temperature sensors, and the corrugated structure of each branch pipeline, which was to evenly feed liquid oxygen from a tank to four engines in parallel. Considering the complexity of the system, the geometry structure needed to be simplified to better study the instability flow of liquid oxygen, and capture the evolution of vortex in the flow field, so the physical object was simplified, including a main pipe, a five-way spherical cavity, and four branch pipes, as shown in Figure 2. The main geometric dimensions remained unchanged, such as the diameter, pipe length, five-way spherical cavity size, elbow size, etc., as shown in Table 1. Moreover, the corrugated pipes were simplified as straight pipes equivalently.

Governing Equations
The phenomenon of flow instability of liquid oxygen in the feeding system with the five-way spherical cavity was governed by the continuity, momentum, and turbulence equations. The Reynolds-averaged Navier-Stokes equations of the incompressible steady state could be expressed here in Cartesian tensor form as follows [39,40]: Continuity equation: Momentum equation: where ρ is the liquid oxygen density; ui is the mass-averaged velocity; ' ' j i u u is the specific Reynolds stress tensor; and ρai is the gravitational body force:

Governing Equations
The phenomenon of flow instability of liquid oxygen in the feeding system with the five-way spherical cavity was governed by the continuity, momentum, and turbulence equations. The Reynolds-averaged Navier-Stokes equations of the incompressible steady state could be expressed here in Cartesian tensor form as follows [39,40]: Momentum equation: where ρ is the liquid oxygen density; u i is the mass-averaged velocity; u i u j is the specific Reynolds stress tensor; and ρa i is the gravitational body force: Energies 2020, 13, 926 6 of 21 where µ t is the turbulent viscosity, k is the turbulence kinetic energy, and ω is the turbulence frequency.

Turbulence Model (SST k-ω)
The k-ω based on the shear stress transport (SST) model takes into account the transport effects of the turbulent shear stress into the formulation of the eddy viscosity, so it provides highly accurate predictions at the onset or flow separation under adverse pressure gradients [41,42]. For this reason, the SST k-ω model was used in the present paper to accurately predict the formation and evolution of vortices in the feedline with the five-way spherical cavity.

Boundary Conditions
During the flight of cryogenic rockets, in order to ensure that the inlet pressure of the turbopump is greater than the saturation pressure corresponding to the local temperature, the ullage pressure of the tank is kept constant by pressurization. The outlet of each branch pipeline is connected to the inlet of the turbopump, which maintains a constant flow. Therefore, considering actual working conditions, the boundary condition of the total pressure was specified at the main pipeline inlet for the feedline with the five-way spherical cavity, the boundary condition of the mass flow rate was set in each branch pipeline outlet, and the rest was set as the boundary condition of the wall, as shown in Figure 2.

Solution Strategy
Considering the flow process of the complex structure for liquid oxygen, an element-based finite volume method was used, which adopted a mesh to discretize the spatial domain first. A coupled solver ANSYS CFX was selected, which consisted of Navier-Stokes equations (for u, v, w, p). In the solution approach, a fully implicit discretization of the equations was used at any given time step. The advection scheme was set as the high-resolution scheme without introducing new extrema. The timescale control was set as the physical timescale, which was about 0.01. For the initial conditions of domain, the Cartesian velocity components, u, v, and w, were given as 1 × 10 −5 m/s, 1 × 10 −5 m/s, and 1 × 10 −5 m/s, respectively. The relative static pressure was designed as 0 Pa, and The turbulence option was set as medium (intensity = 5%). The convergence was judged by the root-mean-square residual value, which was 10 −6 , satisfying the requirements of the convergence criteria.
Considering that an automatic near-wall treatment method provides a gradual switch between low-Reynolds number grids and wall functions without a loss in accuracy for the SST k-ω mode, it was used in the present paper.
In practical working conditions, there are varying degrees of disturbance upstream or downstream of the five-way spherical cavity structure, which leads to the fact that although the geometry is symmetrical, the flow field might be asymmetric. The most typical case is that the flow field of liquid oxygen flowing out from the bottom of the tank may not be symmetrical. Especially under the condition of high pressure and high Reynolds number in the cryogenic rockets, the turbulent flow field is more easily disturbed. Obviously, the disturbance is random. Therefore, in order to study the effect of disturbance on the flow field in the feeding pipelines, the numerical technique applying an artificial disturbance at the inlet of feeding pipelines was adopted.

Grid Generation
The mesh quality has a notable influence on the accurate numerical solution of the liquid oxygen flow and calculating time. In this study, all hexahedral mesh was generated by ANSYS ICEM CFD's mesh generation tools. A 3D block topology model was created to generate Ogrids to ensure good quality meshes for very complex geometry, as shown in Figure 3. Moreover, the wall grid near the Energies 2020, 13, 926 7 of 21 pipeline surface should be densified to accurately resolve the velocity without a loss of efficiency. To better capture the vortex, the computational domain of the five-way spherical cavity was densified.
Energies 2020, 13, x FOR PEER REVIEW 7 of 21 mesh generation tools. A 3D block topology model was created to generate Ogrids to ensure good quality meshes for very complex geometry, as shown in Figure 3. Moreover, the wall grid near the pipeline surface should be densified to accurately resolve the velocity without a loss of efficiency. To better capture the vortex, the computational domain of the five-way spherical cavity was densified.

Y+ Independence Study
The near-wall treatments for wall-bounded turbulent flows are very important for the fidelity of the numerical solution. A fine mesh with y+ around 1 for highly accurate simulations was recommended in some literature. Nevertheless, in complex practical geometry, a too small Y+ could bring fatal damage to the stability and convergence of the numerical calculation, and the calculating time could be greatly increased. Taking the case of the research object in this paper, in order to achieve the requirement of a Y+ about 1 in typical conditions, the distance between the first and second mesh pointing off the wall needed to be reduced to 1 × 10 −7 mm. The convergence of the calculation might be very unstable, and it is difficult to reduce the residual to below the expected value, even resulting in the calculated divergence. In order to verify the influence of Y+ on the results of this paper, the controlling method of Y+ was adopted by changing the distance of the first grid layer to conduct the independent verification. The distance of the first grid layer was set as 1 × 10 −3 mm, 1 × 10 −4 mm, 5 × 10 −5 mm, 1 × 10 −5 mm, 5 × 10 −6 mm, 1 × 10 −6 mm, and 1 × 10 −7 mm, and the corresponding maximum Y+ value in the typical condition was 8626, 653, 320, 65, 33, 9, and 1, respectively. It can be clearly seen from Figure 4 that with the decrease of Y+, the values of the four outlet pressures were basically unchanged, and the maximum relative error was within ±1%, indicating that Y+ was insensitive to the calculation results in this paper. Therefore, the distance close to the wall was set as 1 × 10 −4 mm, and the corresponding average Y+ was 160, which was used as the wall grid treatment scheme for subsequent calculation.

Y+ Independence Study
The near-wall treatments for wall-bounded turbulent flows are very important for the fidelity of the numerical solution. A fine mesh with y+ around 1 for highly accurate simulations was recommended in some literature. Nevertheless, in complex practical geometry, a too small Y+ could bring fatal damage to the stability and convergence of the numerical calculation, and the calculating time could be greatly increased. Taking the case of the research object in this paper, in order to achieve the requirement of a Y+ about 1 in typical conditions, the distance between the first and second mesh pointing off the wall needed to be reduced to 1 × 10 −7 mm. The convergence of the calculation might be very unstable, and it is difficult to reduce the residual to below the expected value, even resulting in the calculated divergence. In order to verify the influence of Y+ on the results of this paper, the controlling method of Y+ was adopted by changing the distance of the first grid layer to conduct the independent verification. The distance of the first grid layer was set as 1 × 10 −3 mm, 1 × 10 −4 mm, 5 × 10 −5 mm, 1 × 10 −5 mm, 5 × 10 −6 mm, 1 × 10 −6 mm, and 1 × 10 −7 mm, and the corresponding maximum Y+ value in the typical condition was 8626, 653, 320, 65, 33, 9, and 1, respectively. It can be clearly seen from Figure 4 that with the decrease of Y+, the values of the four outlet pressures were basically unchanged, and the maximum relative error was within ±1%, indicating that Y+ was insensitive to the calculation results in this paper. Therefore, the distance close to the wall was set as 1 × 10 −4 mm, and the corresponding average Y+ was 160, which was used as the wall grid treatment scheme for subsequent calculation.

Mesh Independence Study
To check the mesh independence and the self-consistency of the present problem, the grid nodes were gradually increased before the numerical calculation, keeping Y+ unchanged, Y+ about 153, and the number of grids was generated from 0.42 million to 17.2 million. The numerical results with the change of grid number are shown in Figure 5. It is clearly seen from this figure that the calculating results remain unchanged basically after the mesh increases to 9.3 million. The maximum relative deviation between 9.3 million and 17.2 million was less than 0.05%. Taking into account the requirements of both mesh independence and computation time, the scheme of 9.3 million was used to calculate the flow instability into the five-way spherical cavity in the present paper.

Validation of Results
According to the numerical model of the liquid oxygen instability flow established above, numerical simulation was carried out under the experimental conditions, and the predicted numerical calculation results were compared with the calculated results of traditional empirical correlations with the increase of the Reynolds number, as shown in Figure 6. The inlet pressure of the

Mesh Independence Study
To check the mesh independence and the self-consistency of the present problem, the grid nodes were gradually increased before the numerical calculation, keeping Y+ unchanged, Y+ about 153, and the number of grids was generated from 0.42 million to 17.2 million. The numerical results with the change of grid number are shown in Figure 5. It is clearly seen from this figure that the calculating results remain unchanged basically after the mesh increases to 9.3 million. The maximum relative deviation between 9.3 million and 17.2 million was less than 0.05%. Taking into account the requirements of both mesh independence and computation time, the scheme of 9.3 million was used to calculate the flow instability into the five-way spherical cavity in the present paper.

Mesh Independence Study
To check the mesh independence and the self-consistency of the present problem, the grid nodes were gradually increased before the numerical calculation, keeping Y+ unchanged, Y+ about 153, and the number of grids was generated from 0.42 million to 17.2 million. The numerical results with the change of grid number are shown in Figure 5. It is clearly seen from this figure that the calculating results remain unchanged basically after the mesh increases to 9.3 million. The maximum relative deviation between 9.3 million and 17.2 million was less than 0.05%. Taking into account the requirements of both mesh independence and computation time, the scheme of 9.3 million was used to calculate the flow instability into the five-way spherical cavity in the present paper.

Validation of Results
According to the numerical model of the liquid oxygen instability flow established above, numerical simulation was carried out under the experimental conditions, and the predicted numerical calculation results were compared with the calculated results of traditional empirical correlations with the increase of the Reynolds number, as shown in Figure 6. The inlet pressure of the

Validation of Results
According to the numerical model of the liquid oxygen instability flow established above, numerical simulation was carried out under the experimental conditions, and the predicted numerical calculation results were compared with the calculated results of traditional empirical correlations Energies 2020, 13, 926 9 of 21 with the increase of the Reynolds number, as shown in Figure 6. The inlet pressure of the pipeline remained unchanged, which was 0.375 MPa. The inlet flow rates were 150, 200, 250, 300, 350, and 400 L/s, respectively, and the Reynolds number ranged from 3.53 × 10 6 to 1.28 × 10 7 . It can be seen from Figure 6 that the calculated values between the numerical calculation and the empirical correlation were in good agreement, and the maximum error between them was less than 1%.
At the same time, the results of the numerical calculation and empirical correlation calculation had good regularity and consistency with the change of the Reynolds number, which indicated that the numerical model developed in this paper, including the geometry simplification, turbulence model selection, setting of boundary conditions and initial conditions, solver strategy, and so on, was reasonable and reliable.
The structure of the feeding pipeline with the five-way spherical cavity could be divided into five parts, as shown in Figure 6. The first part, the third part, and the fifth part were regarded as circular straight pipes. The calculation correlation of the pressure fall of the straight pipe was given as follows: where λ is the friction factor, l is the length of the pipeline, d is the inner diameter of the pipe, ρ is the liquid oxygen density, and u is the average velocity in the pipe. The friction factor, λ, was calculated by the Karman-Prandtl formula [40]: The five-way spherical cavity could be simplified as a combination of a sudden expansion section and a sudden reduction section. The pressure fall of this section was calculated as follows [40]: where ξ 1 and ξ 2 are the local loss coefficient of the sudden expansion and the sudden reduction of the five-way spherical cavity, respectively; A 1 , A 2 , and A 3 are the cross-sectional areas of the first straight pipe, corresponding to the diameter of the five-way spherical cavity and the third straight pipe, respectively; and ξ 3 is the correction coefficient of the five-way spherical cavity, and the value of ξ 3 is considered as 0.61. The calculation correlation for the local loss of the 90 • elbow was written as follows [40]: Therefore, the total pressure loss of the pipeline might be predicted through a set of empirical correlations developed in this paper, as follows:

Re
Numerical results Figure 6. Verification of the numerical simulation.

Evolution Process of Flow Field Inside the Feedline
Through a large number of numerical calculations, it was found that there are three typical flow fields, mirror-symmetric four-vortices structure, mirror-symmetric two-vortices structure, and spindle-like vortex structure, under different perturbations inside the five-way spherical cavity, as shown in Figure 7. It is clearly seen from this figure that the velocity in the flow field could increase with the vortex merging, so the energy loss would rise accordingly. In order to describe the physical field in the feedline with the five-way spherical cavity more clearly, some sections were defined as shown in Figure 8. Section I/III and section II/IV consisted of the central sections of pipe I and III, and central sections of pipe II and IV, respectively. Section I/II and section I/IV were the faces formed by section I/III rotating 45 degrees to the direction of pipe II, and section I/III rotating 45 degrees to the direction of pipe IV, respectively. The velocity field of the mirror-symmetric four-vortices structure and mirror-symmetric two-vortices structure are shown in Figure 9.

Evolution Process of Flow Field Inside the Feedline
Through a large number of numerical calculations, it was found that there are three typical flow fields, mirror-symmetric four-vortices structure, mirror-symmetric two-vortices structure, and spindle-like vortex structure, under different perturbations inside the five-way spherical cavity, as shown in Figure 7. It is clearly seen from this figure that the velocity in the flow field could increase with the vortex merging, so the energy loss would rise accordingly.

Re
Numerical results Figure 6. Verification of the numerical simulation.

Evolution Process of Flow Field Inside the Feedline
Through a large number of numerical calculations, it was found that there are three typical flow fields, mirror-symmetric four-vortices structure, mirror-symmetric two-vortices structure, and spindle-like vortex structure, under different perturbations inside the five-way spherical cavity, as shown in Figure 7. It is clearly seen from this figure that the velocity in the flow field could increase with the vortex merging, so the energy loss would rise accordingly. In order to describe the physical field in the feedline with the five-way spherical cavity more clearly, some sections were defined as shown in Figure 8. Section I/III and section II/IV consisted of the central sections of pipe I and III, and central sections of pipe II and IV, respectively. Section I/II and section I/IV were the faces formed by section I/III rotating 45 degrees to the direction of pipe II, and section I/III rotating 45 degrees to the direction of pipe IV, respectively. The velocity field of the mirror-symmetric four-vortices structure and mirror-symmetric two-vortices structure are shown in Figure 9. In order to describe the physical field in the feedline with the five-way spherical cavity more clearly, some sections were defined as shown in Figure 8. Section I/III and section II/IV consisted of the central sections of pipe I and III, and central sections of pipe II and IV, respectively. Section I/II and section I/IV were the faces formed by section I/III rotating 45 degrees to the direction of pipe II, and section I/III rotating 45 degrees to the direction of pipe IV, respectively. The velocity field of the mirror-symmetric four-vortices structure and mirror-symmetric two-vortices structure are shown in Figure 9. When the incoming fluid was in a uniform and undisturbed condition, the flow field reached a stable state. The four-vortices structure was formed inside the five-way spherical cavity, which presented a mirror-like symmetry. The most significant feature of the four-vortices flow field was the existence of four complete vortex cores in the five-way spherical cavity. It can be observed from Figure 7 and Figure 9 that two sets of spindle-shaped vortices formed in the direction of 45° of each two branch pipes in the five-way spherical cavity structure, which have the same shape and size. At this time, the outlet pressures of four branch pipes were the same, and there was no pressure difference between the four branch pipes. When the incoming fluid was in a uniform and undisturbed condition, the flow field reached a stable state. The four-vortices structure was formed inside the five-way spherical cavity, which presented a mirror-like symmetry. The most significant feature of the four-vortices flow field was the existence of four complete vortex cores in the five-way spherical cavity. It can be observed from Figure 7 and Figure 9 that two sets of spindle-shaped vortices formed in the direction of 45° of each two branch pipes in the five-way spherical cavity structure, which have the same shape and size. At this time, the outlet pressures of four branch pipes were the same, and there was no pressure difference between the four branch pipes.  When the incoming fluid was in a uniform and undisturbed condition, the flow field reached a stable state. The four-vortices structure was formed inside the five-way spherical cavity, which presented a mirror-like symmetry. The most significant feature of the four-vortices flow field was the existence of four complete vortex cores in the five-way spherical cavity. It can be observed from Figures 7 and 9 that two sets of spindle-shaped vortices formed in the direction of 45 • of each two branch pipes in the five-way spherical cavity structure, which have the same shape and size. At this Energies 2020, 13, 926 12 of 21 time, the outlet pressures of four branch pipes were the same, and there was no pressure difference between the four branch pipes.
When there was a disturbance of the incoming flow, a set of spindle-shaped vortices that were mirror symmetrical formed in the spherical cavity in the direction of rotating 45 • for the central section of the branch pipes. The characteristic of this type of flow field was similar to that of the mirror-symmetric four-vortices structure, and the difference was the number of vortices. It can be clearly observed from Figure 9b that the four vortices evolve into a two-vortices structure, which is still mirror symmetrical on section I/IV. Nevertheless, after the merging of four vortices, the velocity field around a pair of vortices cores increased obviously. The velocity field at the bottom of the five-way spherical cavity was also affected, and enhanced sharply.
When the disturbance further increased, and reached a certain level, the flow field of the mirror-symmetric two-vortices structure will be transformed into the spindle-like vortex structure, causing an obvious phenomenon of an unstable pressure decline. In the state of the spindle-like vortex, the pressure effect was different for each branch pipe outlet. The outlet pressures of two branch pipes corresponding to the vortex core was reduced abnormally, since the spindle-like vortex moved to two branch pipes in a spiral way, and the static pressure was converted into dynamic pressure.

Pressure Change Along Flow Direction Under Different Disturbances
In order to describe the change of pressure in the process of flow along the feedline, several sections were defined in the feeding line, as shown in Figure 10. Table 2 shows the specific coordinates of each section. Taking the flow from the main pipe to branch pipe I as an example, 14 sections were defined from the inlet to the outlet, of which the main pipe section, including four sections, the horizontal straight pipe section and the vertical straight pipe section of the branch pipeline included four sections, respectively, and the elbow section of the branch pipe had two sections.
Energies 2020, 13, x FOR PEER REVIEW 12 of 21 When there was a disturbance of the incoming flow, a set of spindle-shaped vortices that were mirror symmetrical formed in the spherical cavity in the direction of rotating 45° for the central section of the branch pipes. The characteristic of this type of flow field was similar to that of the mirror-symmetric four-vortices structure, and the difference was the number of vortices. It can be clearly observed from Figure 9b that the four vortices evolve into a two-vortices structure, which is still mirror symmetrical on section I/IV. Nevertheless, after the merging of four vortices, the velocity field around a pair of vortices cores increased obviously. The velocity field at the bottom of the fiveway spherical cavity was also affected, and enhanced sharply.
When the disturbance further increased, and reached a certain level, the flow field of the mirrorsymmetric two-vortices structure will be transformed into the spindle-like vortex structure, causing an obvious phenomenon of an unstable pressure decline. In the state of the spindle-like vortex, the pressure effect was different for each branch pipe outlet. The outlet pressures of two branch pipes corresponding to the vortex core was reduced abnormally, since the spindle-like vortex moved to two branch pipes in a spiral way, and the static pressure was converted into dynamic pressure.

Pressure Change Along Flow Direction Under Different Disturbances
In order to describe the change of pressure in the process of flow along the feedline, several sections were defined in the feeding line, as shown in Figure 10. Table 2 shows the specific coordinates of each section. Taking the flow from the main pipe to branch pipe I as an example, 14 sections were defined from the inlet to the outlet, of which the main pipe section, including four sections, the horizontal straight pipe section and the vertical straight pipe section of the branch pipeline included four sections, respectively, and the elbow section of the branch pipe had two sections.    The pressure changes along the feedline are shown in Figure 11 when there was no disturbance. It can be seen obviously from the figure that the pressure fall of each branch pipe is basically consistent, about 16 kPa. The pressure drops between section-3 and section-1 of the branch pipes reached 19 kPa, indicating that the maximum pressure fall occurs in the five-way spherical cavity. When the liquid oxygen flowed through section-1 of the branch pipe, the average pressure started to recover and reached about 285 kPa. This is because there were two spiral vortices at the inlet of the branch pipe, which disappeared gradually along the flow direction, and the circumferential velocity of the flow was converted into static pressure energy. After the spiral vortices disappeared completely, the pressure decreased gradually due to the resistance in the flow process. Meanwhile, it can be clearly observed from this figure that the pressure fall slope along the main pipeline and the branch pipeline is almost identical, which means that the flow field in the main pipeline is almost coincident with that in the branch pipeline. The pressure changes along the feedline are shown in Figure 11 when there was no disturbance. It can be seen obviously from the figure that the pressure fall of each branch pipe is basically consistent, about 16 kPa. The pressure drops between section-3 and section-1 of the branch pipes reached 19 kPa, indicating that the maximum pressure fall occurs in the five-way spherical cavity. When the liquid oxygen flowed through section-1 of the branch pipe, the average pressure started to recover and reached about 285 kPa. This is because there were two spiral vortices at the inlet of the branch pipe, which disappeared gradually along the flow direction, and the circumferential velocity of the flow was converted into static pressure energy. After the spiral vortices disappeared completely, the pressure decreased gradually due to the resistance in the flow process. Meanwhile, it can be clearly observed from this figure that the pressure fall slope along the main pipeline and the branch pipeline is almost identical, which means that the flow field in the main pipeline is almost coincident with that in the branch pipeline. In order to further study the influence of disturbance on the pressure fall, the angle, θ, of disturbance was defined on the XOY plane and the YOZ plane, respectively. Figures 12 and 13 display the change of the outlet pressure with the angle of disturbance, which rotates in the counterclockwise direction with reference to the negative direction of the Y axis. Figure 12 shows that before the angle, θ, less than 31 degrees, there was only a very small change of four outlet pressures. When the angle, θ, was greater than 31 degrees, the outlet pressure of pipe I and pipe III had a significant pressure fall compared with that of pipe II and pipe IV, since the spindle-like vortex was formed. Additionally, once the spindle-like vortex was generated, the pressure fall increased with the rise of the disturbance angle. Figure 13 exhibits a similar trend compared with Figure 12. Nevertheless, the anomalous pressure fall distribution was different at the outlet of four branch pipes. The anomalous pressure fall in Figure 12 was in pipe I and pipe III, and that in Figure 13 was in pipe II/pipe IV. It is indicated that In order to further study the influence of disturbance on the pressure fall, the angle, θ, of disturbance was defined on the XOY plane and the YOZ plane, respectively. Figures 12 and 13 display the change of the outlet pressure with the angle of disturbance, which rotates in the counterclockwise direction with reference to the negative direction of the Y axis. Figure 12 shows that before the angle, θ, less than 31 degrees, there was only a very small change of four outlet pressures. When the angle, θ, was greater than 31 degrees, the outlet pressure of pipe I and pipe III had a significant pressure fall compared with that of pipe II and pipe IV, since the spindle-like vortex was formed. Additionally, once the spindle-like vortex was generated, the pressure fall increased with the rise of the disturbance angle. Figure 13 exhibits a similar trend compared with Figure 12. Nevertheless, the anomalous pressure fall distribution was different at the outlet of four branch pipes. The anomalous pressure fall in Figure 12 was in pipe I and pipe III, and that in Figure 13 was in pipe II/pipe IV. It is indicated that the position of the pressure jump was related to the direction of the inlet disturbance, since the direction of the disturbance was always in the same direction as the spindle-like vortex core.

Inlet
the position of the pressure jump was related to the direction of the inlet disturbance, since the direction of the disturbance was always in the same direction as the spindle-like vortex core.
It can be also observed from the two figures that the numerical simulation results are in good agreement with the experimental data when there was a disturbance of about 45 degrees. The results further indicated that the numerical model developed in this paper was reliable, and could accurately predict the instability flow features of liquid oxygen in the feeding pipe with the five-way spherical cavity. Before the anomalous pressure fall, the flow pressure decreased linearly with the enhancement of the disturbance. Once the anomalous pressure fall occurred, the pressure declined exponentially with the enhancement of the disturbance. This is because the flow field structure in the five-way spherical cavity changed from the mirror-symmetrical four-vortices structure to the spindlelike vortex structure. Exp.  Num.  Energies 2020, 13, x FOR PEER REVIEW 14 of 21 the position of the pressure jump was related to the direction of the inlet disturbance, since the direction of the disturbance was always in the same direction as the spindle-like vortex core. It can be also observed from the two figures that the numerical simulation results are in good agreement with the experimental data when there was a disturbance of about 45 degrees. The results further indicated that the numerical model developed in this paper was reliable, and could accurately predict the instability flow features of liquid oxygen in the feeding pipe with the five-way spherical cavity. Before the anomalous pressure fall, the flow pressure decreased linearly with the enhancement of the disturbance. Once the anomalous pressure fall occurred, the pressure declined exponentially with the enhancement of the disturbance. This is because the flow field structure in the five-way spherical cavity changed from the mirror-symmetrical four-vortices structure to the spindlelike vortex structure. Exp.

Figure 12.
Outlet pressure change with the angle, θ, of disturbance on the YZX plane. Num.  It can be also observed from the two figures that the numerical simulation results are in good agreement with the experimental data when there was a disturbance of about 45 degrees. The results further indicated that the numerical model developed in this paper was reliable, and could accurately predict the instability flow features of liquid oxygen in the feeding pipe with the five-way spherical cavity. Before the anomalous pressure fall, the flow pressure decreased linearly with the enhancement of the disturbance. Once the anomalous pressure fall occurred, the pressure declined exponentially with the enhancement of the disturbance. This is because the flow field structure in the five-way spherical cavity changed from the mirror-symmetrical four-vortices structure to the spindle-like vortex structure. Figures 14 and 15 show the velocity distribution and streamline structure of the spindle-like vortex, respectively. The disturbance of 45 degrees was added to the YZX plane. It can be seen from the two figures that the liquid oxygen velocity gradually increases first and then decreases to zero from the inlet of the spherical cavity to the vortex core, which is similar to the characteristics of the Rankine vortex. The spindle-like vortex revolved around the vortex core into the five-way spherical cavity, part of the liquid was thrown out from the top of the branch pipe II inlet, and the other liquid flowed out from the bottom of the branch pipe IV inlet due to the influence of the incoming flow, except for the outflow along the branch pipe I/III corresponding to the vortex core.

Contours of Spindle-Like Vortex
Energies 2020, 13, x FOR PEER REVIEW 15 of 21 Figures 14 and 15 show the velocity distribution and streamline structure of the spindle-like vortex, respectively. The disturbance of 45 degrees was added to the YZX plane. It can be seen from the two figures that the liquid oxygen velocity gradually increases first and then decreases to zero from the inlet of the spherical cavity to the vortex core, which is similar to the characteristics of the Rankine vortex. The spindle-like vortex revolved around the vortex core into the five-way spherical cavity, part of the liquid was thrown out from the top of the branch pipe II inlet, and the other liquid flowed out from the bottom of the branch pipe IV inlet due to the influence of the incoming flow, except for the outflow along the branch pipe I/III corresponding to the vortex core.   Figures 14 and 15 show the velocity distribution and streamline structure of the spindle-like vortex, respectively. The disturbance of 45 degrees was added to the YZX plane. It can be seen from the two figures that the liquid oxygen velocity gradually increases first and then decreases to zero from the inlet of the spherical cavity to the vortex core, which is similar to the characteristics of the Rankine vortex. The spindle-like vortex revolved around the vortex core into the five-way spherical cavity, part of the liquid was thrown out from the top of the branch pipe II inlet, and the other liquid flowed out from the bottom of the branch pipe IV inlet due to the influence of the incoming flow, except for the outflow along the branch pipe I/III corresponding to the vortex core.    Figures 16 and 17 display the distribution of the pressure on the front side and on the side of the spindle-like vortex, respectively. Meanwhile, the disturbance of 45 degrees was added to the YZX plane. When the spindle-like vortex structure was generated in the five-way spherical cavity, the pressure gradually decreased from the outside to the inside. There was a low-pressure region in the center of the vortex core of the spindle-like vortex, which was distributed along the axis of the vortex core. This was also the cause of the anomalous pressure fall in the branch pipelines on the symmetrical side. In addition to the five-way spherical cavity, the pressure fall of other parts conformed to the characteristics of the basic flow dynamics. That is, the pressure declined gradually along the flow direction in the feeding pipeline.

Contours of Spindle-Like Vortex
Energies 2020, 13, x FOR PEER REVIEW 16 of 21 Figures 16 and 17 display the distribution of the pressure on the front side and on the side of the spindle-like vortex, respectively. Meanwhile, the disturbance of 45 degrees was added to the YZX plane. When the spindle-like vortex structure was generated in the five-way spherical cavity, the pressure gradually decreased from the outside to the inside. There was a low-pressure region in the center of the vortex core of the spindle-like vortex, which was distributed along the axis of the vortex core. This was also the cause of the anomalous pressure fall in the branch pipelines on the symmetrical side. In addition to the five-way spherical cavity, the pressure fall of other parts conformed to the characteristics of the basic flow dynamics. That is, the pressure declined gradually along the flow direction in the feeding pipeline.

Stability Analysis of the Flow Field Structure
To further observe the flow stability of the spindle-like vortex, the disturbance of the incoming flow was removed in the numerical calculation, and the outlet pressures of four branch pipelines were monitored, as shown in Figure 18. It is observed evidently that although the outlet pressures of the four branch pipelines recover, the anomalous pressure fall of pipe I/III still exists. The results indicated that the spindle-like vortex still existed, and Figure 19 also confirmed this.  Figures 16 and 17 display the distribution of the pressure on the front side and on the side of the spindle-like vortex, respectively. Meanwhile, the disturbance of 45 degrees was added to the YZX plane. When the spindle-like vortex structure was generated in the five-way spherical cavity, the pressure gradually decreased from the outside to the inside. There was a low-pressure region in the center of the vortex core of the spindle-like vortex, which was distributed along the axis of the vortex core. This was also the cause of the anomalous pressure fall in the branch pipelines on the symmetrical side. In addition to the five-way spherical cavity, the pressure fall of other parts conformed to the characteristics of the basic flow dynamics. That is, the pressure declined gradually along the flow direction in the feeding pipeline.

Stability Analysis of the Flow Field Structure
To further observe the flow stability of the spindle-like vortex, the disturbance of the incoming flow was removed in the numerical calculation, and the outlet pressures of four branch pipelines were monitored, as shown in Figure 18. It is observed evidently that although the outlet pressures of the four branch pipelines recover, the anomalous pressure fall of pipe I/III still exists. The results indicated that the spindle-like vortex still existed, and Figure 19 also confirmed this.  Figure 17. Distribution of the pressure field on the side of the spindle-like vortex.

Stability Analysis of the Flow Field Structure
To further observe the flow stability of the spindle-like vortex, the disturbance of the incoming flow was removed in the numerical calculation, and the outlet pressures of four branch pipelines were monitored, as shown in Figure 18. It is observed evidently that although the outlet pressures of the four branch pipelines recover, the anomalous pressure fall of pipe I/III still exists. The results indicated that the spindle-like vortex still existed, and Figure 19 also confirmed this.  Figure 18. Outlet pressure change of each branch pipe before and after removing disturbance. Figure 19. Flow field structure in the pipeline before and after removing disturbance.
The pressure fall remained unchanged basically after removing different disturbances, as shown in Figure 20. It is indicated that the disturbance of the initial field had little effect on the pressure fall, and once the spindle-like vortex was formed, it is hard to disappear. Nevertheless, there was a certain error between the numerical results and the experimental results. The maximum outlet pressure calculated by the numerical model was 238 kPa, and the maximum outlet pressure measured by the test was 210 kPa. We suggest that this was due to the removal of the disturbance. In practical flight, the disturbance might exist objectively. For instance, the flow from the bottom of the tank was usually swirling. Although it was not observed in the experiment, there are many related reports in the literature [43,44].  Figure 18. Outlet pressure change of each branch pipe before and after removing disturbance. Figure 19. Flow field structure in the pipeline before and after removing disturbance.

LOX
The pressure fall remained unchanged basically after removing different disturbances, as shown in Figure 20. It is indicated that the disturbance of the initial field had little effect on the pressure fall, and once the spindle-like vortex was formed, it is hard to disappear. Nevertheless, there was a certain error between the numerical results and the experimental results. The maximum outlet pressure calculated by the numerical model was 238 kPa, and the maximum outlet pressure measured by the test was 210 kPa. We suggest that this was due to the removal of the disturbance. In practical flight, the disturbance might exist objectively. For instance, the flow from the bottom of the tank was usually swirling. Although it was not observed in the experiment, there are many related reports in the literature [43,44].  Figure 19. Flow field structure in the pipeline before and after removing disturbance.

LOX
The pressure fall remained unchanged basically after removing different disturbances, as shown in Figure 20. It is indicated that the disturbance of the initial field had little effect on the pressure fall, and once the spindle-like vortex was formed, it is hard to disappear. Nevertheless, there was a certain error between the numerical results and the experimental results. The maximum outlet pressure calculated by the numerical model was 238 kPa, and the maximum outlet pressure measured by the test was 210 kPa. We suggest that this was due to the removal of the disturbance. In practical flight, the disturbance might exist objectively. For instance, the flow from the bottom of the tank was usually swirling. Although it was not observed in the experiment, there are many related reports in the literature [43,44].  Figure 20. Outlet pressure change after removing different disturbances.

Effect of the Five-Way Spherical Cavity Structure on Pressure Loss
To study the effect of five-way spherical cavity structure for the pressure loss, the numerical calculation considering the change of the diameter of the five-way spherical cavity was conducted. In the numerical calculation, the scaling ratio was defined as the ratio of the current sizes to the original sizes, which were the spherical cavity with a diameter of 360 mm and the branch pipeline with a diameter of 120 mm. The scaling ratio was varied from a value of 0.01 to 100. That is to say, the two geometric parameters were enlarged or reduced at the same time based on the original size. The angle of disturbance was given as 45 degrees on the YZX plane. The average pressure fall for four branch pipelines in different sizes is shown in Figure 21. The numerical calculation showed that when the size of the five-way spherical cavity was too small (the scaling ratio is less than 0.01), the spindle-like vortex structure in the five-way spherical cavity would disappear. When the scaling ratio was greater than 0.01, the pressure fall increased with the enhancement of the scaling ratio of the fiveway spherical cavity structure. When the scaling ratio was greater than 10, the pressure fall remained unchanged substantially. The results showed that the spindle-like vortex entered a fully developed stage, which makes the low-pressure region in the center of the spindle-like vortex core attain a stable state.

Effect of the Five-Way Spherical Cavity Structure on Pressure Loss
To study the effect of five-way spherical cavity structure for the pressure loss, the numerical calculation considering the change of the diameter of the five-way spherical cavity was conducted. In the numerical calculation, the scaling ratio was defined as the ratio of the current sizes to the original sizes, which were the spherical cavity with a diameter of 360 mm and the branch pipeline with a diameter of 120 mm. The scaling ratio was varied from a value of 0.01 to 100. That is to say, the two geometric parameters were enlarged or reduced at the same time based on the original size. The angle of disturbance was given as 45 degrees on the YZX plane. The average pressure fall for four branch pipelines in different sizes is shown in Figure 21. The numerical calculation showed that when the size of the five-way spherical cavity was too small (the scaling ratio is less than 0.01), the spindle-like vortex structure in the five-way spherical cavity would disappear. When the scaling ratio was greater than 0.01, the pressure fall increased with the enhancement of the scaling ratio of the five-way spherical cavity structure. When the scaling ratio was greater than 10, the pressure fall remained unchanged substantially. The results showed that the spindle-like vortex entered a fully developed stage, which makes the low-pressure region in the center of the spindle-like vortex core attain a stable state.  Figure 20. Outlet pressure change after removing different disturbances.

Effect of the Five-Way Spherical Cavity Structure on Pressure Loss
To study the effect of five-way spherical cavity structure for the pressure loss, the numerical calculation considering the change of the diameter of the five-way spherical cavity was conducted. In the numerical calculation, the scaling ratio was defined as the ratio of the current sizes to the original sizes, which were the spherical cavity with a diameter of 360 mm and the branch pipeline with a diameter of 120 mm. The scaling ratio was varied from a value of 0.01 to 100. That is to say, the two geometric parameters were enlarged or reduced at the same time based on the original size. The angle of disturbance was given as 45 degrees on the YZX plane. The average pressure fall for four branch pipelines in different sizes is shown in Figure 21. The numerical calculation showed that when the size of the five-way spherical cavity was too small (the scaling ratio is less than 0.01), the spindle-like vortex structure in the five-way spherical cavity would disappear. When the scaling ratio was greater than 0.01, the pressure fall increased with the enhancement of the scaling ratio of the fiveway spherical cavity structure. When the scaling ratio was greater than 10, the pressure fall remained unchanged substantially. The results showed that the spindle-like vortex entered a fully developed stage, which makes the low-pressure region in the center of the spindle-like vortex core attain a stable state.

Conclusions
For the anomalous pressure fall phenomenon in the feeding system of four parallel engines of a cryogenic vehicle, a 3D numerical model considering a full-scale structure was developed to accurately capture the formation and evolution of vortices in the multi-branch pipeline with a five-way spherical cavity. Based on the numerical results presented in this paper, the following conclusions could be drawn: (1) The numerical results developed in the present paper were consistent with the experiment results. It was indicated that the established model describes the instability flow features of liquid oxygen reliably and robustly, including geometry simplification, initial condition setting, boundary condition setting, solution strategy, mesh generation, turbulence model selection, and near wall treatment.
(2) It was found that the primary reason of the anomalous pressure was that the spindle-like vortex structure formed due to disturbance of the inlet boundary condition, and there were minimum pressure and velocity values at the vortex core.
(3) With the enhancement of the disturbance, the flow field structure in the five-way spherical cavity evolved from the mirror-symmetrical four-vortex structure to the spindle-like vortex. There was a mirror-symmetrical four-vortex structure in the absence of disturbance, and the outlet pressure of the four branch pipelines was equal. The pressure loss in the spherical cavity was about 19 kPa. When disturbance was added at the entrance of the main pipeline, the mirror symmetrical four vortices gradually merged into the mirror symmetrical two vortices until the disturbance was enhanced to a certain extent, and the mirror-symmetric two-vortex structure was rapidly transformed into the spindle-like vortex. The pressure loss in the feedline on the side corresponding to the vortex core was about 90 kPa when there was about a 45-degree disturbance, and on the other side was about 20 kPa. The numerical results were in agreement with the experimental data.
(4) The change of the flow field was an irreversible behavior. The spindle-like vortex flow field was more stable than the symmetric vortex flow field, and could not change with the disturbance once it was formed. When removing the disturbance, the spindle-like vortex structure still remained stable. Although the outlet pressure of four branch pipeline picked up, the pressure loss on the side corresponding to the spindle-like vortex was still maximum.