Numerical Study on the Cooling Method of Phase Change Heat Exchange Unit with Layered Porous Media

Abstract

The implementation of heat sinks in high-power pulse electronic devices within hypersonic aircraft cabins has been facilitated by the emergence of innovative phase change materials (PCMs) characterized by excellent thermal conductivity and high latent heat. In this study, a representative material, layered porous media filled with paraffin wax, was utilized, and a three-dimensional numerical model based on the enthalpy-porosity approach was employed. A thermal response research was conducted on the Phase Change Heat Exchange Unit with Layered Porous Media (PCHEU-LPM) with different cooling methods. The results indicate that water cooling proved to be suitable for the PCHEU-LPM with a heat flux of 50,000 W/m². Additionally, parametric studies were performed to determine the optimal cooling conditions, considering the inlet temperature and velocity of the cooling flow. The results revealed that the most suitable conditions were strongly influenced by the coolant inlet parameters, along with the position of the PCM interface. Finally, the identification of the parameter combination that minimizes temperature fluctuations was achieved through the Response Surface Analysis method (RSA). Subsequent verification through simulation further reinforced the reliability of the proposed optimal parameters.

1. Introduction

Hypersonic aircrafts are currently regarded as a focal point in the field of aviation development. With speeds exceeding five times the speed of sound, the ability to conduct global flights in remarkably brief durations is demonstrated by these aircraft, emphasizing their efficiency and widespread applicability [1]. In the multifaceted realm of hypersonic aircraft research, thermal management consistently emerges as a paramount consideration. As the aircraft’s velocity increases, substantial aerodynamic heating emanates from the leading edge, resulting in local peak temperatures that may surpass 1800 K [2,3]. The internal structure of the aircraft is infiltrated by significant aerodynamic heating through thermal protection systems, posing challenges in dissipating heat from onboard equipment. Notably, the heat flux density from high-energy devices, such as radar and processors, is suggested by literature [4] to range from 100 to 500 (W/cm²). Due to the fact that hypersonic aircraft typically rely solely on fuel as the primary cooling source, their cooling capacity is limited while simultaneously requiring cooling for various components such as the outer shell structure, propulsion system, equipment compartment, and passenger cabin. This places extremely high demands on the thermal management system. Taking the EU-designed MR2 hypersonic aircraft as an example [5], as depicted in Figure 1, low-temperature fuel flows out from the tank system. A portion of it, within the blue section of the diagram, exchanges heat and evaporates within the low-temperature tank, which penetrates the aircraft shell to address the thermal load. Another portion exchanges heat through a heat exchanger with gases flowing in from the cabin and equipment compartment. After both sections of heat absorption, the mixed fuel, now at temperatures exceeding 1300 K, enters the propulsion unit for combustion. Throughout this thermodynamic cycle, for the equipment compartment and passenger cabin, typically low-temperature levels and narrow temperature adjustment ranges are applicable. Hence, the thermal management system necessitates heat exchange units responsible for fuel and cabin air heat exchange to possess stable temperature control ranges, recyclability, and high heat exchange rates. Upon scrutinizing the characteristics of thermal sources and mission utilization, it becomes evident that the internal heat sources of hypersonic aircraft exhibit extensive spatial distribution, non-uniformity, and significant temporal variations. Thus, effective thermal control, specifically the mitigation of temperature fluctuations and the maintenance of stable temperature levels over extended periods, are recognized as an urgent issue in the study of transient high-heat-flow-density electronic devices.

The conventional thermal management systems for electronic devices in cabin employ single-phase or two-phase liquid circuits [6,7], with the heat sink playing a crucial role in temperature regulation by absorbing heat from the electronic components to prevent temperature elevation and fluctuations. However, pulsed devices such as radar and processor chips discussed in this paper, present unique characteristics including short peak power durations, significant heat loss during transients, and less heat dissipation requirements during the remaining operating periods. Consequently, the utilization of traditional metal or ceramic materials as heat sink substrates would lead to excessive mass and volume of the overall structure to accommodate the peak heat dissipation during transients, resulting in design redundancies. In this context, it is contended that distinct advantages are offered by PCMs, such as high thermal capacity, isothermal or nearly isothermal behavior during phase changes, excellent stability, and reusability. By leveraging these properties, the absorption of thermal shock from the heat source and the suppression of temperature fluctuations in pulsed equipment operating in intermittent mode can be effectively achieved by using PCMs.

In the application research focusing on utilizing PCMs to mitigate temperature fluctuations induced by pulsed heat flow, Waster et al. [8] proposed a heat storage unit incorporating different types of hydrate salt PCMs, along with a fast calculation model validated through experiments. According to their findings, the key factors influencing the temperature control effectiveness of PCMs are thermal conductivity and latent heat. The use of high thermal conductivity hydrate salt materials resulted in a more pronounced suppression of temperature fluctuations. Maxa et al. [9] designed a PCM film to alleviate temperature rise when chip heating exceeds the design threshold. Results indicated a reduction of the temperature peak from 170 °C to 100 °C, significantly expanding the chip’s maximum operating power and the duty cycle under high load conditions.

