Abstract

There is currently a lack of efficient heat transfer analysis methodologies for spiral channel regenerative cooling that has been increasingly applied in liquid rocket engines. To figure out the heat transfer characteristics of the spiral channel regenerative cooling thrust chamber, a simple 1D method based on the traditional semi-empirical formula after correcting the flow velocity is proposed. The accuracy of this approach is verified by the 3D numerical simulation. The verified method is further used to analyze the distribution of inner wall temperature in the test case and optimize the channel’s parameters. The research shows that the maximum inner wall temperature cooled by the spiral channel is 8.5% lower than that of the straight channel under the same channel size and boundary condition, indicating that the application of the spiral channel significantly improves the cooling effect. In addition, the 1D model combined with the second-order response surface model is employed to optimize the channel width, channel height, pitch, and inner wall thickness aiming for the best cooling effect. The calculated maximum temperature of the inner wall after the structure optimization is 586.6 K, which is 29.8% lower than the initial structure before optimization.

1. Introduction

Liquid rocket engine has been widely used in space propulsion system because of its advantages of high thrust, high reliability, and low cost [1]. As the key subassembly of a rocket engine, the thrust chamber provides thrust by accelerating the hot gas to supersonic velocity across the throat. During its operation, the temperature and pressure of a modern high-performance thrust chamber are as high as 3500 K and over 20 MPa, respectively. Hereby, an efficient cooling design is a prerequisite for ensuring structural integrity under such extreme thermomechanical loads [2]. As an efficient and economical cooling method, regenerative cooling is mostly used in the current thrust chamber of a liquid rocket engine by flowing a liquid propellant into the cooling channels machined in the inner liner [3, 4].

Limited by traditional processing methods, the regenerative cooling channels generally adopt straight channels. While by increasing the coolant velocity, the spiral cooling channels are believed to achieve a better cooling effect. As early as the 1930s and 1940s, the Soviet ORM-65 and RD-1-150 engines both exploited the spiral channel structure [5]. Milling is the main manufacturing method for these costly and difficult-to-machine spiral channels [68]. With the development of additive manufacturing technology, the spiral channel is expected to be applied on a larger scale. For example, NASA 40 K-IBf nozzle with spiral channels is manufactured by a typical additive manufacturing technology, namely, LP-DED [9]. The advantage of additive manufacturing is that complex spirals can be machined more easily on the one hand. On the other hand, joint methods such as welding and electroforming can be avoided, so that it is more difficult for connecting parts to fail. Besides the research on the manufacturing technology of the spiral channel, a comprehensive understanding of heat transfer characteristics within the chamber is of utmost importance to quantitatively analyze thermomechanical loads and improve the performance of propulsion systems [10, 11].

At present, there are few related studies on the spiral channel regenerative cooling analysis of the thrust chamber [12]. Heat transfer analysis is basically concentrated on the straight channel, which has formed two categories: the 1D semi-empirical analytical method and 2D/3D numerical simulation based on CFD/FEM. The traditional thrust chamber heat transfer calculation during the preliminary design essentially adopts the 1D method, which is based on the Nusselt-type empirical relations for both hot gas and coolant sides [13]. Bartz first proposes an empirical formula for evaluating the convection heat transfer of hot gas in the thrust chamber [14]. From then on, the Bartz formula combined with the inner wall heat conduction and the tube flow forced convection heat transfer of coolant is widely applied in the thrust chamber heat transfer analysis [15, 16]. With the development of computational fluid dynamics and finite element methods, numerical simulation is more likely to be used to solve the heat transfer problems because of its accuracy and more comprehensive information [17]. Betti et al. [18] employ a numerical procedure to predict the wall temperature and heat flux as well as coolant flow characteristics in a regeneratively cooling thrust chamber. Their results show a reasonable agreement in terms of coolant pressure drop and temperature gain with the published results for the space shuttle main engine (SSME). Furthermore, numerical simulation can also couple the thrust chamber heat transfer and propellant combustion processes together to obtain more accurate thermal load information [19]. A methodology for the prediction of coupled heat transfer of combustion and regenerative cooling is proposed by Song and Sun [2] in a typical LOX/methane thrust chamber. They have shown that hot-gas-side wall heat flux and temperature are nonuniform and vary periodically in the circumferential direction. Rahn et al. [20] conduct the heat transfer simulation based on the frozen nonadiabatic flamelet combustion model. The good agreement with experimental heat load data confirms the accuracy of this conjugate heat transfer simulation.

