Embedded One-Dimensional Orifice Elements for Slosh Load Calculations in Volume-Of-Fluid CFD

Abstract

For CFD liquid sloshing simulations, fine computational mesh resolutions are typically required to model the flow within small flow passages or orifices found in fuel tanks. This work presents a method of replacing the fine computational mesh elements within orifices with large one-dimensional mesh elements that integrate seamlessly with standard finite volume computational elements with the intended advantage of reducing the overall computational cost of CFD simulations. These one-dimensional elements conserve mass and momentum for two-phase flow in incompressible Volume-Of-Fluid CFD. Instead of fully resolving the momentum diffusion term, empirical correlations are used to account for the viscous losses within the orifices for both two- and three-dimensional simulations. The one-dimensional orifice elements are developed and validated against analytical and experimental results using the finite volume CFD code Elemental^®. Furthermore, these elements are tested in a violent sloshing simulation and compared with full-resolution numerical results as well as experimental results. The elements are shown to decrease computational cost significantly by reducing the number of computational elements as well as increasing the simulation time step sizes (due to an increase in element sizes).

1. Introduction

Slosh is a term that often refers to a seemingly chaotic flow of two fluids, usually a liquid and a gas, within a containment structure such as a tank. Common occurrences of slosh are within the fuel tanks of passenger vehicles [1], commercial aircraft [2,3,4], spacecraft [5,6] or in tanker vessels transporting liquefied natural gas [7,8]. The resulting slosh forces are an important consideration in the structural design of tanks as well as the dynamics and control of the vehicles carrying them [4,9].

Due to the inherent non-linearities involved, numerical methods are widely used to study slosh loads. These include Smoothed Particle Hydrodynamics (SPH) [10], Finite Difference Method (FDM) [11], Finite Element Method (FEM) [12] and Boundary Element Method (BEM) [13]. In this work, we focus on the Finite Volume method [14] with a Volume-Of-Fluid (VOF) [15] approach to model two-phase, immiscible fluid flow in the fuel tank slosh loads problem [3].

Baffle structures are often found inside tanks and serve to suppress the fuel slosh loads [16], to act as structural elements [2] or both. Openings in these structures, referred to generally as orifices in this work, can significantly increase computational cost in numerical simulations due to the small geometrical scales [17,18,19,20,21,22]. In addition to small mesh elements, orifices are often the location of high-velocity flow, as fluid is accelerated from one side of the tank to another through a constricting passage. Therefore, the computational cost is increased by both the requirement for more mesh elements and more time steps, since the allowable time step size is decreased in orifices by stability constraints such as the Courant–Friedrichs–Lewy (CFL) condition [23] or viscous diffusion in the case of explicit time integration schemes [24].

An approach to overcome this cost is to model orifice or channel flow instead of directly discretizing and solving the governing Navier–Stokes equations. Ozhan [25] developed a subgrid pressure jump model for use in the development of automotive catalytic converters, replacing small monolithic substrates with a subgrid model coupled to the full solution of the Navier–Stokes equations in the diffuser and outlet nozzle. Porter [26] added to the work of Ozhan by resolving the full Navier–Stokes equations within the inlet section of the monolithic substrate channels while using the pressure loss model for the downstream section (where flow is fully developed). This allows the numerical model to capture the additional losses due to developing flow in the monolithic substrate channels.

Another approach is the coupling of 3D CFD meshes with 1D network codes [27], in which the 1D code is used to model flow in long tubes while the 3D CFD code models the remainder of the flow domain. These methods require a coupling mechanism to communicate between the network code and 3D CFD code.

In this work, we further the methodology of Ozhan [25] and Porter [26] for constant cross-sectional orifices. In our method, we propose to embed the orifice element as part of the mesh. This allows the user to simply modify the mesh in orifice areas by drastically reducing the amount of mesh elements in these regions. A single equation is solved for momentum conservation in the entire domain with different treatment of certain terms in the bulk flow and 1D flow within the orifice. We further expand the 1D model to account for two-phase loss effects by applying existing semi-empirical correlations for flow losses. Another novel aspect of this work is giving careful attention to the VOF advection method to ensure sharp tracking of the liquid–gas interface within the 1D orifice elements. The development and validation of the above is performed in the multi-physics CFD code Elemental^® [19,28,29,30], which is developed by the Industrial CFD (InCFD) Research Group at the University of Cape Town.

2. Materials and Methods

In this section, we will briefly introduce orifice flow loss physics (Section 2.1) and the semi-empirical correlations we use to model this. We will then present the governing equations of our problem (Section 2.2) and the numerical method we use to solve it (Section 2.3).

2.1. Orifice Flow Losses