In theoretical research, Krishnan et al. [10] investigated the feasibility of PCMs in addressing pulse heat dissipation issues in electronic component thermal management through numerical methods. The study employed discrete pulsed heating on one side of the PCM phase change system and air-cooled heat dissipation on the other side. Response analyses were conducted at different pulse frequencies, determining the optimal heat source arrangement and studying the heat transfer efficiency for different aspect ratios at the same volume. Shamberger et al. [11] numerically studied the dynamic response characteristics of PCMs with varying external frequencies, revealing its effective absorption and release of heat and buffering of transient heat pulses at different frequencies. The study found that the internal temperature distribution in the single-phase region was strongly perturbed by oscillating heat boundary conditions. However, in a specific range of the two-phase region, temperature changes were inhibited, resulting in a phase lag (∆φ) and a reduction in the peak temperature change at the heat source (∆T). The magnitude and frequency dependence of this anti-resonance depended on the characteristics of the periodic heating function, material thermophysical properties, slab thickness, and the nature of the applied cooling boundary conditions.

In the aforementioned studies, the application of PCMs were limited by their low thermal conductivity, typically below 1 (W/m·K). The duration of each cycle for pulsed temperature or heat flow often exceeded 1 h, and the thermal load that electronic devices could withstand was relatively small [12,13]. However, recent technological developments, such as low-temperature alloys, hydrate salts, foam metals, microcapsules, etc., have greatly expanded the thermal conductivity of PCMs, even reaching 150 (W/m·K) [14]. The increased thermal conductivity enables PCMs to handle stronger transient features and higher heat consumption, and with the reduced internal thermal resistance of phase change materials, the external heat dissipation intensity can also be correspondingly increased. However, there is currently a lack of research on the applicable scenarios of high thermal conductivity phase change materials, especially the appropriate heat dissipation methods.

In this study, a layered porous medium phase change material was utilized as the foundation for the heat sink. Favorable morphological traits are exhibited by this material, integrating a substantial latent heat capacity with an impressive thermal conductivity of up to 12 (W/m·K). Initially, a three-dimensional numerical model was constructed that incorporates composite boundary conditions, encompassing pulsed heating and forced convective heat dissipation. These conditions were implemented at both the upper and lower extremities of the PCHEU-LPM. Next, a comparison for the effects of different thermal boundary conditions on temperature and phase fields was conducted, and the most suitable heat dissipation method was to be selected. Then, through parameter research, the relationship between the changes in parameters of heat dissipation conditions and temperature fluctuations and uniformity was studied. Finally, through RSA, the optimal parameter values for heat dissipation were found, and a reference for the design of future PCHEU-LPM is provided.

2. Methodology

2.1. Physical Model

A rectangular PCM plate studied in this paper is shown in Figure 2a, and the x, y, z directions are a length of $L=50 \mathrm{m}\mathrm{m}$, a width of $=20 \mathrm{m}\mathrm{m}$, and a height of $H=10 \mathrm{m}\mathrm{m}$. The cross-sectional heat transfer process along with the heat flow direction is shown in Figure 2b, with a uniformly distributed pulsed heat flow ${\tilde{q}}_w$ loaded at the top, and at the bottom is the incoming flow temperature $T_c$ and the incoming flow velocity $U_c$. The PCM heat transfer structure absorbs the transient pulsed heat flow on one side, while forcing convection heat dissipation on the other side. The porous phase change material was studied as a whole, and its physical properties were calculated by combining the metal skeleton and paraffin. The physical properties of the metal skeleton and paraffin are shown in Table 1.

2.2. Control Equations and Thermal Boundary Conditions

The following applies to the object shown in Figure 2, along with the following reasonable assumptions:

(1) The isotropy of PCM material.

(2) The neglect of the heat capacity and thermal resistance of the package housing itself on the heat transfer process.

(3) The disregard of the PCM heat transfer structure dissipating heat to the environment through the surrounding casing.

(4) The thermal conductivity, heat capacity, latent heat, and other thermal physical parameters of PCM materials are solely related to the phase region and are independent of temperature.

(5) The change in density of the molten state liquid PCM is unrelated to viscosity and is only accounted for in the body-force term of the momentum equation, obeying the Boussinesq assumption.

On this basis, the enthalpy-porosity model is used to mathematically describe the PCM before and after the phase change and the melting and solidification process [15,16]. The model uses the porosity $\epsilon$ to represent the phase change state, which is solid at $\epsilon =0$, liquid at $\epsilon =1$, and molten at $0<\epsilon <1$, where the phase change material is considered as a porous media material.

Its N-S control equation is as follows:

(1)$\frac{\partial \rho }{\partial \tau }+\frac{\partial (\rho u)}{\partial X}+\frac{\partial (\rho v)}{\partial Y}+\frac{\partial (\rho w)}{\partial Z}=0$

(2)$\rho \left[\frac{\partial u}{\partial \tau }+\frac{\partial (uu)}{\partial X}+\frac{\partial (uv)}{\partial Y}+\frac{\partial (uw)}{\partial Z}\right]=\left[\frac{\partial }{\partial X}\left(\mu \frac{\partial u}{\partial X}\right)+\frac{\partial }{\partial Y}\left(\mu \frac{\partial u}{\partial Y}\right)+\frac{\partial }{\partial Z}\left(\mu \frac{\partial u}{\partial W}\right)\right]+S_u$

(3)$\rho \left[\frac{\partial v}{\partial \tau }+\frac{\partial (uv)}{\partial X}+\frac{\partial (vv)}{\partial Y}+\frac{\partial (vw)}{\partial Z}\right]=\left[\frac{\partial }{\partial X}\left(\mu \frac{\partial v}{\partial X}\right)+\frac{\partial }{\partial Y}\left(\mu \frac{\partial v}{\partial Y}\right)+\frac{\partial }{\partial Z}\left(\mu \frac{\partial v}{\partial W}\right)\right]-\rho g{\alpha }_v\left(T-T_m\right)+S_v$