In the process of conducting the straight channel heat transfer analysis, researchers find that the parameters of the cooling channel affect the thrust chamber thermal load significantly [21, 22]. Riccius et al. [23] discover that decreasing the inner wall thickness, width, and height of the cooling channel will improve the cooling effect. In order to minimize the thermal load of the regenerative cooling thrust chamber, it is necessary to design the cooling channel parameters through various optimization algorithms [24, 25]. Ji and Sun [26] use three different optimization methods, namely, the response surface method, combined with the line search method and the conjugate gradient method, to optimize the size of the regenerative cooling channel. The optimized structures have the same characteristics of a higher aspect ratio and thinner wall thickness. Xu et al. [27] establish an optimal design model based on the backpropagation (BP) neural network method and the global search algorithm to attain the optimal design results of the regenerative cooling structural parameters. The optimized SSME-MCC regenerative cooling channel reduces both the maximum wall temperature and the coolant pressure loss. Judging from the existing literature, whether it is heat transfer analysis or optimization, researchers pay attention to straight channel regenerative cooling. However, related research on the spiral channel regenerative cooling is quite scarce.

In short, there is a lack of research on the heat transfer and optimization design of spiral cooling channels at present, which directly leads to the inability to quantitatively evaluate the improvement of the cooling effect brought by the spiral channel structure, and the inability to complete reasonable design of channel parameters. In response to the above problems, this paper proposes a 1D calculation method for spiral channels based on the simplified semi-empirical formula. Moreover, a 3D numerical simulation is applied to verify the accuracy of the proposed 1D method, in addition to obtaining the temperature distribution over the entire inner wall domain. At last, to obtain the cooling channel parameters with the best cooling effect, a second-order response surface model is established for optimization. The optimized maximum temperature of the inner wall is reduced by 29.8%.

This paper is arranged as follows: a 1D heat transfer calculation model is proposed in Section 2. In addition, the basic theory of 3D numerical heat transfer simulation that will be used is also introduced. In Section 3, after designing the spiral channel regeneration cooling thrust chamber, we use the 1D method to analyze the heat transfer characteristics and verify its accuracy by 3D numerical simulation. Afterwards, the optimal configuration of the spiral channel parameters is designed by the second-order response surface model in Section 4. Lastly, the main conclusions are stressed in Section 5.

2. Methods

Because of its unique geometric configuration, the flow characteristics of the coolant in the spiral channel are different from those in the straight channel, so the heat transfer analysis is also distinguished from that of the straight channel. This section will propose a 1D heat transfer calculation method for the spiral channel. In addition, 3D numerical models for method validation will be presented.

2.1. 1D Heat Transfer Analysis

A comprehensive thermal model for a regeneratively cooling thrust chamber must account for convection from the hot gas to the inner wall, heat conductance through the wall, and convection from the inner wall to the coolant, as shown in Figure 1 [28]. The 1D coupled flow and heat transfer model based on the criterion relation is a simple and efficient approach that tends to obtain reasonable results compared to 2D/3D numerical simulation. Due to the thin inner wall and high thermal conductivity, this system can quickly establish thermal equilibrium, so it can be regarded as a steady-state analysis without considering time effects. The following assumptions are made: (i)Treat the heat flow as 1D, and only propagate through the inner wall of the thrust chamber in the radial direction(ii)Regard the inner wall of the thrust chamber as a flat plate, the area on both sides of the inner and outer walls of the thrust chamber is considered the same(iii)Ignore the radiation heat transfer process between the gas and the wall(iv)Since the heat transferred to the outside through the outer wall is much smaller than the heat absorbed by the coolant, it is assumed that the outer wall is an adiabatic boundary(v)The convective heat transfer process between the inner wall surface and the coolant side does not affect the combustion flow process of the propellant in the thrust chamber

2.1.1. Convection Heat Transfer of Hot Gas