Figure 1 shows a schematic representation of the contracting and expanding flow within a long orifice, defined as orifices in which the length L to hydraulic diameter $D_H$ ratio is greater than 1.4, i.e., $\frac{L}{D_H}>1.4$. The hydraulic diameter is defined as

$$D_H=\frac{4A_o}{P_w},$$

(1)

where $A_o$ is the orifice cross-sectional area and $P_w$ is the wetted perimeter. The point $x_u$ is at the inlet to the orifice and $x_d$ is at the exit to the orifice where the flow is reattached to the orifice wall, as shown. The vena contracta is the location at which the main fluid stream cross-sectional area is smallest, which is at position $x_c$. Flow through orifices up to the vena contracta is nearly reversible [31]. Flow losses mainly occur downstream of the vena contracta with the majority being downstream of the orifice exit in the downstream jet. In the case of short orifices ($\frac{L}{D_H}<1.4$), the losses within the orifice are negligible [32,33]. For long orifices ($\frac{L}{D_H}>1.4$), the flow reattaches to the orifice wall downstream of the vena contracta and, even though the majority of losses still occur downstream of the orifice, significant flow losses occur within the orifice [34].

This work focuses on modeling the viscous flow losses and corresponding orifice mass flow rate within long orifices for two-phase flow by using one-dimensional, finite volume cells in orifice cavities. The flow outside the orifice is meshed, discretized and solved with the usual numerical discretization. Therefore, we place the onus to account for the losses occurring downstream of the orifice on the existing 2D/3D CFD model, as would be the case in a fully resolved simulation.

The average velocity and pressure over the orifice cross-section (normal to the flow direction) are denoted $\tilde{u}$ and $\tilde{p}$, respectively. To these, we add subscripts that refer to the point where the cross-sectional average velocity or pressure is taken (upstream, downstream or at the vena contracta). Note that the orifices considered in this work all have constant cross-sectional areas; therefore, we simply refer to the orifice cross-sectional area, $A_o=A_u=A_d$. $A_c$ is the cross-sectional area of the vena contracta.

The pressure loss in one-dimensional flow analyses are generally written in the form [35]

$$\Delta \tilde{p}=C\rho {\tilde{u}}^2,$$

(2)

where $\rho$ is the density of the fluid and C is some loss coefficient, determined as a function of the flow conditions and boundary geometry. The effects that contribute to the vast majority of pressure losses through orifices can be placed in three categories. These are viscous losses, flow losses (venturi such as contraction and expansion) and fluid mixing (two-phase flow). These effects are discussed in the following subsections, with a general equation for the loss coefficient given in Section 2.1.4. Note that the entrance losses due to developing flow (as in the case of pipes and ducts) in long orifices can be considered negligible due to the relatively small length-to-diameter ratios [36] and are excluded in this work.

2.1.1. Viscous Losses

The pressure loss along the length L of a long orifice due to viscous effects (including the wall friction) is written in a generic form, assuming isothermal flow, as

$$\Delta p_{visc}=\frac{K_{f}L}{D_H}\frac{\tilde{\rho }{\tilde{u}}^2}{2},$$

(3)

where $K_f$ is the dimensionless Darcy–Weisbach Friction Factor, L is the length of the orifice section and $D_H$ is the hydraulic diameter [37]. The empirical equation by Avci [38] will be used in this work to determine $K_f$ since it is valid over all Reynolds numbers up to ${10}^8$ and is suitable “for the entire range of the Moody chart, that is valid for laminar, transitional and fully turbulent regions” [38]. Furthermore, a smooth surface will be assumed, as is the case in the problems modeled in this work. Avci has a general empirical equation for problems with non-zero surface roughness, which can be implemented and used similarly as the smooth relation used here. The empirical equation by Avci [38] for smooth surfaces is given by

$$K_f=\frac{6.4}{{\left(ln\left(\frac{1}{Re}\right)\right)}^{2.4}}+\left(\frac{64}{Re}-\frac{6.4}{{\left(ln\left(\frac{1}{Re}\right)\right)}^{2.4}}\right)e^{-{\left(\frac{Re}{2560}\right)}^8}.$$

(4)

2.1.2. Flow Contraction and Expansion Losses

The pressure differential due to venturi such as contracting and expanding flow within the orifice is given as [34]

$$\Delta p_{cont}+\Delta p_{exp}=K_{cont}\frac{{\tilde{\rho }}_u{\tilde{u}}^2}{2}+K_{exp}\frac{{\tilde{\rho }}_d{\tilde{u}}^2}{2},$$

(5)

with

$$K_{cont}=\left(\frac{1}{{\sigma }_c^2}-1\right),$$

(6)