(4)$\rho \left[\frac{\partial w}{\partial \tau }+\frac{\partial (uw)}{\partial X}+\frac{\partial (vw)}{\partial Y}+\frac{\partial (ww)}{\partial Z}\right]=\left[\frac{\partial }{\partial X}\left(\mu \frac{\partial w}{\partial X}\right)+\frac{\partial }{\partial Y}\left(\mu \frac{\partial w}{\partial Y}\right)+\frac{\partial }{\partial Z}\left(\mu \frac{\partial w}{\partial W}\right)\right]+S_w$

(5)$\rho c\left[\frac{\partial T}{\partial \tau }+\frac{\partial (uT)}{\partial X}+\frac{\partial (vT)}{\partial Y}+\frac{\partial (wT)}{\partial Z}\right]=\left[\frac{\partial }{\partial X}\left(k\frac{\partial T}{\partial X}\right)+\frac{\partial }{\partial Y}\left(k\frac{\partial T}{\partial Y}\right)+\frac{\partial }{\partial Z}\left(k\frac{\partial T}{\partial Z}\right)\right]+\rho \gamma \frac{\partial \epsilon }{\partial t}$

In the above equation, Equation (1) is the equation of conservation of mass, Equations (1)–(4) are the equations of conservation of momentum in different directions, and Equation (5) is the equation of conservation of energy where $S_u,S_v$, respectively, are the thru-body terms of the momentum equation, which are defined as follows:

(6)$S_u=\frac{{\left(1-\epsilon \right)}^2}{\left({\epsilon }^3+B\right)}{A}_{\mathrm{mushy} }u,S_v=\frac{{\left(1-\epsilon \right)}^2}{\left({\epsilon }^3+B\right)}{A}_{\mathrm{mushy} }v,S_w=\frac{{\left(1-\epsilon \right)}^2}{\left({\epsilon }^3+B\right)}{A}_{\mathrm{mushy} }w$

In Equation (6), the correction parameter $B$ is used so that the defining equation is still meaningful at the PCM solid phase, and is generally taken as a very small value, e.g., ${10}^{-10}$, and $A_{\mathrm{m}\mathrm{u}\mathrm{s}\mathrm{h}\mathrm{y}}$ is the hysteresis parameter, which is intended to ensure that the equation for the Cheer momentum remains continuous in the two-phase region, with the liquid phase rate $\epsilon$ from 0 to 1, and its value varies from ${10}^7~{10}^4$. The value varies uniformly.

The initial condition is a uniform temperature field, as shown in Equation (7), and the boundary condition is a pulsed heat flow boundary on one side, where $Square\left(t\right)$ is a square wave pulse function, and a thermal convection boundary on the other side, with a constant convective heat transfer coefficient of $h$ and an incoming coolant temperature of $T_f$. The remaining four sides are adiabatic boundaries, as shown in Equations (8)–(10):

(7)$T=T_0\hspace{1em}t=0$

(8)${\tilde{q}}_w(Square(t))=-k\frac{\partial T}{\partial Y}\left|{ }_{Y=0,Z=\delta }\right.\hspace{1em}t>0$

(9)$-k\frac{\partial T}{\partial X}{|}_{X=0}=0, k\frac{\partial T}{\partial X}{|}_{X=L}=0, -k\frac{\partial T}{\partial Z}{|}_{Z=0}=0, k\frac{\partial T}{\partial Z}{|}_{\mathrm{Z}=\delta }=0\hspace{1em}t\ge 0$

(10)$h(T-T_f)=-k\frac{\partial T}{\partial Y}{|}_{Y=H,Z=\delta }\hspace{1em}t>0$

2.3. Solving Methods

Equations (1)–(5) are the conservation control equations for the problem shown in Figure 1, and Equations (8)–(10) are the boundary conditions for the problem shown in Figure 1. It should be noted that Equation (5) is not continuous at the solid–liquid phase interface, and therefore Equations (11) and (12) need to be constrained to define the temperature and location on either side of the phase interface as follows:

(11)$k_s\frac{\partial T_s}{\partial n}-k_l\frac{\partial T_l}{\partial n}=\rho \Delta h\frac{ds(x,y,z,t)}{dt}$

(12)$T_s(x,y,z,t)=T_l(x,y,z,t)=T_m$

In the above equation, $\frac{ds(x,y,z,t)}{dt}$ characterizes the differential component of the phase interface position with time. After substituting the constraint Equations (11) and (12) into Equations (1)–(5), they are discretized using the second-order backward difference method, respectively, and solved implicitly in the time and geometric domains using the alternating direction (ADI method) and using the finite element method. The grid cells are divided according to the finite volume method, and the boundary layer grid near the wall is locally encrypted. The mesh division result is shown in Figure 3, which is composed of a structured tetrahedral mesh. There are 5 boundary layers with a thickness of y+ = 0.105 mm on the coupling surface between the PCM and coolant. The results of the grid independence check are shown in Figure 4, and it can be seen that the error is less than 3% when the minimum size of the grid is less than 2.5 mm, which meets the calculation requirements.

2.4. Model Validation

