Simulation of Supersonic Flow in an Ejector Diffuser Using the JPVM

Abstract

The ejectors are used commonly to extract gases in the petroleum industry where it is not possible to use an electric bomb due the explosion risk because the gases are flammable. The steam ejector is important in creating and holding a vacuum system. The goal of this job is to develop an object oriented parallel numerical code to investigate the unsteady behavior of the supersonic flow in the ejector diffuser to have an efficient computational tool that allows modeling different diffuser designs. The first step is the construction of a proper transformation of the solution space to generate a computational regular space to apply an explicit scheme. The second step, consists in developing the numerical code with an-object-oriented parallel methodology. Finally, the results obtained about the flux are satisfactory compared with the physical sensors, and the parallel paradigm used not only reduces the computational time but also shows a better maintainability, reusability, and extensibility accuracy of the code.

1. Introduction

Steam jet ejectors offer a simple, reliable, low-cost way to produce vacuum. They are especially effective in the chemical industry where an on-site supply of the high-pressure motive gas is available. The ejector operation consists of a high-pressure motive gas that enters the steam chest at low velocity and expands through the converging-diverging nozzle. These results show a decrease in pressure and an increase in velocity. Meanwhile, the fluid enters at the suction inlet. The motive fluid, which is now at high velocity, enters and combines with the suction fluid [1, 2].

Ejector operations are classified as critical or noncritical flow mode; in the critical mode; the incoming flow is supersonic; other related experimental work and also work with CFD are in [3, 4]. This paper investigates the behavior of the supersonic flow for value of Mach number of 2.0 in the diffuser. The form of the ejector diffuser is not rectangular; therefore, it is necessary to transform the physical plane to a computational plane where the grid is rectangular.

The implementation was done using JPVM to run the program in a collection of computing system interconnected by one or more networks as a single logical computational resource to reduce the computational time. Figure 1 shows the ejector parts commonly used in the petroleum industry. Ejectors have no moving parts and operate by the action of high-pressure incoming stream like air and other vapors at a lower pressure into the moving stream and thereby removing them from the process system at intermediate pressure.

Figure 1 

Ejector parts.

2. The Governing Equation

We consider the two-dimensional compressible nonlinear Navier-Stokes equation written in the conservation form as

𝜕𝑈+𝜕𝑡𝜕𝐹+𝜕𝑥𝜕𝐺𝜕𝑦=0, (2.1) where 𝑈,𝐹, and 𝐺 are column vectors defined as,

𝑈=𝜌,𝜌𝑢,𝜌𝑣,𝜌𝑒,𝐹=𝜌𝑢,𝜌𝑢𝑣+𝑝−𝜏_𝑥𝑥,𝜌𝑢𝑣−𝜏_𝑥𝑦,𝑢𝜌𝑢(𝑒+²+𝑣²2)+𝑝𝑢−𝑢𝜏_𝑥𝑥−𝑣𝜏_𝑥𝑦+𝑞_𝑥,𝐺=𝜌𝑣,𝜌𝑢𝑣−𝜏_𝑥𝑦,𝜌𝑣²+𝑝−𝜏_𝑦𝑦,𝑢𝜌𝑣(𝑒+²+𝑣²2)+𝑝𝑣−𝑢𝜏_𝑥𝑦−𝑣𝜏_𝑦𝑦+𝑞_𝑦. (2.2) For a calorically perfect gas, it is possibly, to eliminate the energy 𝑒 in favor of 𝑢, 𝑣, 𝑝, and 𝜌 as follows: 𝑒=𝑢²+𝑣²/2+𝑝/(𝛾−1)𝜌. For clarity in calculations, the elements of the vector 𝑈 could be denoted as

𝑈₁𝑈=𝜌,₂𝑈=𝜌𝑢,₃𝑈=𝜌𝑣,₄=1𝑢𝛾−1𝑝+²+𝑣²2𝜌, (2.3) the elements of the column vector 𝐹 are denoted by

𝐹₁𝐹=𝜌𝑢,₂=𝜌𝑢²+𝑝−𝜏_𝑥𝑥,𝐹₃=𝜌𝑢𝑣−𝜏_𝑥𝑦,𝐹₄=𝛾𝑢𝛾−1𝑝𝑢+𝜌𝑢²+𝑣²2−𝑢𝜏_𝑥𝑥−𝑣𝜏_𝑥𝑦+𝑞_𝑥, (2.4) also, the elements of the column vector 𝐺 are denoted by