$$K_{exp}=\left(2-\frac{2}{{\sigma }_c}\right),$$

(7)

and

$${\sigma }_c=\frac{A_c}{A_o}\approx \frac{\pi }{\pi +2}=0.611,$$

(8)

where ${\tilde{\rho }}_u$ and ${\tilde{\rho }}_d$ are the average fluid density within the orifice—upstream and downstream, respectively—of the vena contracta and ${\sigma }_c$ is the contraction coefficient, which in this work is approximated using Milne-Thompson’s method [39] and further motivated by [35,40,41]. Note that, since the orifices considered in this work are all of constant cross-sectional area, the average velocity $\tilde{u}={\tilde{u}}_u={\tilde{u}}_d$. Equation (6) is derived from the energy flux equation assuming steady state and constant enthalpy [31], while Equation (7) is derived from the steady state momentum flux equation [34,35]. Using the approximate value of ${\sigma }_c$ in Equation (8), $K_{cont}=1.679$ and $K_{exp}=-1.273$. The relatively larger positive value of $K_{cont}$ shows how the fluid in the orifice is accelerated (and, therefore, the pressure reduced) without any losses up to the vena contracta. Similarly, the relatively smaller negative value of $K_{exp}$ shows how only some of the pressure is recovered during expansion.

2.1.3. Two-Phase Flow Losses

For two-phase flow, additional losses occur through the orifice. We use the separated flow model [35,42], which takes the total pressure drop over a pipe section or orifice for a single-phase flow and adds a two-phase multiplier ${\varphi }^2$ (determined experimentally). This multiplier accounts for mixing effects and is calculated at every time step according to the flow conditions.

A multiplier coefficient is calculated for viscous losses as well as for contracting and expanding losses. The Baroczy–Chisolm method [35] is used to determine the two-phase multiplier for viscous losses, and the method of Kojasoy [34] (based on the Baroczy–Chisolm method) is used to determine the two-phase losses due to contracting and expanding flow.

For the Baroczy–Chisolm method [35], the dryness fraction is defined as

$$\mathrm{x}=\frac{{\rho }_g{\tilde{u}}_g{A}_g}{{\rho }_l{\tilde{u}}_l{A}_l+{\rho }_g{\tilde{u}}_g{A}_g},$$

(9)

indicating the ratio of the mass flow rate of the gas to the liquid–gas mixture, where the subscripts l and g are used to distinguish between the liquid and gas phase, respectively. Further, $A_l$ and $A_g$ are the cross-sectional areas of the orifice occupied by the liquid and gas, respectively. The void fraction $\xi$ is the ratio between the volume of the orifice occupied by the gas $V_g$ to the total orifice volume $V_{lg}$:

$$\xi =\frac{V_g}{V_{lg}}.$$

(10)

In general, we employ the subscript $lg$ to indicate the two-phase mixture. The mass velocity G is used extensively in one-dimensional two-phase flow literature and is given as

$$G={\rho }_l{\tilde{u}}_l\frac{1-\xi }{1-\mathrm{x}}={\rho }_g{\tilde{u}}_g\frac{\xi }{\mathrm{x}}.$$

(11)

Finally, the physical property coefficient is given as

$${\Gamma }^2=\frac{\Delta p_{g0}}{\Delta p_{l0}},$$

(12)

where subscripts $g0$ and $l0$ refer to the gas or liquid phase, respectively, if it were to flow through the orifice at the mixture’s mass flow rate.

The total pressure differential due to viscous effects of the two-phase mixture $\Delta p_{lg}$ is given as

$$\Delta p_{lg}={\varphi }_{l0}^2\Delta p_{l0}.$$

(13)

The two-phase multiplier, using the Baroczy–Chisolm method, is defined as

$${\varphi }_{l0}^2=1+({\Gamma }^2-1)(B\mathrm{x}(1-\mathrm{x})+{\mathrm{x}}^2),$$

(14)

where the coefficient B can be found in

Table 1. Values of B in the Baroczy–Chisolm method.

$\mathbf{\Gamma }$	$\mathit{G}\left(\frac{\mathit{kg}}{{\mathit{m}}^2\mathit{s}}\right)$	B
$\Gamma \le 9.5$	$G\le 500$	$4.8$
	$500<G<1900$	$\frac{2400}{G}$
	$G\ge 1900$	$\frac{55}{G^{0.5}}$
$9.5<\Gamma <28$	$G\le 600$	$\frac{520}{\Gamma G^{0.5}}$
	$G>600$	$\frac{21}{\Gamma }$
$\Gamma \ge 28$		$\frac{15000}{{\Gamma }^2{G}^{0.5}}$