The experimental setup utilized to validate the aforementioned PCHEU-LPM simulation model is illustrated in Figure 5 and Figure 6. In Figure 5, the PCHEU-LPM (depicted in Figure 6b) is thermally interfaced within a water-cooled compression cycle system (depicted in Figure 6a). At the top of the PCHEU-LPM, an aluminum alloy plate, partially embedded with copper wires, serves as the heating plate (depicted in Figure 6e) to simulate a 12 W heat source. This heating plate is connected to a direct current power supply. Surrounding the PCHEU-LPM is a 3D-printed framework (depicted in Figure 6c) designed for insulation, composed of cured silicone rubber material. The framework features 1 mm apertures on one side, allowing for the insertion of nine K-type thermocouples (depicted in Figure 6d). The bottom of the framework is hollowed out and directly connected to a copper water-cooling head through screws (depicted in Figure 6f).

The experiment is set as two heating modes: constant heat flux heating and pulsed heat flux heating. In the constant heat flux heating mode, the heating power is set at 3.22 W, and the heating duration is 600 s. In the pulsed heat flux heating mode, the peak heating power is 6.15 W, the valley heating power is 0 W, the heating duration is 30 s, and the cooling duration is 110 s. The inlet temperature of the water-cooling system is set at 12 °C with an inlet flow velocity of 0.3 m/s. Three temperature measurement points are arranged along the y-direction in the center of the PCHEU-LPM surface, corresponding to y = 0 mm, y = 3.5 mm, and y = 7 mm, respectively. The experimentally measured temperatures are compared with the numerical results, and the comparative results are depicted in Figure 7a,b. From Figure 7, the temperature difference between the experiment and simulation is less than 5%. The temperature discrepancy is minimal on the top heating surface and maximal on the bottom cooling surface, attributed to the instability introduced by the cooling capacity of the water-cooling system being controlled by the compressor speed.

3. Simulation Results and Discussion

3.1. Influence of Thermal Boundary Conditions on Temperature and Velocity Fields

The loading form of the combination boundary conditions have a direct impact on the distribution of temperature and phase fields in the PCHEU-LPM. On the top of the PCHEU-LPM, a pulsed heat flux is uniformly applied with a heat flux density of 50,000 W/m², while four different cooling methods are applied at the bottom: 1. Natural cooling (NC), 2. Fan cooling (FC), 3. Water cooling (WC), and 4. Microchannel cooling (MC), with corresponding convective heat transfer coefficients of 20 W/m²·K [17], 100 W/m²·K [18], 1500 W/m²·K [19], and 10,000 W/m²·K [19], respectively. The corresponding interface temperature results are plotted in Figure 8. From Figure 8, it can be observed that the use of natural cooling and fan cooling results in poor heat dissipation, with the interface temperature rapidly rising beyond the permissible temperature range of the device. Water cooling exhibits better heat dissipation performance, with lower interface temperature and smaller temperature fluctuations. For microchannel cooling, due to its excessive heat dissipation coefficient, when the pulsed heat flux at the top is loaded to a valley point, the temperature of the PCM rapidly decreases below the phase transition temperature zone, resulting in an increase in interface temperature fluctuations. This is because the PCM’s thermal conductivity plays a dominant role in the internal total thermal resistance, and the excessively high convective heat transfer coefficients are meaningless when the thermal conductivity of the PCM is limited. For the PCM with a thermal conductivity of approximately 12 W/m·K, water cooling proves to be the most suitable cooling method.

Furthermore, the non-steady-state pulsed heating and uneven water cooling will impact the temperature and phase field distribution of the PCHEU-LPM. When subjecting the PCM to a uniformly distributed constant heat flux at the top hot end and a uniformly distributed constant cold flux at the bottom cold end (with a negative surface heat flux density), the resulting temperature field is depicted in Figure 9a, while the corresponding velocity field is illustrated in Figure 10a. Alternatively, when a uniformly distributed constant heat flux is applied to the top of the PCM, and a cooling fluid undergoes cooling from the left to the right at the bottom, the resulting temperature and velocity fields are presented in Figure 9b and Figure 10b, respectively. In the case of a uniformly distributed pulse heat flux at the top of the PCM and simultaneous cooling of a fluid from left to right at the bottom, the resulting temperature and velocity fields are shown in Figure 9c and Figure 10c. The analysis of Figure 9 and Figure 10 reveals that under the scenario where both the hot and cold ends experience uniformly distributed constant heat flux, the temperature field exhibits a layered structure with uniform variation along the direction of the heat flux. However, when the cold end transitions into a cooling fluid flowing from left to right, heat exchange between the cold fluid and the PCM occurs. This leads to the absorption of heat by the cold fluid throughout the process, resulting in a gradual increase in its temperature. Consequently, the temperature difference between the PCM and the cold fluid decreases, manifesting as a temperature difference between the inlet and outlet and a larger cross-sectional velocity field along the process. In the situation where the hot end experiences a pulse heat flow and the cold end becomes a cooling fluid flowing from left to right, the combined effect of frequently changing heat flow and cold flow induces a transformation in the temperature field of the PCM into an uneven surface. Simultaneously, the size of the velocity field rapidly decreases. Although no apparent flow vortices are generated, localized small vortices with minimal changes in velocity gradient are observed.

3.2. Thermal Response of Incoming Flow Temperature

