Rapid design method and new contours for a class of three dimensional supersonic nozzle of arbitrary exit cross section

. 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 ﬂ ow solutions in axisymmetric nozzle. The application is made for nozzle giving a uniform and parallel ﬂ ow at the exit section. The determination of the points of the axisymmetric nozzle contours in a non-dimensional 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 non-dimensional 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 ﬂ ow properties of each point are the same as those of the points of the axisymmetric nozzle. The temperature and the deviation of the ﬂ ow 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.


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][2][3][4][5][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][2][3][4][5][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][9][10][11]. The number of information on the flow parameters (x, y, P, T, r, M, u and c) 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][2][3][4][5][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][2][3][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 second-degree 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][2][3][4][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][2][3][4][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][14][15][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 two-dimensional 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][2][3][4][5][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][2][3][4][5][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][2][3][4][5][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][12][13][14][15][16][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.

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, u, P and r) 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 Du and D 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 e depend essentially on Du, D, 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: With The M and T at HT are connected by [18,19]: The ratios r/r 0 and P/P 0 are calculated by [18,19]: 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]: To validate the numerical results of design, one can calculate the relative error committed the ratio of the critical sections by the following relation: 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 g = 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][19][20] so as not to have dissociation of the molecules.

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 non-dimensional 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 non-dimensional 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: The flow properties (M, P/P 0 , T/T 0 , r/r 0 and u) 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: 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 b) as well as (M, T/T 0 , P/P 0 , r/r 0 and u) 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: 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 b as shown in

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 , r/r 0 and u) 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: The coordinates (y Q , z Q ) of this point can be determined by: The value of b 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: 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.
With 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 , r/r 0 and u) 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 b by computer from relation (17) generally gives the angle in the interval [Àp/2, p/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)  For the second and third points we will have the same problem. The computation by computer will confuse the calculation of the angle b for these two points by the relation (17). Then the angle b of the second quadrant point (y, Àz) is equal to p + the angle b 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 b = 0.0 (no problem for this point) and the second is b = p. The computation by computer of the angle b corresponding to the point of the negative axis by the relation (17) still gives b = 0.0. Then we must correct this result by adding p for this point.
The third problem encountered is the computation of the angle b 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 b = p/2 for the first point on the positive axis and b = Àp/2 for the second point on the negative axis. In this way it is possible to calculate exactly the corresponding angle b 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 u by the following relations: With The ratio T/T 0 at the point Q is equal to T Q /T 0 . The values of M, r/r 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.

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 b. By repeating the process for all selected sections of the nozzle, it can be confirmed that the nozzle is of circular cross-section. 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 b. Then it is necessary to find for each section that the flow parameters do not depend on the angle b. 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 non-circular 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.  Table 1 shows the effect of Du of the Kernel region on the boundary and internal mesh parameters N F , N S , N P and the precision e for the axisymmetric nozzle when using the MOC with D = 0.1, M E = 3.00 and T 0 = 2000 K. Table 2 shows the effect of D of the transition region on boundary and internal mesh parameters N F , N S , N P and e obtained for the axisymmetric nozzle by using MOC when Du = 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 e for the axisymmetric nozzle using the MOC when Du = 0.1 deg, D = 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 e for the axisymmetric nozzle using the MOC when Du = 0.1 deg, D = 0.1 and M E = 3.00. Table 5   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 Du, D, 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][2][3][4][5][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][2][3][4][5][6].
According to table 5, when Du = 0.01 deg and D = 0.01, the accuracy of the results of the axisymmetric nozzle is equal to e = 0.027% which is an acceptable error for calculation. But in parallel the methods presented in the references [1][2][3][4][5][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][2][3][4][5][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 Du and D 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][2][3][4][5][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), 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, u, T/T 0 , r/r 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 u from u * to u Max at the point of inflection and then decreases at u = 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 Table 5. Unfavorable case on the axisymmetric nozzle mesh for Du = 0.01 (deg), D = 0.01, M E = 5.00 and T 0 = 3500 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 cross-sectional area as the axisymmetric nozzle in addition to the preservation of M E . Figures 13, 14  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).

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 non-powerful 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 b, while the flow in the nozzle 3D depends on three variables (x, r and b).
As a perspective, and in order to see the improvement of the performances of the nozzle 3D relative to the axisymmetric nozzle, the same cross-section 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.