[35] and the physical property coefficient $\Gamma$ can be determined using Equation (12). $\Delta p_{l0}$ and $\Delta p_{g0}$ are found using Equation (3) for the liquid

$$\Delta p_{l0}=\frac{K_{fl0}L}{D_H}\frac{G^2}{2{\rho }_l}$$

(15)

and gas phase

$$\Delta p_{g0}=\frac{K_{fg0}L}{D_H}\frac{G^2}{2{\rho }_g},$$

(16)

with the mass velocity G calculated using (11). $K_{fg0}$ and $K_{fl0}$ are found using Equation (4). Note that the Reynolds number in Equation (4) is, for two-phase flow problems, determined as

$$Re={\displaystyle \frac{G{D}_H}{{\tilde{\mu }}_u}},$$

(17)

where ${\tilde{\mu }}_u$ is the viscosity of the respective phase being dealt with. Since the two-phase effects are accounted for separately, the hydraulic diameter $D_H$ is calculated using the full perimeter of the orifice.

Finally, assuming that the liquid and gas velocities are the same, the total two-phase viscous pressure differential is calculated, using Equations (13)–(15), as

$$\Delta p_{visc}={\varphi }_{l0visc}^2\Delta p_{l0visc}={\varphi }_{l0visc}^2\frac{K_{fl0}L}{D_H}\frac{G^2}{2{\rho }_l}.$$

(18)

Chisolm [35] and Kojasoy [34] developed methods to expand the application of two-phase multipliers to obstructions such as orifices in pipes and ducts. These methods are used only to determine the two-phase multiplier for the pressure difference due to the contraction and expansion caused by orifices. The two-phase multiplier due to the effects of an obstruction such as an orifice is given as [35]

$${\varphi }_{orif}^2=\left(1+\left(\frac{{\rho }_l}{{\rho }_g}-1\right)(B_{orif}\mathrm{x}(1-\mathrm{x})+{\mathrm{x}}^2)\right),$$

(19)

where the dryness fraction $\mathrm{x}$ is as defined previously;

$$B_{orif}=\frac{1}{S^n},$$

(20)

with n determined experimentally; and the velocity ratio or slip ratio S is defined as

$$S={\left(1+\mathrm{x}\left(\frac{{\rho }_l}{{\rho }_g}-1\right)\right)}^{\frac{1}{2}}.$$

(21)

Kojasoy [34] approximated $n=0$ downstream of the vena contracta since the expanding flow is well mixed and $n=0.15$ upstream of the vena contracta for long orifices. This approximation is further motivated by [34,43].

Finally, in this work, the total two-phase pressure differential due to contracting and expanding flow, respectively, are calculated using Equations (13) and (19)–(21), as

$$\Delta p_{cont}={\varphi }_{l0cont}^2\Delta p_{l0cont}=\left(1+\left(\frac{{\rho }_l}{{\rho }_g}-1\right)\left(\frac{1}{S^{0.15}}\mathrm{x}(1-\mathrm{x})+{\mathrm{x}}^2\right)\right)K_{cont}\frac{G^2}{2{\rho }_l}$$

(22)

and

$$\Delta p_{exp}={\varphi }_{l0exp}^2\Delta p_{l0exp}=\left(1+\left(\frac{{\rho }_l}{{\rho }_g}-1\right)(\mathrm{x}(1-\mathrm{x})+{\mathrm{x}}^2)\right)K_{exp}\frac{G^2}{2{\rho }_l},$$

(23)

where $K_{cont}$ and $K_{exp}$ are defined in Equations (6) and (7).

2.1.4. Total Two-Phase Losses within Orifices

The total pressure differential over a long orifice accounting for two-phase flow in this work is, therefore, modeled by adding loss terms in Equations (18), (22) and (23) such that

$$\begin{array}{cc}\hfill \Delta p_{total}& =\Delta p_{visc}+\Delta p_{cont}+\Delta p_{exp}\hfill \\ \hfill & ={\varphi }_{l0visc}^2\frac{K_{fl0}L}{D_H}\frac{G^2}{2{\rho }_l}+{\varphi }_{l0cont}^2{K}_{cont}\frac{G^2}{2{\rho }_l}+{\varphi }_{l0exp}^2{K}_{exp}\frac{G^2}{2{\rho }_l}\hfill \\ \hfill & =C\frac{G^2}{{\rho }_l},\hfill \end{array}$$

(24)

where C represents a global loss coefficient

$$C=\frac{{\varphi }_{l0visc}^2\frac{K_{fl0}L}{D_H}+{\varphi }_{l0cont}^2{K}_{cont}+{\varphi }_{l0exp}^2{K}_{exp}}{2}.$$

(25)

The form of Equation (25) allows users to input alternative models for the loss coefficient, if required.

