**Appendix B. Staggered-Grid Finite-Volume Approach**

As described by Torregrosa et al. [34], the quasi-3D model uses a staggered grid with two different basic elements: volumes and connectors. Volumes have associated scalar information such as pressure, density or temperature, as well as the cell volume. Connectors contain vector information, like flow velocity, momentum or the orientation of the connector in space, together with its area. With this configuration, a volume might have attached as many connectors as needed, but a connector will always connect only two volumes. In Figure A2 two volumes and a connector are shown schematically (volumes and connectors actually have an undefined shape).

**Figure A2.** Basic mesh elements, definition of velocity projections and notation of volumes and connectors.

The method is based on the solution of the 3D Euler conservation equations without source term:

$$
\partial \rho / \partial t + \nabla \cdot (\rho \mathbf{u}) = 0 \tag{A4}
$$

$$
\partial(\rho \mathbf{u})/\partial t + \nabla \cdot (\rho \mathbf{u} \mathbf{u}) = -\nabla p \tag{A5}
$$

$$\left[\partial(\rho\epsilon\mathbf{o})/\partial\mathbf{t} + \nabla \cdot [(\rho\epsilon\mathbf{o}) + p)\mathbf{u}\right] = 0\tag{A6}$$

together with the perfect gas equation of state. Here, *ρ* is the density, *p* is the pressure, **u** is the velocity vector, ∇· indicates the divergence, and *e*0 is the specific stagnation internal energy, whose expression for a perfect gas is:

$$\mathfrak{e}\_0 = \mathfrak{c}\_\upsilon T + \mathfrak{u}^2/2 \tag{A7}$$

where *cv* is the specific heat capacity at constant volume and *T* is the fluid static temperature.

However, in the context of a finite volume method on a staggered grid, the key issue is how and where these equations are solved. The mass conservation equation, being scalar, is solved in the volumes. Upon discretization it becomes:

$$
\rho^{n+1} = \rho^n + \frac{\Delta t}{V} \sum\_{c=1}^{N\_c} \rho\_c^n u\_c^n A\_c \tag{A8}
$$

where *u* is the flow velocity. Superscript *n* indicates the time step, Δ*t* is the time interval, *V* the volume of the cell, *Nc* the number of connectors, *A* is the cross section, and subscript *c* indicates that the variable is taken at a connector (otherwise it is assumed that the variable is taken at the volumes).

Following a similar procedure, the discretized energy equation is written as:

$$(\rho \varepsilon\_0)^{n+1} = (\rho \varepsilon\_0)^n + \frac{\Delta t}{V} \sum\_{c=1}^{N\_c} \rho\_c^n e\_0^n u\_c^n A\_c + \frac{\Delta t}{V} \sum\_{c=1}^{N\_c} p\_c^n u\_c^n A\_c \tag{A9}$$

The momentum equation in a 3D case consists of three coupled equations. In this context, the momentum equation is solved in the connectors and only in the direction orthogonal to its surface. This is achieved by projecting the flow velocity of the two adjacent volumes into that direction, as shown in Figure A2, where the velocity *uc* in the connector, and the projections of the volume flow velocity from its left and right, *uLn* and *uRn*, are shown. As a result, the momentum in the connector follows a one-dimensional momentum equation, whose discretization gives:

$$(\rho\_{\varepsilon}u\_{\varepsilon}A\_{\varepsilon})^{n+1} = (\rho\_{\varepsilon}u\_{\varepsilon}A\_{\varepsilon})^{n} + (\Delta t/\Delta L) \left[ \left(\rho u\_{n}^{2} + p\right)\_{L} + \left(\rho u\_{n}^{2} + p\right)\_{R} \right] A\_{\varepsilon} \tag{A10}$$

Here, *un* is the velocity projection onto the direction orthogonal to the connector surface and subscripts *R* and *L* denote the variables taken from the volumes at the right and left of the connector, respectively. With this simplification, a single equation is solved for each connector instead of three coupled equation for each volume. While usually the number of connectors is higher than the number of volumes in a 3D mesh, the fact that the momentum equations for each direction are not coupled reduces drastically the computation time when compared to a regular 3D method.

Once the momentum is calculated, its value is used to compute the mass and energy conservation equations in the next time step. It is worth noticing that some terms in the momentum and energy equations, like density or pressure, must be evaluated in the connectors, but are only calculated in the volumes. To solve this, an upwind approach is adopted, the required values being taken from the right or left volumes, depending on the flow direction.

Finally, the momentum in the volumes is computed by distributing the momentum calculated in the connectors among the volumes they connect, taking into account their size. In a uniform 1D mesh, half the momentum is thus assigned to each volume. As the orientation of the connectors in space is known, the resulting momentum vector of a volume can be calculated from the vector sum of the momentum in the connectors:

$$\left(\rho\_{\varepsilon}\mathbf{u}V\right)\_{\upsilon}^{n+1} = \frac{1}{2}\sum\_{c=1}^{N\_{\varepsilon}} \left(\rho\mathbf{u}\_{\varepsilon}A\_{c}\Delta L\right)\_{\upsilon}^{n+1} \tag{A11}$$

The resulting method is a second-order accuracy method based on an explicit scheme with a staggered grid, as shown in Figure A3.

**Figure A3.** Scheme of the staggered mesh and the associated time marching.

The fact that the resulting scheme offers second-order accuracy, together with the simplifications adopted in the momentum equation, results in nonphysical oscillations, especially in the vicinity of significant pressure gradients. In order to mitigate such overshoots, two different flux limiters have been used here: the momentum diffusion term (MDT) proposed in [35] and the flux-corrected transport (FCT) methodology proposed in [34].

In the case of the MDT, the main goal is to add a diffusion term to the momentum equation so that the mass flux computed at that connector is conveniently limited. With this purpose, the momentum

flux density tensor used in the momentum Equation (A5) can be modified, in a way similar to that used for incorporating viscosity effects, as follows:

$$
\partial(\rho \mathbf{u})/\partial t + \nabla \cdot (\rho \mathbf{u} \mathbf{u} + \mathbf{D}) = -\nabla p \tag{A12}
$$

where the tensor **D** is assumed to depend linearly on the local momentum gradients, i.e.,

$$\mathbf{D} = \epsilon \nabla(\rho \mathbf{u}) \tag{A13}$$

The scalar quantity  has the dimensions of a kinematic viscosity and can thus be interpreted as a momentum diffusion coefficient. With this prescription, the contribution of the diffusion term ∇·**D** will only be relevant if significant gradients exist, and any resulting spurious oscillations will be damped.

Adding the projection of Equation (A12) onto the direction of a connector to the discretized momentum equation one gets:

$$\left(\rho\_{\varepsilon}u\_{\varepsilon}A\_{\varepsilon}\right)^{n+1} = \left(\rho\_{\varepsilon}u\_{\varepsilon}A\_{\varepsilon}\right)^{n} + \left(\frac{\Delta t}{\Delta L}\right)\left\{\left[\left(\rho u\_{n}^{2} + p\right)\_{L} + \left(\rho u\_{n}^{2} + p\right)\_{R}\right]A\_{\varepsilon} + \left(\vec{D}\_{\rm Ln} - \vec{D}\_{\rm Rn}\right)\right\} \tag{A14}$$

where *D L* and *D R* are the projections onto the connector direction of tensor **D** = ∇(*ρ***u***<sup>A</sup>*) computed in the two adjacent volumes. Following [35], the momentum diffusion coefficient is computed by relating the mesh size and the time step to the local flow velocity at the volume, as:

$$
\epsilon = \frac{|\mathbf{u}|}{2} (\Delta L - |\mathbf{u}| \Delta t) \tag{A15}
$$

and the gradient of mass flow ∇(*ρ***u***<sup>A</sup>*) is also computed from the projections of the mass flows of the adjacent connectors onto each direction.

Considering now the FCT technique, it was initially conceived to be used in a finite differences scheme, but it was adapted to a finite volume staggered grid in [34]. In the FCT technique three stages can be identified: a transport stage based on the scheme considered, a diffusion stage where the numerical dispersion is reduced, and an anti-diffusion stage in which the accuracy of the scheme where the solution is smooth is restored. The diffusion operator is defined as:

$$D\_j(\mathcal{W}) = \theta\left(\mathcal{W}\_{\dot{j}+1/2}\right) - \theta\left(\mathcal{W}\_{\dot{j}-1/2}\right) \tag{A16}$$

with,

$$
\theta\left(W\_{j+1/2}\right) = \theta / 4\left(W\_{j+1} - W\_j\right) \tag{A17}
$$

where *Wj* is the conservative variable computed at cell *j* in the transport stage, *j* ± 1/2 indicates that the conservative variable is evaluated in the midpoint between *j* and *j* ± 1, and the factor *ϑ* is a positive real number with *ϑ* ≥ 1/2. Calculating the diffusion via damping, as suggested in [34], the new variable *Wj* is obtained as:

$$
\overline{\mathcal{W}}\_{j}^{n+1} = \mathcal{W}\_{j}^{n+1} + D\_{j}(\mathcal{W}^{n}) \tag{A18}
$$

When considering a staggered mesh finite volume method, the FCT is only applied to the momentum equation, since this is the main source of oscillations. Therefore, all the calculations in the volumes are made using *Wnj* = (*ρucAc*)*nj* as the conservative variable, whereas the calculations in the midpoints make use of the variables evaluated at the volumes adjacent to the connector.

Finally, the non-linear anti-diffusion operator *Aj* is defined as:

$$A\_j(\mathcal{W}) = \Psi\left(\mathcal{W}\_{j-1/2}\right) - \Psi\left(\mathcal{W}\_{j+1/2}\right) \tag{A19}$$

Making use of the anti-diffusive limited flow defined in [37] gives:

$$\Psi\left(\mathcal{W}\_{\hat{\jmath}+1/2}\right) = \text{smax}\left[0, \min\left((5/8)s\,\Delta\mathcal{W}\_{\hat{\jmath}-1/2}, (1/8)\left|\Delta\mathcal{W}\_{\hat{\jmath}+1/2}\right|, (5/8)s\,\Delta\mathcal{W}\_{\hat{\jmath}+3/2}\right)\right] \tag{A20}$$

where s = sign -Δ *Wj*−1/2 , Δ *Wj*−1/2 = *Wj* − *Wj*−1, Δ *Wj*+1/2 = *Wj*+<sup>1</sup> − *Wj* and Δ *Wj*+3/2 = *Wj*+<sup>2</sup> − *Wj*+1. Then, according to [34], the phoenical form should be used, so that:

$$
\tilde{\mathcal{W}}\_{\dot{j}}^{n+1} = \overline{\mathcal{W}}\_{\dot{j}}^{n+1} + A\_{\dot{j}} \left( \mathcal{W}^{n+1} \right) \tag{A21}
$$

The anti-diffusion stage equations can be adapted to the staggered grid mesh finite volume method in a way similar to that used with the diffusion stage.

### **Appendix C. 1D Method with Pressure Loss-Based Junction Model**

In this case, a collocated one-dimensional finite volume method is used for all the calculations inside the pipes. The Euler equations of fluid dynamics simplified for one-dimensional flow in a straight uniform duct can be expressed as:

$$\frac{d\overline{\mathbf{W}}\_{i}}{dt} = \frac{d}{dt} \begin{pmatrix} \rho \\ \rho u \\ \rho \varepsilon\_{0} \end{pmatrix}\_{i} = \frac{A(F\_{i-1,i} - F\_{i,i+1})}{V\_{i}} \tag{A22}$$

Here, *Wi* is the cell-averaged state vector of cell *i*, *<sup>F</sup>i*−1,*<sup>i</sup>* and *<sup>F</sup>i*,*i*+<sup>1</sup> are the inter-cell fluxes between cells *i* − 1 and *i* and between *i* and *i* + 1, respectively, and the other symbols refer to the same magnitudes as in Appendix B, with *u* being now the axial velocity of the flow. The inter-cell fluxes are computed by an approximate solution of the Riemann problem as described by Toro et al. [38]. The state vector is extrapolated to the cell boundaries to compute the fluxes by means of a Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) approach as described in [39], while the solution is propagated in time using Heun's method, leading to a second order in time and space, total variation diminishing scheme.

While the main one-dimensional flow inside the ducts is simulated, the effects of the geometry of the junction are modelled. The connections of the ducts to the junction are solved using an auxiliary small zero-dimensional element. Each of the one-dimensional branches is connected to that zero-dimensional element making use of the Riemann variables to compute the fluxes at their corresponding boundary. At each connection, it is assumed that a certain amount of stagnation pressure is lost, depending on the angle of the junction and of the ratio between the outflow mass flow . *mout* passing through the branch of interest and the inflow mass flow . *min*. The pressure loss coefficients are computed following the expressions given in [15,16], and are defined as the ratio of the difference in stagnation pressure between the outflow branch and the inflow branch to the dynamic pressure (*ρu*2/2) of the inflow branch. Finally, the total pressure loss coefficient *K* for a three-branch junction with the same section in all the branches can be estimated as:

$$K = \left(\frac{\dot{m}\_{out}}{\dot{m}\_{in}}\right)^2 - \frac{3}{2}\frac{\dot{m}\_{out}}{\dot{m}\_{in}} + \frac{1}{2} \tag{A23}$$

when the branch is collinear with the inflow branch, and,

$$K = \left(\frac{\dot{m}\_{out}}{\dot{m}\_{in}}\right)^2 + 1 - 2\frac{\dot{m}\_{out}}{\dot{m}\_{in}}\cos\left(\frac{3}{4}\theta\right) \tag{A24}$$

for the lateral branch, when the flow is split between a collinear and a lateral branch. In this case, *θ* is the angle between the lateral branch and the axial outflow branch, so that 0◦ means that both outflow branches are parallel. The same expression applies when the inflow branch is not parallel to any of the outflow branches: in that case, the angle *θ* is measured between the inflow branch and the other outflow branch.

In the auxiliary zero-dimensional element, the gas is again considered as a perfect gas, and the mass and energy equations are solved:

$$\frac{\text{dm}}{\text{dt}} = \sum \dot{m} \tag{A25}$$

$$\frac{\mathbf{d}(mc\_{\upsilon}T)}{\mathbf{d}t} = \sum \dot{m}h\_0 \tag{A26}$$

where *m* is the mass trapped in the zero-dimensional element, .*m* is the mass flow, positive when it enters the element, and *h*0 is the specific stagnation enthalpy associated with the mass moving inside or outside of the element. These two equations set an additional limitation to the maximum possible time step.


© 2017 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
