Theoretical and Experimental Heat Transfer in Solid Propellant Rocket Engine

ABSTRACT:

Accurate determination of heat flux is an important task not only in the designing aspect, but also in the performance analysis of rocket engines. In this purpose, this work deals with the heat flux determination in a combustion chamber through the inverse method. In this approach, the transient heat flux is determined from the experimental temperature data measured at the outer sidewall of the rocket engine. In this work the physical phenomenon was modeled by the transient one-dimensional heat equation in cylindrical coordinates and the material properties of the chamber were considered constant. Furthermore, the model is solved using the inverse heat conduction problem with least squares modified by the addition of Tikhonov regularization term of zero-order. Moreover, the sensitivity coefficients were obtained by Duhamel’s theorem. Through the regularization parameter, it was able to generate acceptable results even when using data with considerable experimental errors.

KEYWORDS:
Combustion chambers; Heat flux; Heat conduction; Ill-posed problems

INTRODUCTION

In a thrust chamber (nozzle and combustion chamber), the amount of energy transferred as heat to the chamber walls is between 0.5% and 5% of the total generated energy (Sutton 1992). Nevertheless, this amount could be enough to cause structural failure. Thus, to prevent the chamber and nozzle walls from failing, it is necessary to predict the heat flux accurately. Furthermore, accurate determination of heat flux is also important in the calculation of rocket engine performance and for the cooling system design.

Convective and radiative heat transfer must be determined for characterization of total heat flux. The convective heat transfer coefficient typically depends on many fluid physical properties. The computation of the radiation heat transfer must be accounted by the surfaces emissivity and the absorption and scattering coefficients of the fluid mixture. However, information about these parameters is not always easily found. Considering the propellant applied in this work (potassium nitrate with sucrose, KNSu), the combustion products that affects radiation heat transfer are predominantly water vapor, carbon dioxide and carbon monoxide, which are strong radiation emitters and absorbers (Incropera et al. 2007). In addition, the presence of liquid phase particles promotes radiation scattering (Incropera et al. 2007). Therefore, calculation of the radiation phenomenon is very complicated. Owing to the difficulties previously reported, the heat flux prediction has large uncertainties.

In order to minimize these uncertainties, the Inverse Heat Conduction Problem (IHCP), or Inverse Problem (IP), can be applied. The standard heat conduction problem (Direct Problems or DP) is concerned with the determination of the temperature distribution in the interior of the solid material. In contrast, the IP aims to obtain the boundary condition by using measured temperature history at one or further locations in a solid material (Beck et al. 1985). The outstanding advantage in using inverse problem is to avoid the computation of physical properties of gas combustion in order to obtain the heat flux from combustion (Kimura 1987).

The primary difficulty in solving inverse problems is associated to the fact that they are ill-posed. According to Hadamard, DP are well posed since the solution exists, is unique and stable under slight changes in input data. The problems that do not satisfy one or more of these conditions are referred to as ill-posed (Ozisik and Orlande 2000). In order to overcome the difficulty in solving ill-posed problems, a variety of analytic and numerical approaches have been proposed.

Stoltz (1960) and Burgraff (1964) proposed analytical solutions. Stoltz (1960) employed the Duhamel theorem to solve the heat equation for slab and spheres. Burgraff (1964) presented a solution based on series. Analytical solutions are restricted to linear problems, thus IP the technique was extended to nonlinear problems by the application of numerical methods. The first researcher to apply a nonlinear problem in IP was Beck (1970). In his work, the Fourier equation was solved using temperature dependent properties (specific heat and conductivity). Williams and Curry (1977) calculated the heat flux for the nonlinear case in materials of high and low thermal conductivity using the least squares method. Mehta (1981) employed the Finite Difference Method (FDM) to determine the heat flux in a rocket engine using inverse problem approach.

Tikhonov and Arsenin (1977) proposed the modified least squares method for a regularization term. The authors showed that this method is efficient in solving systems of ill-conditioned equations resulting from data with considerable experimental errors. The authors’ emphasis was performing an analysis of the stability of existing methods and the proposed method. The authors also discussed methods for determining optimal values of the regularization parameter. Alifanov and Egorov (1985) presented the conjugate gradient method. In this minimization method, the regularization is made during the iterative procedure in an implicit way.

Simpler techniques require less computational effort; however, they may not achieve satisfactory stability. On the other hand, sophisticated techniques can generate further accurate results, although the computational cost is usually higher. Currently, several methods are combined to take advantage of two or various approaches in order to obtain a competitive algorithm (hybrid).

Much work has been done to obtain heat flux in the nozzle, especially in the throat region (Mehta 1981; Kimura 1987; Patire Jr. 2010). Therefore, the purpose of this work is to predict the heat flux in a rocket combustion chamber using IP and least square method modified by the Tikhonov regularization parameter (Tikhonov and Arsenin 1997). The goal is to estimate the data of heat flux values in a rocket combustion chamber by using IP.

PROBLEM DEFINITION

The combustion chamber utilized in this work has a cylindrical geometry. Thus, the physical model can be described as a hollow cylinder of internal and external radius R_i and R_o, respectively. In Fig. 1 the heat flux q”(t) from the combustion is applied on the boundary surface at the internal radius R_i. It is reasonable to consider the external surface, at radius R_o, insulated. The experimentally determined temperature Y is measured by thermocouples located on the outer surface, as exposed in Fig. 1.