2.2. Fluid Governing Equations

2.2.1. Conservation Equations in a Non-Inertial Reference Frame

Consider a fluid volume $V\left(t\right)$ shown in Figure 2. $V\left(t\right)$ is enclosed by a surface $A\left(t\right)$ with outward pointing unit normal vector $\mathit{n}$. The volume is located in a non-inertial reference frame (x-y-z); contains a fluid with density $\rho$ and viscosity $\mu$; and is subject to a force, which is divided into a surface force ${\mathit{F}|}_A$ and a body force ${\mathit{F}|}_V$. The non-inertial reference frame (x-y-z) is free to move relative to some inertial reference frame (X-Y-Z). In this work, we will only consider the linear translation of the non-inertial reference frame.

The general momentum conservation equation for a fluid in the non-inertial reference frame system described above can be written as

$$\begin{array}{cc}\hfill \frac{\partial \rho \mathit{u}}{\partial t}& +\nabla ·(\rho \mathit{u}\otimes \mathit{u})+\nabla p-\nabla ·\mu (\nabla \mathit{u}+\nabla {\mathit{u}}^T)\hfill \\ \hfill & =\rho \left(\mathit{g}-{\mathit{a}}_B\right),\hfill \end{array}$$

(26)

where p is the static pressure, $\mathit{u}$ is the velocity of the fluid relative to the non-inertial reference frame (x-y-z) and $\mathit{g}$ is the body force acceleration due to gravity expressed as $-9.81$ m/s${}^2$ in the positive Y-direction. ${\mathit{a}}_B$ is the linear acceleration of the non-inertial reference frame relative to X-Y-Z. The conservation of mass for an incompressible fluid is given by

$$\nabla ·\mathit{u}=0.$$

(27)

The mass and momentum conservation equations are known as the incompressible Navier–Stokes equations. For CFD simulations in this work, the tank is fixed to the non-inertial reference frame x-y-z, which undergoes a time-varying acceleration ${\mathit{a}}_B$.

2.2.2. Orifice Flow Loss Model

We now seek to expand the momentum conservation Equation (26) to describe both the bulk flow as well as orifice flow losses as per Equation (24). Consider an orifice with diameter $D_H$ and length L, as shown in Figure 3. The orifice is fixed in the x-y-z reference frame, which is subjected to an acceleration ${\mathit{a}}_B$ while the pressures at the inflow and outflow coordinates ${\mathit{x}}_1$ and ${\mathit{x}}_2$ are denoted by $p_1$ and $p_2$, respectively. In this case, $\Delta p_{total}=p_2-p_1$. Unit vector ${\mathit{n}}_O$ indicates the 1D direction vector of the flow within the orifice in the x-y-z reference frame. The viscous term in the above momentum equation is now expanded to include the orifice loss model as follows:

$$\nabla ·\mu (\nabla \mathit{u}+\nabla {\mathit{u}}^T)=\left(\nabla ·\mu (\nabla \mathit{u}+\nabla {\mathit{u}}^T)\right)\left(1-H_O\right)-\frac{C\rho {\tilde{u}}_O^2}{L}{H}_O{\mathit{n}}_O,$$

(28)

where ${\mathit{n}}_O$ is the unit normal vector along the orifice in the direction of the flow and $H_O$ is a Heaviside function with value 1 within the orifice and 0 elsewhere.

Combining Equations (26) and (28) gives the momentum conservation equation for two-phase incompressible flow in both bulk and orifice regions as follows:

$$\begin{array}{cc}\hfill \frac{\partial \rho \mathit{u}}{\partial t}+\nabla ·(\rho \mathit{u}\otimes \mathit{u})+\nabla p-\left(\nabla ·\mu (\nabla \mathit{u}+\nabla {\mathit{u}}^T)\right)(1-H_O)& \hfill +\frac{C\rho {\tilde{u}}_O^2}{L}{H}_O{\mathit{n}}_O\hfill \\ \hfill & \hfill \approx \rho \mathit{g}-\rho {\mathit{a}}_B,\hfill \end{array}$$

(29)

where the nomenclature is as previously defined.

2.3. Numerical Method

The numerical method applied to solve the fluid governing Equations (27) and (29) is described next.

2.3.1. Brief Overview of Vertex-Centered VOF Method

Elemental^® uses a vertex-centered finite volume [44] fractional-step scheme [19] such that the control volumes are formed as a geometric dual to the cells. In the case of two-phase flow, the fluids are assumed immiscible. A one-fluid formulation [45] is used so that the same set of equations is solved everywhere in the computational domain, regardless of the phase or phases present in the computational cell. To track the phases, the VOF equation is based on

