Effect of geometrical factors interactions on design optimization process of a natural gas ejector

Abstract

Enhancing the ejector entrainment ratio plays an important role in the ejector performance. In this article, a surrogate-based optimization approach along with computational fluid dynamics technique has been employed to optimize the entrainment ratio of a single-phase ejector working with natural gas. Nine ejector geometrical factors have been varied to maximize the ejector entrainment ratio. Validation results of the presented computational fluid dynamics model were in a good agreement with the experimental data from the literature with an average error of 0.6% in the critical mode. Reported results showed that the optimum design achieves entrainment ratio of 19.45% at 12, 2, and 5.2 MPa motive pressure, induced pressure, and discharge pressure, respectively. Moreover, the primary nozzle convergent angle and throat length are insignificant factors. Furthermore, secondary nozzle inclination angle has a minor effect on the entrainment ratio of the optimum design.

Introduction

Natural gas treatment and production plants often choke the natural gas coming from high-pressure wells or compression station in order to meet pipelines or plant pressure. Recovery of the wasted energy in choking process has the potential to make the natural gas production more energy-efficient. In addition, low-pressure gas from the later stages of the process system is to be flared for being uneconomical to recycle. The flaring of low-pressure gas is an energy inefficient process and detrimental to the environment. The high-pressure gas coming from the high-pressure wells or compression station could be used to boost gas production from aging wells or to recompress and recycle the low-pressure gas from the later process stages.

Ejectors have no moving parts, which subsequently increases its reliability and reduce maintenance cost significantly. Moreover, ejectors are distinguished by construction and operation simplicity. Ejector has been introduced by Sarshar ¹ as an energy-efficient and environment-friendly solution to recover the wasted energy in choking process and to reduce gas-to-flare amount. Ejectors are mechanical devices by which a high-pressure stream could be used to induce and accelerate a low-pressure stream.

Because of these advantages, many research works have been conducted on the application of ejectors in natural gas production industry. For example, Andreussi et al. ² presented the design method and result of field tests for a multiphase ejector installed to boost the production from a depleted well in the Gulf of Mexico. Melancon and Nunn ³ discussed the benefits of using the ejector vapor recovery unit in the United States, Reddell field to recover a high-Btu content vent gas and reported the system economics and operational modes. Also, Wu et al. ⁴ carried out a field application economic analysis of the natural gas ejector, being used to boost the production of low-pressure gas well and investigated its performance at different operating conditions. The above-mentioned research works all agreed on the potency of ejectors as a natural gas production boosting device and as a waste gas recovery tool.

Despite the advantages of ejectors, they are still not widely used in natural gas production industry due to its low performance specially at a relatively high back pressure and low induced flow pressure. As a result, numerous investigations have been conducted to improve their performance. Chong et al. ⁵ optimized numerically the geometrical parameters of a natural gas ejector for boosting low-pressure natural gas production. A maximum pressure ratio of 60% has been obtained at 12 MPa motive pressure, 3 MPa induced pressure, and 5.2 MPa discharge pressure. The optimal diameter ratio of mixing tube-to-primary nozzle throat, length-to-diameter ratio of mixing tube, and inclination angle of mixing chamber are 1.6, 4.0, and 28°, respectively.

The entrainment ratio (ω) is used as a key parameter in many research works such as Chen et al.^6,7 and Xue et al. ⁸ The entrainment ratio is defined as the ratio between the secondary flow rate, m_s, and the primary flow rate, $m_p,(\omega =m_s/m_p)$ . Ejectors can operate at three different modes. First, if the back pressure is below the critical back pressure, the primary flow and the secondary flow are both choked, causing constant mass flow rate. As a result, the entrainment ratio stays constantly as the back-pressure changes in this range, this is called the critical mode. Second, if the back-pressure increases beyond the critical pressure, entering into subcritical mode, the entrainment ratio will drop sharply due to absence of the secondary flow choking. Finally, any further increase in back-pressure will lead to back flow or malfunction mode, where the entrainment ratio is zero and, a reverse flow might occur. ⁹

Computational fluid dynamics (CFD) investigation was carried out by Chen et al. ⁶ to optimize some structural parameters of natural gas ejector by maximizing the entrainment ratio. The simulations revealed that, the optimal inclination angle of the mixing chamber, diameter ratio of the mixing tube to primary nozzle throat, length-to-diameter ratio of the mixing tube and inclination angle of the diffuser are 14°, 1.7, 5.0, and 1.43°, respectively. Based on their previous studies, Chen et al. ⁷ conducted an optimization study to maximize both the entrainment ratio and pressure ratio. The results showed that, the maximum entrainment ratio achieved when the primary nozzle exit position was varied from 3.6 to 7.2 mm and the mixing tube length-to-diameter ratio varied from 2 to 8, while the maximum pressure ratio achieved when the primary nozzle exit position ranged from 1.2 to 7.2 mm and the mixing tube length-to-diameter ratio ranged from 3 to 7.

Xue et al. ⁸ compared the performance of supersonic air ejector with three different nozzles shapes, conical nozzle, petalage nozzle, and crenation nozzle. Compared to the conical nozzle, the entrainment ratio of petalage nozzle was slightly smaller, but the critical back pressure increased by 5.2%. The entrainment ratio of crenation nozzle was slightly higher, and the critical back pressure decreases by 2.1%. Chen et al. ¹⁰ proposed another design enhancement to maximize ejector capacity. The proposed enhancement was to install bypass at ejector low-pressure district. The results showed that, the optimal length ratio is about 1.1, and the ejector performance with the optimal bypass has relatively higher entrainment ratio and lower critical back-pressure compared to conventional the ejector.

The above-mentioned investigations reported a significant effect of ejector geometrical factors on its performance. However, no structural optimization on the primary nozzle of natural gas ejectors has been reported. Moreover, there is no study considered the interactions between the natural gas ejector geometrical factors, since the one factor at a time approach is dominant in all studies. In this study, the surrogate modeling along with CFD are used to obtain the optimum value of nine natural gas ejector geometrical factors that maximize the entrainment ratio. According to this approach, the interactions between ejector factors are considered in the optimization process.