Figure 1
Geometry of the problem.

The assumptions adopted for this problem are one-dimensional transient heat transfer (in radial direction) and combustion chamber with constant physical properties.

The 1-D transient heat conduction equation in cylindrical coordinates is written as Eq. 1 (Ozisik 1993):

(1) $\mathit{\rho }{\mathit{c}}_{\mathit{p}}\frac{\partial \mathit{T}}{\partial \mathit{t}}=\frac{\mathit{k}}{\mathit{r}}\frac{\partial }{\partial \mathit{r}}\left(\mathit{r}\frac{\partial \mathit{T}}{\partial \mathit{r}}\right),$

where k is the thermal conductivity [W/(mK)]; ρ is density (kg/m³); c is specific heat [J/(kg K)]; T is the temperature (K), the dependent variable; and r (m) and t (s) are the independent variables.

The boundary condition on the outer surface is given by Eq. 2 (Ozisik 1993):

(2) ${\left.\frac{\partial \mathit{T}}{\partial \mathit{r}}\right|}_{\mathit{r}=\mathit{Ro}}=0$

The hypothesis of insulation on outer surface is appropriate since the heat transferred by natural convection is negligible compared to heat transferred from combustion gases.

The boundary condition on the inner surface is given by Eq. 3 (Ozisik 1993):

(3) ${\left.-\mathit{k}\frac{\partial \mathit{T}}{\partial \mathit{r}}\right|}_{\mathit{r}=\mathit{Ri}}={\mathit{q}}^{″}\left(\mathit{t}\right)$

where k is the thermal diffusivity (m²/s), defined as the ratio between the conductivity and density (ρ) versus specific heat (c).

The initial condition is the initial temperature T₀ (Eq. 4):

(4) $\mathit{T}={\mathit{T}}_0.$

DIRECT PROBLEM

As previously commented, DP is the standard procedure which is concerned with determination of temperature distribution inside a solid. All properties, boundary conditions and initial condition are prescribed, including the boundary condition on the inner surface q”(t). Thus, the DP can be readily solved by classical solution techniques.

The DP was solved numerically using the Finite Volume Method (FVM) (Versteeg and Malalasekera 1995). The physical domain was discretized using a uniform mesh (all elements has length Δr), as shown in Fig. 2:

Figure 2
Uniform one-dimensional mesh.

Figure 2 also illustrates both the internal and the boundary volumes for N elements. Each internal volume P has a neighbor to the left, W, and to the right, E. Similarly, each volume has faces to the left, w, and to the right, e.

The 1-D transient heat conduction equation in cylindrical coordinates, Eq. 1, was discretized using FVM, which provides a general equation:

(5) ${\mathit{a}}_{\mathit{P}}{\mathit{T}}_{\mathit{P}}={\mathit{a}}_{\mathit{w}}{\mathit{T}}_{\mathit{W}}+{\mathit{a}}_{\mathit{e}}{\mathit{T}}_{\mathit{E}}+{\mathit{b}}_{\mathit{P}}.$

The approximation implemented to obtain the coefficients and the source term, Eq. 5, to the internal volumes was Central Differencing Scheme (CDS) (Maliska 2004; Versteeg and Malalasekera 1995). The coefficients a_P, a_w and a_e and the source term b_P are defined as follows (Eqs. 6 to 9):

(6) ${\mathit{a}}_{\mathit{P}}=\frac{∆{\mathit{r}}^2{\mathit{r}}_{\mathit{P}}}{\mathit{\alpha }∆\mathit{t}}+{\mathit{r}}_{\mathit{e}}\mathit{\theta }+{\mathit{r}}_{\mathit{w}}\mathit{\theta },$

(7) ${\mathit{a}}_{\mathit{e}}={\mathit{r}}_{\mathit{e}}\mathit{\theta },$

(8) ${\mathit{a}}_{\mathit{w}}={\mathit{r}}_{\mathit{w}}\mathit{\theta },$

(9) ${\mathit{b}}_{\mathit{P}}={\mathit{T}}_{\mathit{E}}^0{\mathit{r}}_{\mathit{e}}\left(1-\mathit{\theta }\right)+{\mathit{T}}_{\mathit{W}}^0{\mathit{r}}_{\mathit{w}}\left(1-\mathit{\theta }\right)+{\mathit{T}}_{\mathit{P}}^0\left(\mathit{\theta }-1\right)\left[{\mathit{r}}_{\mathit{e}}+{\mathit{r}}_{\mathit{w}}\right]+\frac{∆{\mathit{r}}^2{\mathit{r}}_{\mathit{P}}}{\mathit{\alpha }∆\mathit{t}}{\mathit{T}}_{\mathit{P}}^0,$

where θ is the parameter that defines the formulation implemented in time. In this work, θ = 0.5 (Crank-Nicholson) was chosen.

The approximation implemented to obtain the coefficients and source term of the first volume, N = 1, was the Downstream Differencing Scheme (DDS) (Maliska 2004; Versteeg and Malalasekera 1995). The coefficients and source term are defined as follows (Eqs. 10 to 13):

(10) ${\mathit{a}}_{\mathit{P}}=\frac{∆{\mathit{r}}^2{\mathit{r}}_{\mathit{P}}}{{\mathit{r}}_{\mathit{e}}\mathit{\alpha }∆\mathit{t}}+\mathit{\theta },$