$$\frac{\partial H_{\alpha }}{\partial t}+\nabla ·\mathit{u}{H}_{\alpha }=0,$$

(30)

where $H_{\alpha }$ is a Heaviside function having a value of 1 in the liquid phase and 0 in the gas phase. The volume fraction or Volume-Of-Fluid, ${\alpha }_i$, is the average value of $H_{\alpha }$ in computational cell $V_i$ such that

$${\alpha }_i=\frac{1}{V_i}{\int }_V{H}_{\alpha }dV.$$

(31)

This means that $\alpha =1$ for a cell completely filled with liquid and $\alpha =0$ for a cell completely filled with gas. The advection of $\alpha$, the second term in (30), is calculated using the CICSAM interface capturing scheme [46]. The density and viscosity of the fluid in each cell are calculated using an arithmetic average and are given by

$$\rho =\alpha {\rho }_1+(1-\alpha ){\rho }_0$$

(32)

and

$$\mu =\alpha {\mu }_1+(1-\alpha ){\mu }_0,$$

(33)

where subscripts 1 and 0 denote the liquid and gas, respectively. The mass and momentum Equations (26) and (27) are solved using a split projection method [19].

2.3.2. Spatial Discretization of Diffusive Term

Figure 4 shows 2D meshes containing orifice elements: a structured mesh on the left and an unstructured mesh on the right. The dual-cell control volumes are demarcated with dotted lines. Nodes that do not form part of the boundary nor boundary elements are termed interior nodes (e.g., D and G). The orifice boundary nodes are denoted as $O_{Ci}$ and the internal orifice nodes as $O_{Ai}$ and $O_{Bi}$ for left- and right-hand side orifice nodes, respectively.

The diffusive term for an interior node, e.g., D, is discretized as usual:

$${\left.\nabla ·\mu \left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)\right|}_D=\frac{1}{V_D}\sum _{fϵA_D}{\left.\left(\mu \left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)·\Gamma \right)\right|}_f,$$

where $V_D$ is the volume of the dual cell at node D, $A_D$ is the surface (in dotted lines in the figure) enclosing $V_D$ and the subscript f indicates a specific face section of $A_D$. $\Gamma$ is the outward pointing normal vector of the face with magnitude equal to the area of the face. The gradient at the face is calculated using the compact stencil formulation [47].

In the case where the discretization takes place at the orifice nodes, Equation (28) is used to calculate the viscous losses through the orifice. Here, the diffusion contribution at any orifice node (e.g., $O_A$, $O_B$ or $O_{C1}$) is replaced with the empirical loss term by setting $H_O=1$. As an example, the diffusion term at node $O_A$ on the right-hand side of Figure 4 would be computed as follows:

$$\begin{array}{cc}\hfill \nabla ·\mu \left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)=\frac{1}{V_{O_A}}(& {\left.\mu \left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)·\Gamma \right|}_{f_{G-O_A}}\hfill \\ \hfill & +{\left.\mu \left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)·\Gamma \right|}_{f_{D-O_A}}\hfill \\ \hfill & +{\left.\mu \left(\nabla \mathit{u}+\nabla {\mathit{u}}^T\right)·\Gamma \right|}_{f_{P-O_A}}\hfill \\ \hfill & -V_{O_A}\frac{C\rho {\tilde{u}}_O^2}{L}{\mathit{n}}_O).\hfill \end{array}$$

(34)

Here, the subscript $f_{G-O_A}$ refers to the face connecting nodes G and $O_A$. Note that the empirical loss term accounts for the contributions of the faces connecting $O_{C1}$, $O_{C2}$ and $O_B$ to $O_A$.

2.3.3. Numerical Calculation of Orifice Loss Coefficient

The loss coefficient C in the discretized diffusion equation can be determined using Equation (25), repeated here for readability:

$$C=\frac{{\varphi }_{l0visc}^2\frac{K_{fl0}L}{D_H}+{\varphi }_{l0cont}^2{K}_{cont}+{\varphi }_{l0exp}^2{K}_{exp}}{2}.$$

Since only nodal values of the fluid and flow properties are known, the equations from Section 2.1 must be discretized appropriately. This section gives an overview of how the loss coefficient is determined numerically.

The volume flow rate of each phase (e.g., air and water) at orifice nodes i, ${\dot{V}}_{g_i}$ and ${\dot{V}}_{l_i}$, are, respectively, calculated at every time step as

$${\dot{V}}_{g_i}=\sum _{A_i}{\left.\left((1-{\alpha }_i)A_i(\mathit{u}·{\mathit{n}}_O)\right)\right|}_{A_i}$$

(35)

and

