# Study on Heat Exchange Structure Design and Propulsion Performance of Solar Thermal Thruster

Authors: Haoran Zhang; Minchao Huang; Xiaoping Hu
Journal: International Journal of Aerospace Engineering
 (2022)

---

## 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.


---

## Body

## 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 Figure1. 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.Figure 1
The structure diagram of solar thermal propulsion.Figure 2
The configuration diagram of platelet heat exchanger.
(a)
Front view(b)
Sectional view(c)
Upper viewUnder 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.Table 1
Thermal parameter of the platelet material.
MaterialDensityg/cm3Thermal conductivityW/m⋅KSpecific heatJ/kg⋅KMelting point °CMolybdenum10.093.73282622Alumina fiber0.20.012002500


## 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 Figure3, 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 3.8×10−8m. According to Knudsen number calculation formula, Kn=7.6×10−5<0.001, which meant the continuum model was applicable.Figure 3
The unit and meshing.
(a)
Unit(b)
MeshingActually, 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]:
(1)−kAr∂Tw∂rr=r0=IηCAc−σεArTw4−T4,where Tw is the wall temperature, k is the thermal conductivity of platelet, Ar is the radiation area of platelet, I is the radiation intensity of solar, C is the concentrated ratio, Ac is the open area of concentrator, η is the concentrator efficiency, σ is the Stefan-Boltzmann constant, ε is the radiant emissivity of heat exchanger surface, and T 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.Table 2
Boundary conditions of unsteady-state simulation.
ZoneBoundary typeValueInletMass flow rateTemperature1.14×10−6kg/s300 KPressure0.8 (MPa)OutletPressure0Radiation wallHeat flux1.2228×10−6−5.67×10−8×0.65×Tw4W/m2Outer wallRadiation5.67×10−8×0.65×Tw4InterfaceCoupled\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.
(2)P=IηCAc,ς=ṁheṁhi+P,where P is the solar power, ς is the system thermal efficiency, ṁ is the mass flow rate, and hi,he is the gas specific enthalpy of inlet and outlet, respectively.Figure 4
3D unsteady-state CFD simulation results.
(a)
Variation of gas outlet temperature and radiation wall temperature with time, whereDeltaT=Tw−Tout(b)
Variation of wall temperature with time at different radial positions, whereDeltaTw=Twr=0.016−Twr=0.025(c)
Variation of heat exchange efficiency with timeFigure4(a) suggests that the outlet temperature of the gas reached the maximum 2387 K at time t=0.1s, and decreased to 1935 K at time t=48s, 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 (wherer=0.016m, 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 t=3.2s 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 Table3.Table 3
Boundary conditions of 3D steady-state CFD simulation.
ZoneBoundary typeValueInletMass flow rate1.14~4.38×10−6kg/sPressure0.8 (MPa)OutletPressure0Radiation wallTemperature2400 KConvection wallConvective heat transfer coefficient10 (W/m2K)InterfaceCoupled\The temperature distribution of the gas passing through the heat exchanger at different mass flow rates is shown in Figures5(a)–5(f). The outlet temperature decreased with the increase of mass flow rate. When the mass flow rate was 1.14×10−6kg/s, 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.Figure 5
Temperature distribution of propellant at different mass flow rate.
(a)
ṁ=1.14×10−6kg/s(b)
ṁ=1.725×10−6kg/s(c)
ṁ=2.31×10−6kg/s(d)
ṁ=3.00×10−6kg/s(e)
ṁ=3.69×10−6kg/s(f)
ṁ=4.38×10−6kg/sThe fluid domain was calculated by infinitesimal method, and the division of infinitesimal in distributed passage and metering passage is shown in Figures6(a) and 6(b), respectively. The section was established in a direction perpendicular to the flow direction.Figure 6
Division of fluid domain in different passage.
(a)
Division of fluid domain in distributed passage, whereH=0.001m, dr=0.0002m, θ=22.5°, and l=θ×r(b)
Division of fluid domain in metering passage, wherer̅=0.0025m and dr=0.0002m


### 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.(3)Tf=∫ρCpuTdA∫ρCpudA,where Tf is the section average temperature of propellant, ρ is the propellant density, Cp is the specific heat of propellant at constant pressure, u is the velocity component normal to the section, and T is the section temperature.The wall temperature took the average temperature of the part in contact with the gas.(4)Tw=∫Tsdl∫dl,where Tw is the section average wall temperature, Ts is the wall temperature at boundary, and l 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 Figures7(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.Figure 7
Variation of section average temperature along radial direction between gas and wall in different mass flow rate.
(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 rateThe heat flux at the wall was determined by the temperature gradient:(5)qr=−k∇T,where q is the heat flux, k is the thermal conductivity of gas, and ∇ is the gradient operator. The convective heat transfer coefficient hr at the section could be expressed as:
(6)hr=qrTf−Tw.Thus, the Nusselt number and Reynolds number at the section are calculated by:(7)Nur=hrLk,Rer=ρuLμ,where Nur is the Nusselt number at the section, k is the thermal conductivity of fluid,L is the characteristic length of fluid domain, Rer 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.Figure 8
Variation of section average Nu and Re along radial direction in different mass flow rate.
(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 rateAs Figure8 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.Figure 9
Contour of propellant temperature along radial direction in different mass flow rate.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:(8)A1=θr22−r12,dA1=θr+dr2−r2,Nu¯1=∫NurdA1A1=∑r=r1r=r22θrNurdrA1,Re¯1=∫Rerdrr2−r1=∑r=r1r=r2Rerdrr2−r1.In metering passage:(9)Nu¯2=∫NurdA2A2=∑r=r1r=r2Nurdrr1−r0,Re¯2=∫RerdAA2=∑r=r1r=r2Rerdrr1−r0,where Nu¯,Re¯ is the average Nusselt number and Reynolds number (subscripts 1 and 2 indicate different passages), r0=0.016m, r1=0.019m, r2=0.025m, θ=π/8, and A 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:(10)Nu¯1=0.9381Re¯1+41.950.493.In metering passage:(11)Nu¯2=3.1×10−8Re¯22.56+5.04.The results of nonlinear fitting are shown in Figure10, 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.Figure 10
Correlation of average Nusselt number in different passages.
(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 environmentAt the same time, since the temperature difference on the inner wall of the heat exchanger was less than 100 K (as shown in Figure7(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:
(12)m·=ρuA.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]:
(13)Ma=ua,u=ṁρA,a=kR0T,ρ=PRT,where Ma is the Mach number, u is the gas velocity, ṁ is the mass flow rate of propellant, a is the local acoustic velocity, A is the sectional area of the metering passage, R0 is the gas constant,ρ is the propellant density, and k is the specific heat ratio. By substituting the maximum mass flow rate 4.38×10−6kg/s into formula (13), the maximum Mach number can be estimated as:
(14)Mamax=ṁPAR0Tk=4.38×10−68×105×π×0.00252/48.314×23000.028×1.4=7.79×10‐4<<0.3.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:(15)Δp=ξṁ22ρA2,where A 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:(16)S1r=θrL12,Re1r=m·L1μS1r=2m·μθr,Re¯1=∫Re1rdrr2−r1=2m·μ1θr2−r1lnr2r1,Nu¯1=0.9381Re¯1+41.950.493.In metering passage:(17)S2=πr¯24,Re2=m·L2μ2S2=8m·μ2πr¯,Re¯2=∫Re2drr1−r0=Re2=8m·μ2πr¯,Nu¯2=3.1×10−8Re¯22.56+5.04,where S is the cross-sectional area normal to the flow direction, L is the characteristic length of fluid domain, μ is the dynamic viscosity at qualitative temperature,ṁ is the mass flow rate, r0=0.016m, r1=0.019m, r2=0.025m, and r¯=0.00025m. The energy equation in the heat exchanger is as follows:For gas:In distributed passage:(18)h¯1=Nu¯1k1L1,A1=θr22−r12,Tf1=Tin+Tmid2,ṁCp1Tmid−Tin=h¯1A1Tw−Tf1.In metering passage:(19)h¯2=Nu¯2k2L2,A2=πr¯r1−r02,Tf2=Tmid+Tout2,ṁCp2Tout−Tmid=h¯2A2Tw−Tf2.For solid:(20)Ew=ICAcη−σεAwTw4−T04,MCpwdTwdt=Ew−h¯1A1Tw−Tf1−h¯2A2Tw−Tf2,where h¯ is the average convective heat transfer coefficient; Tin, Tmid, and Tout are the section average temperature at inlet, channel junction, and outlet, respectively; Cpw, Cp1, and Cp2 are the specific heat at constant pressure of heat exchanger, gas in distributed passage, and gas in metering passage, respectively; Ew is the radiation power; I is the radiation intensity of light source; C is the concentrated ratio; Ac is the open area of concentrator; η is the concentrator efficiency; σ is the Stefan-Boltzmann constant; ε is the radiant emissivity of heat exchanger surface; Aw is the surface area of heat exchanger; Tf is the qualitative temperature; M is the heat exchanger mass; Tw is the heat exchanger temperature; and T0 is the ambient temperature.


### 4.4. Comparison with 3D Results

As shown in Figure11, 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.Figure 11
Comparison of 3D simulation and 1D simulation results under different mass flow rates.
(a)
Comparison of 3D simulation and 1D simulation results when the mass flow rate was1.14×10−6kg/s(b)
Comparison of 3D simulation and 1D simulation results when the mass flow rate was3.00×10−6kg/s(c)
Comparison of 3D simulation and 1D simulation results when the mass flow rate was4.38×10−6kg/sAs 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:(21)AeAt=ΓPe/Pc1/k2k/k−11−Pe/Pck−1/k,Γ=k2k+1k+1/2k−1,where Ae is the nozzle outlet cross-sectional area, At is the throat cross-sectional area, k is the adiabatic index of nitrogen, Pc is the chamber pressure, and Pe is the nozzle outlet pressure.
(22)ue=R⋅Tc2kk−11−pepck−1/k1/2,(23)Isp=∫0tFtdt∫0tm·dt=ue+AepeRTcAtpcΓ,(24)F=ṁ⋅ue⋅η1+Aepe⋅η2−p0,where R=R0/Mnitrogen is the gas constant, Tc is the temperature in chamber, ue is the exhaust velocity of nozzle, Isp is the specific impulse, F is the vacuum thrust,η1,η2 is the estimates of velocity loss and pressure loss in the nozzle part [32, 33], and P0 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.Figure 12
The framework of the 1D model in Simulink.


### 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 Figure13, and the schematic diagram of STP system is shown in Figure 14.Figure 13
The configuration of each component of STP system.
(a)
Propellant supply system(b)
Vacuum chamber(c)
Xenon lamp light source(d)
Thruster(e)
Platelet heat exchanger(f)
Thrust measuring deviceFigure 14
The schematic diagram of STP system.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.Figure 15
Variation of mass flow rate and thrust with time in the experiment.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 Table4, and the comparison between the experimental results and the simulation results of 1D simplified model is shown in Figure 16.Table 4
Parameters in simulation.
ParametersValueInlet pressure8×105paInlet temperature293 KConcentrator diameter0.1 mLight source power2000 WNozzle throat diameter0.0013 mNozzle outlet diameter0.0136 mExpansion angle of nozzle expansion section15°Expansion ratio of nozzle expansion section109.4430Mass flow rate2.8017×10−4kg/sConcentrator efficiency85%Nozzle efficiency95%Figure 16
Comparison between the experimental results and the simulation results of 1D simplified model.
(a)
Comparison between experimental and simulation values ofF(b)
Comparison between experimental and simulated values ofIspIt can be concluded from Figure16 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.


## 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 Table3.Table 3
Boundary conditions of 3D steady-state CFD simulation.
ZoneBoundary typeValueInletMass flow rate1.14~4.38×10−6kg/sPressure0.8 (MPa)OutletPressure0Radiation wallTemperature2400 KConvection wallConvective heat transfer coefficient10 (W/m2K)InterfaceCoupled\The temperature distribution of the gas passing through the heat exchanger at different mass flow rates is shown in Figures5(a)–5(f). The outlet temperature decreased with the increase of mass flow rate. When the mass flow rate was 1.14×10−6kg/s, 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.Figure 5
Temperature distribution of propellant at different mass flow rate.
(a)
ṁ=1.14×10−6kg/s(b)
ṁ=1.725×10−6kg/s(c)
ṁ=2.31×10−6kg/s(d)
ṁ=3.00×10−6kg/s(e)
ṁ=3.69×10−6kg/s(f)
ṁ=4.38×10−6kg/sThe fluid domain was calculated by infinitesimal method, and the division of infinitesimal in distributed passage and metering passage is shown in Figures6(a) and 6(b), respectively. The section was established in a direction perpendicular to the flow direction.Figure 6
Division of fluid domain in different passage.
(a)
Division of fluid domain in distributed passage, whereH=0.001m, dr=0.0002m, θ=22.5°, and l=θ×r(b)
Division of fluid domain in metering passage, wherer̅=0.0025m and dr=0.0002m


## 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.(3)Tf=∫ρCpuTdA∫ρCpudA,where Tf is the section average temperature of propellant, ρ is the propellant density, Cp is the specific heat of propellant at constant pressure, u is the velocity component normal to the section, and T is the section temperature.The wall temperature took the average temperature of the part in contact with the gas.(4)Tw=∫Tsdl∫dl,where Tw is the section average wall temperature, Ts is the wall temperature at boundary, and l 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 Figures7(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.Figure 7
Variation of section average temperature along radial direction between gas and wall in different mass flow rate.
(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 rateThe heat flux at the wall was determined by the temperature gradient:(5)qr=−k∇T,where q is the heat flux, k is the thermal conductivity of gas, and ∇ is the gradient operator. The convective heat transfer coefficient hr at the section could be expressed as:
(6)hr=qrTf−Tw.Thus, the Nusselt number and Reynolds number at the section are calculated by:(7)Nur=hrLk,Rer=ρuLμ,where Nur is the Nusselt number at the section, k is the thermal conductivity of fluid,L is the characteristic length of fluid domain, Rer 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.Figure 8
Variation of section average Nu and Re along radial direction in different mass flow rate.
(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 rateAs Figure8 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.Figure 9
Contour of propellant temperature along radial direction in different mass flow rate.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:(8)A1=θr22−r12,dA1=θr+dr2−r2,Nu¯1=∫NurdA1A1=∑r=r1r=r22θrNurdrA1,Re¯1=∫Rerdrr2−r1=∑r=r1r=r2Rerdrr2−r1.In metering passage:(9)Nu¯2=∫NurdA2A2=∑r=r1r=r2Nurdrr1−r0,Re¯2=∫RerdAA2=∑r=r1r=r2Rerdrr1−r0,where Nu¯,Re¯ is the average Nusselt number and Reynolds number (subscripts 1 and 2 indicate different passages), r0=0.016m, r1=0.019m, r2=0.025m, θ=π/8, and A 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:(10)Nu¯1=0.9381Re¯1+41.950.493.In metering passage:(11)Nu¯2=3.1×10−8Re¯22.56+5.04.The results of nonlinear fitting are shown in Figure10, 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.Figure 10
Correlation of average Nusselt number in different passages.
(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 environmentAt the same time, since the temperature difference on the inner wall of the heat exchanger was less than 100 K (as shown in Figure7(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:
(12)m·=ρuA.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]:
(13)Ma=ua,u=ṁρA,a=kR0T,ρ=PRT,where Ma is the Mach number, u is the gas velocity, ṁ is the mass flow rate of propellant, a is the local acoustic velocity, A is the sectional area of the metering passage, R0 is the gas constant,ρ is the propellant density, and k is the specific heat ratio. By substituting the maximum mass flow rate 4.38×10−6kg/s into formula (13), the maximum Mach number can be estimated as:
(14)Mamax=ṁPAR0Tk=4.38×10−68×105×π×0.00252/48.314×23000.028×1.4=7.79×10‐4<<0.3.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:(15)Δp=ξṁ22ρA2,where A 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:(16)S1r=θrL12,Re1r=m·L1μS1r=2m·μθr,Re¯1=∫Re1rdrr2−r1=2m·μ1θr2−r1lnr2r1,Nu¯1=0.9381Re¯1+41.950.493.In metering passage:(17)S2=πr¯24,Re2=m·L2μ2S2=8m·μ2πr¯,Re¯2=∫Re2drr1−r0=Re2=8m·μ2πr¯,Nu¯2=3.1×10−8Re¯22.56+5.04,where S is the cross-sectional area normal to the flow direction, L is the characteristic length of fluid domain, μ is the dynamic viscosity at qualitative temperature,ṁ is the mass flow rate, r0=0.016m, r1=0.019m, r2=0.025m, and r¯=0.00025m. The energy equation in the heat exchanger is as follows:For gas:In distributed passage:(18)h¯1=Nu¯1k1L1,A1=θr22−r12,Tf1=Tin+Tmid2,ṁCp1Tmid−Tin=h¯1A1Tw−Tf1.In metering passage:(19)h¯2=Nu¯2k2L2,A2=πr¯r1−r02,Tf2=Tmid+Tout2,ṁCp2Tout−Tmid=h¯2A2Tw−Tf2.For solid:(20)Ew=ICAcη−σεAwTw4−T04,MCpwdTwdt=Ew−h¯1A1Tw−Tf1−h¯2A2Tw−Tf2,where h¯ is the average convective heat transfer coefficient; Tin, Tmid, and Tout are the section average temperature at inlet, channel junction, and outlet, respectively; Cpw, Cp1, and Cp2 are the specific heat at constant pressure of heat exchanger, gas in distributed passage, and gas in metering passage, respectively; Ew is the radiation power; I is the radiation intensity of light source; C is the concentrated ratio; Ac is the open area of concentrator; η is the concentrator efficiency; σ is the Stefan-Boltzmann constant; ε is the radiant emissivity of heat exchanger surface; Aw is the surface area of heat exchanger; Tf is the qualitative temperature; M is the heat exchanger mass; Tw is the heat exchanger temperature; and T0 is the ambient temperature.


## 4.4. Comparison with 3D Results

As shown in Figure11, 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.Figure 11
Comparison of 3D simulation and 1D simulation results under different mass flow rates.
(a)
Comparison of 3D simulation and 1D simulation results when the mass flow rate was1.14×10−6kg/s(b)
Comparison of 3D simulation and 1D simulation results when the mass flow rate was3.00×10−6kg/s(c)
Comparison of 3D simulation and 1D simulation results when the mass flow rate was4.38×10−6kg/sAs 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:(21)AeAt=ΓPe/Pc1/k2k/k−11−Pe/Pck−1/k,Γ=k2k+1k+1/2k−1,where Ae is the nozzle outlet cross-sectional area, At is the throat cross-sectional area, k is the adiabatic index of nitrogen, Pc is the chamber pressure, and Pe is the nozzle outlet pressure.
(22)ue=R⋅Tc2kk−11−pepck−1/k1/2,(23)Isp=∫0tFtdt∫0tm·dt=ue+AepeRTcAtpcΓ,(24)F=ṁ⋅ue⋅η1+Aepe⋅η2−p0,where R=R0/Mnitrogen is the gas constant, Tc is the temperature in chamber, ue is the exhaust velocity of nozzle, Isp is the specific impulse, F is the vacuum thrust,η1,η2 is the estimates of velocity loss and pressure loss in the nozzle part [32, 33], and P0 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.Figure 12
The framework of the 1D model in Simulink.


## 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 Figure13, and the schematic diagram of STP system is shown in Figure 14.Figure 13
The configuration of each component of STP system.
(a)
Propellant supply system(b)
Vacuum chamber(c)
Xenon lamp light source(d)
Thruster(e)
Platelet heat exchanger(f)
Thrust measuring deviceFigure 14
The schematic diagram of STP system.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.Figure 15
Variation of mass flow rate and thrust with time in the experiment.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 Table4, and the comparison between the experimental results and the simulation results of 1D simplified model is shown in Figure 16.Table 4
Parameters in simulation.
ParametersValueInlet pressure8×105paInlet temperature293 KConcentrator diameter0.1 mLight source power2000 WNozzle throat diameter0.0013 mNozzle outlet diameter0.0136 mExpansion angle of nozzle expansion section15°Expansion ratio of nozzle expansion section109.4430Mass flow rate2.8017×10−4kg/sConcentrator efficiency85%Nozzle efficiency95%Figure 16
Comparison between the experimental results and the simulation results of 1D simplified model.
(a)
Comparison between experimental and simulation values ofF(b)
Comparison between experimental and simulated values ofIspIt can be concluded from Figure16 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 chamberFurthermore, 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 was1.14×10−6kg/s, 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 to1.14×10−6kg/s. 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 instrumentThe 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.



---

*Source: 1348289-2022-07-06.xml*