Mathematical model

In this work, the commercial software ANSYS Fluent 16.2 is employed as a numerical solver. ¹¹ The flow in the ejector is considered to be compressible, steady-state, and turbulent flow. The two-dimensional (2D) axisymmetric simplification is used in this study. The conservation equations are given by

Mass conservation

$$\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x}\left(\rho v_x\right)+\frac{\partial }{\partial r}\left(\rho v_r\right)+\frac{\rho v_r}{r}=0$$

(1)

Momentum conservation

$$\begin{array}{l}\frac{\partial \rho v_x}{\partial t}+\frac{1}{r}\frac{\partial }{\partial x}\left(r\rho v_x{v}_x\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r\rho v_r{v}_x\right)=-\frac{\partial p}{\partial x}+\frac{1}{r}\frac{\partial }{\partial x}\\ \left[r\eta \left(2\frac{\partial v_x}{\partial x}-\frac{2}{3}\left(\frac{\partial v_x}{\partial x}+\frac{\partial v_r}{\partial r}+\frac{v_r}{r}\right)\right)\right]+\frac{1}{r}\frac{\partial }{\partial r}\left[r\eta \left(\frac{\partial v_x}{\partial r}+\frac{\partial v_r}{\partial x}\right)\right]+F_x\end{array}$$

(2)

$$\begin{array}{c}\frac{\partial \rho v_r}{\partial t}+\frac{1}{r}\frac{\partial }{\partial x}\left(r\rho v_x{v}_r\right)+\frac{1}{r}\frac{\partial }{\partial r}\left(r\rho v_r{v}_r\right)\hfill \\ =-\frac{\partial \rho }{\partial r}+\frac{1}{r}\frac{\partial }{\partial x}\left[r\eta (\frac{\partial v_x}{\partial r}+\frac{\partial v_r}{\partial x})\right]\hfill \\ +\frac{1}{r}\frac{\partial }{\partial r}\left[r\eta (2\frac{\partial v_r}{\partial r}-\frac{2}{3}(\frac{\partial v_x}{\partial x}+\frac{\partial v_r}{\partial r}+\frac{v_r}{r}))\right]+F_r\hfill \end{array}$$

(3)

where x and r are the axial and radial coordinates, respectively. v_x and v_r are the axial and radial velocities, respectively. While F_x and F_r are the external forces in axial and radial directions, respectively. ρ is gas density and $\eta$ is the gas dynamic visocity.

Renormalization group (RNG) $k-\epsilon$ model is recommended to be adopted to simulate the turbulence phenomena inside the ejector, that is, because it is suitable to simulate complex shear flows involving rapid strain which is the case in this study.^6,7,12 The RNG $k-\epsilon$ model is based on two extra transport equations for the turbulent kinetic energy, $k$ , and turbulent dissipation rate, ε ¹¹

$$\frac{\partial \rho k}{\partial t}+\frac{\partial \rho k{u}_i}{\partial x_i}=\frac{\partial }{\partial x_j}\left[(\eta +\frac{{\eta }_{\mathrm{tu}}}{{\alpha }_k})\frac{\partial k}{\partial x_j}\right]+\mathrm{Pk}-\rho \epsilon -2\rho \epsilon M_{\mathrm{tu}}^2$$

(4)

$$\frac{\partial \rho \epsilon }{\partial t}+\frac{\partial \rho \epsilon u_i}{\partial x_i}=\frac{\partial }{\partial x_j}\left[(\eta +\frac{{\eta }_{\mathrm{tu}}}{{\alpha }_{\epsilon }})\frac{\partial \epsilon }{\partial x_j}\right]+C_2\frac{\epsilon }{k}\mathrm{Pk}-C_{\epsilon 2}^{*}\rho \frac{{\epsilon }^2}{k}$$

(5)

where $C_{\epsilon 2}^{*}=C_3+(\left(C_1{\beta }^3(1-\beta /C_4)\right)/(1+C_5{\beta }^3))$ , $\beta =\mathrm{\Omega }k/\epsilon$ and $\mathrm{\Omega }=\sqrt{2{\mathrm{\Omega }}_{\mathrm{ij}}{\mathrm{\Omega }}_{\mathrm{ij}}}$ . $M_{\mathrm{tu}}$ is the turbulent Mach number that accounts for compressibility effect and the model constants have the following default values: C₁ = 0.0845, C₂ = 1.42, C₃ = 1.68, C₄ = 4.38, and C₅ = 0.012.

In this study, the entrainment ratio, used as a key parameter, is not affected by the flow boundary layer according to these conditions. This conclusion is obtained from a comparison between these CFD results and Chong et al. ¹² results. Chong solved the flow boundary layer and the entrainment ratio that are in a good agreement with these results. Therefore, the scalable wall function is the best selection to solve the near wall cells in this study. The scalable wall functions is used to force the usage of the log law in conjunction with the standard wall functions approach since it is so hard to adjust ejector walls $y^{+}$ as required because of the huge variations in gas speed inside the supersonic ejector. The average wall $y^{+}$ values during the simulations are almost 42.

The second-order upwind discretization scheme is adopted to solve the convection terms. While, the coupled pressure-based algorithm is selected to obtain ejector pressure field at low dominant Mach number. ¹³ The boundary conditions used to solve the governing equations are pressure inlets for the motive and induced flow inlets, and pressure outlet for ejector outlet. The flow condition of the natural gas ejector used in this study is supersonic with zero velocity at the ejector inlets. As a result, the ejector walls are considered to be adiabatic due to the small residence time of gas. The temperature is kept at 300 K. Turbulent intensity and hydraulic diameter, listed in Table 1, are used for turbulent boundary conditions.

Table 1.

Boundary condition values used during the validation and the optimization simulations.