(11) ${\mathit{a}}_{\mathit{e}}=\mathit{\theta },$

(12) ${\mathit{a}}_{\mathit{w}}=0,$

(13) ${\mathit{b}}_{\mathit{P}}={\mathit{T}}_{\mathit{E}}^0\left(1-\mathit{\theta }\right)+{\mathit{T}}_{\mathit{P}}^0\left(\mathit{\theta }-1\right)+{\mathit{T}}_{\mathit{P}}^0\frac{∆{\mathit{r}}^2{\mathit{r}}_{\mathit{P}}}{{\mathit{r}}_{\mathit{e}}\mathit{\alpha }∆\mathit{t}}+\frac{{\mathit{q}}^{″}\left(\mathit{t}\right){\mathit{r}}_{\mathit{w}}∆\mathit{r}}{{\mathit{r}}_{\mathit{e}}\mathit{k}}.$

Finally, the approximation implemented to obtain the coefficients and the source term of the last volume, N, was the Upstream Differencing Scheme (UDS) (Maliska 2004; Versteeg and Malalasekera 1995). The coefficients and source term of the last volume are defined as (Eqs. 14 to 17):

(14) ${\mathit{a}}_{\mathit{P}}=\frac{∆{\mathit{r}}^2{\mathit{r}}_{\mathit{P}}}{{\mathit{r}}_{\mathit{w}}\mathit{\alpha }∆\mathit{t}}+\mathit{\theta },$

(15) ${\mathit{a}}_{\mathit{e}}=0,$

(16) ${\mathit{a}}_{\mathit{w}}=\mathit{\theta },$

(17) ${\mathit{b}}_{\mathit{P}}={\mathit{T}}_{\mathit{W}}^0\left(1-\mathit{\theta }\right)+{\mathit{T}}_{\mathit{P}}^0\left(\mathit{\theta }-1\right)+{\mathit{T}}_{\mathit{P}}^0\frac{∆{\mathit{r}}^2{\mathit{r}}_{\mathit{P}}}{{\mathit{r}}_{\mathit{w}}\mathit{\alpha }∆\mathit{t}}.$

The system of linear equations of Eq. 5 was solved using the TriDiagonal Matrix Algorithm (TDMA) (Maliska 2004). The DP was implemented using the FORTRAN 90 language with double precision. Ten volumes were employed to achieve the solution.

INVERSE PROBLEM

The mathematical formulation of IP is identical to the DP, except that the boundary condition in Eq. 3, q”(t), is unknown. The IP consists in the determination of the heat flux q”(t) taking into account the known physical properties and geometry data.

Ill-posed problems do not satisfy one or more Hadamard’s conditions of existence, uniqueness and stability. The existence and uniqueness of an IP are guaranteed by requiring that the inverse solution minimize the least squares norm rather than make it necessarily zero.

The measured temperature on the outer surface is represented by Eq. 18:

(18) $\begin{array}{cc}{\mathit{Y}}_{\mathit{j}}\left({\mathit{t}}_{\mathit{j}}\right),& \mathit{j}=1,2,\dots \mathit{M},\end{array}$

where t_j is the instant in which the temperature is measured; and M is the total number of measurements.

The estimated temperature on the outer surface is represented by Eq. 19:

(19) $\begin{array}{cc}{\hat{\mathit{T}}}_{\mathit{j}}\left({{\hat{\mathit{q}}}^{″}}_{\mathit{j}}\right),& \mathit{j}=1,2,\dots \mathit{M},\end{array}$

the symbol ^ denotes the estimated values.

In order to solve the IP, it is required that the estimated temperature on the outer wall (computed from the solution of the direct problem by the Fourier Equation) should match the measured temperatures as closely as possible over a specified tolerance. In order to perform that, the least square norm is minimized. Here, the least square norm is modified by addition of a zero-order regularization term (Ozisik 1993; Arsenin and Tikhonov 1977; Beck et al. 1985):

(20) $\mathit{S}\left({\hat{\mathbf{q}}}^{″}\right)={\sum _{\mathit{j}=1}^{\mathit{M}}\left[{\mathit{Y}}_{\mathit{j}}-{\hat{\mathit{T}}}_{\mathit{j}}\left({\hat{\mathbf{q}}}^{″}\right)\right]}^2+\mathit{\alpha }\sum _{\mathit{j}=1}^{\mathit{M}}{{\hat{\mathit{q}}}^{″}}_{\mathit{j}}^2,$

where S is the sum of squares; and α is the regularization parameter.

The traditional least squares on the right-hand side of Eq. 20 are added to a zero-order regularization term in order to reduce instability or oscillations that are inherent to the solution of ill-posed problems. The instability problem is solved by using the regularization term. The next step in the analysis is the minimization of the least square Eq. 20 by differencing it regarding each of the unknown heat flux components (Ozisik 1993). After some algebra, it is possible to write the solution for the heat flux in the matrix form:

(21) ${\mathbf{q}}^{″}={\left({\mathbf{X}}^{\mathit{T}}\mathbf{X}+\mathit{\alpha }\mathbf{I}\right)}^{-1}{\mathbf{X}}^{\mathit{T}}\left(\mathbf{Y}-{\mathbf{T}}_0\right),$