To investigate the impact of the incoming coolant flow temperature on the temperature response of the PCHEU-LPM heating surface, the coolant’s incoming temperature was varied from 283 K to 303 K while keeping the other parameters constant. Figure 11 depicts the response curves of the heating surface interface temperature $T_h$ for different incoming flow temperatures. Figure 12 presents the maximum temperature values $T_{h,max}$, and the temperature fluctuations ${∆T}_h$ at the interface for different incoming flow temperatures $T_f$. Figure 13 shows the maximum gradient of interface temperature ${\nabla T}_{h,max}$ at the corresponding peak and valley values of pulse thermal excitation with different $T_f$. Figure 14 shows the three-dimensional velocity profiles and phase field contours of the PCHEU-LPM at the same time instances. The results depicted in Figure 11 and Figure 12 indicate that a higher incoming flow temperature leads to a reduced temperature difference between the bottom of the PCHEU-LPM and the cooling fluid, resulting in an elevated overall temperature level and an increase in the interface temperature. However, the temperature fluctuations of the heating surface exhibit a decreasing trend followed by an increasing trend as the incoming flow temperature rises. This observation suggests the existence of an optimal coolant temperature. When the incoming flow temperature is below 283 K, the convective heat transfer temperature difference becomes excessively large, causing the interface temperature of the PCHEU-LPM to rapidly decline by more than 15 K during the pulse cycle when the heat stimulus is zero, dropping below the phase transition temperature. Conversely, when the incoming flow temperature exceeds 283 K, the convective heat transfer temperature difference becomes too small, resulting in a rapid increase in the interface temperature during the pulse cycle when the heat stimulus is 50,000 W/m². In the absence of heat stimuli, effective heat dissipation becomes challenging, leading to heat accumulation. According to Figure 13, the inlet temperature of the coolant mainly affects the uniformity of the interface temperature at the peak moment of thermal excitation. The higher the inlet temperature, the larger the corresponding temperature gradient. Based on the analysis of the results in Figure 12, Figure 13 and Figure 14, it is found that the distribution of the phase zone is the reason for the interface temperature response and temperature uniformity changes. When the interface is in the two-phase region instead of being undercooled or overheated, corresponding to $T_f$ = 283 K in Figure 14, the temperature fluctuation is the lowest. When the flow intensity in the liquid phase region of Figure 14 increases with $T_f$, the temperature non-uniformity in Figure 13 increases with the increase in natural convection intensity.

In summary, the interface temperature fluctuations and temperature uniformity of the PCHEU-LPM are influenced by the motion mode of the phase interface position. The position of the phase interface is related to the inlet temperature of the coolant. When the inlet temperature is appropriate, the phase interface always maintains oscillatory motion near the end face, as shown in Figure 15.

3.3. Thermal Response of Incoming Flow Velocity

To investigate the impact of the incoming coolant flow velocity on the temperature response of the PCHEU-LPM heating surface, the incoming velocity $U_f$ was varied from 0.1 m/s to 2.5 m/s while keeping other parameters constant. Figure 16 presents the response curves of the heating surface interface temperature $T_{\mathrm{h}}$ for different $U_f$, while Figure 17 displays the corresponding maximum temperature values $T_{\mathrm{h},\mathrm{m}\mathrm{a}\mathrm{x}}$ and temperature fluctuations ${∆T}_{\mathrm{h}}$ at the interface. Additionally, Figure 18 shows the maximum gradient of interface temperature ${\nabla T}_{h,max}$ at the corresponding peak and valley values of pulse thermal excitation with different $U_f$. Moreover, Figure 19 provides three-dimensional velocity profiles and phase field contour plots of the PCHEU-LPM at the valley time (t = 400 s) and peak time (t = 450 s) of the pulse cycle for different incoming flow velocities. The results depicted in Figure 16 and Figure 17 indicate that a higher incoming flow velocity leads to a larger convective heat transfer coefficient, resulting in the PCHEU-LPM reaching a temperature that closely aligns with and remains constant at the cooling fluid temperature. Consequently, the peak temperature at the surface decreases, and temperature fluctuations diminish. However, Figure 17 reveals that the suppression effect on temperature fluctuations follows a parabolic curve, gradually diminishing as the incoming flow velocity increases. As shown in Figure 18, ${\nabla T}_{h,max}$ decreases as $U_f$ increases, and the decreasing trend gradually flattens. Although the larger the $U_f$ is, the better the heat dissipation effect; for the PCHEU-LPM, there is a limit to its heat dissipation capacity, which depends on the thermal conductivity of the PCM. The reason is consistent with the analysis in Section 3.1, because the thermal conductivity accounts for the dominant factor in the total thermal resistance of the heat transfer system. Figure 19 illustrates that the distribution of the phase zone affects the temperature fluctuations and uniformity at the interface. When the phase interface approaches the PCM interface as $U_f$ increases, the interface is in a phase transition process and temperature fluctuations are suppressed. The lower the velocity field that characterizes the intensity of natural convection, the smaller the interface temperature gradient, and the better the uniformity.

Overall, the interface position of the PCHEU-LPM is inversely correlated with the inlet velocity of the coolant. The position of the phase interface affects the temperature fluctuations and uniformity of the interface. The higher the inlet flow velocity, the closer the phase interface is to the PCM end face, but there is a certain extreme value, as shown in Figure 20.

3.4. Method for Solving Optimal Values

Existing literature [20] suggests that under simple harmonic boundary condition excitation, the phase interface position X(t) correlates with parameters such as k, L, density ρ, and latent heat Δh. By using the method of separating variables, the interface temperature can be decomposed into T1 and T2, as shown in Equations (13) and (14):

(13)$T_1(x,t)=T_0+A(t)[x-X(t)]+B(t){[x-X(t)]}^2-\frac{\left(T_{\mathrm{c}}-T_0\right)\{x-X(t)\}}{H\{1-X(t)/H+1/Bi\}}$