Boundary	Pressure (Mpa)	Turbulent intensity (%)	Hydraulic diameter (m)
Motive flow (validation)	1	10	0.02
Induced flow (validation)	0.5	10	0.0025
back flow (validation)	0.4–0.7	5	0.02
Motive flow (optimization)	14	10	0.02
Induced flow(optimization)	2	10	0.076
Back flow (optimization)	5.2	5	0.02

Due to the symmetry of the ejector, the 2D axisymmetric simplification is used for calculation time reduction purpose. The 2D full domain and the 2D axisymmetric domain are demonstrated in Figure 1(a) and (b). A structured mesh with hexahedral elements was generated using ANSYS meshing tool as shown in Figure 1(c). Methane is selected as the working fluid because it represents more than 85% of natural gas volume. ¹⁴ The density is modeled using the real gas equation of Soave Redlich Kwong, equation (6) ¹¹

$$P=\frac{\mathrm{RT}}{V-b}-\frac{a}{V(V+b)}$$

(6)

$$b=0.08664\frac{R{T}_C}{P_C}$$

(7)

$$\begin{array}{c}a=0.42747\frac{R^2{T}_C^2}{P_C}\hfill \\ {[1+(0.48+1.574w+0.176w^2)(1-\sqrt{\frac{T}{T_C}})]}^2\hfill \end{array}$$

(8)

where P, V, T, and R are the gas pressure, specific volume, temperature, and universal constant, respectively. $T_C$ and $P_C$ are the gas critical temperature and pressure and $w$ is gas acentric factor. Fluid viscosity, specific heat capacity, and thermal conductivity are derived from NIST database ¹⁵ and given by