where X is the matrix of sensitivity.

Equation 21 is the solution of the inverse heat conduction problem. One can notice that the unknown surface heat flux depends on the regularization term α, the matrix of sensitivity X and the measured temperature data Y, only. No information a priori about the behavior of the heat flux on outer surface is needed.

The matrix of sensitivity is by definition the first derivatives of the dependent variable (i.e., temperature) with respect to the unknown parameter (i.e., heat flux), as exposed in Eq. 22:

(22) $\mathbf{X}\equiv \frac{\partial \mathbf{T}}{\partial {{\mathbf{q}}^{″}}^{\mathit{T}}}=\left[\begin{array}{cccc}\frac{\partial {\mathit{T}}_1}{\partial {{\mathit{q}}^{″}}_1}& \frac{\partial {\mathit{T}}_1}{\partial {{\mathit{q}}^{″}}_2}& \cdots & \frac{\partial {\mathit{T}}_1}{\partial {{\mathit{q}}^{″}}_{\mathit{M}}}\\ \frac{\partial {\mathit{T}}_2}{\partial {{\mathit{q}}^{″}}_1}& \frac{\partial {\mathit{T}}_2}{\partial {{\mathit{q}}^{″}}_2}& & \frac{\partial {\mathit{T}}_2}{\partial {{\mathit{q}}^{″}}_{\mathit{M}}}\\ ⋮& & & \\ \frac{\partial {\mathit{T}}_{\mathit{M}}}{\partial {{\mathit{q}}^{″}}_1}& \frac{\partial {\mathit{T}}_{\mathit{M}}}{\partial {{\mathit{q}}^{″}}_2}& \cdots & \frac{\partial {\mathit{T}}_{\mathit{M}}}{\partial {{\mathit{q}}^{″}}_{\mathit{M}}}\end{array}\right],$

The matrix of sensitivity represents the changes of dependent variable with respect to the changes of the unknown parameter. Once the direct problem can be promptly solved by classical solution techniques, the sensitivities are easily determined. For this study, the direct problem was solved numerically through the Duhamel’s theorem (Ozisik 1993).

STABILITY ANALYSIS

To examine the accuracy of the predictions by the inverse analysis, a study case is proposed. Consider that a known heat flux having a triangular shape, Eq. 23, is applied on the inner surface of hollow cylinder of Fig. 1.