(14)$T_2(x,t)=-2\frac{\Delta h}{\pi \mathrm{c}}\times {\displaystyle \sum _{n=1}^{\infty }\frac{{(-1)}^n-1}{n}}\times \sin \left(\frac{n\pi }{L}x\right){\displaystyle {\int }_0^t\frac{dx(\tau )}{d\tau }}{e}^{-\frac{a{n}^2{\pi }^2(t-\tau )}{L^2}}d\tau$

In the above equation, the term A(t) and B(t) present, respectively:

(15)$A(t)=\frac{L}{2cX(t)}\left[-1+\sqrt{1-\frac{4c{q}_w^{\prime \prime }X(t)}{kL}}\right]$

(16)$B(t)=\frac{L}{8c{\left\{X(t)\right\}}^2}{\left[-1+\sqrt{1-\frac{4c{q}_w^{\prime \prime }X(t)}{kL}}\right]}^2$

The term $H$ in the above equation needs to be supplemented with the transcendental equation for the position of the phase interface as a constraint:

(17)$\frac{T_{\mathrm{ss} }-T_0}{q_w}=\frac{H}{k_{\mathrm{s}}}\left(\frac{1}{Bi}+\frac{H-x}{H}\right)$

Recognizing that Equation (13) encompasses both non-steady-state boundary conditions and a transcendental term H, rendering a direct solution challenging, the Response Surface Analysis (RSA) method is employed. Initially, parameters associated with temperature fluctuations, as expressed in Equation (18), are identified. Subsequently, a target function for minimizing temperature fluctuations, as depicted in Equation (19), is formulated. Finally, a Response Surface Model (RSM) is established to represent the relationship between the target function and the set of parameters. Analysis results based on Equation (18) indicate that, in addition to the inlet temperature and inlet velocity parameters of the cooling fluid at the heat dissipation boundary, the interface temperature fluctuation is influenced by the thermal conductivity, thickness, density, and latent heat of the phase change material (PCM). To simplify the model, sensitivity analysis is conducted on these four sets of material properties with respect to temperature fluctuations, yielding the results presented in Figure 21. According to the results depicted in Figure 21, parameters ρ and λ are identified as highly sensitive, while α and Δ belong to the category of low-sensitivity parameters. Consequently, ρ and λ are treated as fixed values, and an RSM model is exclusively established for α and Δ, as illustrated in Figure 22. Ultimately, by identifying the parameter set corresponding to the lowest point in Figure 22, the optimal values for the rectangular PCHEU are determined and documented in Table 2. By comparing with the unoptimized interface temperature curve, the temperature fluctuation value decreased to 77.9% after adopting the optimal parameter value, as shown in Figure 23.

(18)$\Delta T=f_1(k,L,\rho ,\Delta h)*f_2(T_f,U_f)$

(19)$OBJ=min(\Delta T)$

4. Conclusions

In this paper, we present a three-dimensional model of a heat transfer system, which is coupled with pulsed heating and forced convective heat dissipation boundary conditions, to investigate the temperature response of the heated interface. The key findings are as follows:

(1) The emergence of high thermal conductivity phase change materials has made it possible for PCMs to be used as a substrate for heat exchangers rather than heat accumulators.

(2) It is important to choose the appropriate heat dissipation method for phase change material heat exchangers. This article selects the PCHEU-LPM with a thermal conductivity of 12 W/m·K, and the most suitable heat dissipation method is water cooling.

(3) The temperature fluctuation and uniformity of the PCM interface are directly related to the position of the phase interface, and the temperature fluctuation is minimized when the oscillation motion of the phase interface is always maintained near the interface under pulse thermal excitation. The position of the phase interface is related to the inlet temperature and inlet velocity of the coolant. The inlet temperature cannot be too high or too low, and there is an optimal value. The higher the inlet speed, the better, but there is a limit, beyond which the temperature control effect no longer changes.

(4) The optimal values of coolant inlet temperature and inlet velocity are also related to the thermal conductivity and thickness of PCMs. Through response surface methodology, the optimal parameter set can be obtained with the goal of minimizing temperature fluctuations.

By investigating these aspects, our study enhances the understanding of composite heat transfer systems and provides valuable insights for optimizing the design and performance of phase change heat exchangers.

Author Contributions

Conceptualization, J.-Y.Z.; Methodology, R.-J.Z.; Software, R.-J.Z.; Validation, R.-J.Z.; Formal Analysis, R.-J.Z. and J.-Y.Z.; Investigation, R.-J.Z. and J.-Y.Z.; Resources, J.-Y.Z.; Data Curation, J.-Y.Z.; Writing—Original Draft Preparation, R.-J.Z.; Writing—Review & Editing, J.-Z.Z.; Project Administration, J.-Z.Z. Funding Acquisition, J.-Z.Z. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