$${\dot{V}}_{l_i}=\sum _{A_i}{\left.\left({\alpha }_i{A}_i(\mathit{u}·{\mathit{n}}_O)\right)\right|}_{A_i},$$

(36)

where $A_i$ is the area of each of the faces connecting orifice nodes from opposite sides. With reference to Figure 4, these are the faces connecting $O_{A1}$ and $O_{B1}$, $O_{A2}$ and $O_{B2}$, etc. Further, ${\alpha }_i$ is the VOF value at the cell face as calculated by CICSAM and ${\mathit{n}}_O$ is a unit normal along the orifice in the direction of the average flow. The dryness fraction and mass velocity from Equations (9) and (11) are calculated numerically as

$${\mathrm{x}}_i=\frac{{\rho }_g{\dot{V}}_{g_i}}{{\rho }_l{\dot{V}}_{l_i}+{\rho }_g{\dot{V}}_{g_i}}$$

(37)

and

$$G_i=\frac{{\rho }_l{\dot{V}}_{l_i}+{\rho }_g{\dot{V}}_{g_i}}{A_i},$$

(38)

respectively.

The friction loss coefficient $K_{fl0}$ is calculated using Equation (4), where the Reynolds number is determined from the mixture mass velocity $G_i$ and the liquid viscosity ${\mu }_{l_i}$. The expansion and contraction loss coefficients are given in Equations (6) to (8). Finally, the two-phase multipliers in Equation (25) are calculated using Equations (14), (22) and (23).

2.3.4. CICSAM Correction for Orifice Elements

A key assumption inherent to the CICSAM method is that the vertex is located at the volume centroid. On stretched meshes, this is not the case when using the vertex-centered FV scheme. Since there is a significant amount of stretching introduced in the mesh at orifice cells in this work, the aforementioned assumption introduces a particularly large error. Figure 5 shows a zoomed image of the left of Figure 4, with orifice element cells/volumes in dashed lines and nodes in black dots. Notice the large disparity between the nodes positioned at $x_1$ and $x_2$ and the centroid locations of the corresponding cells, which are indicated with the ×’s ($x_1^{*}$ and $x_2^{*}$). Therefore, for CICSAM to calculate VOF fluxes at the surfaces that bound these cells, the positions of the centroids of the orifice element cells ($x_1^{*}$ and $x_2^{*}$) are used instead of the nodal positions of those cells ($x_1$ and $x_2$). The efficacy of this correction is demonstrated next.

3. Results

The one-dimensional orifice elements have been implemented into the Elemental^® code. The basic implementation was validated using various test cases.

3.1. VOF Advection Test Case

The first test-case assesses if the liquid–gas interface is propagated sharply in 1D orifice elements via the CICSAM scheme (Section 2.3.4). For this purpose, a liquid packet (Fluid 1) is propagated in a channel at a constant velocity of 1 m/s, as shown in Figure 6. The flow is inviscid with slip boundary conditions at the top and bottom boundaries, while the channel contains embedded 1D elements, as depicted in Figure 7. The results in this figure show the fluid interface before and after passing through the one-dimensional orifice elements, both with and without the proposed CICSAM modification. In the case where the centroid locations are not corrected, it can be seen that the interface smears considerably while passing through the one-dimensional orifice elements. However, with the modification applied, the fluid interface remains sharp and matches the expected position (white dotted lines), thus validating the proposed method.

3.2. Violent Slosh Test Case Results

The remainder of the validation work will be focused on a baffled tank undergoing violent slosh as detailed in [48]. As shown in Figure 8, the baffles are placed centrally and contain three long orifices. The tank is subjected to sinusoidal acceleration in the x-direction, which is ramped up over the first $25\%$ of the simulation. The resulting slosh-induced pressures at points $p_1$ and $p_2$ were measured and reported. The $75\%$ fill case is modeled via 2D CFD simulations using Elemental^®. To evaluate the efficacy of the developed orifice element model, the simulations were conducted both with and without the one-dimensional orifice elements followed by comparing the predicted pressure difference between the pressure probes to the measured values. The experimental results are confidential; therefore, time and pressure difference values are normalized with their respective maximum values.

3.2.1. Full-Resolution Simulation Results

A numerical simulation of the experiment is performed, where the 2D governing equations are discretized over the full fluid volume (including the volume inside the orifices). The computational mesh is shown in Figure 9 and is similar to an earlier publication [19]. Note the relatively small computational element sizes required in and around the orifices. Figure 10 shows the local Reynolds numbers of the flow within the respective orifices over time. Note that the Reynolds number is calculated with the velocity being the average velocity in the orifice and the length scale being the height of the orifice. Comparison of the flow pressure differentials will be presented next along with 1D orifice mesh results.