(23) ${\mathit{q}}^{″}\left(\mathit{t}\right)=\left\{\begin{array}{cc}5\mathit{t}\hfill & 0\le t\le 1\hfill \\ -5\mathit{t}+10\hfill & 1<t\le 2\hfill \end{array}\right..$

The physical properties of the cylinder material are listed in Table 1.

Table 1
Physical properties hypothetic problem.

Parameter	Property
Type of material	Aluminum ABNT 6061-T5
Initial temperature	300 K
Specific mass	2700.0 (kg/m3)
Specific heat	896.0 (J/kg.K)
Thermal conductivity	180.0 (W/m.K)
External radius	40.0 (mm)
Internal radius	30.0 (mm)

The measured data is simulated by solving the direct problem using the known heat flux. Thus, by solving the direct problem the temperature over the whole domain can be achieved, including the temperature on the outer surface. Then the measured data is generated by introducing a random error ω, as follows (Eq. 24):

(24) $\mathit{Y}={\mathit{T}}_{\mathit{exact}}+\mathit{\omega },$

where ω has values between -5.0 and 5.0. The exact value of the surface heat flux is shown in Eq. 23.

Figure 3 shows the temperature generated by introducing random errors, Y, and the temperature in the inner wall, T_w. The temperature simulated Y, i.e. the outer wall, is utilized in the code program developed to estimate the heat flux q”.

Figure 3
Temperature in the outer and inner wall in the study case.

Figure 4 presents the heat flux q” estimated from Eq. 21 using different values of regularization parameters α. Applying values for regularization parameter such as 1.10^-12, the solution exhibits oscillatory behavior, whereas for α = 1.10^-10 the solution is damped and deviates from the exact result. The regularization parameter that approximates the solution satisfactorily to exact result is 1.10^-11.

Figure 4
Effect of the regularization parameter α on the stability of inverse solution.

Figure 5 shows a diagram of the computational procedure. The first step is to read the data of temperature in relation to time, Y(t_j). Temperature data is fed just once. Other data are: physical properties, internal and external diameter of combustion chamber, and the sensitivity coefficient.

Figure 5
Diagram of the computational procedure.

After that reading data, the computer code will calculate the sensitivity matrix, the transposed sensitivity matrix and obtain identity matrix.

The heat flux is calculated using Eq. 21. The heat flux calculation consists in a single step, therefore convergence procedures are not necessary.

LITERATURE COMPARISON

In this subsection, analytical solutions from Kacynski et al. (1987) and Burgraff (1964) are compared to the methodology presented in this work.

Burgraff (1964) presents a solution for the heat flux based in expansion series. His correlation is available for slab, solid cylinder and solid sphere. Thus, to accomplish the comparison, the solutions for the current model and Kacynski (1987) were computed in rectangular coordinates. The same material and physical properties listed in Table 2 were considered. The temperature data considered was the T_exact, obtained by solving the direct problem. The exact values for temperature must be chosen due to the analytical solutions being ill-posed and do not deal with measurement errors.

Table 2
Diameter of MTP engines.

Engine	Le (mm)	Di (mm)	De (mm)
NRB25/26/27	180.0	38.1	44.45
U25/26	280.0	44.46	50.8

The correlation developed by Kacynski et al. (1987) is an analytical solution significantly simplified by the hypothesis that the temperature derivative over time is independent of the spatial derivative over time (Eq. 25). This correlation was developed for a slab considering the same assumptions:

• Constant properties;

• Linear rate of temperature rise;

• Adiabatic outer wall;

• No axial or circumferential conduction;

derivative of temperature with respect to time independent of the position.

(25) ${\mathit{q}}^{″}=\frac{\mathit{k}}{\mathit{\alpha }}\frac{\mathit{dY}}{\mathit{dt}}{\mathit{x}}_{\mathit{i}}\left(\frac{{\mathit{x}}_{\mathit{o}}}{{\mathit{x}}_{\mathit{i}}}-1\right)$

where x_o (mm) corresponds to the outer wall of combustion chamber, R_o, and x_i (mm) corresponds to the inner wall, R_i.

The derivative of temperature with respect to time was computed using DDS (Eq. 26) (Maliska 2004).

(26) $\frac{\mathit{dY}}{\mathit{dt}}=\frac{{\mathit{Y}}_{\mathit{j}}-{\mathit{Y}}_{\mathit{j}-1}}{∆\mathit{t}}$

Figure 6 presents the heat fluxes obtained by Kacynski et al. (1987) correlation by the current model, and Burggraf (1964) analytical solution.

Figure 6
Comparison between Kacynski et al. (1987), Burggraf (1964) analytical correlation, and the solution for IP of this work

One can note good agreement between the current work and Burgraff (1964) solution. The Kacynski et al. (1987) solution presents a delay, possibly due to the assumption of temperature derivative with respect to time and not with respect to the position.

GRID INDEPENDENCY

Consider the same hypothetic problem approached in the Stability Analysis section. To evaluate the grid independency, the hypothetic problem was solved with different quantity of volumes in radial direction (N = 10, 100 and 1000) and a fixed number of time advances (M = 100).

Figure 7 shows the results for the different quantity of volumes. It is possible to note that the result is independent on the quantity of volumes.

Figure 7
Effect of quantity of volumes for hypothetic problem with M = 100 (fixed).

To evaluate the effect of time advance, the hypothetic problem was solved using different time advances (M) and a fixed number of volumes (N = 10). Figure 8 shows the results have no influence in the results.

Figure 8
Effect of the quantity of time advances for hypothetic problem with N = 10 (fixed).

EXPERIMENTAL FACILITIES

The static tests were performed at the Hydraulic Machines Laboratory (LMH) ate Federal University of Paraná (UFPR). For the tests, the engine was fixed horizontally on a test bench, which allows the measurement of the thrust through strain gauges. Outer wall temperatures were measured by type K thermocouples fixed on the outer wall of engine. According to the manufacturer, the thermocouple response time is under 0.3 s. In order to evaluate the effect of acquisition frequency in the results, tests were performed employing 200 Hz, 2 Hz and 1 Hz. The data acquisition system employed was HBM Spider 8 and Catman 4.5 software.

Space Activities Laboratory (LAE) of UFPR developed Netuno-RB (NRB) and Urano (U) with the purpose of studying propellant performance parameters. The engines are made of Aluminum 6063-T6, cylindrical in shape. The noteworthy dimensions of NRB and Urano engines are shown in the Table 2. These dimensions are indicated in Fig. 9.

Figure 9
Essential parts of NRB and Urano engine and thermocouples location.

The physical properties of 6063-T6 employed in this work (ASM 2016) are listed in Table 3.

Table 3
Physical properties of Netuno-RB and Urano engine.

NRB and U	Values
Thermal conductivity (W/m.K)	209.0
Specific heat (J/kg.K)	900.0
Density (kg/m3)	2690.0

The total length without nozzle and cap for NRB and Urano is 180 and 280 mm, respectively. Figure 9 illustrates the NRB or Urano rocket engine and its parts. Smectite is a non-reactive substance deposited between cap and casing to seal the propellant grain in order to avoid gas losses during burning.

Powder of smectite (approximately 0.1 kg) was placed in the chamber casing and conformed by mechanical press. Afterwards, powder propellant was placed in the engine and similarly conformed by mechanical press. The propellant apllied was potassium nitrate (KNO₃)/sucrose (C₁₂H₂₂O₁₁), known as KNSu. The composition was 65 wt% of KNO₃ and 35 wt% of C₁₂H₂₂O₁₁.

The tests were performed without nozzle. The position of thermocouples on the surface of the combustion chamber is presented in Fig. 9. In Table 4, the distance T4 and T5 are showed in mm. Thermocouple T4 placed on the outer surface of engine is positioned next to propellant grain and thermocouple T5 is positioned on the outer surface next to nozzle entrance.

Table 4
Position of thermocouples during the static test.

Engine	T5 (mm)	T4 (mm)
NRB25/26/27	15	20
U25/26	20	40

Five tests were performed. In all tests the thermocouples were positioned in T4 and T5 position. NBR25, NBR26 and NBR27 were tested using 1 Hz, 2 Hz and 200 Hz acquisition frequency, respectively. U25 and U26 were tested using 1 Hz and 200 Hz. The temperature history of thermocouples T4 and T5 are presented in Fig. 10. The behavior of temperature history is very similar, therefore only the NBR25 temperature history is shown.

Figure 10
Experimental temperature data of a) T4 and b) T5 for NRB25.