Parameter name
$L$	Length	mm
$H$	Width	mm
$\delta$	Height	mm
$q$	Heat flow	W/m2
$T$	Temperature	K
$U$	X-direction speed	m/s
$V$	Y-direction speed	m/s
$W$	Z-direction speed	m/s
$S$	Body-force	N/m3
S	Phase interface location	m
$h$	Convective heat transfer coefficient	W/m2-K
Δh	Latent heat of phase change	KJ/kg
Greek alphabet
$\tau$	Time	s
$\epsilon$	Porosity or liquid phase ratio
$\mu$	Sport viscosity	Pa-s
$\rho$	Density	Kg/m3
$k$	Thermal conductivity	W/m-k
${\rm Y}$	Latent heat of phase change	KJ/kg
$\alpha$	Thermal diffusion coefficient	m2/s
$\phi$	Liquid fraction
Subscript
L	Liquid phase
$s$	Solid phase
$n$	Directional vector at the phase interface
W	External thermal excitation
M	Phase change state
C	Cold fluids
0	Initial state
Ss	Steady state, time-averaged value of transient temperature fluctuations
Abbreviation
PCHEU-LPM	Phase Change Heat Exchange Unit with Layered Porous Media
NC	Natural cooling
FC	Fan cooling
WC	Water cooling
MC	Microchannel cooling

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Thermodynamic Cycle Layout of MR2 Aircraft.

Figure 2. Schematic diagram of the heat transfer process with phase change. (a) Three-dimensional heat transfer process; (b) Two-dimensional cross-section diagram in the direction of heat transfer process.

Figure 3. Structured grid division for PCHEU-LPM and water channels.

Figure 4. Grid independence test results.

Figure 5. Experimental setup for the validation of PCHEU-LPM model effectiveness.

Figure 6. Photographs of components inside the test section.

Figure 7. Comparison between experimental temperature results and numerical results. (a) Constant heating mode; (b) Pulse heating mode.

Figure 8. Influence of different heat dissipation methods on the hot end temperature of PCHEU-LPM.

Figure 9. Temperature field distribution under different boundary condition combinations. (a) Constant heat flux with constant cold flow; (b) Constant heat flux with coolant flow; (c) Pulsed heat flux with coolant flow.

Figure 10. Velocity field distribution under different boundary condition combinations. (a) Constant heat flux with constant cold flow; (b) Constant heat flux with coolant flow; (c) Pulsed heat flux with coolant flow.

Figure 11. Response curves of the heating surface interface temperature [Forumla omitted. See PDF.] for different incoming flow temperatures.

Figure 12. Maximum temperature values [Forumla omitted. See PDF.] and temperature fluctuations [Forumla omitted. See PDF.] at the interface for different incoming flow temperatures.

Figure 13. Maximum gradient of interface temperature [Forumla omitted. See PDF.] with different [Forumla omitted. See PDF.].

Figure 14. Velocity profiles and phase field contours at 450 s of the pulse cycle, considering various incoming flow temperatures.

Figure 15. The variation of non-dimensional phase interface position with different Tf.

Figure 16. Response curves of the heating surface interface temperature Th for different incoming flow velocities.

Figure 17. Maximum temperature values [Forumla omitted. See PDF.] with different [Forumla omitted. See PDF.] [Forumla omitted. See PDF.] and temperature fluctuations [Forumla omitted. See PDF.] at the interface for different incoming flow velocities.

Figure 18. Maximum gradient of interface temperature [Forumla omitted. See PDF.] with different Uf.

Figure 19. Velocity profiles and phase field contours at 450 s of the pulse cycle, considering various incoming flow velocities.

Figure 20. The variation of non-dimensional phase interface position with different [Forumla omitted. See PDF.].

Figure 21. Parameter sensitivity analysis chart.

Figure 22. Response Surface Analysis results.

Figure 23. Comparison of interface temperature results between optimizing parameters with benchmark parameters.

Table 1

Thermal physical parameters of materials.

Name of Material	Thermal Conductivity (W/m-K)	Density(kg/m3)	Specific Heat Capacity(J/kg-K)	Latent Heat(kJ/kg)	Melting Temperature(K)
Paraffin	0.558	900	2170	220	323
Copper foam	385	7900	3900	/	/

Table 2

Calculation results of optimal parameters.

Parameter Name	Value	Maximum Temperature	Temperature Fluctuation
k	15.53 (W/mK)	313.2 K	11.5 K
L	0.042 m
Tf	296 K
Uf	2.5 m/s

References