The heat transfer from the combustion gases to the chamber walls of the thrust chamber occurs mainly through forced convection, which includes the heat transfer of the laminar part close to the wall and the heat transfer of the turbulent part in the boundary layer. The basic correlation for this complicated convective heat transfer can be expressed by where is the convective heat transfer heat flux, is the gas-side convective heat transfer coefficient, and is the adiabatic wall temperature of the gas, which can be obtained by [29] where is the stagnation temperature of hot gas, is the local Mach number, is the local recovery coefficient, is the effective recovery coefficient (0.90~0.98), and is the specific heat ratio of hot gas. The heat transfer coefficient in equation (1) is calculated using the Bartz formula [14]: where , , and denote diameter at the throat, area at the throat, and local cross-sectional area, respectively. , , and denote dimensionless factor, characteristic velocity, and throat radius of curvature, respectively. is the chamber pressure. The effect of the boundary layer is considered in the derivation of the Bartz formula.

2.1.2. Heat Conduction on the Inner Wall

The heat flux transferred from the gas wall surface to the liquid wall surface through the inner wall can be expressed as follows: where is the thermal conductivity of the material of the chamber wall and is the thickness of the inner wall.

2.1.3. Convection Heat Transfer of the Coolant

In the thermal boundary layer around the coolant, the heat exchange between the coolant and the inner wall of the thrust chamber is performed by forced convection. The heat flux can be calculated by the following formula: where is the convective heat transfer coefficient between the coolant and the inner wall surface, which can be calculated by formula (6). where is the heat transfer coefficient between the coolant and the wall surface without ribs, which can be calculated by equation (7), and is the rib coefficient, which can be calculated by equation (8) [30].

In the formula, is the channel width of the cooling channel, is the rib width of the cooling channel, and is the channel height of the cooling channel. is the Nusselt number, which can be calculated by [30] where and are the Prandtl number of the coolant and the Prandtl number of the coolant on the wall, respectively.

2.1.4. Analysis of Coolant Flow in Spiral Channel

Since the thrust chamber is segmented into relatively small sections along the axis direction in the 1D analysis, every section of the spiral channel can be treated as an equivalent straight channel following the law of mass conservation as shown in Figure 2. The equivalent flow per axial length of the spiral channel can be expressed as where θ is the angle between the spiral channel and the thrust chamber axis. and are the total flow mass of coolant in the spiral channel and straight channel, respectively.

Within the same time , the total coolant flow mass through the spiral channel and the straight channel can be expressed as where is the cross-sectional area of the cooling channel and and are the flow velocities of coolant in the spiral channel and straight channel, respectively. According to equations (10) and (11a) and (11b),

Therefore, the equivalent flow velocity in the spiral channel is faster than that of the straight channel, which is also a manifestation of the improved cooling effect for the spiral channel.

2.1.5. Solution Methodology

Based on the analysis of the above three heat transfer processes and the flow characteristics in the spiral channel, a calculation program is developed in the Fortran language. As shown in the flow chart of Figure 3, the main processes are summarized as follows: (i)Input the thrust chamber design parameters, and divide the chamber into multiple sections(ii)Calculate the thermodynamic parameters of hot gas per section such as the Prandtl number, specific heat capacity, and the Mach number by the Rocket Propulsion Analysis (RPA) Lite Edition [31]. Besides, the coolant parameters under different temperatures are obtained by establishing the interface of REFPROP (NIST) [32](iii)An initial hot-gas-side wall temperature is assumed to calculate heat transferred from the gas to the inner wall, as well as the convection heat transfer from the inner wall to the coolant(iv)According to the inner wall heat conduction formula, the hot-gas-side wall temperature is updated. The calculation will iterate until the error between and is within the acceptable range(v)After completing the solution of the above section, update the temperature and pressure of the coolant to solve the following sections. Finally, the distribution of inner wall temperature is acquired

2.2. 3D Numerical Heat Transfer Analysis

The main purpose of the 3D numerical heat transfer simulation in this paper is to verify the accuracy of the proposed 1D steady-state method; therefore, only the steady-state flow and heat transfer in the thrust chamber are considered. In the numerical simulation of heat transfer, the hot gas at the inlet is assumed to be completely combusted, and the atomization, vaporization, and combustion processes of propellant are not considered. The continuity, momentum, and energy equations governing the compressible and stable 3D flow of hot gas and coolant are given by [33]:

The shear stress in the preceding equations can be expressed as

The standard -epsilon equation is used for the turbulent flow of gas and coolant:

In the process of fluid-solid coupling heat transfer, the fluid domain and the solid domain follow different governing equations. In order to realize the exchange of data in different domains, the fluid and solid interface must meet the continuous boundary condition, namely, [34]

The above-governing equations of the fluid domain are discretized by the finite volume method. The viscous term and convective term are calculated by the central difference scheme and the second-order upwind scheme, respectively. The pressure-based coupled algorithm is applied to solve the heat transfer process.

3. Analysis of Heat Transfer Results

This paper takes a cryogenic liquid oxygen/methane engine thrust chamber as the research object, and its main design parameters are shown in Table 1. The engine is preliminarily designed for hot run testing, and liquid nitrogen is the best choice acting as the coolant considering our laboratory conditions. In this section, the spiral channel cooling configuration of the thrust chamber is designed according to the spiral equation. The 1D simplified algorithm is applied to analyze the results of spiral channel regenerative cooling, mainly the temperature distribution of the inner wall and the coolant. Also, the 3D temperature distribution of the thrust chamber of the spiral channel is analyzed. Finally, the accuracy of the proposed 1D simplified calculation method is verified by the results of the 3D method.

3.1. Geometric Model Design

As to the design of the thrust chamber spiral channel, the trajectory of the spiral must be determined first. The schematic diagram of the spiral is shown in Figure 4(a) with the center of the nozzle inlet as the origin and the axial direction as the -axis. The spiral line can be regarded as the spiral motion trajectory of point on the thrust chamber. The motion of this point can be divided into three directions, namely, the linear motion along the -axis and the circular motions and perpendicular to the -axis. is a linear motion with uniform velocity, while the motion perpendicular to the -axis is a circular motion with uniform angular velocity. Under the circumstances, the trajectory of point is the spiral, which can be expressed by equation (17) in the Cartesian coordinate system:

In the formula, is the pitch of the spiral, that is, the distance of point moving along the -axis after making a circular motion. Double arc nozzle is adopted in this article. This kind of nozzle configuration is suitable for the thrust chamber in this paper with relatively small shrinkage. The function of the double arc nozzle profile is shown in

The parameters of the cooling channel refer to the regenerative cooling design results of similar thrust chambers. The number of cooling channels is determined to be 30. According to the restriction relation between cooling channel dimensions and preliminary heat transfer calculation, the inner wall thickness and cooling channel width and height are 1 mm, 1.8 mm, and 3.5 mm, respectively. This size will be used for subsequent 1D and 3D heat transfer analyses. Based on the above-designed spiral and channel parameters, a 3D model of the thrust chamber of the spiral channel regenerative cooling is established by using a 3D modeling software, as shown in Figure 4(b).

3.2. 1D Heat Transfer Result

The heat transfer characteristics of the designed spiral channel thrust chamber are analyzed using the 1D simplified algorithm. The designed thrust chamber is divided into 100 sections along the axial direction, and the thermodynamic parameters of the hot gas in each section are firstly calculated by the Rocket Propulsion Analysis (RPA) Lite Edition. Table 2 lists the main parameter values of typical sections.

Figure 5 shows the calculation results of the intercepted convergence and expansion sections. Since there is no sudden shrinkage and expansion in the cooling channel, the changes in temperature, heat flux, and pressure are very smooth. As shown in Figure 5(a), the heat flux increases first, reaching the maximum value upstream of the throat, and then decreases along the axial direction, which is in line with the general distribution law of thrust chamber heat flux. Figure 5(b) demonstrates the temperature evolution of the hot-gas-side inner wall () and the coolant inner wall (). Both share the same trend and reach the maximum value of 835.5 K and 730.8 K upstream of the throat, which is consistent with the location of the maximum heat flux.

The coolant continuously absorbs heat in the process of flowing from the nozzle outlet to the inlet and gradually increases its temperature by 71.44 K. As illustrated in Figure 5(c), the fastest heating rate is located in the throat region. At the axial position of about 31 mm, the coolant is in a supercritical state, which persists until the coolant exits. In addition, due to the existence of flow resistance, the coolant endures a pressure loss of 0.071 MPa, as shown in Figure 5(d).