𝐺₁𝐺=𝜌𝑣,₂𝐺=𝜌𝑢𝑣,₃=𝜌𝑣²𝐺+𝑝,₄=𝛾𝑢𝛾−1𝑝𝑣+𝜌𝑣²+𝑣²2−𝑢𝜏_𝑥𝑦−𝑣𝜏_𝑦𝑦+𝑞_𝑦, (2.5) the viscous stress terms are written in terms of velocity gradients as

𝜏_𝑥𝑦=𝜇(𝜕𝑢+𝜕𝑦𝜕𝑣𝜏𝜕𝑥),_𝑥𝑥=𝜆(∇⋅𝐕)+2𝜇(𝜕𝑢𝜏𝜕𝑥),_𝑦𝑦=𝜆(∇⋅𝐕)+2𝜇(𝜕𝑣𝜕𝑦). (2.6) Likewise, the components of the heat flux vector (from Fourier's law of heat) are defined as

𝑞_𝑥=−𝑘𝜕𝑇,𝑞𝜕𝑥_𝑦=−𝑘𝜕𝑇.𝜕𝑦 (2.7)

The variations of the properties as dynamic viscosity and thermal conductivity are considered temperature-dependent and they are computed by Sutherland's law as

𝜇(𝑇)=𝜇₀(𝑇𝑇₀)^3/2𝑇₀+𝑆,𝑇+𝑆𝑘(𝑇)=𝜇(𝑇)𝛾𝑅.(𝛾−1)Pr (2.8)

3. The Numerical Scheme

An explicit predictor-corrector scheme for (2.1) is formulated as follows.

Predictor:

𝑈^∗_𝑖,𝑗=𝑈^𝑡_𝑖,𝑗−Δ𝑡Δ𝑥(𝐹^𝑡_𝑖+1,𝑗−𝐹^𝑡_𝑖,𝑗)−Δ𝑡Δ𝑦(𝐺^𝑡_𝑖+1,𝑗−𝐺^𝑡_𝑖,𝑗)+𝑆_𝑖,𝑗(𝑈_𝑖,𝑗), (3.1) where the artificial viscosity 𝑆 is given by

𝑆^𝑡_𝑖,𝑗=𝐶_𝑥|𝑝^𝑡_𝑖+1,𝑗−2𝑝^𝑡_𝑖,𝑗+𝑝^𝑡_{𝑖−1,𝑗}|𝑝^𝜏_𝑖+1,𝑗+2𝑝^𝜏_𝑖,𝑗+𝑝^𝜏_{𝑖−1,𝑗}(𝑈^𝑡_𝑖+1,𝑗−2𝑈^𝑡_𝑖,𝑗+𝑈^𝑡_{𝑖−1,𝑗})+𝐶_𝑦|𝑝^𝑡_𝑖+1,𝑗−2𝑝^𝑡_𝑖,𝑗+𝑝^𝑡_{𝑖−1,𝑗}|𝑝^𝑡_𝑖+1,𝑗+2𝑝^𝑡_𝑖,𝑗+𝑝^𝑡_{𝑖−1,𝑗}(𝑈^𝑡_𝑖,𝑗+1−2𝑈^𝑡_𝑖,𝑗+𝑈^𝑡_{𝑖,𝑗−1}), (3.2) where the 𝐶_𝑥 and 𝐶_𝑦 are two parameters; typical values of 𝐶_𝑥 and 𝐶_𝑦 range from 0.01 to 0.3, for this application 𝐶_𝑥 = 𝐶_𝑦 = 0.1.

Corrector:

𝑈^𝑘+1_𝑖,𝑗=12{𝑈^𝑡_𝑖,𝑗+𝑈^∗_𝑖,𝑗−Δ𝑡Δ𝑥(𝐹^∗_𝑖+1,𝑗−𝐹^∗_𝑖,𝑗)−Δ𝑡Δ𝑦(𝐺^∗_𝑖+1,𝑗−𝐺^∗_𝑖,𝑗)+𝑆_𝑖,𝑗(𝑈^∗_𝑖,𝑗)}. (3.3)

4. Validating the Numerical Code

Before managing the ejector's more complex geometry, it is possible to test the code with a simpler case with analytical solution in order to compare against the result of the numerical solution. The case is an expansion corner as sketched in Figure 2; for this problem an exact analytical solution exists in order to obtain a reasonable feeling for the accuracy of the numerical simulations result for the ejector shape; in this case in (2.1) the vector 𝐺=0 and the viscous stress terms are suppressed [5].

Figure 2 

Physical plane.

It is necessary to establish some details of the particular problem to be solved. The physical plane drawn in Figure 3 is considered. The flow at the upstream boundary is at Mach 2.0 with a pressure, density, and temperature equal to 1×10⁵ N/m², 1.23 Kg/m³ and 286 K, respectively. The supersonic flow is expanded at an angle of 5.352^∘; the reason for this choice is to have analytical solution. As shown in Figure 3, the calculations will be made in the domain from 𝑥=0 to 𝑥=66 in and from the wall from 𝑦=0 to 𝑦=40 in.

Figure 3 

Expansion corner initial parameters.

The location of the expansion corner is at 10 in, for this geometry; the variation of ℎ=ℎ(𝑥) is given by

ℎ(𝑥)=40,for𝑥≤10in,40+(𝑥−10)tan𝜃,for𝑥≥10in. (4.1)

We can construct a proper transformation as follows (Figure 4): let ℎ denote the local height from the lower to the upper boundary in the physical plane, clearly ℎ=ℎ(𝑥). Denote the 𝑦 location of the solid surface (the lower boundary in the physical plane) by 𝑦_𝑠, where 𝑦_𝑠=𝑦_𝑠(𝑥).

Figure 4 

Computational plane.

We define the transformation as

𝜉=𝑥,𝜂=𝑦−𝑦_𝑠(𝑥).ℎ(𝑥) (4.2)

With this transformation, in the computational plane 𝜉 varies from 0 to 𝐿 and 𝜂 varies from 0.0 to 1.0; 𝜂=0.0 corresponds to the surface in the physical plane, and 𝜂=1.0 corresponds to the upper boundary. The lines of constant 𝜉 and 𝜂 form a rectangular grid in the computational plane (Figure 4). The lines of constant 𝜉 and 𝜂 are also sketched in the physical plane; they form a rectangular grid upstream of the corner and a networks of divergent lines downstream of the corner.

The partial differential equations for the flow are numerically solved in the rectangular space and therefore must be appropriately transformed for use in the computational plane. That is, the governing equation must be transformed into terms dealing with 𝜉 and 𝜂. The derivative transformation is given by the following equations:

𝜕=𝜕𝜕𝑥(𝜕𝜉𝜕𝜉𝜕𝜕𝑥)+(𝜕𝜂𝜕𝜂𝜕𝜕𝑥),=𝜕𝜕𝑦(𝜕𝜉𝜕𝜉𝜕𝜕𝑦)+(𝜕𝜂𝜕𝜂𝜕𝑦). (4.3)

The metrics (𝜕𝜉/𝜕𝑥), (𝜕𝜂/𝜕𝑥), (𝜕𝜉/𝜕𝑦), and (𝜕𝜂/𝜕𝑦), in (4.3), are obtained from the transformation given by (4.2), that is, 𝜕𝜉/𝜕𝑥=1 and 𝜕𝜉/𝜕𝑦=0:

𝜕𝜂=𝜕𝑥𝜕[𝑦−𝑦_𝑠(𝑥)/ℎ(𝑥)]=𝜕𝜕𝑥[𝑦𝜕𝑥𝜕ℎ(𝑥)]−[𝑦𝜕𝑥^𝑠(𝑥)]ℎ(𝑥)=−𝑦ℎ^′(𝑥)[ℎ(𝑥)]²𝑦−[^′𝑠(𝑥)ℎ(𝑥)−𝑦_𝑠(𝑥)ℎ^′(𝑥)[ℎ(𝑥)]²]=ℎ^′(𝑥)[𝑦_𝑠(𝑥)−𝑦][ℎ(𝑥)]²−𝑦^′_𝑠(𝑥)ℎ(𝑥)[ℎ(𝑥)]²; (4.4) if 𝑦=ℎ(𝑥)𝜂+𝑦_𝑠(𝑥), then

𝜕𝜂𝜕𝑥=ℎ^′(𝑥)[𝑦_𝑠(𝑥)−[ℎ(𝑥)𝜂+𝑦_𝑠(𝑥)]][ℎ(𝑥)]²−𝑦^′_𝑠(𝑥)=ℎℎ(𝑥)_′(𝑥)[−ℎ(𝑥)𝜂][ℎ(𝑥)]²−𝑦^′_𝑠(𝑥)𝜂ℎ(𝑥)=−ℎℎ(𝑥)^′1(𝑥)−𝑦ℎ(𝑥)^′𝑠(𝑥); (4.5) if 𝑥<10, then 𝑦^′_𝑠(𝑥)=0 and ℎ^′(𝑥); if 𝑥≥10=𝑦^′_𝑠(𝑥)−tan𝜃 and ℎ^′(𝑥)=tan𝜃, so the metrics could be finally expressed as 𝜕𝜉/𝜕𝑥=1,𝜕𝜉/𝜕𝑦=0: 𝜕𝜂=𝜕𝑥0,if𝑥<10in,(1−𝜂)tan𝜃ℎ(𝑥),if𝑥>10in,𝜕𝜂=1𝜕𝑦.ℎ(𝑥) (4.6)

The mesh for the expansion corner and for the ejector is automatically generated to avoid stability numerical problems using the CFL criterion where the values of Δ𝜉_𝑖 are given by

Δ𝜉_𝑖=𝐶Δ𝑦|tan𝜃±𝜇|_max, (4.7) where 𝜃=𝑎sin(1/𝑀_𝑖,𝑗) and 𝜇=𝑎tan(𝑣_𝑖,𝑗/𝑢_𝑖,𝑗).

The mesh obtained is shown in Figure 5. The size of the ejector grid is 79×41 discrete points.

Figure 5 

Mesh generated for the expansion corner.

In the results of the Mach number, pressure and density for the expansion corner problems are shown in the results. The analytical solution of the Mach Number for the expansion corner could be obtained with the Prandtl-Meyer function (4.8):

𝑓(𝑀)=𝛾+1𝛾−1tan⁻¹𝛾−1𝛾+1(𝑀²−1)−tan⁻¹𝑀²−1. (4.8)

The analytical solution of the expansion corner is shown in Figure 6.

Figure 6 

Analytical solution of the expansion corner.

The leading edge of the expansion fan makes an angle 𝜇₁ with respect to the upstream flow direction, and the trailing edge of the wave makes an angle 𝜇₂ with respect to the downstream flow direction. The angles 𝜇₁ and 𝜇₂ are mach angles, defined as 𝜇₁=sin⁻¹(1/𝑀₁) and 𝜇₂=sin⁻¹(1/𝑀₂).

The numerical solution of the expansion corner is depicted in Figure 7.

Figure 7 

Numerical solution of the expansion corner.

The numerical solution matches the analytical with an error of less than 1% percent at all points and the expansion Mach fan is well formed in the numerical solution.

5. The Transformation of the Physical Ejector Diffuser

The appropriate transformation to generate a boundary-fitted coordinate system for the ejector diffuser is defined as follows:

𝜉=𝑥,𝜂=𝑦−𝑦_𝑠(𝑥)𝑦_𝑧(𝑥)−𝑦_𝑠(𝑥), (5.1)(5.2) where 𝑦_𝑠(𝑥) and 𝑦_𝑧(𝑥) are defined as

𝑦_𝑠8(𝑥)=0.0,for𝑥≤10,−8167.875(𝑥−10),for10<𝑥≤177.875,8.0,for177.875<𝑥≤299.125,𝑦89(𝑥−299.125)+8,for𝑥>299.125,_𝑧−8(𝑥)=42.0,for𝑥≤10,8167.875(𝑥−10)+42,for10<𝑥≤177.875,34.0,for177.875<𝑥≤299.125,89(𝑥−299.125)+34,for𝑥>299.125. (5.3)

The metrics in (4.3) given by (5.1) and (5.2) are calculated as follows: 𝜕𝜉/𝜕𝑥=1, 𝜕𝜉/𝜕𝑦=0, 𝜕𝜂/𝜕𝑦=1/(𝑦_𝑧(𝑥)−𝑦_𝑠(𝑥)), by (5.2), 𝑦=𝜂[𝑦_𝑧(𝑥)−𝑦_𝑠(𝑥)]+𝑦_𝑠(𝑥), so 𝜕𝜂/𝜕𝑥=(𝜂[𝑦^′_𝑠(𝑥)−𝑦^′_𝑧(𝑥)]−𝑦^′_𝑠(𝑥))/(𝑦_𝑧(𝑥)−𝑦_𝑠(𝑥)). The metric 𝜕𝜂/𝜕𝑥 can be obtained differentiating 𝑦_𝑠(𝑥) and 𝑦_𝑧(𝑥):

𝑦^′_𝑠8(𝑥)=0,for𝑥≤10,−8167.875,for10<𝑥≤177.875,0,for177.875<𝑥≤299.125,𝑦89,for𝑥>299.125,_′𝑧−8(𝑥)=0,for𝑥≤10,8167.875,for10<𝑥≤177.875,0,for177.875<𝑥≤299.125,89,for𝑥>299.125. (5.4)

Therefore

𝜕𝜂=−𝜕𝑥0.0,for𝑥≤10,9.53×10⁻²𝜂−4.76×10⁻²−42.95+9.53×10⁻²𝑥−,for10<𝑥≤177.875,0.0,for177.87<𝑥≤299.125,−0.17𝜂+8.98×10⁻²27.77−0.17𝑥,for𝑥>299.125. (5.5) In Figures 8 and 9 is depicted the distribution of the metrics 𝜕𝜂/𝜕𝑥 and 𝜕𝜂/𝜕𝑦, respectively. The size of the ejector grid is 1101×41 discrete points and it was generated using the CFL marching criterium.

Figure 8 

Distribution of 𝜕𝜂/𝜕𝑥.

Figure 9 

Distribution of 𝜕𝜂/𝜕𝑦.

The boundary condition in the outlet is Neumann and in the walls is applied the nonslipping condition. The flux field is initialized with a pressure, density, and temperature equal to 1×10⁵ N/m², 1.23 Kg/m³, and 286 K, respectively; the velocity field is initialized with 𝑢=𝑣=0.0 m/s. In the upstream boundary a velocity of Mach 2.0 (𝑢=678 m/s and𝑣=0.0 m/s) is injected. The simulation is performed for 5,000,000 time steps, with a Δ𝑡 calculated using the CFL criterium (see Figure 10).

Figure 10 

Ejector plane with the initial parameters.

5.1. Viscous Stress Terms

Of course is necessary to use the metrics transformation for the viscous stress terms. Once applied the metrics the terms are expanded as follows:

𝜏_𝑥𝑦=𝜇(𝜕𝑢(𝜕𝜂𝜕𝜂𝜕𝑦)+𝜕𝑣+𝜕𝜉𝜕𝑣(𝜕𝜉𝜕𝜂𝜏𝜕𝑥)),_𝑥𝑥=𝜆(𝜕𝑢+𝜕𝜉𝜕𝑢(𝜕𝜂𝜕𝜂𝜕𝑥)+𝜕𝑣(𝜕𝜂𝜕𝜂𝜕𝑦))+2𝜇(𝜕𝑢+𝜕𝜉𝜕𝑢(𝜕𝜂𝜕𝜂𝜏𝜕𝑥)),_𝑦𝑦=𝜆(𝜕𝑢+𝜕𝜉𝜕𝑢(𝜕𝜂𝜕𝜂𝜕𝑥)+𝜕𝑣(𝜕𝜂𝜕𝜂𝜕𝑦))+2𝜇(𝜕𝑣(𝜕𝜂𝜕𝜂𝜕𝑦)). (5.6)

For the heat flux vector the transformation is

𝑞_𝑥=−𝑘(𝜕𝑇+𝜕𝜉𝜕𝑇(𝜕𝜉𝜕𝜂𝑞𝜕𝑥)),_𝑦=−𝑘(𝜕𝑇𝜕𝜂+(𝜕𝜂𝜕𝑦)). (5.7) for this application the values of 𝑇₀=273𝐾,𝑆=110.5,𝜇₀=1.68×10⁻⁵,𝛾=1.4, and Pr=0.71 are used.

6. Parallelizing the Numerical Scheme

The JPVM (Java Parallel Virtual Machine) library is a software system for explicit message passing-based distributed memory MIMD parallel programming in Java. The library supports an interface similar to the C and FORTRAN interfaces provided by the Parallel Virtual Machine (PVM) system, but with syntax and semantics enhancements afforded by Java and better matched to Java programming styles. The similarity between JPVM and the widely used PVM system supports a quick learning curve for experienced PVM programmers, thus making the JPVM system an accessible, low-investment target for migrating parallel applications to the Java platform. At the same time, JPVM offers novel features not found in standard PVM such as thread safety, multiple communication end-points per task, and default-case direct message routing. JPVM is implemented entirely in Java and is thus highly portable among platforms supporting some version of the Java Virtual Machine. This feature opens up the possibility of utilizing resources commonly excluded from network parallel computing systems such as Mac-, Windows- and Linux-based systems [6, 7].

The method used is an explicit finite-difference technique which is second-order accurate in both space and time with artificial viscosity.

In the predictor step the governing equation is written in terms of forward differences (predictor), later on is written in backward differences (corrector). Following the scheme we can divide the computation into tasks that perform over a part of the solution space only in the axis 𝜂; so if we have 𝑚 discrete points that divide the space over 𝜂, then the computation could be performed into 𝑘 tasks, where the number of points in which each task operates is 𝑚/𝑘; if the division is not exact, then the number of remaining points will be distributed among the tasks. For example, if we have 40 points on 𝜂, and we want to divide it into 6 tasks, then there will be 4 tasks that take charge of working on 7 points and two tasks on 6 points. If the calculation is wanted to divide among 4 tasks, then each task will work on 10 points.

If 𝑚=40 and 𝑘=4, then the iterators 𝑗₁,𝑗₂,𝑗₃,𝑗₄ of each task will work on the following points: 𝑗₁=1,…,10, 𝑗₂=11,…,20, 𝑗₃=21,…,30, 𝑗₄=31,…,40.

Figure 11 shows the set-up of four tasks working over the same solution space, the arrows that arise of the tasks indicate the messages, and the arrows ahead of each task indicate the direction of the calculation.

Figure 11 

Fine grained.

The main issue in a traditional parallel algorithm for finite differences problem, is the intensive message passing, obtaining a poor performance in distributed memory systems, where the latency time is higher. Suppose that we want to advance on 300 discrete points over the axis 𝜉 and we have 4 tasks; so in every iteration it is necessary to send 6 messages; therefore it will be 300×6=1800 messages that holds just one value, and this is just for the case of the variable 𝐹₁. Additionally, also we have to create 3 tasks for the other flux variables 𝐹₂,𝐹₃,𝐹₄, that give a total of 1800×4=7200 messages.

To reduce the quantity of messages it is necessary to increase the tasks granularity, assigning a tasks for flux variable. Figure 12 depicted four tasks working in their own flux variable.

(a)

(b)

(c)

(d)

Figure 12 

Coarse grained.

Creating four tasks by every term of the flux 𝐹₁, 𝐹₂, 𝐹₃, and 𝐹₄, the task 𝑞(𝑞=1,2,3,4) carries out the necessary calculations to obtain the new value of 𝐹_𝑞 at 𝑖+1. The procedure is listed as follows.

(1)Calculate the forward differences of the predictor step.(2)Calculate the artificial viscosity of the predictor step.(3)Calculate the predicted value of 𝐹_𝑞.(4)Calculate the rearward differences of the corrector step, using the predicted values of 𝐹_𝑞.(5)Calculate the artificial viscosity for corrector step.(6)Calculate the average derivative.(7)Calculate the value of the flux variable 𝐹_𝑞 at 𝑖+1.

Additionally to the four tasks (slave tasks), we need the master task to control the execution of the slave tasks and achieve the following calculations.

(1)Calculate Δ𝜉 at 𝑖 to obtain the predicted values of the flux variables 𝐹_𝑞.(2)Apply the boundary conditions.(3)Adjust the flux variables 𝐹 at boundary.

The former procedure indicates that the work of every task has to perform over the points collection (𝑖,𝑗) that conform the bidimensional mesh. Every task 𝑞(𝑞=1,2,3,4) will calculate the new values of the flux variables 𝐹_𝑞 at localizations Δ𝜉₁,Δ𝜉₁+Δ𝜉₂,…,∑^𝑛_𝑖=1Δ𝜉_𝑖.

We follow a similar reasoning if we want to increase the granularity of tasks.

Figure 13 depicted the setup of four tasks working over its own flow field, in this case every task 𝑞(𝑞=1,2,3,4) carried out the calculations to find the new values at each grid point (𝑖,𝑗) of the variables 𝑈_𝑞 in the time 𝜏₀+1, 𝜏₀+2, …, 𝜏₀+𝑘.

(a)

(b)

(c)

(d)

Figure 13 

Coarse-grained 3D.

The advantage of this paradigm for the parallelism is that we can made quickly hybridizations between a coarse and fine granularity, and the parallel tasks are created recursively easily.

7. Simulation Results

In this section we show the results of the simulations about the expansion corner and the ejector.

7.1. Expansion Corner Results

The expansion corner results for the Mach number, density, and temperature are shown in Figures 14 and 15 in steady state.

Figure 14 

Mach number in the expansion corner.

Figure 15 

Density in the expansion corner.

7.2. Ejector Results

Figures 16, 18, 20, and 22 show the contour graphs of the density and the Figures 17, 19, 21, and 23 of the Mach number, after 0.2 s, 0.9 s, 1.6 s, and 5.0 s of real simulation when the flux is stable.

Figure 16 

Density through the diffuser section at 0.2 s

Figure 17 

Mach number through the diffuser section at 0.2 s.

Figure 18 

Density through the diffuser section at 0.9 s.

Figure 19 

Mach number through the diffuser section at 0.9 s.

Figure 20 

Density through the diffuser section at 1.65 s.

Figure 21 

Mach number through the diffuser section at 1.65 s.

Figure 22 

Density through the diffuser section at 5.0 s.

Figure 23 

Mach number through the diffuser section at 5.0 s.

In Figure 24, is depicted the profile of the Mach number in the different sections of the diffuser; when the compression, transfer; and expansion occur, three numerical visors are taken at 5 inches, 21 inches, and 37 inches in the vertical, and compared with the experimental processed results when the flux is stable [1].

Figure 24 

Mach profile of fluid through the diffuser section when the flux is stable.

8. Performance Results

To gain a better perspective of the performance options, the JPVM was compiled using a native compiler gcj (gnu java compiler) to avoid the overload of the virtual machine for Linux and Windows to generate a native version and carry out comparatives with the byte-code version. The creation and communication is under the same machine; a slow Pentium IV processor was used because is easier to measure the time. Table 1 shows the creation time of 1, 2, 4, 8, and 16 tasks.

Table 1 

Required memory.

In Table 1 we can see that the time of creation is practically the same in both versions.

In Table 2 we can observe that there is not significative difference between the native version and the byte-code version in the communication time of the tasks. Nevertheless the advantage of the native version is the saving memory storage. About the execution time, gcj can produce programs that are faster than the byte code version around 30%.

Table 2 

Communication time.

The cluster configuration used to execute the simulation is

(i)14 nodes, every node with two-processor Xeon DP 2.8 Ghz, 4 GB RAM, in total 28 cores.

For this application is supposed 10% of the serial part, the part that could not be parallelized (when the primitive variables are decoded); so the max speed-up for 16 tasks is 14.5; Table 3 shows the performance metrics for 16 tasks.

Table 3 

Performance Metrics using JPVM over the cluster.

9. Conclusions

A numerical parallel code has been developed to simulate the supersonic flow in the ejector diffuser in a parallel-distributed system with an object-oriented methodology. The model was validated with cases where there is an exact solution. The design of a parallel program allows to reduce the execution time; nevertheless the design and build of a parallel program is complex due to the no deterministic execution and in a scientific computing ambient is more difficult build and debug because it is necessary to manage several data, and so it is necessary to use a good methodology to minimize the risk of bugs in the program construction.

We have discussed implementations of object-oriented design using Java in computational fluid dynamics simulations. This also provides the benefits of better maintainability, reusability, and extensibility of the code. For examples, is possible to create recursively new parallel tasks to control easily the granularity; for this reason Java is a serious language suitable for demanding applications in science and engineering and the JPVM is a good MPI-based tool for parallel cluster computing.

Finally, this program is a good tool to investigate the behavior of the flux in the ejector diffuser and shows how to transform the solution space in an easy way.

Nomenclature

𝑢:	Component velocity in the 𝑥 direction (m/s)
𝑣:	Component velocity in the 𝑦 direction (m/s)
𝑇:	Temperature (K)
𝑃:	Pressure (N/m3)
𝜌:	Density (Kg/m3)
𝑒:	Energy
𝑀:	Mach number
𝑐𝑝:	Specific heat capacity at constant pressure ((J/(kg K)))
𝛾:	Ratio of heat capacities
Pr:	Prandtl number
𝑞𝑖:	Heat flux along 𝑥𝑖 direction (W/m2)
𝑅:	Ideal gas constant ((J/(kg K)))
𝑆:	Sutherland’s constant (K)
𝑘:	Thermal diffusivity (W/(m K))
𝑣:	Kinematic viscosity (m2/s).

Acknowledgment

This research was supported by the Mexican Petroleum Institute by Grant 204005.