3.2.2. Results Using 1D Orifice Element Model

The results for cases where the computational mesh in the orifices are replaced with one-dimensional orifice elements are presented now. Three different computational meshes with different numbers of orifice elements (see

Table 2. Number of one-dimensional orifice elements in each computational mesh.

	No. of Orifice Elements
	Top Orifice	Middle Orifice	Bottom Orifice
Mesh A	3	7	8
Mesh B	3	4	4
Mesh C	3	3	3

) are used and shown in Figure 11. Zoomed images of the mesh around the orifice areas are shown in Figure 12. Note that, for the top orifice, all three meshes have three orifice elements in the vertical direction since the flow passage is very small.

The measured experimental as well as the numerically calculated full-resolution simulation pressure difference, $p_1-p_2$, are compared with that of the one-dimensional orifice element model simulations in Figure 13 and Figure 14. The flow rate over time through each orifice in the one-dimensional orifice element model is plotted and compared with the full-resolution simulation results in Figure 15, Figure 16 and Figure 17. The total mass in the left-hand side of the tank over time for each simulation is plotted in Figure 18. Additionally, the relative mass difference in the left-hand side of the tank of the one-dimensional orifice element model simulations compared with the full-resolution simulation are shown in Figure 19 for Meshes A to C. Here, the relative mass difference is given by $\frac{m_F-m_O}{m_F}$, where $m_F$ and $m_O$ are, respectively, the masses in the left-hand side of the tank calculated with full resolution and the one-dimensional orifice element model simulations.

4. Discussion

Figure 13 and Figure 14 show the pressure difference over the tank walls calculated using the one-dimensional orifice element model with each computational mesh. It can be seen that the results are in good agreement with both the experimental and full-resolution numerical results. Another metric to consider when evaluating the one-dimensional orifice element model is the flow rates through the orifices.

The flow rates through the orifices (Figure 15, Figure 16 and Figure 17) using the one-dimensional orifice element model for both Mesh A and Mesh B are in reasonable agreement with the full-resolution model for the first $25\%$ of the simulation (Mesh C results are slightly lower). This also correlates with the similarity in the position of liquid–gas interface between the meshes during this period, as shown on the left of Figure 20. From this point onward, the slosh becomes particularly chaotic and the differences in interface position between the different mesh coarseness levels grow more pronounced, as shown on the right of Figure 20. This is reflected in the larger differences in flow rate predicted through the orifices in the latter part of the simulation. Despite this, the overall amount of liquid contained in the two compartments is still predicted with reasonable accuracy. This is evident from the total mass in the left tank (Figure 18 and Figure 19) being similar in the one-dimensional element model simulations and the full-resolution simulation.

The computational cost of the full-resolution CFD simulation and the one-dimensional element model simulations are presented in

Table 3. Comparison of computational cost of full-resolution and one-dimensional element model simulations of violent slosh test case.

	Full Resolution	1D Element Model
		Mesh A	Mesh B	Mesh C
No of Nodes	3757	3458	3173	2263
Time Steps	154,606	129,170	131,062	130,958
CPU Time	67,786 s	58,540 s	46,892 s	33,286 s

. Note that serial calculations were performed in all cases due to the small mesh sizes. As per the table, the one-dimensional element model requires less mesh nodes. This also requires less time steps to complete the same simulation due to the larger allowable time-step sizes (larger elements).

It is interesting to note that Meshes B and C require slightly more time steps than Mesh A, even though the majority of the mesh element sizes are much larger. This can be attributed to the high flow velocities through the top orifice. The average velocity over time through the top orifice is six times the velocity through the middle and bottom orifices. All three meshes have similarly sized elements around the top orifice; therefore, the time step sizes are limited by those elements. Importantly, the one-dimensional element model simulations all require less CPU time, without any significant loss in the accuracy of the pressure differential over the tank walls or the relative mass difference in the tanks. Mesh C achieves the largest speed-up factor of 2 compared with the full-resolution model.

5. Conclusions

In this work, a method to replace fine computational mesh elements within long orifices with one-dimensional elements is presented; this is to gain computational speed when modeling violent slosh in baffled tanks. The one-dimensional orifice element model was implemented into the finite volume VOF code Elemental^®, and integrated seamlessly in the computational mesh with standard finite volume elements for two-phase flow. To achieve this, special consideration had to be taken with regard to discretization of the CICSAM-based VOF advection scheme. Two-dimensional simulations using full-resolution CFD and different meshes using the 1D orifice elements were compared with experimental results of pressure differences across an accelerated baffled tank. The one-dimensional orifice element model was found to decrease the simulation CPU time by up to $50\%$ while still achieving high accuracy in predicted tank pressures.