3.3. 3D Numerical Result

The 3D numerical simulation is employed to calculate the heat transfer result of the spiral channel. In addition to verifying the accuracy of the proposed 1D heat transfer calculation method for the spiral channel, a more detailed inner wall temperature distribution in 3D space can also be obtained. After that, the temperature distribution of the whole solid domain is studied, which provides a more comprehensive understanding of the heat transfer characteristics.

3.3.1. Grid Generation and Grid Independent Test

Different from the straight channel regenerative cooling thrust chamber, the spiral channel regenerative cooling thrust chamber is not circumferentially symmetrical, which requires a full 3D simulation model. In addition, the unstructured tetrahedral volume grid is mainly applied because of the complex shapes of the hot gas domain, coolant domain, and inner wall domain. To obtain a high-quality heat transfer solution, the boundary layer grid is set at the interface of the fluid and solid domains. The meshing result is shown in Figure 6.

In the simulation, both the parameters of the cooling channel and the boundary conditions such as the coolant mass inlet flow rate are the same as those of the 1D calculation. The inner wall thickness and cooling channel width and height are 1 mm, 1.8 mm, and 3.5 mm, respectively. The boundary conditions are set as shown in Figure 6. The inlet pressure of hot gas is 3 MPa, and the total temperature is 3449 K. The mass fractions of CO2, H2O, and CH4 are 0.5381, 0.4402, and 0.0217, respectively. The inlet mass flow rate of coolant is 0.05 kg/s, and the temperature is 70 K. The physical parameters of liquid nitrogen, such as the density, specific heat, viscosity, and thermal conductivity, are imported from the National Institute of Standards and Technology (NIST) data library. The no-slip boundary condition is specified for the thrust walls, and the fluid-solid interfaces are set as coupling conditions [32].

It is necessary to perform the grid independence test after grid discretization. In this paper, three different levels of fineness of meshes are used, namely, coarse grid, base grid, and fine grid. They are different in the element size of the three computational domains and the number of mesh layers in the boundary layer. Taking the inner wall as an example, the sizes in these three grids are 1 mm, 0.5 mm, and 0.3 mm, respectively. Figure 7 presents the calculated maximum gas-side wall temperature along the axial direction for various grid levels. It can be seen that the coarse grid will overestimate the inner wall temperature, while the calculation result of the fine grid is consistent with the basic grid with an error below 1%, which proves that the basic grid solution is independent of the grid. Therefore, the base grid shown in Figure 6 is adopted for subsequent studies. There are 8,450,544 elements and 2,349,038 nodes in the base grid.

3.3.2. Results of 3D Numerical Simulation

Figure 8(a) depicts the temperature distribution of the hot-gas-side inner wall surface. Along the spiral trajectory of the coolant, the wall temperature exhibits a regular feature; that is, the temperature is the lowest at the coolant inlet and reaches the maximum upstream of the throat. In addition, the temperature results of these 30 cooling trajectories are consistent, and they can be considered to follow a spiral symmetry law. These laws are similar to the general regenerative cooling thrust chamber wall temperature in Reference [2].

We also analyze the temperature distribution of the entire inner and outer wall domains. The 3D structure is sectioned along the radial direction to the 2D result shown in Figure 8(b). During the stable operation of the thrust chamber, the outer wall maintains a relatively low temperature, whose maximum value is only about 320 K. Regardless of the axial position, the maximum temperature of the inner wall is located in the center of the cooling channel. Since the tangent plane is not in the direction of the cooling channel spiral, the temperature change in the axial direction is not completely continuous. Also, it can be seen from the subplot that the temperature gradient of the inner wall and the rib along the radial direction is relatively large, especially at the connection between the inner wall and the rib.

3.4. Result Analysis and Verification

The distributions of gas-side wall temperature between the 1D simplified calculation method and the 3D simulation are compared for method validation. Moreover, to validate the improvement of the cooling effect of the spiral channel structure, the 1D and 3D thrust chamber models with straight channels are constructed. The model shares consistent coolant temperature, mass flow, cooling tank parameters, etc. with that of the spiral channel.