RESULTS AND DISCUSSION

Table 5 summarizes the tests performed. Applying the methodology discussed in the previous sections, the heat fluxes were obtained for thermocouple T4 and T5 for each test.

Table 5
Position of thermocouples during the static test.

Test	Engine	T5 (mm)	T4 (mm)	Data acquisition frequency (Hz)
1	NRB25	15	20	1
2	NRB26	15	20	2
3	NRB27	15	20	200
4	U25	20	40	1
5	U26	20	40	200

Figures 11 to 14 show the heat fluxes obtained for all tests performed. The heat fluxes acquired from the thermocouple T4 for tests 1, 2 and 3 are presented in Fig. 11. The heat fluxes acquired from the thermocouple T5 for tests 1, 2 and 3 are presented in Fig. 12. In the same vein, the heat fluxes from T4, for U25 and U26 are presented in Fig. 13 and the heat fluxes for U25 and U26 engines are presented in Fig. 14.

Figure 11
Heat flux from temperature history of thermocouple T4, of NBR25, 26 and 27.

Figure 12
Heat flux from temperature history of thermocouple T5, of NBR25, 26 and 27.

Figure 13
Heat flux from temperature history of thermocouple T4, U25 and U26.

Figure 14
Heat flux from temperature history of thermocouple T5, of U25 and U26.

The firing duration of NRB engines was on average 3 s. The firing duration of Urano engine was on average 4 s. The heat flux peaks of the tests were superposed.

One can note that the results are very similar despite different data acquisition frequency of the tests. This can be regard to the fact that current methodology is based in derivative of temperatures, that is, the temperature variation in respect to the heat flux, or sensitivities. A peak of heat in Fig. 14 is probably related to iteration error due to a greater amount of math operations.

ERROR ANALYSIS

The results of the present work are subjected to errors from: the mathematical model, experimental measurement and numerical solution. Errors from the mathematical model arise from hypotheses that simplify the physical phenomenon; i.e., the assumption of heat transferring in one-dimension, constant properties, and insulated outer wall. The error from the mathematical model is very difficult to obtain since the actual heat flux value is unknown. An approximation to the actual heat flux could be obtained if heat flux was measured precisely. In this case, specifically, the experimental errors should be estimated.

Experimental errors in this work come from temperature measurement. It is expected to contribute toward the majority of total error (sum of errors from mathematical model, experimental measurement and numerical solution). The calibration reports gave an approximated idea of range of measurement errors. The thermocouples employed in this work have measurement errors almost linear, of 0.4 K at 273 K and 1.4 K at 393 K.

Furthermore, numerical errors are also present in the solution. However, it is expected to present a minor contribution to the total error.

In order to estimate the numerical error, the measured temperature should not be included in simulations. Thus, only the direct problem was considered in this analysis. Table 6 shows an estimative of discretization error U of the direct problem, obtained through Richadson estimator (Roache 1998) based in the apparent order Pu. These results were obtained by simulation of the hypothetic problem, then, the boundary conditions and the physical properties are known as described in Stability Analysis. In order to achieve an estimative to discretization error, the direct problem was considered. Therefore no experimental errors were embedded in the results. In addition, the results present no iteration error due to the linear system of equations is solved by TDMA.

Table 6
Estimative of discretization error (U) in the direct problem (D.P.)

ϕ	∆t	M	Numerical solution (K)	Pu	U
1	0.2	10
2	0.1	20	343.833821320819
3	0.05	40	337.220558396801
4	0.025	80	333.901761686173	0.994701998687	-3.34330689
5	0.0125	160	332.238022579065	0.996231040865	-1.67246615
6	0.00625	320	331.404873509277	0.997782664935	-0.83571599
7	0.003125	640	330.987952798396	0.998801609461	-0.41761421
8	0.001563	1280	330.779402491387	0.999377604961	-0.20873037
9	0.000781	2560	330.675104416838	0.999682911583	-0.10434394
10	0.000391	5120	330.622949594699	0.999839971494	-0.05216639
11	0.000195	10240	330.596870730558	0.999919613163	-0.02608177
12	9.77E-05	20480	330.583830934388	0.999959716249	-0.01304052

As expected, apparent order converges to 1 due to approximations for derivatives are of first order. Apparent order is obtained from Eq. 27:

(27) $\mathit{Pu}=\frac{\log \left(\frac{{\varphi }_2-{\varphi }_1}{{\varphi }_3-{\varphi }_3}\right)}{\log \left(\mathit{r}\right)}$

where ϕ₁ is the variable in the most coarse mesh; ϕ₂ is the variable in the medium mesh; and ϕ₃ is the variable in the fine mesh. The variable ϕ implemented in error determination was the temperature in the outer surface in the instant 1.0 s. Grid refinement was applied in time variable by increasing M volumes. The number of volumes in space discretization was always the same in all simulations (10 volumes).

Factor r in Eq. 28 is the ratio between the coarse grid and medium mesh or medium mesh and fine mesh. The ratio r in this analysis was 2.