$$\eta =\{\begin{array}{c}0.01726268+0.00077648T-1.54\times {10}^{-5}{T}^2-1.71\times {10}^{-7}{T}^3+1.15\times {10}^{-9}{T}^4\hfill \\ \multicolumn{1}{c}{\mathrm{for}100<T<179.9\mathrm{K}}\\ \multicolumn{1}{c}{-0.1156432+0.0029433T-3.09\times {10}^{-5}{T}^2+1.72\times {10}^{-7}{T}^3-5.33\times {10}^{-10}{T}^4}\\ \multicolumn{1}{c}{\mathrm{for}179.9<T<300\mathrm{K}}\end{array}$$

(9)

$$C_p=\{\begin{array}{c}-658,833.1+23,960.32T-346.62T^2+2.48T^3-0.0088T^4\mathrm{for}100<T<205.2\mathrm{K}\hfill \\ \multicolumn{1}{c}{20,900,000-399,555.3T+3049.36T^2-11.62T^3+0.022T^4\mathrm{for}205.2<T<300\mathrm{K}}\end{array}$$

(10)

$$K_{\mathrm{cond}}=-10.759+0.4545-0.00782T^2+7.26\times {10}^{-5}{T}^3-3.94\times {10}^{-7}{T}^4\mathrm{for}100<T<300\mathrm{K}$$

(11)

Figure 1.

The ejector flow domain. (a) Ejector 2D wholeflow domain, (b) ejector 2D axisymmetric half domain, and (c) part of the generated mesh for ejector 2D axisymmetric flow domain.

Mesh independency and validation study

Roache^16–18 suggested the grid convergence index (GCI) as a quantitative measure for the grid convergence. The GCI is based upon a grid refinement error estimator derived from the theory of the generalized Richardson extrapolation. The GCI, equation (12), is a measure of how far the computed value is away from the value of the asymptotic numerical value. Consequently, it indicates how much the solution would change with a further refinement of the grid. A small value of GCI indicates that the computation is within the asymptotic range

$$\mathrm{GC}{I}^{F\mathrm{i}\mathrm{ne}}=\frac{F_s\left|ϵ\right|}{r^p-1}$$

(12)

where $F_s$ is the safety factor which equals to 1.25 for three grids comparison, $ϵ$ is the relative error measure of the key variable $f$ , $r$ is the grid refinement ratio, and p is the order of the discretization method. r and $ϵ$ are mathematically defined as

$$r={\left(\frac{N_1}{N_2}\right)}^{\frac{1}{\mathrm{Di}}}$$

(13)

$$ϵ=\frac{f_2-f_1}{f_2}$$

(14)

where $N_1$ and $N_2$ are the fine and coarse grid number of cells, respectively, and $\mathrm{Di}$ equals to 2 for 2D geometry. $f_1$ and $f_2$ are respectively, the fine and coarse grid numerical solutions. In three meshes grid refinement study, where the grid refinement ratio does not equal $(r_{12}\ne r_{23})$ , the order of the discretization method $p$ can be calculated iteratively as

$$p=0.5p_0+0.5\frac{\ln \left(\frac{(r_{12}^p-1)e_{23}}{(r_{23}^p-1)e_{12}}\right)}{\ln \left(r_{12}\right)}$$

(15)

where $p_0$ is the previous iteration of $p$ . The following equation is used as a first guess of $p$ , the iteration will stop if $|p-p_0|<1\times {10}^{-5}$

$$r={\left(\frac{N_1}{N_2}\right)}^{\frac{1}{\mathrm{Di}}}$$

(16)

where $e_{23}=f_3-f_2$ and $e_{12}=f_2-f_1$ .

The solution is in the asymptotic range when α, which defined in equation (17), approaches unity $(\alpha \approx 1)$ . The value of $f$ when the grid spacing approaches zero $(h\to 0)$ can be estimated using the Richardson extrapolation as in equation (18). The ratio $R$ which stated in equation (19) classified the possible convergence conditions into three groups as following:

• Oscillatory convergence, $R<0$ .

• Monotonic convergence, $0<R<1$ .

• Divergence, $R>1$

$$\alpha =\frac{r_{12}^p{\mathrm{GCI}}_{12}^{\mathrm{fine}}}{{\mathrm{GCI}}_{23}^{\mathrm{fine}}}$$

(17)

$$f_{\mathrm{exact}}=f_1+\frac{f_1-f_2}{r_{12}^p-1}$$

(18)

$$R=\frac{{ϵ}_{12}}{{ϵ}_{23}}$$

(19)

In this study, over 90 geometrically different ejectors are needed for the optimization process. It is too much time consuming to study the mesh independency on each ejector. So, the authors proposed the average cell area $(A_{C_{\mathrm{avg}}})$ method to estimate the appropriate mesh size for each ejector. $A_{C_{\mathrm{avg}}}$ is the sum of each cell area in the domain divided by the number of cells, $A_{C_{\mathrm{avg}}}$ defined as

$$A_{C_{\mathrm{avg}}}=\frac{\sum _{\stackrel{·}{l}=1}^N{A}_i}{N}=\frac{A_D}{N}$$

(20)

where $A_i$ , $A_D$ , and $N$ are the cell area, domain area, and the number of cells in the domain, respectively. According to $A_{C_{\mathrm{avg}}}$ calculations, the appropriate number of cells can be easily calculated. Therefore, a sample ejector geometry picked from the initial sampling plan for the grid convergency index study and the average cell area $A_{C_{\mathrm{avg}}}$ is calculated. $A_{C_{\mathrm{avg}}}$ value for this study is used to estimate the appropriate number of cells for the out of study ejectors. ANSYS meshing tool is used to generate a structured hexahedral elements mesh. Three grid levels are used for the $\mathrm{GCI}$ analysis with a minimum 10% mesh refinement $(r=1.1)$ as recommended by Roache. ¹⁶ Table 2 summarizes the grid convergency calculations using GCI method and three grid levels.

Table 2.

Grid convergence index calculations results.

$A_D(\mathrm{m}{\mathrm{m}}^2)$	i	Ni	$A_{c_{\mathrm{avg}}}(\mathrm{m}{\mathrm{m}}^2)$	ωi (%)	${\mathrm{GCI}}_{i,i+1}^{\mathrm{fine}}(\mathrm{\%})$	α	R
2333	$0$ a			−5.68		1.0079	0.7293
	1	294,360	0.0079	−5.766
					1.7825
	2	235,190	0.0099	−5.812
					2.7548
	3	191,875	0.0121	−5.875

a

The value at zero grid space (f_exact equation (18)).

The following conclusions have been obtained from the GCI analysis:

• The obtained values for α are close to 1, so the results are in the asymptotic range.

• The ratio $R$ is in the monotonic convergence range $(0<R<1)$ .

• There is a reduction in the $\mathrm{GCI}$ value for the successive grid refinements $({\mathrm{GCI}}_{12}^{\mathrm{fine}}<{\mathrm{GCI}}_{23}^{\mathrm{fine}})$ . This indicates that the dependency of the numerical results on the cell size has been reduced.

• The extrapolated value (f_exact) is only slightly changed from the finest grid solution. This confirms that the grid independency solution has been achieved.

• The minimum appropriate value for $A_{C_{\mathrm{avg}}}$ is 0.0079. However, the value of 0.007 will be used for the future grid generations to insure the mesh independency for different geometry configurations.

The accuracy and the stability of the numerical solution are highly correlated to mesh quality. Therefore, the mesh quality characteristic is carefully monitored during mesh generation. The minimum orthogonal quality, maximum skewness, and maximum aspect ratio for the constructed grids throughout this study are 0.85, 0.3, and 30, respectively. Meshes are constructed manually using ANSYS meshing tool for every simulated ejector during the optimization process with $A_{C_{\mathrm{avg}}}$ equal to 0.007, and the maximum wall $y^{+}$ value is 190. Refining the mesh to get smaller $y^{+}$ value will dramatically increase mesh number and its effect on solution accuracy will be negligible.

The CFD model is validated with the experimental data from Chong et al. ¹² In Chong et al.’s study, air is used as a working fluid in the experiments. Therefore, air is used during the validation process and the real gas equation of Soave Redlich Kwong is used for density calculations with the temperature-dependent polynomial functions provided by ANSYS Fluent database for air. The simulation is performed at 1 MPa motive pressure, 0.2 MPa secondary pressure, and discharge pressure ranging from 0.4 to 0.7 MPa. The mesh independency study is conducted first using the GCI analysis before the validation process at 0.4 MPa discharge pressure. According to the GCI analysis summarized in Table 3, the mesh independence is achieved.

Table 3.

Grid convergence index calculations results for ejector model validation.

$A_D(\mathrm{m}{\mathrm{m}}^2)$	i	Ni	$A_{c_{\mathrm{avg}}}(\mathrm{m}{\mathrm{m}}^2)$	ωi (%)	${\mathrm{GCI}}_{i,i+1}^{\mathrm{fine}}(\mathrm{\%})$	α	R
2341	$0$ a			45.79		0.9995	0.3415
	1	257,550	0.009	45.766
					0.06
	2	206,056	0.0113	45.745
					0.1182
	3	153,600	0.0152	45.683

a

The value at zero grid space (f_exact equation (18)).

Figure 2 shows that predictions of the entrainment ratio are in a good agreement with the experimental data with an average error of 0.6% in the critical mode. In addition, the model is able to predict ejector operational modes. Therefore, the CFD model is used in the optimization process to predict with considerable accuracy the entrainment ratio.

Figure 2.

Comparison of entrainment ratio experimental data and CFD results for different discharge pressures at 1 MPa motive pressure and 0.2 MPa induced pressure.

Optimization methodology

The most popular approaches used for the optimization are the one factor at a time approach and the factorial approach in which the objective function affected by more than one factor. The one factor at a time method depends on selecting a starting point, or baseline and set of levels for each factor, then successively varying each factor over its range with the other factors held constant at the baseline level. The major disadvantage of this method is that, it does not take into consideration any possible interaction between the factors. This leads to the failure of the one factor at a time method to produce the same effect on the response at different levels of another factor. While, factors are varied together in the factorial approach. Hence, the interactions between factors can be investigated in which better results can be reached. ¹⁹ In this phase, the regression Kriging surrogate model is used to utilize the factorial approach in natural gas ejector geometrical optimization process. The optimization process passes through four main steps: defining the optimization goal and factor ranges, constructing the sampling plan, building the surrogate model, and finally defining the infill process and convergence criterion. The coupling between the CFD results and the optimization tool is performed manually.

Optimization goal and factors

The optimization goal is to maximize the entrainment ratio, ω. In order to achieve such objective, the geometrical factors of the natural gas ejector, shown in Figure 3, such as the primary nozzle convergent angle, ${\theta }_{\mathrm{pc}}$ , primary nozzle divergent angle, ${\theta }_{\mathrm{pd}}$ , primary nozzle exit diameter, $D_p$ , primary nozzle exit position, $\mathrm{NXP}$ , secondary nozzle inclination angle, ${\theta }_s$ , mixing tube diameter, $D_{\mathrm{mt}}$ , diffuser inclination angle, ${\theta }_D$ , mixing tube length, $L_{\mathrm{mt}}$ , and primary nozzle throat length, $L_t$ have been varied. All lengths and diameters are normalized by the primary nozzle throat diameter, $D_t$ , and expressed as dimensionless geometric ratios. The objective function becomes as follows

$$\omega =f({\theta }_{\mathrm{pc}},{\theta }_{\mathrm{pd}},\frac{D_p}{D_t},\frac{\mathrm{NXP}}{D_t},{\theta }_s,\frac{D_{\mathrm{mt}}}{D_t},{\theta }_D,\frac{L_{\mathrm{mt}}}{D_t},\frac{L_t}{D_t})$$

(21)

Figure 3.

Natural gas ejector geometrical factors.

Ejector primary nozzle throat diameter kept constant at 4.6 mm, and the entire ranges at which these factors studied are summarized in Table 4.

Table 4.

Ejector geometrical factors ranges under study.

Factor	${\theta }_{\mathrm{pc}}$	${\theta }_{\mathrm{pd}}$	$D_p/D_t$	$\mathrm{NXP}/D_t$	${\theta }_s$	$D_{\mathrm{mt}}/D_t$	${\theta }_D$	$L_{\mathrm{mt}}/D_t$	$L_t/D_t$
Range	4–14	1–15	1.1–1.45	0.2174–8.4783	13.5–38	1.55–2.97	1–15	1.55–60	0–1

Kriging and initial sampling plan construction

Kriging is an interpolating technique which considers both the distance and the degree of the variation between known data points when predicted values in unknown areas. ²⁰ A kriged estimate is a weighted linear combination of the known sample values around the point to be estimated. Due to the great number of the Kriging model parameters, it is well suited for a highly nonlinear function with multi extremes. The maximum likelihood estimation (MLE) method is used to tune regression Kriging parameters. The MLE is a method of estimating the parameters of a statistical model given observations, by finding the parameter values that maximize the likelihood of making the observations, given the parameters. The concentrated ln-likelihood function is given by ²⁰

$$\ln \left(L\right)\approx -\frac{n}{2}\ln \left({\sigma }^2\right)-\frac{1}{2}\ln |\psi +\lambda I|$$

(22)

$${\sigma }^2=\frac{{(y-1\mu )}^T{(\psi +\lambda I)}^{-1}(y-1\mu )}{n}$$

(23)

$$\mu =\frac{1^T{(\psi +\lambda I)}^{-1}y}{1^T{(\psi +\lambda I)}^{-1}1}$$

(24)

where μ and ${\sigma }^2$ are the mean and the standard deviation, respectively; $n$ is the number of the sample points; $y$ is a column vector of the observations; $1$ is an n × 1 column vector of ones; $T$ refers to matrix transpose; and $\psi$ is the correlation matrix which correlate the observed data, with each other, using the basis function expression as following

$$\psi =\begin{array}{ccc}\mathrm{cor}[Y\left(x^1\right),Y\left(x^1\right)]\hfill & \dots \hfill & \mathrm{cor}[Y\left(x^1\right),Y\left(x^n\right)]\hfill \\ \multicolumn{1}{c}{⋮}& \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{\mathrm{cor}[Y\left(x^n\right),Y\left(x^1\right)]}& \dots \hfill & \mathrm{cor}[Y\left(x^n\right),Y\left(x^n\right)]\hfill \end{array}$$

(25)

where

$$\mathrm{cor}[Y\left(x^{(i)}\right),Y\left(x^{(l)}\right)]=\exp (-\sum _{j=1}^K{\Theta }_j|x_j^{(i)}-x_j^{(l)}{|}^{P_j})$$

(26)

The correlation matrix depends on the absolute distance between the sample points $|x_j^{\left(i\right)}-x_j^{\left(l\right)}|$ and the unknown parameters ${\mathrm{\Theta }}_j$ and $P_j$ that will be optimized to maximize the concentrated ln-likelihood function equation (22). Kriging basis function uses different $\mathrm{\Theta }$ and $P$ for each factors, since $K$ represents the number of factors that allows the width and the smoothness of the function to vary from factor to factor. The regression term $(\lambda I)$ which appears in the concentrated ln-likelihood function is used to add the regression constant $(\lambda )$ to the leading diagonal of the correlation matrix $\Psi$ using the $n\times n$ identity matrix (I). The regression constant is also unknown and must be optimized.

In order to tune Kriging model parameters, Latin hypercubes technique has been used to construct a sampling plan of 90 points. The constructed sampling plan responses were evaluated using CFD and then used to find the optimum $\mathrm{\Theta }$ , $P$ , and $\lambda$ that maximizes the concentrated ln-likelihood function. The optimized parameters are used by Kriging predictor to predict new values based on the observed data as follows

$$\hat{y}=\mu +{\psi }^T{(\Psi +\lambda I)}^{-1}(y-1\mu )$$

(27)

Infill and convergency criterion

The expected improvement infill strategy has been used to improve the accuracy in the region of the predicted optimum by Kriging to obtain an accurate optimal value quickly, and enhance the general accuracy of the model. The expected improvement approach computes how much improvement can be expected if decided to sample at a given point (x) as follows ²⁰

$$\begin{array}{c}E\left[I\left(x\right)\right]=(y_{\min }-\hat{y}\left(x\right))\frac{1}{2}\hfill \\ +[\frac{1}{2}\mathrm{erf}\left(\frac{y_{\min }-\hat{y}\left(x\right)}{s\sqrt{2}}\right)+s\frac{1}{\sqrt{2\pi }}\exp \left[\frac{-{(y_{\min }-\hat{y}\left(x\right))}^2}{2s^2}\right]]\hfill \end{array}$$

(28)

where $\mathrm{erf}(·)$ is the Gauss error function encountered in integrating the normal distribution

$$\mathrm{erf}\left(\frac{y_{\min }-\hat{y}\left(x\right)}{s\sqrt{2}}\right)=\frac{2}{\sqrt{\pi }}\sum _{n=0}^{\infty }\frac{{(-1)}^n{\left(\frac{y_{\min }-\hat{y}\left(x\right)}{s\sqrt{2}}\right)}^{2n+1}}{n!(2n+1)}$$

(29)

$s^2$ is the re-interpolation mean squared error in a Gaussian process–based prediction and can be calculated as

$$s^2={\sigma }_{\mathrm{ri}}^2[1-{\psi }^T{\Psi }^{-1}\psi +\frac{1-1^T{\Psi }^{-1}\psi }{1^T{\Psi }^{-1}1}]$$

(30)

${\sigma }_{\mathrm{ri}}^2$ is the re-interpolation standard deviation which defined as

$${\sigma }_{\mathrm{ri}}^2=\frac{{(y-1\mu )}^T{(\Psi +\lambda I)}^{-1}\Psi {(\Psi +\lambda I)}^{-1}(y-1\mu )}{n}$$

(31)

The non-bolded “1” is a single one. $y_{\min }$ is the minimum response value in the sampling plan. The expected improvement for the prediction, calculated by equation (28) is for a minimization problem. Therefore, the negative of the observed data has been taken so that the problem could be treated as minimization problem and the data mapped into $[0,1]$ space. Genetic algorithm in MATLAB global optimization toolbox ²¹ has been used as a global search algorithm to optimize Kriging model parameters using equation (22), and to search the constructed model using equation (28) for the point which maximizes the expected improvement for the prediction. Genetic algorithm properties have been optimized to give the best optimization result for the same initial population for both the concentrated ln-likelihood function, and the expected improvement function. Genetic algorithm operators and parameters used in this study have been summarized in Table 5. The ranges at which Kriging parameters $\mathrm{\Theta }$ , $P$ , and $\lambda$ are searched are ${10}^3\mathrm{to}{10}^{-6}$ , $1-2$ , and $1\mathrm{to}{10}^{-6}$ , respectively.

Table 5.

Genetic algorithm operators and parameters.

Parameter	Concentrated ln-likelihood function	Expected improvement function
Population type	Double vector	Double vector
Population size	700	500
Fitness scaling	Rank	Rank
Selection operation	Remainder	Roulette wheel
Elite count	35	25
Crossover operation	Scattered	Arithmetic
Crossover fraction	0.8	0.8
Mutation operation	Adaptive feasible	Adaptive feasible

Figure 4 demonstrates the procedures of optimization process in order to obtain the optimum ejector geometrical factors.A suitable convergence criterion is adopted to stop the surrogate infill process which is a hard decision for the expectation of improvement infill strategy, since the estimated mean squared error of Gaussian process-based models is often an underestimate and the search may be stopped prematurely. ²⁰ The surrogate infill process is stopped if the expectation of improvement is smaller than 1% of the range of observed objective function values for ten successive infill processes. After the convergency criterion for the infill process reached which is based on maximizing the expectation of improvement, the infill process based on maximizing Kriging prediction starts. The infill process based on maximizing Kriging prediction is terminated as the maximum value of the objective function founded through the infill process, based on maximizing the expectation of improvement, does not change for 20 successive infill processes.

Figure 4.

Flowchart of ejector optimization procedure.

The construction of the required ejector geometries and corresponding grids was performed manually for each ejector. The setting inputs to simulation software such as boundary conditions, physical model, and solver settings were applied though ANSYS Fluent journal commands. The entrainment ratio evaluation from ANSYS Fluent software and input to MATLAB optimization code were performed manually. All simulations have been performed using a laptop with 8 GB of RAM and 2.40 GHz CPU. In average, the calculation time is around 8–10 h for each run.

Results and discussions

The initial regression Kriging model has been trained with a set of 90 points. Most of the tested ejectors in the initial 90 points sampling plan are in the back-flow mode, since the entrainment ratio is negative and only 11 ejectors from the data set work with a positive entrainment ratio. As mentioned above, the infill process was stopped if the convergency criteria reached. The infill process based on maximizing the expectation of improvement was stopped after 106 infill processes, and the maximum value of the entrainment ratio was not changed for the 20 infill processes based on maximizing Kriging prediction. The final sampling plan of 216 points has been reached. The simulated 126 infill points have 103 ejectors working with a positive entrainment ratio. The total number of the simulated ejectors working with a positive entrainment ratio was 114. In contrast to the initial sampling plan, the median of the entrainment ratios is positive which means that the infill process is concentrated in the area of the maximum entrainment ratio as demonstrated in Table 6.

Table 6.

Sampling plans of statistical data of the entrainment ratio.

Sampling plan	Maximum (%)	Minimum (%)	Median (%)	Mean (%)	Standard deviation (%)
Initial	12.67	−403.03	−86.16	−102.04	±92.37
Final	19.45	−403.03	2.81	−40.48	±80.81

The model proved its ability to fit the data set successfully since the training root mean square error for the data set is smaller than $6.6\times {10}^{-3}$ with a coefficient of correlation close to $1$ $(R^2>0.999)$ for both sampling plans. This indicates the Kriging model superior capability of data fitting compared with traditional models used in curve fitting. The tuned Kriging parameters for the initial and the final sampling plan are summarized in Table 7. Kriging parameter, $\mathrm{\Theta }$ , could be used to establish the order of importance of variables. According to tuned $\mathrm{\Theta }$ , the initial and converged models both emphasize on the insignificancy of the primary nozzle convergent angle and its throat length $({\theta }_{\mathrm{pc}},L_t)$ on the entrainment ratio and the significance of the other factors except the mixing tube length, $L_{\mathrm{mt}}$ . The significant effect of $L_{\mathrm{mt}}$ has been tested using CFD by varying $L_{\mathrm{mt}}$ while the other factors kept constant at their optimum values. The mixing tube length effect on the entrainment ratio has been tested at the maximum and minimum values indicated in Table 4. The ejector with the maximum $L_{\mathrm{mt}}$ value gives entrainment ratio of $3.21\mathrm{\%}$ while the ejector with the minimum $L_{\mathrm{mt}}$ value gives entrainment ratio of $-13.65\mathrm{\%}$ .

Table 7.

Converged regression Kriging model tuned parameters.

Factor	${\Theta }_{\mathrm{initial}}$	$P_{\mathrm{initial}}$	${\lambda }_{\mathrm{initial}}$	${\Theta }_{\mathrm{final}}$	$P_{\mathrm{final}}$	${\lambda }_{\mathrm{final}}$
${\theta }_{\mathrm{pc}}$	${10}^{-6}$	1.92	${10}^{-6}$	${10}^{-6}$	1.955	${10}^{-3}$
${\theta }_{\mathrm{pd}}$	0.0665	2		0.3157	1.748
$D_p$	0.2197	2		0.8361	2
$\mathrm{NXP}$	1.0525	1.33		2.7496	1.978
${\theta }_s$	1.4845	2		2.11	2
$D_{\mathrm{mt}}$	0.8376	2		5.4884	1.93
${\theta }_D$	0.1578	2		0.7729	2
$L_{\mathrm{mt}}$	$0.1589$	2		${10}^{-6}$	1.01
$L_t$	${10}^{-6}$	1.91		${10}^{-6}$	1.05

Furthermore, the optimum value of $L_{\mathrm{mt}}$ gives an ejector with entrainment ratio of $19.45\mathrm{\%}$ . Simulation results confirmed the initial model result on the significance of $L_{\mathrm{mt}}$ and the failure of the final model to estimate the significant effect of $L_{\mathrm{mt}}$ .

Figure 5 illustrated the main effect of the significant factors according to the converged Kriging model prediction on the entrainment ratio. The plot depicted by varying one factor at a time in the range stated in Table 4, while the other factors were kept constant at the values which gives the maximum entrainment ratio of (19.35%) which stated in Table 7.

Figure 5.

Main effect plot of six ejector geometrical factor according to converged Kriging model.

The negative values of the entrainment ratio presented reveal the sensitivity of the studied factors. The demonstrated trends in Figure 5 showed an agreement with the trends found in the literature works.^6,7 The effect of the primary nozzle divergent angle, ${\theta }_{\mathrm{pd}}$ , on the entrainment ratio is very small at smaller angle values. However, increasing ${\theta }_{\mathrm{pd}}$ beyond some point will reduce the expansion time inside the primary nozzle divergent section and a normal shock wave is formed at the nozzle exit. This leads to increase in the pressure and the primary flow exit from the secondary nozzle inlet. Increasing the primary nozzle exit diameter, $D_p$ , leads to increase in the primary nozzle divergent section length which increases the expansion time of the primary flow inside the primary nozzle divergent section and also increases the jet diameter and subsequently increases the entrainment ratio. Nevertheless, at some point increasing $D_p$ causes the entrainment ratio to decrease again. The entrainment ratio also decreases by increasing the primary nozzle exit position, $\mathrm{NXP}$ , so the optimum value is near the minimum value of $\mathrm{NXP}$ in the investigated range. The secondary nozzle inclination angle, ${\theta }_s$ , has no significant effect on the entrainment ratio compared to the other factors. However, the optimum value is at the minimum value.

Kriging model ranked the mixing tube diameter, $D_{\mathrm{mt}}$ , as the most effective factor on the entrainment ratio. The mixing tube diameter, $D_{\mathrm{mt}}$ , with the jet diameter forms the annuals diameter in which the induced flow accelerates, larger $D_{\mathrm{mt}}$ causes the annuals diameter to become larger and the primary jet flow effect on the near wall flow is small which decreases the induced flow speed. In addition, much larger $D_{\mathrm{mt}}$ causes back flow to take place in the mixing tube and the flow exits from the secondary nozzle inlet. However, too much small $D_{\mathrm{mt}}$ decreases both annuals diameter and the entrainment ratio. The optimum $D_{\mathrm{mt}}$ is near the minimum value of the investigated range. Diffuser inclination angle, ${\theta }_D$ , is the parameter used to determine diffuser length and how far the back pressure from mixing tube exit. In large diffuser angle, the flow is separated from the wall and locally recirculated inside the diffuser. Furthermore, large diffuser angle produces a higher pressure at the mixing tube end which causes the entrainment ratio to decrease. Small angle diffuser solves the recirculation problem. However, a small angle leads to longer diffuser and consequently, increases energy loss due to friction. Thus, the optimum ${\theta }_D$ is also near the minimum value of the investigated range.

Figure 5 reveals that, the one factor at a time plot comprises no information about the interaction between the factors, since it just demonstrates the main effect of the factors under study. In order to build a mental image of the constructed Kriging model, two-factor slices have been extracted from the model and arranged as tiles in a matrix as shown in Figure 6. The matrix of the contour plots in Figure 6 is constructed by plotting diffuser inclination angle, ${\theta }_D$ , and the primary nozzle exit diameter, $D_p$ , against each other in each tile, the x-axis in each tile represents ${\theta }_D$ while the y-axis represents $D_p$ . The rows represent different values of the primary nozzle exit position, $\mathrm{NXP}$ , and the columns represent the mixing tube diameter, $D_{\mathrm{mt}}$ , values. As shown in Figure 6, at small $D_{\mathrm{mt}}$ , and $\mathrm{NXP}$ values, the optimum values of ${\theta }_D$ and $D_p$ are always restricted at small ${\theta }_D$ value, and large $D_p$ value. Increasing $D_{\mathrm{mt}}$ at small NXP values will increase the optimum ${\theta }_D$ value. Moreover, increasing $D_{\mathrm{mt}}$ and NXP decrease the optimum $D_p$ value according to Kriging prediction. This behavior demonstrates the interaction between the factors. The most trusted area in Kriging prediction is the area at which the infill process concentrated which is the area of small $D_{\mathrm{mt}}$ and $\mathrm{NXP}$ values (the area near the optimum values).

Figure 6.

Kriging nested tile plot for the top four factors based on converged sampling plan, in each tile the x-axis represents ${\theta }_D$ and the y-axis represent $D_p$ .

The optimum ejector

The optimum values for the nine ejector geometrical factors are summarized in Table 8. The optimum ejector is simulated to validate the Kriging result. The grid independency study for the optimum ejector design is performed to be sure that, the obtained results are grid independence. The GCI analysis results in Table 9 demonstrate that, the fine grid is sufficient to represent the ejector flow characteristics with a high resolution.

Table 8.

The optimum values of the ejector geometrical factors based on the final 216 observations.

Factor	${\theta }_{\mathrm{pc}}$	${\theta }_{\mathrm{pd}}$	$D_p/D_t$	$\mathrm{NXP}/D_t$	${\theta }_s$	$D_{\mathrm{mt}}/D_t$	${\theta }_D$	$L_{\mathrm{mt}}/D_t$	$L_t/D_t$
Optimum value	11.2071	7.0631	1.4069	1.9178	13.5	1.8045	5.5723	15.7391	0.6387

Table 9.

Grid convergence index calculations results for the optimum ejector design.

$A_s(\mathrm{m}{\mathrm{m}}^2)$	i	Ni	$A_{c_{\mathrm{avg}}}(\mathrm{m}{\mathrm{m}}^2)$	ωi (%)	${\mathrm{GCI}}_{i,i+1}^{\mathrm{fine}}(\mathrm{\%})$	α	R
2938.4	$0$ a			19.95		0.9975	0.9728
	1	501,900	0.0059	19.45
					3.2413
	2	385,800	0.0076	19.40
					3.5672
	3	301,525	0.0097	19.35

a

The value at zero grid space (f_exact equation (18)).

The optimum ejector design with the fine mesh size is simulated and the flow field is illustrated in Figure 7. Figure 7(a) and (b) represents Mach number and pressure contour, respectively. Figure 7(d) illustrates the pressure distribution along ejector center line, and Figure 7(d) demonstrates the sonic lines inside the ejector. As shown in Figure 7(a) and (b), the primary flow expanded in the primary nozzle and discharged in the mixing section before the mixing tube. The flow continually expands in the mixing section because its pressure is still higher than the pressure in the mixing section as revealed in Figure 7(c). After that the diamond waves are induced in the jet core. The high velocity difference between the primary and the induced flow forms the shear stress layer. The shear stress layer along with the mixing tube wall forms the converging duct for the secondary flow in which the induced flow accelerated to the sonic speed and choked since the sonic line is very close to the wall as shown in the sonic line plot in Figure 7(d). The sonic line plotted by connecting the center of the cells which have Mach number value of one.

Figure 7.

The flow characteristics of the optimum ejector design. (a) Mach number flow field of the optimum ejector design, (b) pressure flow field of the optimum ejector design, (c) static pressure distribution along the optimum ejector design center line, and (d) sonic lines plot of the optimum ejector design.

After the choking of the induced flow, strong shock has formed in the mixing tube and the second series of shock waves induced. At 0.04 m, after the beginning of the mixing tube, the second series of shock waves vanished and the static pressure gradually recovers to the back-pressure at the end of the mixing tube and the diffuser. Because of the choking of the induced flow, the entrainment ratio of the optimum ejector design is no longer affected by the reduction in the back-pressure of the ejector.

Conclusion

In this study, the effect of nine ejector geometrical factors on the entrainment ratio is investigated and the optimum ejector design is obtained at 12 MPa motive pressure, 2 MPa induced pressure, and 5.2 MPa discharge pressure using CFD technique, and surrogate-based optimization approach. The CFD model is validated against experimental data from the literature and the validation showed a good agreement with the experimental data with an average error of $0.6\mathrm{\%}$ in the critical mode. According to the CFD and Kriging results, the following conclusions can be drawn:

• The primary nozzle convergent angle and throat length have no effect on the entrainment ratio.

• Kriging model based on the final data set failed to predict the significance of the mixing tube length.

• Mixing tube diameter, primary nozzle exit position, and secondary nozzle incitation angle are the most effective factors on the entrainment ratio based on Kriging prediction. When all factors are optimum, small effect of the secondary nozzle incitation angle has been obtained.

• The optimum primary nozzle convergent and divergent angles are 11.2° and 7.06°, respectively. While the optimum primary nozzle exit diameter, throat length, and primary nozzle exit position and nozzle incitation angle are 6.47, 2.9, and 8.82 mm and 13.5°, respectively. In addition, the optimum mixing tube diameter and length are 8.3 and 72.4 mm, respectively; and the optimum diffuser angle is 5.57°. The optimum design achieves entrainment ratio of $19.45\mathrm{\%}$ .

• Surrogate-based optimization approach using Kriging as a surrogate model along with CFD represents as an effective way to optimize the ejector geometrical factors, and to investigate their effects on the entrainment ratio.