The 1D calculation method of the straight channel regenerative cooling is similar to that of the spiral channel, except that the processing of the equivalent velocity in equations (10)–(12) is omitted. The 3D numerical solution of straight channel regenerative cooling also undergoes a grid-independent check. The relevant solution settings are the same as the spiral channel regenerative cooling. Figure 9 shows the temperature distribution of the hot-gas-side wall and the 2D solid domain obtained by straight channel regenerative cooling. It can be found that the temperature distribution along the axial direction is similar to that of the spiral channel regenerative cooling; that is, the temperature reaches a maximum value upstream of the throat. The biggest difference between the two is that the temperature distribution of the straight channel cooling is circumferentially symmetrical.

The temperature data of the entire hot-gas-side surface in spiral and straight channel regenerative cooling obtained by 3D simulation are compared in Figure 10. We extract the temperature distribution data along the axial direction of the entire hot-gas-side wall, forming a wide band of data. The corresponding results calculated from the 1D methods are also plotted in Figure 10. The two methods share the same temperature evolution and the locations of the dangerous cross-section. However, it is clear that the 1D method using the Bartz formula overestimates the temperature of the inner wall, which is consistent with the conclusions in previous literatures [2, 35, 36]. A possible reason for the overestimation is that the Bartz formula does not take into account the effect of thickness variation of the boundary layer on the gas heat flux [29, 37]. In the 1D method, the temperature of the dangerous cross-section upstream of the throat decreased from 913.9 K to 835.5 K, a decrease of 78.4 K, after using the spiral channel relative to the straight channel regenerative cooling. In the 3D numerical simulation, the corresponding dangerous cross-section temperature decreased by 76.9 K from 686.6 K to 609.7 K. The deviation of the inner wall temperature reduction calculated by these two methods is only 1.95%, which verifies the accuracy of the 1D method proposed in this paper.

The cooling effect of the spiral channel is better than that of the straight channel in the entire axial range of the thrust chamber. For example, at the dangerous cross-section, the inner wall temperature using the spiral channel is 78.4 K (8.6%) lower than that of the straight channel. Furthermore, the thermal and mechanical loads of the thrust chamber can be significantly reduced by using the spiral channel. It should also be noted that the improvement of the cooling effect by adopting the spiral channel is at the cost of a greater coolant pressure drop. In this case, the coolant pressure drops of the spiral channel and straight channel are 0.071 MPa and 0.044 MPa, respectively. Therefore, the advantages and disadvantages of spiral channels should be considered comprehensively in practice.

4. Optimization Design

The structural parameters of the cooling channel have a significant influence on the regenerative cooling effect. It is necessary to optimize the spiral channel regenerative cooling structure to minimize the thermal load of the inner wall. This section begins with an introduction to the methods and principles used in optimizing the design. After that, the proposed method in Section 2.1 is employed to optimize the main parameters of the spiral channel.

4.1. Theory of Optimal Design

In order to obtain the structure of the spiral channel under the best cooling effect, the optimization design research is carried out by taking the channel width, channel height, pitch, and inner wall thickness as design variables and inner wall temperature as objective functions. For the cooling channel optimization, the response surface method is commonly used, because the relationship between design variables and design objective is highly nonlinear and cannot be represented by an explicit function. As an optimization method that combines experimental technology and applied mathematical statistics, the response surface method can obtain an accurate approximation function relationship in a local range through fewer experiments. A hypersurface response surface function is constructed by the least squares method. The typical second-order polynomial response surface approximate function is expressed as follows: where is the polynomial coefficient, is the number of design variables, and is the design variable.

For the purpose of obtaining a reliable response surface function, it is vital to select suitable sample points. The Latin hypercube experimental design method (LHD) is employed here to determine the initial sample points. This method can effectively reflect the overall distribution of random variables through sampling values. In addition, the method can better fit nonlinear relationships while possessing effective space-filling ability. The basic principle of LHD is that for mutually independent random variables , each random variable is divided into intervals. The sample extracted from the interval segment in random variable is as follows:

Referring to the optimization algorithm, the nonlinear programming by quadratic Lagrangian (NLPQL) is adopted because of the nonlinear and continuous characteristics of the cooling channel optimization problem. The algorithm assumes that the objective function is continuously differentiable, thereby expanding the objective function into a second-order Laplace equation and linearizing the constraints. Finally, the optimization calculation is transformed into a nonlinear quadratic programming problem to acquire the minimum value of the objective function. The second-order equation is improved by the quasi-Newton formula, and the stability of the algorithm is improved by linear search. The mathematical model of the optimization problem in this paper is written as where represents the objective function, represents the constraints, and represents the design variables.

4.2. Optimization Results and Analysis

An optimization design study is carried out in order to obtain the structural parameters of the spiral channel with the best cooling effect. A second-order response surface model is employed for fitting. The ranges of the design variables take into account the constraints among dimensions, the processing capacity of our laboratory, basic strength requirements, and other constraint conditions. Table 3 gives the specific value range of each parameter.

The heat transfer analysis method uses the proven 1D method proposed in Section 2.1. It is simple and efficient, which can improve the efficiency of the optimal design. In order to obtain a reliable response surface function, the Latin hypercube experimental design method is used to design the initial sample points. The number of collection points is 200, and the specific parameters of each point are shown in Table 4.

Figure 11(a) shows the effect of channel width and channel height on the maximum inner wall temperature when the pitch and the wall thickness of the thrust chamber are constant. It can be seen from the figure that the larger the channel width and the smaller the channel height, the lower the maximum inner wall temperature. The results could be explained that both the increase in channel width and the decrease in channel height reduce the flow area with other parameters unchanged. As a result, the coolant flowing velocity increases. According to equation (9), the convective heat transfer coefficient on the coolant side will increase with the increase of the flow velocity. In addition, it can be found that the effect of the channel height on the maximum temperature of the inner wall is less than that of the channel width. This is owing to the reduction of the channel height that will also generate an opposite effect; that is, the reduction of the channel effect coefficient in equation (8) leads to a decrease in heat dissipation.

Figure 11(b) depicts the effect of pitch and wall thickness on the maximum inner wall temperature when the channel width and channel height are kept unchanged. It is clear that the smaller the pitch, the smaller the wall thickness, and the lower the maximum temperature of the inner wall. When the wall thickness is reduced, on the one hand, the heat flux conducted by the inner wall is larger, and on the other hand, the temperature of the side wall surface of the coolant increases. Both factors lead to an increase in the convective heat transfer intensity of the coolant, which improves the cooling effect. When the pitch is reduced, the angle θ in equation (12) increases and the coolant velocity accelerates, resulting in a better cooling effect. In addition, we can find that the influence of the inner wall thickness on the maximum temperature of the inner wall is greater than that of the pitch.

The influence of the four independent variables of the spiral channel regenerative cooling on the maximum temperature of the inner wall is analyzed, as shown in Figure 12. Since the second-order response surface is used in this paper, there are 14 parameter items in total. In the figure, blue indicates a positive correlation, and red indicates a negative correlation. Most of them are positively correlated with temperature; that is, the higher the value, the higher the temperature. Obviously, the height and width of the cooling channel have the greatest influence on the maximum temperature of the inner wall, while the influence of the pitch is not obvious.

The optimization design is carried out using the nonlinear programming by quadratic Lagrangian (NLPQL) algorithm introduced in Section 4.1. The final optimization results are as follows: the channel height is 5.0 mm, the channel width is 0.8 mm, the pitch is 180.0 mm, and the inner wall thickness is 0.8 mm. Under this set of design parameters, the maximum temperature of the inner wall is 586.6 K. Compared with the calculated result of 835.5 K before optimization, the temperature has dropped by 29.8%. It should be noted that these optimization laws and results are only applicable to the thrust chamber in this paper. The prerequisites have a great influence on the results, such as the specific mass flow rate, equivalence ratio, chamber dimensions, the type of propellants, and the material of the thrust chamber. In conclusion, by carrying out the optimal design of the spiral channel, the correlation between the inner wall temperature and the main design parameters has been clarified. What is more, the cooling effect of the optimized spiral channel has been greatly improved, which is of great significance to its design improvement.

5. Conclusions

