Issue 
Mechanics & Industry
Volume 19, Number 4, 2018



Article Number  403  
Number of page(s)  11  
DOI  https://doi.org/10.1051/meca/2018003  
Published online  12 November 2018 
Regular Article
Rapid design method and new contours for a class of three dimensional supersonic nozzle of arbitrary exit cross section
Institute of Aeronautics and Space Studies, University of Blida 1,
BP 270,
Blida
09000, Algeria
^{*} email: z_toufik270169@yahoo.fr
Received:
2
June
2017
Accepted:
12
December
2017
The aim of this work is to develop a new and rapid numerical method for designing a new contour of the supersonic nozzle with arbitrary exit cross sections as a class of three dimensional nozzle extracted from the calculation of stationary flow solutions in axisymmetric nozzle. The application is made for nozzle giving a uniform and parallel flow at the exit section. The determination of the points of the axisymmetric nozzle contours in a nondimensional manner with respect to the arbitrary throat radius is necessary. The exit section of the nozzle must be discretized at several points. The radii of the points of the throat section are determined by equalizing the ratio of the axisymmetric critical sections corresponding to each selected point of the exit section. Each point passes a contour of the nozzle to the throat where their position is determined by the multiplication of the nondimensional positions of the axisymmetric nozzle point by the throat radius of this contour. A uniform portion is added at the end of each contour to compensate the decrease in its length. Longitudinal discretization of the nozzle is necessary. The flow properties of each point are the same as those of the points of the axisymmetric nozzle. The temperature and the deviation of the flow at each point of any section are determined by their interpolations between two successive points of each contour. The Mach number, the pressure and the density are determined accordingly. The application is made for air at high temperature lower than the dissociation threshold of the molecules.
Key words: Supersonic flow / supersonic axisymmetric minimim length nozzle / high temperature / calorically imperfect gas / method of characteristics / specific heat at constant pressure / error
© AFM, EDP Sciences 2018
1 Introduction
Numerical computing obliged us to think of developing efficient methods that meet the demands of convergence with high precision and in a smallest possible computer execution time. Generally the execution time for design problems is quite high despite the power of calculator used since there is repetition of calculation. The methods developed for determining the contours of the three dimensional nozzles are based on the determination of the internal stream lines using the definition of the stream function [1–6]. As the flow is considered as perfect and nonviscous, the stream lines can represent rigid surfaces and hence a wall shape of obstacles [7]. In our case these stream lines can represent the wall of the nozzle [1–6]. The use of these methods requires a very high computational time since the search of the points of the stream lines is based on the use of the segments in upward and downward characteristics from the developed internal mesh in the axisymmetric nozzle, which becomes impossible for fine meshes, in a time when a fine mesh give a high precision. Theoretically the fine meshes gives a good accuracy for the design, but in parallel, it falls failing for the computer design of the three dimensional nozzles based on the stream lines search.
The segment number of the internal mesh characteristics of the axisymmetric nozzle depends on many of the descending characteristics in the Kernel region and the step used in the transition region. Generally the number of points found on the wall is much greater than the number of descending characteristics in the Kernel region. We are interested in axisymmetric nozzles giving a uniform and parallel flow to the exit section [8–11]. The number of information on the flow parameters (x, y, P, T, ρ, M, θ and ψ) at each internal mesh point becomes very large, which requires a very large storage reservation in the computer memory, and which gives in parallel a very important calculation time for the design of the three dimensional and axisymmetric nozzles. For a small number of N_{S} and N_{P}, the design becomes possible, but with moderate precision. This gives that the methods developed in references [1–6] are limited by a certain maximum of discretization. It is not necessary, in the general case, to know the properties of the flow inside the nozzle. Only the axisymmetric nozzle contour is desired. But if the design is made as a passage for the design of the three dimensional nozzle, it is necessary to determine the internal mesh information according to references [1–4].
According to [1], the design is made only for the PG model case, calorically and thermally perfect from the design of the axisymmetric nozzle using the internal mesh for the search of the stream lines points.
The reference [2] deals with the design of the three dimensional nozzles from a given axisymmetric shape taken as a seconddegree polynomial for the use of the PG model, using the definition of a stream line always in an internal mesh. References [3–4] deal with the problem presented in reference [2] but only for a square section. In addition a truncation is presented in the reference [4].
In reference [5], the author made the design of the nozzles 3D of arbitrary exit section by the method of the characteristics in the case of PG model.
In the reference [6], the authors made the design of the nozzles 3D from the axisymmetric nozzle in the context of the HT model using the axisymmetric nozzle design and the search of stream lines from the internal mesh.
In references [1–5], the PG model is used for the design. Not to mention the effectiveness of the method, this model gives good results if the M_{E} does not exceed about 2.00, and T_{0} does not exceed the 1000 K. approximately [6,8,11,12]. The HT model is used for the same design problem in view of the presentation of the corrections to the results obtained in the references [1–5] in the case where M_{E} exceeds 2.00 and T_{0} exceeds 1000 K. By way of information, the HT model presents a generalization of the PG model.
References [13–16] have studied and analyzed analytically and numerically the supersonic flow by determining the parameters of flow through the Laval nozzle, as well as CFD validation. The study is done for the twodimensional nozzle using PG model already discussed. In reference [16], the authors studied the outside adaptation of the nozzle with consideration of the shock birth and a separation of the flow.
Regarding the design of the nozzle 3D according to references [1–6], for N_{T} and N_{L} points chosen, one needs to treat N_{T} × N_{L} × N_{S} segments. It is true that a technique is normally opted not to go through the N_{S} segments for even point of the nozzle 3D. But as even, the number of segments to be processed is very large. It can be said that the method developed in references [1–6] is practically limited by a maximum of discretization of the axisymmetric nozzle, at a time when the high discretization is preferred in order to obtain very good results.
The aim of this work is to solve the disadvantages found in the methods used for the design of the nozzles 3D presented in references [1–6], by developing a new very fast and efficient numerical method that can design the new contours of the nozzles 3D of arbitrary exit section in the HT model as a special class of Laval nozzle, only from the boundary axisymmetric nozzle wall, without going through the use of the stream function and storage of the internal meshing information, to derive the properties of the flow at all points of the wall of the nozzle 3D. The first problem is to design the axisymmetric nozzle in dimensionless values with respect to the throat radius, that is, determine the positions of the wall in terms of (x/y_{*} and y/y_{*}). The axisymmetric MLN is chosen in our work seen its practical and current use in aerospace applications of missiles and satellites launchers because it has better performances compared to other types of supersonic nozzles in addition to its simplicity of construction [8,11–17]. The application is made for the case of air with HT, lower than the dissociation threshold of the molecules. The flow is considered as adiabatic, irrotational and stationary.
2 Axisymmetric MLN design contour
We are interested in the MLN giving a uniform and parallel flow at the exit section due to elimination of the thrust losses caused by the inclination of the wall and to have a complete expansion through the nozzle. The general shape of the MLN as well as the various zones of the flow are shown in Figure 1. The contour search (x, y) of the nozzle wall as well as the flow properties (M, T, θ, P and ρ) at all flow internal points and the boundary points are done simultaneously and by the use of the HT MOC [8]. The application is made for air, lower than the dissociation threshold of the molecules.
The design results of the axisymmetric nozzle depend on the mesh quality used, which is itself dependent on two important parameters which are Δθ and Δ characterizing respectively the refinement and the mesh quality in the Kernel and transition zones. The internal mesh therefore contains N_{S} total segments of the descending rising characteristics and N_{P} internal points. The wall of the nozzle then contains N_{F} points. The values of N_{F}, N_{S}, N_{P} and ε depend essentially on Δθ, Δ, M_{E} and T_{0}.
The equations for axisymmetric supersonic non viscous flow of perfect gas with HT by the MOC can be summarized by: $$\{\begin{array}{l}\frac{{C}_{P}(T)}{2\text{}H(T)}\text{}\sqrt{{M}^{2}(T)1\text{}}dT+d\theta =\frac{\mathrm{sin}\theta \text{}\mathrm{sin}\mu}{y\text{}\mathrm{cos}(\theta \mu )}\text{}dx\\ \frac{\mathrm{d}y}{\mathrm{d}x}=tg(\theta \mu )\end{array}\text{,}$$(1) $$\{\begin{array}{l}\frac{{C}_{P}(T)}{2\text{}H(T)}\text{}\sqrt{{M}^{2}(T)1\text{}}dTd\theta =\frac{\mathrm{sin}\theta \text{}\mathrm{sin}\mu}{y\text{}\mathrm{sin}(\theta +\mu )}\text{}dx\\ \frac{\mathrm{d}y}{\mathrm{d}x}=tg(\theta +\mu )\end{array}\text{.}$$(2) With $$H\left(T\right)={\int}_{T}^{{T}_{0}}{C}_{P}\left(T\right)\text{}dT\text{.}$$(3)
The M and T at HT are connected by [18,19]: $${M}^{2}\left(T\right)=\frac{2H\left(T\right)}{\gamma \left(T\right)\text{}R\text{}T}\text{.}$$(4) With $$\gamma \left(T\right)=\frac{{C}_{P}\left(T\right)}{{C}_{P}\left(T\right)R}\text{,}$$(5) $$\mu =\mathrm{arcsin}\left(\frac{1}{M}\right)\text{.}$$(6)
The ratios ρ/ρ_{0} and P/P_{0} are calculated by [18,19]: $$\frac{\rho}{{\rho}_{0}}=\text{exp}\left({\int}_{T}^{{T}_{0}}\frac{{C}_{P}\left(T\right)}{\gamma \left(T\right)\text{}R\text{}T}\hspace{0.17em}dT\right)\text{,}$$(7) $$\frac{P}{{P}_{0}}=\frac{T}{{T}_{0}}\frac{\rho}{{\rho}_{0}}\text{.}$$(8)
As the flow through the throat and the exit section is unidirectional, the theoretical ratio of critical sections at HT, given by equation (9) remains valid for the convergence of the numerical found results [18,19]: $${\left(\frac{{y}_{E}}{{y}_{*}}\right)}_{Exact}^{2}=\frac{{A}_{E}}{{A}_{*}}=\text{exp}\left({\int}_{{T}_{E}}^{{T}_{*}}\left[\frac{{C}_{P}\left(T\right)}{\gamma \left(T\right)\text{}R\text{}T}\frac{{C}_{P}\left(T\right)}{2\text{}H\left(T\right)}\right]\hspace{0.17em}dT\right)\text{.}$$(9) To validate the numerical results of design, one can calculate the relative error committed the ratio of the critical sections by the following relation: $${\u03f5}_{({y}_{E}/{y}_{*})}\text{}(\%)\text{}=\text{}\left\text{}1\text{}\text{}\frac{{({y}_{E}/{y}_{*})}_{Computed}}{{({y}_{E}/{y}_{*})}_{Exact}}\text{}\right\text{}\times \text{}100\text{.}$$(10)The value of (y_{E}/y_{*})_{Computed} illustrated in relation (10) represents the dimensionless exit radius determined according to the chosen problem discretization.
For air R = 287.1029 J/(kg K). For the PG model, the value of γ = 1.402. The range of the variation of M_{E} is up to 5.00 and for T_{0} is up to 3500 K according to the references [8,11,18–20] so as not to have dissociation of the molecules.
Fig. 1 Presentation of the axisymmétric MLN flow field. 
3 Presentation of the developed method
A new very rapid and effective method is proposed for designing the novel contours of the nozzle 3D from the shape of the axisymmetric nozzle without going through the search of the stream lines and the properties of the flow in the internal mesh of the axisymmetric nozzle. This time, we are only interested in the contour properties of the axisymmetric nozzle wall. Then during the design of the axisymmetric nozzle, the properties of the internal flow mesh are omitted and not stored, which make in this step the fairly quick calculation process which is a reported advantage in this first step. During the design of the axisymmetric nozzle, the radius of the throat is chosen in the nondimensional form. The coordinates of the contour points of the nozzle wall are calculated by (x/y_{*}, y/y_{*}). Consequently, the length and the exit radius of the nozzle are normalized again by y_{*}. The values L/y_{*} and y_{E}/y_{*} respectively are obtained. The choice of the exit section shape of the nozzle 3D is arbitrary, which must be discretized at several points and placed in the exit section of the axisymmetric nozzle. The radii y_{C} of the points of the critical section, which the latter has the same shape as the exit section, are determined by equalizing the ratio of the corresponding axisymmetric critical sections to each selected point of the exit section. Each point of the exit section passes a contour of the nozzle 3D to the throat.
For a given radius of the throat 0 < y_{C }< y_{*} of the axisymmetric nozzle, it is possible to determine the position of the forming points of a contour which can represent a wall of an intermediate nozzle of the nozzle 3D whose position of these points is determined by multiplying the nondimensional positions of the axisymmetric nozzle by the throat radius of this contour. Then the position of the new point Q, see Figure 2, is determined from the position of point P by the following relations :
${x}_{Q}=\left(\frac{{x}_{P}}{{y}_{{}^{*}}}\right)\text{}{y}_{C}$,${y}_{Q}=\left(\frac{{y}_{P}}{{y}_{{}^{*}}}\right)\text{}{y}_{C}\text{.}$(11)
The flow properties (M, P/P_{0}, T/T_{0}, ρ/ρ_{0} and θ) at this new point Q of this contour are the same as those of the point P of the axisymmetric nozzle.
By way of information, the reference mark chosen on the section of the nozzle is named by Oyz. Then the axis of the nozzle is the Ox.
The intermediate contour which can also represent the wall of the nozzle 3D is included in the main axisymmetric shape as shown in Figure 2. The longitudinal length of this contour is less than the main length of the axisymmetric nozzle. It is given by: $${L}_{2}\left({y}_{C}\right)=\left(\frac{{L}_{1}}{{y}_{{}^{*}}}\right)\text{}{y}_{C}\text{.}$$(12) The value given by relation (12) represents the length of the intermediate contour between points C and D of Figure 2. A uniform portion is added at the end of each contour to compensate the decrease in its length.
By varying the value of y_{C}, we can find other contours we will use to form the nozzle 3D.
In the first place, it is necessary to arbitrarily choose the shape of the exit section of the nozzle 3D which is to be placed in the exit section of the axisymmetric nozzle as shown in Figure 3. The discretization of this shape in an N_{T} number of points and the (y, z, r and β) as well as (M, T/T_{0}, P/P_{0}, ρ/ρ_{0} and θ) in all these points is necessary. The flow properties are the same for all points of the exit section of the axisymmetric nozzle and are equal to the values corresponding to M = M_{E}.
For our method it is not necessary to calculate the value of the stream function at each point as it is used in references [1,6], so a slight reduction in calculation time for our method, given the large number of mathematical operations won.
A particle which runs through the flow field on this contour from point C will arrive with an exit Mach number M_{E} at the point D. The portion of this contour between the point D and the point E of Figure 2 is added since the contour ends at point E because the exit section of the nozzle 3D is placed exactly at the exit section of the axisymmetric nozzle. This contour portion is horizontal since M_{E} is reached at point C. The flow properties on this segment are the same and they are equal to those of the exit section of the nozzle.
So we are only interested in the physical properties at the points of this contour. The problem therefore becomes the determination of the ordinate y_{C} < y_{*} which will identify this contour and which verifies the equality of the ratios of the axisymmetric critical sections between the axisymmetric main nozzle and that of the contour which forms the nozzle 3D. Since the contour passes through the Q points of the exit section as shown in Figure 3, the equality of the ratios of the axisymmetric sections gives the following result: $${y}_{C}^{2}=\frac{{z}_{Q}^{2}+{y}_{Q}^{2}}{\left(\frac{{A}_{E}}{{A}_{*}}\right)}\text{.}$$(13) For each point Q of the exit section of Figure 3, there is a contour which passes through this point and passes through the longitudinal length of the nozzle to the throat and cuts the throat at a point having a radius y_{C} < y_{*} found by the relation (13). The number of points found on each intermediate contour is the same and it is equal to that of the axisymmetric main nozzle plus 1, since point E has been added according to Figure 2. Each contour forms a plane of the nozzle inclined by an angle β as shown in Figure 3. Then all the points of the intermediate contour have the same value of the angle β.
Fig. 2 Intermediates contour of the axisymmetric nozzle for nozzle 3D. 
Fig. 3 Transverse discretization of the exit section of the nozzle 3D. 
4 Flow properties in a vertical section
The points found on each contour are not placed exactly on straight sections. It is highly recommended to determine the shape of the nozzle in each section by determining the positions (y, z) and the properties (M, P/P_{0}, T/T_{0}, ρ/ρ_{0} and θ) on each point of a section.
The discretization of the longitudinal direction of the nozzle in N_{L} section must be presented as illustrated in Figure 4. In this Figure, only the abscissa x_{Q} of the N_{T} points of each section is known.
The point Q to be searched for the nozzle 3D is located on a limited segment between two successive points of corresponding contour of the nozzle 3D according to Figure 4. The segment containing the point Q of the nozzle 3D wall must verify the following condition: $$({x}_{L}{x}_{Q})\text{}\times \text{}({x}_{R}{x}_{Q})\text{}\le \text{}0\text{.}$$(14) The coordinates (y_{Q}, z_{Q}) of this point can be determined by : $${y}_{Q}=r\text{}\mathrm{cos}(\beta )\text{,}$$(15) $${z}_{Q}=r\text{}\mathrm{sin}(\beta )\text{.}$$(16)The value of β is already determined at each point of the exit section chosen in the first step. It is the same at all points of the corresponding contour. It can be determined by the following relation: $$\beta =a\mathrm{tan}\left(\frac{{y}_{Q}}{{z}_{Q}}\right)\text{.}$$(17)The point Q in relation (17) is that of the exit section and not those given by relations (15) and (16). The point Q given by (15) and (16) and that used in relation (17) lie on the same longitudinal plane formed by the contour under consideration but are not merged.
The value of r can be determined by the following relationship, referring to Figure 3. $$r=\frac{({r}_{R}{r}_{L})}{({x}_{R}{x}_{L})}({x}_{Q}{x}_{L})+{r}_{L}\text{.}$$(18) With $${r}_{L}=\sqrt{{y}_{L}^{2}+{z}_{L}^{2}}\text{,}$$(19) $${r}_{R}=\sqrt{{y}_{R}^{2}+{z}_{R}^{2}.}$$(20)The point Q in relation (14) represents all the points of the Figure 3. The properties at points L and R of Figure 4 represent those of intermediate contour already calculated. It can also be considered that the Figure 3 also represents the shape of the section of the nozzle 3D in each longitudinal station.
The coordinates (x, y, z) and the flow properties (M, T/T_{0}, P/P_{0}, ρ/ρ_{0} and θ) at each point of the section can be determined for each contour of the nozzle 3D by starting the calculation from the exit section towards the throat of the nozzle. The abscissa of these points belonging to a section of the nozzle of Figure 4 is given by x_{Q} chosen arbitrarily. The points in bold circles in Figure 4 represent the intermediate contour determined from the main axisymmetric shape and corresponding to the throat radius y_{C}. The points in hollow circles illustrated in Figure 4 represent the points at the selected sections of the nozzle 3D, where it is desired to determine their shapes and the physical properties.
The computation of the angle β by computer from relation (17) generally gives the angle in the interval [−π/2, π/2]. Then we have a cosine and sinus calculation problem of the four coordinates (y, z) and (−y, z), (y, −z) and (−y, −z) from relations (15) and (16).
The first and the fourth points have the same ratio y/z. Then the calculation of the angle by the relation (17) gives the angle corresponding to the point of the first quadrant. Then for the point of the third quadrant (fourth point), one will have a problem of calculation of the corresponding angle. To solve this problem, the angle of the third quadrant must be corrected by 2πβ where β is the angle for the point of the first quadrant.
For the second and third points we will have the same problem. The computation by computer will confuse the calculation of the angle β for these two points by the relation (17). Then the angle β of the second quadrant point (y, −z) is equal to π + the angle β of the point (−y, z) of the fourth quadrant.
The second problem is found for the particular points (y = 0.0, z) and (y = 0.0, −z) of the horizontal axis. The angle for the first point is β = 0.0 (no problem for this point) and the second is β = π. The computation by computer of the angle β corresponding to the point of the negative axis by the relation (17) still gives β = 0.0. Then we must correct this result by adding π for this point.
The third problem encountered is the computation of the angle β by computer for the particular points (y, z = 0.0) and (−y, z = 0.0) of the vertical axis. We will have a division by zero if we use the relation (17). In this case, the calculation is ignored by the relation (17) by assignment β = π/2 for the first point on the positive axis and β = −π/2 for the second point on the negative axis.
In this way it is possible to calculate exactly the corresponding angle β for all points of Figure 3 according to its position and there will be no cosine and sinus calculation problem in relation (15) and (16). Therefore, the position of the corresponding point is exactly calculated.
The flow parameters at this new point Q of the Figure 4 can be determined by first interpolating the parameters T and θ by the following relations: $${T}_{Q}=\frac{{S}_{LQ}}{{S}_{LR}}({T}_{R}{T}_{L})+{T}_{L}\text{,}$$(21) $${\theta}_{Q}=\frac{{S}_{LQ}}{{S}_{LR}}({\theta}_{R}{\theta}_{L})+{\theta}_{L}\text{.}$$(22) With $${S}_{LQ}=\sqrt{{({y}_{L}{y}_{Q})}^{2}+{({z}_{L}{z}_{Q})}^{2}}\text{,}$$(23) $${S}_{LR}=\sqrt{{({y}_{L}{y}_{R})}^{2}+{({z}_{L}{z}_{R})}^{2}}\text{.}$$(24)The ratio T/T_{0} at the point Q is equal to T_{Q}/T_{0}. The values of M, ρ/ρ_{0} and P/P_{0} at the point Q can be determined by relations (4), (7) and (8) respectively, by replacing T by T_{Q} that obtained by the relation (21). The integration of the relation (8) is done by the Gauss Legendre quadrature [21,22] for the speed of calculation with high accuracy.
By doing the same procedure for all the N_{T} points of the selected section of the nozzle 3D of Figure 3 in all the N_{L} stations of Figure 4, it is possible to determine the positions of all the points of the selected sections of the nozzle 3D.
Fig. 4 Longitudinal discretization of the nozzle 3D. 
5 Comparison and validation
The validation of the developed method is controlled by choosing a circular section of the nozzle 3D and calculating the corresponding complete shape of the nozzle 3D at each longitudinal station by applying the developed method and comparing it with the results of the axisymmetric nozzle. The shape of the nozzle in all longitudinal sections must be exactly circular in shape. To confirm the circular shape, we calculate the vector radius of all the points forming a section. This radius must be the same in all respects and independently of the angle β. By repeating the process for all selected sections of the nozzle, it can be confirmed that the nozzle is of circular crosssection. Hence the method is validated. It is possible to widen the validation by calculating the flow properties at eatch point, and to calculate the corresponding angle β. Then it is necessary to find for each section that the flow parameters do not depend on the angle β. In this case, it can be said that the program realized on the basis of the developed method gives good results for any chosen section.
For an arbitrarily chosen noncircular section, the shape of the section at all longitudinal stations must not necessarily be kept on a scale except with the shape of the exit section and the throat.
6 Results and comments
Table 1 shows the effect of Δθ of the Kernel region on the boundary and internal mesh parameters N_{F}, N_{S}, N_{P} and the precision ε for the axisymmetric nozzle when using the MOC with Δ = 0.1, M_{E} = 3.00 and T_{0} = 2000 K.
Table 2 shows the effect of Δ of the transition region on boundary and internal mesh parameters N_{F}, N_{S}, N_{P} and ε obtained for the axisymmetric nozzle by using MOC when Δθ = 0.1 deg, M_{E} = 3.00 and T_{0} = 2000 K.
Table 3 shows the effect of M_{E} on the boundary and internal mesh N_{F}, N_{S}, N_{P} and ε for the axisymmetric nozzle using the MOC when Δθ = 0.1 deg, Δ = 0.1 and T_{0} = 2000 K.
Table 4 shows the effect of T_{0} on the boundary and internal mesh N_{F}, N_{S}, N_{P} and ε for the axisymmetric nozzle using the MOC when Δθ = 0.1 deg, Δ = 0.1 and M_{E} = 3.00.
Table 5 represents the values of N_{F}, N_{S}, N_{P} and ε of a case generally giving a very fine mesh in the axisymmetric nozzle. This mesh is obtained when Δθ = 0.01 deg, Δ = 0.1 and M_{E} = 3.00 and T_{0} = 2000 K.
It is noted that the internal mesh N_{S}, N_{P} characteristics and the boundary mesh N_{F} characteristics of an axisymmetric nozzle depend mainly on Δθ, Δ, M_{E} and T_{0}. The number N_{F} of points obtained on the boundary of the wall is very small compared to the numbers N_{S} of the segments and N_{P} of internal mesh points in an axisymmetric nozzle. Then the methods presented in the references [1–6] based on the search of the internal stream lines from internal mesh require a very high computer execution time and which becomes impossible for the fine discretization as the example of the table 5. While our present method based on the points N_{F} of the boundary of the nozzle requires a very small execution time and that becomes negligible before the time of execution of the methods of the references [1–6].
According to table 5, when Δθ = 0.01 deg and Δ = 0.01, the accuracy of the results of the axisymmetric nozzle is equal to ε = 0.027% which is an acceptable error for calculation. But in parallel the methods presented in the references [1–6] use N_{S} = 13317795 and N_{P} = 6675864 for the search of the profiles of the nozzle 3D. While our method uses only the points of the boundary equal in this case to N_{F} = 5542. From the references [1–6], the N_{S} and N_{P} numbers are used to search for each N_{L} × N_{T} points of the nozzle 3D. This demonstrates that the method used in those references is limited by a maximum of N_{S} and N_{P} corresponding to the precise values of Δθ and Δ of discretization. Given the obtained number N_{F}, the computer execution time of our method is very small and incomparable to that of the method of [1–6].
Figures 5 and 6 represent the intermediate contours of the axisymmetric nozzle wall, respectively without and with the uniform horizontal line BE of Figure 2 for different values of y_{C}. The example presented is for y_{C} = 0.2 (curve 1), y_{C} = 0.5 (curve 2), and y_{C} = 0.8 (curve 3) and y_{C} = 1.0 (curve 4), while the main nozzle is taken for y_{C} = 1.0. The intermediate contours of Figure 6 will be used for determining the contours of the wall of the nozzle 3D after selecting the exit section. The example taken is for M_{E} = 300 and T_{0} = 2000 K.
Figures 7, 8, 9, 10 and 11 represent respectively the variation of M, θ, T/T_{0}, ρ/ρ_{0} and P/P_{0} on the intermediate contours of the nozzles walls of the Figure 6. The axisymmetric nozzle has an inflection point considering the increase of θ from θ^{*} to θ_{Max} at the point of inflection and then decreases at θ = 0.0 at the exit of the nozzle.
Figure 12 shows the shapes of the nozzle 3D for 16 selected sections. The shape 1 has a circular cross section. The shape 2 quadrilateral, the shape 3 square, the shape 4 is a half ellipse. The other 12 sections are arbitrarily chosen to showcase the power of our numerical program. It should be noted that the shape of the throat section is identical to that chosen in the exit section. While the intermediate transverse shapes are totally different from that of the exit section except for the circular shape 1. This difference can be clearly seen on the 11 presented shapes. Our program can do any complex unsymmetrical section. It should be informed that the sections numbers 2, 4, 8, 13, 14, 15 and 16 do not contain the axis of symetry of the parent nozzle.
The nozzles shown in Figure 12 all have the same longitudinal length L = 9.0504 and the same ratio of the crosssections A_{E}/A_{*} = 4.9958 as those of the axisymmetric nozzle, considering that they are designed for M_{E} = 3.00 and T_{0} = 2000 K.
To make the comparison between the nozzle 3D and the axisymmetric from the point of view of the mass obtained and the thrust force delivered by the wall, it is necessary to design 3D shapes having the same critical crosssectional area as the axisymmetric nozzle in addition to the preservation of M_{E}.
Figures 13, 14, 15, 16 and 17 show respectively the variation of M, θ, T/T_{0}, ρ/ρ_{0} and P/P_{0} as a function of the transverse angle β of the nozzle 11 of Figure 12. The nozzle 11 has been taken as an example of an illustration to show that the physical parameters in the nozzle of arbitrary 3D section depend essentially on x, y and z if the Cartesian coordinates are used or on x, r and β if the cylindrical coordinate is used, with the exception of the circular section that the physical parameters do not depend on β. For this reason the axisymmetric nozzle is named by the quasithreedimensional nozzle, while the other sections are cases of threedimensional flows. It is noted that in the throat section and the exit section of the nozzle 3D, the flow parameters do not depend on the angle β since the flow in both sections is unidirectional and even parallel. In the throat section M = M^{*} which represents the value of M just after the expansion. Then through the nozzle there is an abrupt increase in M between M = 1.0 and M = M^{*} at the same point A of Figure 1, then a progressive increase along the wall of the nozzle between M = M^{*} at point A, to M ═ M_{E} at the exit section (point E in Fig. 1).
Figure 14 shows that there is an increase in θ of θ = θ^{*} at the point A to θ = θ_{Max}>θ^{*} at the inflection point which is in the vicinity of the throat, then decrease of θ from θ = θ_{Max} of inflection point to θ = 0.0 at the exit of the nozzle. On an intermediate section, the angle θ depends on the angle β of Figure 3. Only the circular section, the radius r of the crosssections does not depend on the angle β. For the other noncircular sections, the radius r depends essentially on the angle β.
The results shown in Figures 13, 14, 15, 16 and 17 are valid for any shape of the exit section except for the circular section originally placed at a coordinate reference (nozzle 1 of Fig. 12).
Δθ effect on the mesh of the axisymmetric nozzle for Δ = 0.1, M_{E} = 3.00 and T_{0} = 2000 K.
Δ effect on the mesh of the axisymmetric nozzle for Δθ = 0.1 (deg), M_{E} = 3.00 and T_{0} = 2000 K.
M_{E} effect on the mesh of the axisymmetric nozzle for Δθ = 0.1 (deg), Δ = 0.1 and T_{0} = 2000 K.
T_{0} effect on the mesh of the axisymmetric nozzle for Δθ = 0.1 (deg), Δ = 0.1 and M_{E} = 3.00.
Unfavorable case on the axisymmetric nozzle mesh for Δθ = 0.01 (deg), Δ = 0.01, M_{E} = 5.00 and T_{0} = 3500 K.
Fig. 5 Intermediate contours of the axisymmetric nozzle without the uniform part BE. 
Fig. 6 Intermediate contours of the axisymmetric nozzle with the uniform part BE. 
Fig. 12 Transverse contours of nozzle 3D for 16 exit section shapes when M_{E} = 3.00 for T_{0} = 2000 K. 
7 Conclusion
A novel numerical effective and rapid method for the design of nozzle 3D of arbitrary exit section nozzles has been presented with respect to the axisymmetric nozzle applied to any type of nozzle in the case of HT or PG models. The following conclusions can be drawn from this work:
The design of the axisymmetric nozzle in the nondimensional shape is necessary.
The choice of unsymmetrical exit section of the nozzle 3D is arbitrary.
The developed numerical program can process any natural gas. It suffices to add the function C_{P}(T) and R of gas in question and to calculate the enthalpy H(T). The application is made only for air.
An infinity of shape of the nozzles 3D according to the choice of the parameters M_{E}, T_{0}, the shape of the exit section of the nozzle 3D and the gas can be obtained.
The most important of this method is that it is independent of the mesh results of the axisymmetric nozzle, which is not the case for the other methods used until now. This advantage has enabled us to considerably reduce the calculation time compared to the other methods currently used.
The nozzle 3D design by this method depends solely on the position and properties of the boundary points of the axisymmetric nozzle.
The calculation accuracy of this method depends on the design accuracy of the axisymmetric nozzle.
The execution time of this new method is decreased in a very considerable and important way and becomes negligible compared to the calculation time of the old design methods currently used.
The calculation of the current function is not necessary for this method for the determination of the coordinates of the points of the nozzle 3D, which is not the case for the other methods.
New contours of the nozzle 3D are found by this new method.
The design of the nozzle 3D can be done by this method using a simple nonpowerful calculator with a very high discretization without interrupting the execution of numerical program.
The shape of the throat section of the nozzle 3D is identical with that of the exit section. These two sections verified the ratio of the critical sections of the axisymmetric nozzle.
Between the throat and the exit section, the nozzle 3D takes an arbitrary shape which is not identical to the throat and exit.
The axisymmetric form is used to compare the results found by our method.
As a very advanced technique to further improve the performances of the nozzle 3D, the nozzle 3D can be truncated in a section close to the exit section. In this case a gain of a portion of the mass is observed and in parallel a slight loss of thrust is noticed.
The flow in the axisymmetric nozzle depends on two variables (x and y) and does not depend on the angle β, while the flow in the nozzle 3D depends on three variables (x, r and β).
As a perspective, and in order to see the improvement of the performances of the nozzle 3D relative to the axisymmetric nozzle, the same crosssection of the throat between the two nozzles is maintained, in addition to fixing another parameter among the exit Mach number, the mass of the nozzle or the thrust coefficient. The two remaining parameters determine the performance of the new obtained shape.
A second problem in perspective consists in removing the uniform flow portion of the wall from the nozzle found for each contour in order to decrease in addition the mass of the nozzle 3D with the same Mach number performances and the delivered thrust force.
A third problem of improving the performances of the nozzle 3D is to place the exit section of the nozzle 3D in the uniform zone of the axisymmetric nozzle near the end of the Kernel region in order to decrease the maximum possible the undesirable uniform area.
Nomenclature
x: Abscissa of a section of the nozzle
y, z: Position of a point in the section
r: Radius of a point in a section
β: Angle formed by the vector radius between two points of a section
Δθ: Flow deviation step chosen in the expansion center A of Figure 1
Δ: Step progress of the characteristics on the uniform BE line in the transition zone
R: Thermodynamic constant of gas
C_{P}: Specific heat at constant pressure
N_{L}: Number of longitudinal sections of the three dimensional nozzle
N_{T}: Number of points in a crosssection of the three dimensional nozzle
N_{F}: Number of points on the boundary of the axisymmetric nozzle
N_{S}: Number of mesh segments of axisymmetric nozzle
N_{P}: Number of mesh points of axisymmetric nozzle
MOC: Method Of Characteristics
Indices
0: Stagnation condition (combustion chamber)
Acknowledgments
The author acknowledges Khaoula, AbdelGhani Amine, Ritadj and Assil Zebbiche and Mouza Ouahiba for granting time to prepare this manuscript.
References
 J.C. Evvard, S.H. Maslen, Threedimensional supersonic nozzles and inlets of arbitrary exit cross section, NASA TN 2688, 1952 [Google Scholar]
 A. Haddad, J.B. Moss, Aerodynamic design for Supersonic nozzles of arbitrary cross section, J. Propuls. Power 6 (1990) 740–746 [CrossRef] [Google Scholar]
 M.K. Aukin, A.N. Kraiko, D.A. Lubinov, V.E. Makarov, N.I. Tillyaeva, Designing three dimensional nozzles to achieve near uniform supersonic or hypersonic flow in the rectangular exit section, Phys. Astron. 30 (1995) 787–794 [Google Scholar]
 I.E. Beckwith, H.W. Ridyard, N. Cromer, The aerodynamic design of high Mach number nozzles utilizing axisymmetric flow with application to a nozzle of square test section, NACA TN 2711, 1952 [Google Scholar]
 A.I. Rylov, Design of supersonic asymmetric nozzles, Phys. Astron. 12 (1977) 414–420 [Google Scholar]
 O. Abada, T. Zebbiche, A. Abdallah ElHirtsi, Threedimensional supersonic minimum length nozzle design at high temperature for arbitrary exit cross section, Arab J. Sci. Eng. 39 (2014) 8233–8245 [CrossRef] [Google Scholar]
 Jr. J.D. Anderson, Modern compressible flow with historical perspective, 2nd edition, McGrawHill Book Company, New York, USA, 1982 [Google Scholar]
 T. Zebbiche, Stagnation temperature effect on the supersonic axisymmetric minimum length nozzle design with application for air, Adv. Space Res. 48 (2011) 1656–1675 [CrossRef] [Google Scholar]
 B.M. Argrow, G. Emanuel, Comparison of minimum length nozzles, J. Fluid Eng. 110 (1988) 283–288 [Google Scholar]
 L.Z. Dumitrescu, Minimum length axisymmetric laval nozzles, AIAA J. 13 (1975) 520–532 [CrossRef] [Google Scholar]
 C.R. Peterson, P.G. Hill, Mechanics and Thermodynamics of Propulsion, AdditionWesley Publishing Company Inc., New York, USA, 1965 [Google Scholar]
 G.P. Sutton, O. Biblarz, Rocket propulsion elements, 8th edition, John Wiley and Sons, 2010, New York, USA [Google Scholar]
 S.M. Patel, D.S. Mane, M. Raman, Concepts and CFD analysis of Laval nozzle, Int. J. Mech. Eng. Technol. 7 (2016) 221–240 [Google Scholar]
 H.D. Deshpande, S.S. Vidwans, P.R. Mahale, R.S. Joshi, K.R. Jagtap, Theoretical & CFD analysis of laval nozzle, Int. J. Mech. Prod. Eng. 22 (2014) 33–36 [Google Scholar]
 V. Venkatesh, C.J. Reddy, Modelling and simulation of supersonic nozzle using computational fluid dynamics, Int. J. Nov. Res. Interdiscip. Stud. 12 (2015) 16–27 [Google Scholar]
 G.M. Kumar, D.X. Fernando, R. Muthu Kumar, Design and optimization of lavel nozzle to prevent shock induced flow separation, Adv. Aerosp. Sci. Appl. 13 (2013) 119–124 [Google Scholar]
 S.T. Travis, Introduction to rocket science and engineering, CRC Press, Taylor and Francis group, USA, 2006 [Google Scholar]
 T. Zebbiche, Z. Youbi, Effect of stagnation temperature on the supersonic flow parameters with application for air in nozzles, Aeronaut. J. 111 (2007) 31–40 [CrossRef] [Google Scholar]
 T. Zebbiche, Z. Youbi, Supersonic flow parameters at high temperature, Application for Air in nozzles, German Aerospace Congress 2005, DGLR20050256, Friendrichshafen, Germany, 2005 [Google Scholar]
 B.J. McBride, S. Gordon, M.A. Reno, Coeffcients for calculating thermodynamic and transport properties of individual species, NASA TM 4513, 1993 [Google Scholar]
 B. Démidovitch, I. Maron, Eléments de calcul numérique, Editions MIR, Moscow, Russia, 1987 [Google Scholar]
 A. Raltson, A. Rabinowitz, A first course in numerical analysis, McGraw Hill Book Company, 1985, New York, USA [Google Scholar]
Cite this article as: T. Zebbiche, Rapid design method and new contours for a class of three dimensional supersonic nozzle of arbitrary exit cross section, Mechanics & Industry 19, 403 (2018)
All Tables
Δθ effect on the mesh of the axisymmetric nozzle for Δ = 0.1, M_{E} = 3.00 and T_{0} = 2000 K.
Δ effect on the mesh of the axisymmetric nozzle for Δθ = 0.1 (deg), M_{E} = 3.00 and T_{0} = 2000 K.
M_{E} effect on the mesh of the axisymmetric nozzle for Δθ = 0.1 (deg), Δ = 0.1 and T_{0} = 2000 K.
T_{0} effect on the mesh of the axisymmetric nozzle for Δθ = 0.1 (deg), Δ = 0.1 and M_{E} = 3.00.
Unfavorable case on the axisymmetric nozzle mesh for Δθ = 0.01 (deg), Δ = 0.01, M_{E} = 5.00 and T_{0} = 3500 K.
All Figures
Fig. 1 Presentation of the axisymmétric MLN flow field. 

In the text 
Fig. 2 Intermediates contour of the axisymmetric nozzle for nozzle 3D. 

In the text 
Fig. 3 Transverse discretization of the exit section of the nozzle 3D. 

In the text 
Fig. 4 Longitudinal discretization of the nozzle 3D. 

In the text 
Fig. 5 Intermediate contours of the axisymmetric nozzle without the uniform part BE. 

In the text 
Fig. 6 Intermediate contours of the axisymmetric nozzle with the uniform part BE. 

In the text 
Fig. 7 Variation of M on the contours of Figure 6. 

In the text 
Fig. 8 Variation of θ on the contours of Figure 6. 

In the text 
Fig. 9 Variation of T/T_{0} on the contours of Figure 6. 

In the text 
Fig. 10 Variation of ρ/ρ_{0} on the contours of Figure 6. 

In the text 
Fig. 11 Variation of P/P_{0} on the contours of Figure 6. 

In the text 
Fig. 12 Transverse contours of nozzle 3D for 16 exit section shapes when M_{E} = 3.00 for T_{0} = 2000 K. 

In the text 
Fig. 13 Variation of M on the wall of the transverse sections of the nozzle 11 of Figure 12. 

In the text 
Fig. 14 Variation of θ on the wall of the transverse sections of the nozzle 11 of Figure 12. 

In the text 
Fig. 15 Variation of T/T_{0} on the wall of the transverse sections of the nozzle 11 of Figure 12. 

In the text 
Fig. 16 Variation of ρ/ρ_{0} on the wall of the transverse sections of the nozzle 11 of Figure 12. 

In the text 
Fig. 17 Variation of P/P_{0} on the wall of the transverse sections of the nozzle 11 of Figure 12. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.