1. Gou, J.; Xiao, S.; Hu, J.; Gao, G.; Gong, C. Research Progress of Aerodynamic Heat Dissipation, Transport and Conversion Technologies of Hypersonic Vehicles. Yuhang Xuebao/J. Astronaut.; 2022; 43, pp. 983-999. [DOI: https://dx.doi.org/10.3873/j.issn.1000-1328.2022.08.001]

2. Behrens, B.; Müller, M. Technologies for Thermal Protection Systems Applied on Re-Usable Launcher. Acta Astronaut.; 2004; 55, pp. 529-536. [DOI: https://dx.doi.org/10.1016/j.actaastro.2004.05.034]

3. Gori, F.; Corasaniti, S.; Worek, W.M.; Minkowycz, W.J. Theoretical Prediction of Thermal Conductivity for Thermal Protection Systems. Appl. Therm. Eng.; 2012; 49, pp. 124-130. [DOI: https://dx.doi.org/10.1016/j.applthermaleng.2011.07.012]

4. Mehta, J.; Charneski, J.; Wells, P. Unmanned Aerial Systems (UAS) Thermal Management Needs, Current Status, and Future Innovations. Proceedings of the 10th International Energy Conversion Engineering Conference; Atlanta, GA, USA, 30 July–1 August 2012; ISBN 978-1-62410-190-8

5. Balland, S.; Fernandez Villace, V.; Steelant, J. Thermal and Energy Management for Hypersonic Cruise Vehicles—Cycle Analysis. Proceedings of the 20th AIAA International Space Planes and Hypersonic Systems and Technologies Conference; Glasgow, UK, 6–9 July 2015; American Institute of Aeronautics and Astronautics: 2015.

6. Lee, S.H.; Mudawar, I.; Hasan, M.M. Thermal Analysis of Hybrid Single-Phase, Two-Phase and Heat Pump Thermal Control System (TCS) for Future Spacecraft. Appl. Therm. Eng.; 2016; 100, pp. 190-214. [DOI: https://dx.doi.org/10.1016/j.applthermaleng.2016.01.018]

7. Guo, W.; Li, Y.; Li, Y.-Z.; Zhong, M.-L.; Wang, S.-N.; Wang, J.-X.; Li, E.-H. An Integrated Hardware-in-the-Loop Verification Approach for Dual Heat Sink Systems of Aerospace Single Phase Mechanically Pumped Fluid Loop. Appl. Therm. Eng.; 2016; 106, pp. 1403-1414. [DOI: https://dx.doi.org/10.1016/j.applthermaleng.2016.06.140]

8. Waser, R.; Ghani, F.; Maranda, S.; O’Donovan, T.S.; Schuetz, P.; Zaglio, M.; Worlitschek, J. Fast and Experimentally Validated Model of a Latent Thermal Energy Storage Device for System Level Simulations. Appl. Energy; 2018; 231, pp. 116-126. [DOI: https://dx.doi.org/10.1016/j.apenergy.2018.09.061]

9. Maxa, J.; Novikov, A.; Nowottnick, M.; Heimann, M.; Jarchoff, K. Phase Change Materials for Thermal Peak Management Applications with Specific Temperature Ranges. Proceedings of the 2018 17th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm); San Diego, CA, USA, 29 May–1 June 2018; IEEE: 2018; pp. 92-101.

10. Krishnan, S.; Garimella, S.V. Analysis of a Phase Change Energy Storage System for Pulsed Power Dissipation. IEEE Trans. Comp. Packag. Technol.; 2004; 27, pp. 191-199. [DOI: https://dx.doi.org/10.1109/TCAPT.2004.825758]

11. Shamberger, P.J.; Hoe, A.; Deckard, M.; Barako, M.T. Effects of Boundary Conditions on the Dynamic Response of a Phase Change Material. Proceedings of the 2020 19th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm); Orlando, FL, USA, 21–23 July 2020; IEEE: 2020; pp. 985-992.

12. Kim, T.Y.; Hyun, B.-S.; Lee, J.-J.; Rhee, J. Numerical Study of the Spacecraft Thermal Control Hardware Combining Solid–Liquid Phase Change Material and a Heat Pipe. Aerosp. Sci. Technol.; 2013; 27, pp. 10-16. [DOI: https://dx.doi.org/10.1016/j.ast.2012.05.007]

13. Gilmore, D. Spacecraft Thermal Control Handbook; Aerospace Press: Chantilly, VA, USA, 2002; ISBN 978-1-884989-11-7

14. Lin, Y.; Jia, Y.; Alva, G.; Fang, G. Review on Thermal Conductivity Enhancement, Thermal Properties and Applications of Phase Change Materials in Thermal Energy Storage. Renew. Sustain. Energy Rev.; 2018; 82, pp. 2730-2742. [DOI: https://dx.doi.org/10.1016/j.rser.2017.10.002]

15. Gürel, B. A Numerical Investigation of the Melting Heat Transfer Characteristics of Phase Change Materials in Different Plate Heat Exchanger (Latent Heat Thermal Energy Storage) Systems. Int. J. Heat Mass Transf.; 2020; 148, 119117. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2019.119117]

16. Youssef, W.; Ge, Y.T.; Tassou, S.A. CFD Modelling Development and Experimental Validation of a Phase Change Material (PCM) Heat Exchanger with Spiral-Wired Tubes. Energy Convers. Manag.; 2018; 157, pp. 498-510. [DOI: https://dx.doi.org/10.1016/j.enconman.2017.12.036]

17. Huang, Z.; Xie, N.; Zheng, X.; Gao, X.; Fang, X.; Fang, Y.; Zhang, Z. Experimental and Numerical Study on Thermal Performance of Wood’s Alloy/Expanded Graphite Composite Phase Change Material for Temperature Control of Electronic Devices. Int. J. Therm. Sci.; 2019; 135, pp. 375-385. [DOI: https://dx.doi.org/10.1016/j.ijthermalsci.2018.09.031]

18. Fan, L.-W.; Xiao, Y.-Q.; Zeng, Y.; Fang, X.; Wang, X.; Xu, X.; Yu, Z.-T.; Hong, R.-H.; Hu, Y.-C.; Cen, K.-F. Effects of Melting Temperature and the Presence of Internal Fins on the Performance of a Phase Change Material (PCM)-Based Heat Sink. Int. J. Therm. Sci.; 2013; 70, pp. 114-126. [DOI: https://dx.doi.org/10.1016/j.ijthermalsci.2013.03.015]

19. Zhou, J.; Cao, X.; Zhang, N.; Yuan, Y.; Zhao, X.; Hardy, D. Micro-Channel Heat Sink: A Review. J. Therm. Sci.; 2020; 29, pp. 1431-1462. [DOI: https://dx.doi.org/10.1007/s11630-020-1334-y]

20. Mazzeo, D.; Oliveti, G.; De Simone, M.; Arcuri, N. Analytical Model for Solidification and Melting in a Finite PCM in Steady Periodic Regime. Int. J. Heat Mass Transf.; 2015; 88, pp. 844-861. [DOI: https://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.04.109]