Although spiral channel regenerative cooling has been used increasingly, there is no effective heat transfer analysis method. Aiming at this problem, this paper proposes a simplified 1D inner wall temperature prediction method for spiral channel regenerative cooling. Based on the traditional semi-empirical formula, the method corrects the flow velocity of the coolant according to the characteristics of the spiral channel. The proposed method is verified by 3D numerical simulation based on a pressure-based coupled algorithm and a conjugate manner. The two methods are consistent in predicting the magnitude of the temperature reduction at the dangerous cross-section of the inner wall with an error of 1.95% in the test case. The characteristics of spiral channel regenerative cooling are analyzed by the proposed 1D method and the 3D numerical simulation, and the optimal design of the spiral channel is studied by the 1D method.

For the spiral channel regenerative cooling, the inner wall temperature first decreases and then increases along the axial direction, reaching the maximum value upstream of the throat. Then, it gradually decreases to its minimum at the gas outlet. Under the same structural parameters, the maximum temperature of the inner wall surface cooled by the spiral channel is 835.5 K, which is 8.5% lower than the maximum temperature of 913.9 K of the straight channel.

While using the 3D method for method validation, the temperature distribution of the entire wall solid domain is also analyzed. Regardless of the axial position, the maximum temperature of the inner wall is located in the center of the cooling channel. The temperature gradient of the inner wall and the rib along the radial direction is relatively large.

Taking channel width, channel height, pitch, and inner wall thickness as design variables and taking the inner wall temperature as the objective function, the response surface methodology is applied to obtain the optimal structure that achieves the best cooling effect. For the thrust chamber and its prerequisite in this paper, the final optimization results are as follows: the channel height is 5.0 mm, the channel width is 0.8 mm, the pitch is 180.0 mm, and the inner wall thickness is 0.8 mm. The calculated maximum temperature of the inner wall after the structure optimization is 586.6 K, which is reduced by 29.8%.

In short, this paper proposes an efficient heat transfer analysis method for spiral channel regenerative cooling, which can realize an accurate prediction of the maximum inner wall temperature. This method has strong applicability in the early stage of thrust chamber design.

Nomenclature

:Local cross-sectional area of nozzle (m2)
:Cross-sectional area of cooling channel (m2)
:Area at the throat (m2)
:Channel width (m)
:Rib width (m)
:Characteristic velocity (m/s)
:Specific heat (J/(kg∙K))
:Diameter at the throat (m)
:Hydraulic diameter (m)
:Function of the double arc nozzle profile (-)
:Channel height (m)
:Pitch of the spiral (m)
:Convective heat transfer coefficient in the hot-gas side (W/(m2∙K))
:Convective heat transfer coefficient in the coolant side (W/(m2∙K))
:Convective heat transfer coefficient in the coolant side without ribs (W/(m2∙K))
:Specific heat ratio of hot gas (-)
:Different axial lengths of nozzle (m)
:Mach number (-)
:Nusselt number (-)
:Chamber pressure (Pa)
:Prandtl number of hot gas (-)
:Prandtl number of the coolant (-)
:Prandtl number of the coolant on the wall (-)
:Total flow mass of coolant (Kg)
:Heat flux (W/m2)
:Local recovery coefficient (-)
:Effective recovery coefficient; also throat radius of curvature (-)
:Reynolds number (-)
:Curvature at different positions of nozzle (m)
:Adiabatic wall temperature of hot gas (K)
:Stagnation temperature of hot gas (K)
:Hot-gas-side wall temperature (K)
:Coolant side wall temperature (K)
:Coolant temperature (K)
:Time (s)
:Flow velocity of coolant in the spiral channel (m/s)
:Flow velocity of coolant in the straight channel (m/s)
: coordinate (m)
: coordinates of curvature center at different positions of nozzle (m)
: coordinate (m)
: coordinates of curvature center at different positions of nozzle [m]
: coordinate (m)
:Rib effect coefficient (-)
:Angle between the spiral channel and the thrust chamber axis (rad)
:Heat conductivity coefficient of the inner wall (W/(m2 K))
:Thermal conductivity of the coolant (W/(m2 K))
:Viscosity of hot gas (Pa∙s)
:Density of coolant (kg/m3)
:Dimensionless factor (-).

Data Availability

No underlying data was collected or produced in this study.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

This research was funded by the National Natural Science Foundation of China (grant number 11902013).