Abstract
This paper designed a platelet heat exchanger in the solar thermal thruster and analyzed the unsteady-state conjugate heat transfer characteristics between heat exchanger and propellant. The conjugate heat transfer (CHT) computational fluid dynamics (CFD) simulation of the 3D model of the platelet under steady-state conditions was carried out with different mass flow rates to find the empirical correlation between the average Nusselt number and the average Reynolds number. The unsteady-state 1D simplified model of the heat exchanger was established using a loose coupling algorithm based on quasi-steady flow domain and finally verified by experiments. The results show that the platelet structure could heat the working medium to more than 2380 K with the heat transfer efficiency of 87% and produce a peak thrust of 0.57 N and specific impulse of 2200 m/s; in steady state, the outlet temperature and heat transfer efficiency of the heat exchanger were stable at 1950 K and 69%. Moreover, 1D model could accurately reflect the real heat exchange situation to a certain extent, the simulation error was less than 5% compared with the 3D model, and the calculation time was greatly shortened, making it more convenient to adjust the heat exchange strategy. The experimental results were consistent with the simulation results at the initial stage of heat exchange, and the difference was mainly reflected in the steady-state stage, which might be caused by the lack of precision of the experimental equipment.
1. Introduction
In recent years, human beings have higher requirements for the propulsion performance of satellites with the more complex space tasks. The chemical propulsion system can provide large thrust in an instant, but its specific impulse is only 100~400 s and requires a lot of fuel [1]. Electric propulsion is the research hotspot at present, which is used in orbit transfer, formation flight, and on orbit attitude maintenance. Although its specific impulse can reach more than 1000 s, the typical thrust is relatively low, generally 1 μN~400 mN, and it means more time will be taken to complete the task at the mean time [2–4]. The 40 cm ion thruster developed by NEXT (NASA’s evolutionary xenon thruster) had a thrust of only 364 mN at a power of 10.5 kW [5]. A miniature vacuum cathode arc thruster was designed with a specific impulse of 1571 s and a thrust to power ratio of 16.3 μN/W [6]. A very low power cylindrical Hall thruster for nanosatellite “PROITERES-3” under development in Osaka Institute of Technology was designed with a specific impulse of 1570s and a thrust of 2.87mN [7]. Although the technology of cold gas propulsion is mature, it does not make full use of the properties of working medium, so the thrust and specific impulse are correspondingly low, which provided the thrust of 70~1200 mN [8, 9]. By contrast, the solar thermal propulsion system can obtain a large velocity increment due to its high specific impulse and moderate thrust, which can make up for the gap in space transfer tasks [10–12].
The principle of solar thermal propulsion (STP) is to heat the propellant through the sunlight gathered by the concentrator and then make the propellant expand through the Laval nozzle to generate thrust [13, 14]. The heat exchanger is the most important structure in the solar thermal thruster. The wall of the heat exchanger accumulates heat through radiation heat exchange with the concentrator and then transfers the heat to the propellant. Therefore, the efficiency of heat exchanger determines the gas outlet temperature and thrust. Most of the early heat exchangers used spiral channel structure, which has the characteristics of simple structure and easy processing. Markopoulos P. et al. [15] designed a direct endothermic thrust chamber. In the experiment, H2 was used as propellant, the inlet working medium mass flow was 0.25 g/s, the temperature was 333 K, the thrust chamber pressure was 0.1 MPa, and the temperature after heating reached 2264 K. Although the spiral heat exchanger could heat the working medium to the ideal temperature, the length of the spiral channel had reached 6 m. The increase of the length of the chamber meant the increase of the mass and the decrease of the performance of the thruster. Therefore, it is necessary to design a more efficient heat exchanger. The platelet structure is mainly used for transpiration cooling in rocket engine [16, 17]. Xing B.Y. et al. designed a multilayer platelets heat exchanger in STP system based on the transpiration cooling technology of platelet [18]. In their study, the number of platelets was 20 layers, the thickness of a single platelet was 1 mm, and the depth of the control channel was 0.1 mm, and the heat transfer area of the propellant between the platelets was about 5 ~ 10 times larger than that of the spiral channel under the same conditions. The simulation results showed that when the wall temperature of the platelet structure reached 2400 K, the propellant could be heated to about 2300 K, which proved the feasibility of the structure.
In fact, the metal wall temperature cannot be maintained above 2400 K all the time due to the coupling effect. The intermittent heating strategy should be adopted to keep the working medium at a high temperature, hence the need for an analysis of the unsteady heat transfer of the heat exchanger. However, the amount of calculation of 3D CFD simulation is huge, and the calculation time may reach tens or even hundreds of hours according to the difference of calculation accuracy and time steps. Therefore, researchers have developed solutions with different calculation efficiency and accuracy for different problems. The transient tight coupling method can get effective results, but it consumes a lot of computing resources and is difficult to solve practical complex engineering problems [19–21]. The loose coupling method based on the quasi-steady flow field considers that in the whole process of fluid-solid coupling heat transfer, the flow field is in several quasi-steady states, and each quasi-steady flow field is solved by the steady-state Navier-Stokes equations; the research results show that the algorithm can greatly improve the calculation efficiency, but the deviation of the calculation results is large, which is mainly due to the large deviation between the treatment method of completely isolating the flow field from other parts and the actual coupling relationship [22–24]. De Giorgi and Fontanarosa [25] proposed a novel quasi-one-dimensional model for the performance estimation of a vaporizing liquid microthruster (VLM); in their study, the performance of 1D model well agreed with the experimental data, with a maximum estimated error of 7.3% on the thrust and the specific impulse.
To solve these problems, this paper designed a new platelet structure based on the temperature requirements considering the continuity requirements of propellant and the difficulty of manufacturing process. According to the symmetry of the structure, the platelet was simplified to 1D model. A new loose coupling algorithm for global transient tight coupling heat transfer based on quasi-steady flow field was adopted, and the steady-state algorithm was used to update the flow field and solve the energy equations of fluid and solid. The results of 3D CFD simulation and 1D simulation were compared and analyzed, and an experimental device was built to verify the above results.
2. Model of Platelet Heat Exchanger
The structure of STP is shown in Figure 1. It consists of four parts: These are refractive concentrator, platelet heat exchanger, insulation layer, and Laval nozzle, respectively. Distributed passage is used to increase the heat exchange area, and metering passage is used to accelerate the flow of propellant. Figure 2 shows the details of platelet heat exchanger, there are 9 layers of platelets, the thickness of each layer is 2 mm, and the gap between layers is 2 mm. The inner and outer diameters of the platelet heat exchanger are 32 mm and 50 mm, respectively; the vertical height of the heat exchange core is 40 mm; and the length of distributed passage and metering passage is 3 mm and 6 mm, respectively. Each layer has 8 metering passage in evenly distributed around the heat exchange core to accelerate the flow of working medium, and the diameter of single metering passage is 0.5 mm.
(a) Front view
(b) Sectional view
(c) Upper view
Under the high-temperature condition, the heat resistance of thruster materials should be first considered. Because of the difficulty of processing high-temperature-resistant ceramic materials and poor sealing, high-temperature-resistant single crystal molybdenum was selected as the main material of thrust chamber [26], and the outside of the heat exchanger was wrapped with a layer of thermal insulation material in order to protect other parts, which was made of alumina fiber. The main thermal parameters are shown in Table 1.
3. 3D Unsteady-State CFD Simulation of Conjugate Heat Transfer
Owing to the axisymmetric structure of the platelet, a single unit could be taken for calculation during simulation. The unit model and grid division are shown in Figure 3, the circumferential angle of the unit was 22.5°, structured mesh generation was adopted, and densification was carried out at the metering passage, and the total number of grids was 47060. Since the minimum characteristic length of the heat exchanger was 0.5 mm, which belonged to the category of microchannel, the Knudsen number of the propellant must be calculated to judge whether it met the continuity condition. Nitrogen was used as propellant because of the safety, and its molecular average free path is . According to Knudsen number calculation formula, , which meant the continuum model was applicable.
(a) Unit
(b) Meshing
Actually, the heat transfer of the platelet heat exchanger was complex, including the radiation heat transfer between the outer wall and the concentrator, the convective heat transfer between the inner wall and nitrogen, and the heat conduction of itself; thus, it was essential to simplify the model. The effect of solar radiation on platelet can be expressed in the form of flux [27]: where is the wall temperature, is the thermal conductivity of platelet, is the radiation area of platelet, is the radiation intensity of solar, is the concentrated ratio, is the open area of concentrator, is the concentrator efficiency, is the Stefan-Boltzmann constant, is the radiant emissivity of heat exchanger surface, and is the ambient temperature(0 K in space). The flow was thought as laminar flow since the Reynolds number was lower than 2300, and the velocity sliding wall condition was used because of the characteristics of microchannels [28]. The SIMPLE algorithm was adopted for the calculation method, and nitrogen was regarded as an ideal gas by using NIST real gas model. In the initial state, the fluid domain temperature was set to 300 K, while the solid domain temperature was set to 2400 K. The first-order discrete scheme was adopted for the time term, where the time step was 0.001 s and the number of time steps was 50000. Specific boundary conditions are shown in Table 2.
The variation of gas outlet temperature and heat exchanger wall temperature with time were calculated, respectively, and formula (2) was used to evaluate the efficiency of the heat exchanger, and the kinetic energy of gas was ignored in calculation since it was three orders of magnitude smaller than the heat transfer energy, and the results are shown in Figure 4. where is the solar power, is the system thermal efficiency, is the mass flow rate, and is the gas specific enthalpy of inlet and outlet, respectively.
(a) Variation of gas outlet temperature and radiation wall temperature with time, where
(b) Variation of wall temperature with time at different radial positions, where
(c) Variation of heat exchange efficiency with time
Figure 4(a) suggests that the outlet temperature of the gas reached the maximum 2387 K at time , and decreased to 1935 K at time , while the wall temperature decreased from 2400 K to 1963 K, and the temperature difference DeltaT between the wall and propellant developed gradually with time. In Figure 4(b), the variation of wall temperature under different radial length (where, 0.019 m, 0.022 m, and 0.025 m, respectively) is analyzed, where DeltaTw reflected the temperature difference between the inside and outside wall of the heat exchanger and the value of DeltaTw reached the maximum 87 K at time and then gradually decreased and finally stabilized at about 83 K. The maximum heat exchange efficiency of the platlet could reach 87%, and with the weakening of the heat exchange effect, the heat exchange efficiency in the steady state remained about 69% (Figure 4(c)).
4. Model Simplification and Result Comparison
In engineering practice, the focus is on the fluid temperature at the outlet of the heat exchanger rather than the internal temperature field. Therefore, it is necessary to simplify the model to save computational resources. As the energy exchange mode between propellant and working medium is mainly convective heat transfer, it is necessary to solve the convective heat transfer coefficient before simplifying the model. Moreover, since the platelet structure was proposed for the first time, there was no mature empirical formula that could be used before, and it is necessary to fit the heat transfer coefficient according to the structural characteristics.
4.1. 3D Steady-State CFD Simulation of Conjugate Heat Transfer
Considering that there is no time term in the empirical formula, the coupled temperature field was calculated by steady-state method, and the Reynolds number was changed by changing the mass flow rate at the inlet to explore the relationship between Nusselt number and Reynolds number. The boundary conditions of steady-state calculation are shown in Table 3.
The temperature distribution of the gas passing through the heat exchanger at different mass flow rates is shown in Figures 5(a)–5(f). The outlet temperature decreased with the increase of mass flow rate. When the mass flow rate was , the average static temperature at the outlet was 2387 K. Nevertheless, when the mass flow rate was 4.38 × 10-6 kg/s, the average static temperature at the outlet was only 2028 K, about 350 K lower than the maximum outlet temperature.
(a)
(b)
(c)
(d)
(e)
(f)
The fluid domain was calculated by infinitesimal method, and the division of infinitesimal in distributed passage and metering passage is shown in Figures 6(a) and 6(b), respectively. The section was established in a direction perpendicular to the flow direction.
(a) Division of fluid domain in distributed passage, where , , , and
(b) Division of fluid domain in metering passage, where and
4.2. Calculation and Nonlinear Fitting of Nusselt Number
The influence of flow migration term should be considered in the calculation of section average temperature. where is the section average temperature of propellant, is the propellant density, is the specific heat of propellant at constant pressure, is the velocity component normal to the section, and is the section temperature.
The wall temperature took the average temperature of the part in contact with the gas. where is the section average wall temperature, is the wall temperature at boundary, and is the boundary between wall and gas section.
Variation of section average temperature along radial direction between gas and wall in different mass flow rate is shown in Figures 7(a) and 7(b). The figure revealed that the overall temperature of fluid and solid decreased with the increase of mass flow rate due to the influence of migration term.
(a) Variation of gas temperature along radial direction in different mass flow rate
(b) Variation of wall temperature along radial direction in different mass flow rate
The heat flux at the wall was determined by the temperature gradient: where is the heat flux, is the thermal conductivity of gas, and is the gradient operator. The convective heat transfer coefficient at the section could be expressed as:
Thus, the Nusselt number and Reynolds number at the section are calculated by: where is the Nusselt number at the section, is the thermal conductivity of fluid, is the characteristic length of fluid domain, is the Reynolds number at the section, and is the dynamic viscosity of fluid. Variation of section average Nusselt number and Reynolds number along radial direction in different mass flow rate are shown in Figure 8.
(a) Variation of Nusselt number along radial direction in different mass flow rate
(b) Variation of Reynolds number along radial direction in different mass flow rate
As Figure 8 suggests, the Nusselt number in the inlet section is much larger than that in other areas due to the boundary layer effect, and the temperature boundary layer gradually increases with the increase of mass flow rate as was highlighted in Figure 9. Although the Reynolds number increased with the increase of mass flow rate, the change of Nusselt number was not obvious, and the Nusselt number and Reynolds number in different passages varied differently, which meant that the nonlinear fitting needed to be carried out in different zones.
Since the small characteristic length of the platelet heat exchanger, the boundary layer effect should not be ignored. The average Nusselt number of the whole heat exchange surface rather than the local Nusselt number should be considered in the heat transfer calculation. The average Nusselt number and Reynolds number of the whole wall can be obtained from integration method:
In distributed passage:
In metering passage: where is the average Nusselt number and Reynolds number (subscripts 1 and 2 indicate different passages), , , , , and is the infinitesimal area.
In general, the empirical correlation of Nusselt number is usually expressed in the form of power series of Reynolds number and Prandtl number [29]. Since the Prandtl number of gas is approximately constant, the Nusselt number is only related to Reynolds number. The correlation of average Nusselt number was obtained by fitting with nonlinear least square method:
In distributed passage:
In metering passage:
The results of nonlinear fitting are shown in Figure 10, and the sum of squares of residuals in distributed passage and metering passage is 0.000221 and 0.0001514, respectively, which proved that the fitting result was credible.
(a) In distributed passage
(b) In metering passage
4.3. Establishment of 1D Simplified Model
The loose coupling algorithm based on quasi-steady flow field was used to simplify the model. The algorithm assumed that the flow field was in several quasi-steady states, and the energy equations of fluid and solid were solved by tight coupling transient. According to the structural characteristics of the heat exchanger, the following assumptions need to be made before simplification: (i)The flow is a one-dimensional, laminar flow(ii)The flow satisfies the continuity assumption(iii)The thermophysical parameters of the heat exchanger are set as constants(iv)There is no chemical reaction or deposition(v)The nonuniformity of flow properties over the nozzle cross section is ignored(vi)The heat radiation of the heat exchanger exists only on the solid surface(vii)The pressure is 0 Pa and the temperature is 0 K in vacuum environment
At the same time, since the temperature difference on the inner wall of the heat exchanger was less than 100 K (as shown in Figure 7(b)), it could be assumed that the temperature distribution in the heat exchanger was uniform, which meant lumped parameter method could be used. The qualitative temperature of fluid was equal to the average temperature of inlet and outlet, and the thermophysical parameters were determined by the qualitative temperature. According to the continuity assumption, the mass flow rate can be expressed as:
The compressibility of gas is usually measured by Mach number [30]. When the Mach number is less than 0.3, the relative change of gas density caused by the change of gas flow velocity is very small; hence, the gas can be treated as an incompressible fluid. In this paper, the maximum flow rate appeared in the metering passage, and the Mach number could be expressed by the following formulas [31]: where is the Mach number, is the gas velocity, is the mass flow rate of propellant, is the local acoustic velocity, is the sectional area of the metering passage, is the gas constant, is the propellant density, and is the specific heat ratio. By substituting the maximum mass flow rate into formula (13), the maximum Mach number can be estimated as:
Since the maximum Mach number is much less than 0.3, the gas in the heat exchanger can be considered as incompressible gas, and it is considered that the pressure loss is mainly caused by fluid acceleration, reflux, and other factors. The pressure drop is: where is the sectional area of the component and is the local loss coefficient.
Therefore, the average Nusselt number and average Reynolds number in different passages can be calculated as:
In distributed passage:
In metering passage: where is the cross-sectional area normal to the flow direction, is the characteristic length of fluid domain, is the dynamic viscosity at qualitative temperature, is the mass flow rate, , , , and . The energy equation in the heat exchanger is as follows:
For gas:
In distributed passage:
In metering passage:
For solid: where is the average convective heat transfer coefficient; , , and are the section average temperature at inlet, channel junction, and outlet, respectively; , , and are the specific heat at constant pressure of heat exchanger, gas in distributed passage, and gas in metering passage, respectively; is the radiation power; is the radiation intensity of light source; is the concentrated ratio; is the open area of concentrator; is the concentrator efficiency; is the Stefan-Boltzmann constant; is the radiant emissivity of heat exchanger surface; is the surface area of heat exchanger; is the qualitative temperature; is the heat exchanger mass; is the heat exchanger temperature; and is the ambient temperature.
4.4. Comparison with 3D Results
As shown in Figure 11, the results of 1D simulation are consistent with that of 3D simulation. Overall, the gas outlet temperature of 1D simulation is slightly lower than that of 3D simulation, while the solid wall temperature is just the opposite, which might be due to the deviation between qualitative temperature and actual temperature. The maximum estimation error was about 5% for both gas outlet temperature and gas outlet temperature.
(a) Comparison of 3D simulation and 1D simulation results when the mass flow rate was
(b) Comparison of 3D simulation and 1D simulation results when the mass flow rate was
(c) Comparison of 3D simulation and 1D simulation results when the mass flow rate was
As the gas outlet temperature is difficult to measure in experiment, the thrust is usually used to replace the temperature. Taking Laval nozzle as an example, its mathematical model can be expressed as: where is the nozzle outlet cross-sectional area, is the throat cross-sectional area, is the adiabatic index of nitrogen, is the chamber pressure, and is the nozzle outlet pressure. where is the gas constant, is the temperature in chamber, is the exhaust velocity of nozzle, is the specific impulse, is the vacuum thrust, is the estimates of velocity loss and pressure loss in the nozzle part [32, 33], and is the exit pressure, which is 0 MPa in space.
The 1D model of solar thermal thruster was established by combining formulas (12)–(22) with MATLAB Simulink. The framework of the model is shown in Figure 12.
4.5. Experiment of STP System
The experimental system consists of propellant supply system, thruster test platform, vacuum chamber, xenon lamp light source, and thruster. The propellant supply system is mainly composed of stop valve, pressure reducing valve, pressure gauge, gas cylinder, and pipeline, and the thruster test platform is mainly composed of data acquisition card, various sensors, and corresponding acquisition software. The configuration of each component of the system is shown in Figure 13, and the schematic diagram of STP system is shown in Figure 14.
(a) Propellant supply system
(b) Vacuum chamber
(c) Xenon lamp light source
(d) Thruster
(e) Platelet heat exchanger
(f) Thrust measuring device
The parameters that can be directly measured by the experimental system include propellant volume flow rate, wall temperature of thrust chamber, and thrust and pressure in the thrust chamber.
4.6. Comparison with Experimental Results
In the experiment, the inlet pressure of propellant was 0.8 MPa, the ambient pressure was about 50 Pa, the inlet volume flow rate of propellant was maintained at 14.5 nL/min, the power of xenon light source was 2000 W, the solenoid valve of propellant tank would be opened when the platelet temperature reached 2400 K, and as the thruster was wrapped with thermal insulation material, the emissivity ε can be taken as 0.65. The experiment results are shown in Figure 15.
The results suggest that the propellant flow remained basically stable, while the thrust increased rapidly from zero, then decreased gradually, which might be caused by the weakening of heat transfer effect and the decrease of environmental vacuum. The fluctuation of thrust and volume flow rate in the experiment was caused by the measurement accuracy of the instrument.
To verify the accuracy of 1D model, the parameter settings in the simulation are shown in Table 4, and the comparison between the experimental results and the simulation results of 1D simplified model is shown in Figure 16.
(a) Comparison between experimental and simulation values of
(b) Comparison between experimental and simulated values of
It can be concluded from Figure 16 that each variable had a maximum at the beginning, the maximum thrust and specific impulse in the simulation results were 0.56 N and 2100 m/s, and the variation trend of experimental value was similar to that of simulation value in the initial stage of heat exchange, but with the increase of time, the simulation model reached steady state, while the experimental value continued to decrease.
5. Conclusion and Discussion
In this paper, a new platelet heat exchange structure of solar thermal thruster was proposed. The structure could not only increase the heat exchange area, but also divert the flow of propellant. Compared with spiral tube heat exchanger, the platelet could heat the propellant to a higher temperature and improve the heat exchange efficiency. The CHT characteristics between heat exchanger and propellant were analyzed. In 3D unsteady-state simulation, the variation of gas temperature and wall temperature with time under sliding wall and radiation boundary conditions were calculated. In the 3D steady-state simulation, the infinitesimal method was used to divide the flow field, and the empirical correlation of Nusselt number was fitted by nonlinear regression method. The convective heat transfer coefficient was solved in the light of the obtained empirical correlation. A loose coupling algorithm based on quasi-steady flow field was used to establish a 1D simplified model, which was compared with the 3D simulation results and experimental results to confirm the accuracy of the 1D model. However, the results of the simplified model were consistent with the experimental results in the initial stage of heat exchange. Then, the variables in the simulation gradually reached the steady state, while they continued to decrease in the experiment. Though there is a certain difference between the theoretical value and the experimental value, the overall error is within the acceptable range. The reasons for the continuous decrease of variables in the experiment might be as follows: (i)The accuracy limitation of thrust measurement system(ii)Lack of tightness of vacuum chamber(iii)The nitrogen was not pumped away in time after discharged from the thruster, resulting in the continuous decline of vacuum in the vacuum chamber
Furthermore, the following conclusions and innovations can be drawn: (1)In the calculation unit, the propellant was heated to more than 2380 K when the wall temperature reached 2400 K through the platelet heat exchanger, while the inlet pressure was 0.8 MPa and the mass flow rate was , which was 200 K higher than that of the traditional spiral tube heat exchanger, and the heat exchange efficiency could reach 87%. Owing to the fluid-solid coupling effect, the heat exchanger temperature and fluid temperature finally decreased to 1963 K and 1935 K, about 18.7% lower than the maximum temperature, and the heat exchange efficiency decreased to 69%.(2)The propellant inlet mass flow rate had a significant effect on the outlet temperature of the heat exchanger. The average static temperature at the outlet reached 2387 K when the mass flow rate was set to . Consequently, the mass flow rate should be reduced as much as possible to obtain a higher temperature so that the propellant could be fully heated. According to the structural characteristics of the platelet heat exchanger and the flow characteristics of the propellant, the empirical correlations were obtained by nonlinear fitting at the distributed passage and the metering passage, and the residual errors of the fitting results were 0.000221 and 0.0001514, respectively, which confirmed the credibility of the simplified method(3)The error between 1D simulation results and 3D simulation results was less than 5%, but when the simulation time was set to 50s, the time spent in 3D CFD simulation exceeded 10 hours, which increased rapidly with the increase of grid number. However, 1D model only took a few seconds and could obtain the variation of key parameters. Therefore, 1D model had obvious advantages in unsteady-state calculation. The 1D model was consistent with the experimental results at the initial stage of heat exchange, but the variables in the experiment gradually decreased with time and finally failed in reaching the steady state, which might be caused by the insufficient sealing of the experiment and the limitation of the accuracy of the instrument
The above conclusions provide a new design guide for STP system of platelet heat exchange structure. The future work will focus on improving the experimental environment, adopting more accurate measuring instruments and strengthening the sealing of the vacuum chamber. Other more competitive working mediums will be used as propellants, such as hydrogen and ammonia.
Data Availability
The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy or ethical restrictions.
Conflicts of Interest
We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work and there is no professional or other personal interest of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.