(28) $\mathit{r}=\frac{{\mathit{h}}_1}{{\mathit{h}}_2}=\frac{{\mathit{h}}_2}{{\mathit{h}}_3}$

Richadson estimator (Roache 1998) based in the apparent order Pu is defined in Eq. 29:

(29) ${\mathit{U}}_{\mathit{RI}}=\frac{\left({\varphi }_3-{\varphi }_3\right)}{\left({\mathit{r}}^{\mathit{Pu}}-1\right)}$

Figure 15 shows the error estimation (U) in relation to variation of time length. ∆t. The estimated error diminishes as the grid refinement in time variable increases.

Figure 15
Estimative of discretization error of direct problem.

Theoretical analysis can also be worthwhile to compare the theoretical heat flux and the heat flux obtained by the inverse approach. The hypothesis in theoretical heat flux calculation in this analysis is: radiation heat transfer is minimum and turbulent flow.

Convection heat transfer coefficient, h, was calculated through Gnielinski correlation for turbulent flow in a cylindrical duct section as follows (Eq. 30) (Incropera et al. 2007):

(30) $\mathit{Nu}=\frac{\mathit{hD}}{\mathit{k}}=\frac{\left(\raisebox{1ex}{$\mathit{f}$}\!\left/ \!\raisebox{-1ex}{$8$}\right.\right)\left({\mathit{Re}}_{\mathit{D}}-1000\right)\mathrm{Pr}}{1+12.7{\left(\raisebox{1ex}{$\mathit{f}$}\!\left/ \!\raisebox{-1ex}{$8$}\right.\right)}^{1/2}\left({\mathrm{Pr}}^{2/3}-1\right)};\left\{\begin{array}{c}3000<{\mathit{Re}}_{\mathit{D}}<5.{10}^6\hfill \\ 0,5<\mathrm{Pr}<2000\hfill \end{array}\right.$

where Nu is Nusselt number; Re is the Reynolds number; Pr is the Prandtl number; D is the internal duct diameter (m); k is the thermal conductivity of the gases, W/(mK); and f is the Moody friction factor.

Moody friction factor is described by Colebrook-White correlation (Eq. 31):

(31) $\frac{1}{\sqrt{\mathit{f}}}=-2{\log }_{10}\left(\frac{\mathit{e}}{3,7\mathit{D}}+\frac{2,51}{\mathit{Re}\sqrt{\mathit{f}}}\right)$

where e is the equivalent roughness (m), in this case considered 1.10^-4 m.

Reynolds number is defined as follows (Eq. 32):

(32) ${\mathit{Re}}_{\mathit{D}}=\frac{4\dot{\mathit{m}}}{\mathit{\pi }\mathit{D}\mathit{\mu }}$

where ṁ is the mass flow rate (kg/s).

Table 7 shows the data achieved for convection heat flux coefficient calculation. The third column shows the source, that is, how that data was attained. Mass flow rate was calculated through the ratio between the mass of propellant and the burn time in NRB25 test. Through Propep software (ROCKETWORKBENCH 2016) viscosity and Prandtl number from chemical composition of propellant were obtained. The reference utilized to determine the thermal conductivity of combustion gases and the ratio between roughness (e) and tube diameter (D) was Incropera et al. (2007).

Table 7
Data for convection heat flux coefficient calculation.

Data	Value	Source
ṁ	5.0029.10-2 (kg/s)	NRB25 burn data
µ	1.0465.10-4 (Pa.s)	Propep
Pr	0.7	Propep
e/D	0.00262	Reference
k	0.056 [W/(mK)]	Reference

Results for Reynolds and Nusselt number are showed in Table 8. Moreover, the convection heat transfer was calculated through Eq. 33:

(33) $\mathit{q}=\mathit{h}\left({\mathit{T}}_{\mathit{s}}-{\mathit{T}}_{\infty }\right)$

Table 8
Results.

Data	Value
Re	15,976.0
Nu	143.36
h	210.72 [W/(m2K)]
q"	271.4 (kW/m2)

where T_s is the temperature of combustion chamber before test (19 °C); and is the temperature of the hot gases (1302.47 °C), given by Propep software.

Heat flux by convection is listed in Table 8. This value met heat flux obtained trough inverse approach acceptably.

It is reasonable to note that the total error contributes to the unstable behavior inherent of ill-posed problems. Thus, it is important to minimize other sources of errors (not only measurement errors) that could lead to difficulties in the regularization parameter estimation.

Finally, the method applied in this work, least squares, is of maximum likelihood for measurements errors for the following assumptions: uncorrelated errors, constant deviation and normal distribution. These assumptions were approximately met.

CONCLUSION

A methodology for heat flux determination using experimental temperature history was tested. This methodology for solution of the inverse problem is computed from a linear equation; consequently, the solution is given directly. Therefore, the methodology in this work is not an iterative process, resulting in a simple program code. The essential advantage of the method is to avoid the calculation of the physical properties of the gases from combustion. Another point to emphasize is that the methodology presented is not very sensitive to experimental errors.

In order to verify the programming code for solving IP a hypothetic problem was suggested. A program code was developed for solving IP generated results, which were compared to the results accomplished through analytical correlations from the literature. The results agreed satisfactorily. Furthermore, the methodology was applied to achieve heat flux from temperature data of static tests accomplished at UFPR. The maximum heat flux of NRB was estimated in 280 kW/m² from thermocouple T4 or T5. For Urano, the maximum heat flux was estimated in 300 kW/m² from thermocouple T5. The heat flux calculated through theoretical approach showed an acceptable agreement in relation to inverse method approach. In addition, the results indicate that acquisition frequency does not have any effect in heat flux estimation.

• FUNDING
  Conselho Nacional de Desenvolvimento Científico e Tecnológico (CAPES) [https://doi.org/10.13039/501100003593] Grant No: Finance Code 001, Project n. 20 of Edital Pró-Estratégia 50/2011
  Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CNPq) [https://doi.org/10.13039/501100002322] Grant No: 309365/2013-9

ACKNOWLEDGMENTS

The authors would acknowledge the “UNIESPAÇO Program” of The Brazilian Space Agency (AEB), CAPES (Coordination for the Improvement of Higher Education Personnel, Brazil) and CNPq (National Counsel of Technological and Scientific Development, Brazil) by physical and financial support given for this work. The authors would also acknowledge the journal editor and reviewers by their suggestions and corrections.

REFERENCES

• Alifanov OM, Egorov YV (1985) Algorithms and results of solving the inverse heat conduction boundary problem in two-dimensional formulation. Journal of Engineering Physics 48(4):489-496. https://doi.org/10.1007/BF00872080
  » https://doi.org/10.1007/BF00872080
• [ASM] Aerospace Specification Metals (2016) Aluminun 6063-T5. Aerospace Specification Metals; [accessed 2018 June 10]. http://asm.matweb.com/search/SpecificMaterial.asp?bassnum=MA6063T5
  » http://asm.matweb.com/search/SpecificMaterial.asp?bassnum=MA6063T5
• Beck JV (1970) Nonlinear estimation applied to the non linear inverse heat conduction problem. Journal Heat Mass Transfer 13(4):703-716. https://doi.org/10.1016/0017-9310(70)90044-X
  » https://doi.org/10.1016/0017-9310(70)90044-X
• Beck JV, Blackwell B, St. Clair Jr. CR (1985) Inverse heat conduction-Ill-posed Problems. New York: Wiley Interscience.
• Burggraf OR (1964) An exactsolution of the inverse problem in heat conduction theory and applications. Journal of Heat Transfer 86(3):373-380. https://doi.org/10.1115/1.3688700
  » https://doi.org/10.1115/1.3688700
• Incropera FP, Dewitt DP, Bergman TL, Lavine AS (2007) Fundamentals of heat and mass transfer. Hoboken: John Wiley and Sons.
• Kacynski JK, Pavli A, Smith T (1987) Experimental evaluation of heat transfer on a 1030:1 area ratio rocket nozzle. (NASA-TP-2726). NASA Technical Paper.
• Kimura AL (1987) Transferência de calor em motor-foguete (PhD Thesis). São José dos Campos: Instituto Tecnológico de Aeronáutica. In Portuguese.
• Maliska CR (2004) Transferência de Calor e Mecânica dos Fluidos Computacional. 2nd ed. Rio de Janeiro: Livros Técnicos e Científicos.
• Mehta RC (1981) Estimation of heat transfer coefficient in a rocket nozzle. AIAA Journal 19(8):1085-1086. https://doi.org/10.2514/3.7846
  » https://doi.org/10.2514/3.7846
• Ozisik MN (1993) Heat Conduction. 2nd ed. Hoboken: John Wiley & Sons.
• Ozisik MN, Orlande RB (2000) Inverse heat transfer: fundamentals and applications. New York: Taylor & Francis.
• Patire Jr. H (2010) Determinação da carga térmica na parede da câmara de empuxo de propulsores bipropelente usando método inverso (PhD Thesis). São José dos Campos: Instituto Tecnológico de Aeronáutica. In Portuguese.
• Roache PJ (1998) Verification and validation in computational science and engineering. Albuquerque: Hermosa.
• ROCKETWORKBENCH (2016) Cpropep-Web. Rocketworkbench; [accessed 2018 October 30]. http://rocketworkbench.sourceforge.net/equil.phtml
  » http://rocketworkbench.sourceforge.net/equil.phtml
• Stolz G (1960) Numerical solutions to an inverse problem of heat conduction for simple shapes. Journal of Heat Transfer 82(1):20-25. https://doi.org/10.1115/1.3679871
  » https://doi.org/10.1115/1.3679871
• Sutton GP (1992) Rocket Propulsion Elements. 6th ed. New York: John Wiley & Sons.
• Tikhonov AN, Arsenin YV (1977) Solutions of ill-posed problems. Washinghton: Winston & Sons.
• Versteeg H and Malalasekera W (1995) An introduction to computational fluid dynamics: the finite volume method. London: Prentice Hall.
• Williams SD, Curry DM (1977) An Analytical and Experimental Study for Surface Heat Flux Determination. Journal of Spacecraft and Rockets 14(10):632-637. https://doi.org/10.2514/3.27987
  » https://doi.org/10.2514/3.27987

Edited by

Section Editor: T John Tharakan

Publication Dates

• Publication in this collection
  26 Aug 2019
• Date of issue
  2019

History

• Received
  06 Feb 2018
• Accepted
  30 Oct 2018