Numerical studies on thermal performance of a novel edge cooling design in fuel cell stacks for unmanned aerial vehicles

1 Introduction

SUnmanned aerial vehicles (UAVs) have had a broad spectrum of applications in the past decade and have received considerable attention (Gong et al., 2022; Mohsan et al., 2022; Idrissi et al., 2022). Conventional batteries, traditionally used as UAV power sources, cannot provide long endurance without a substantial increase in system mass (Gadalla and Zafar, 2016). As an alternative power source for UAVs, proton exchange membrane fuel cells (PEMFCs) are seen to be more promising owing to their zero emissions, high electrical conversion efficiencies (Kim et al., 2011), as well as high specific energy compared to existing batteries (Gong et al., 2017). Especially, PEMFCs can greatly extend the endurance of UAVs as a result of rapid startup and low operating temperature (El-Sharkh et al., 2004). Statistics indicate that the endurance of UAVs powered by PEMFCs increases by 4–5 times compared with existing UAVs based on batteries (Rottmayer and Miller, 2011). The optimum operating temperature for current PEMFCs ranges from 60 °C to 100 °C (Ravishankar and Arul Prakash, 2014). Either higher temperature or lower temperature has an adverse effect on the performance of fuel cells (Sasmito et al., 2011). To ensure a uniform reaction during the PEMFC operation, the temperature difference of the active zone should be less than 5 °C–10 °C (Barbir, 2013). Therefore, it is important to take measures (e.g.,: optimize the flowfield of bipolar plates (Toghyani et al., 2021) and improve the cooling strategies) to effectively remove the heat generated in a fuel cell and provide a uniform temperature distribution for the active area.

Effective cooling of PEMFC stacks is challenging, especially for UAV applications, in which the power system requires a simple structure and light weight. For UAV applications, the cooling system must maintain the stack within an allowable temperature range while adding as little mass, volume, and auxiliary power consumption as possible. Because pumps, fans, heat exchangers, and flow passages contribute to both balance-of-plant mass and parasitic power, the cooling architecture directly affects system efficiency, cost, reliability, and achievable UAV endurance.

Several PEMFC cooling strategies have been investigated, including increased cathode airflow (Yim et al., 2008; Sohn et al., 2005), separate cooling channels (Wu et al., 2009; Jung et al., 2008; Santarosa et al., 2007), fins or edge cooling (Flückiger et al., 2007), liquid cooling using deionized water or antifreeze (Eckl et al., 2004; Rodatz et al., 2004; Laurencelle et al., 2015), and phase-change cooling (Faghri and Guo, 2005; Zhang and Kandlikar, 2012). Figure 1 compares a liquid-cooling system with an air-cooled heat-spreader system. The latter can reduce system mass and complexity while improving reliability (Burke and Colozza, 2010; Kösters and Gao, 2021).

Three air-cooling strategies are particularly relevant to UAV PEMFC stacks: (i) area cooling, in which cooling air passes through dedicated stack passages; (ii) the use of excess cathode air as both reactant and coolant; and (iii) edge cooling, in which air removes heat from the exposed edges of thermally conductive plates, as shown in Figure 2a. Dedicated cooling passages increase stack complexity, whereas excess cathode airflow may promote membrane drying and require substantial pumping power (Sohn et al., 2005; Flückiger et al., 2007). Edge cooling, also described as cooling with heat spreaders, conducts heat from the active region toward the plate edges, where it is removed by air (Flückiger et al., 2007; Zhang and Kandlikar, 2012). Its comparatively simple balance-of-plant architecture makes it attractive for lightweight UAV applications.

Flückiger et al. (2007) analyzed and optimized a 500 W edge-air-cooled PEMFC stack using graphite separator plates as heat spreaders, while Sasmito et al. (2011) numerically evaluated several stack thermal-management strategies. These studies established the feasibility of edge cooling but primarily considered airflow supplied from terminal regions. As stack length and thermal load increase, longer in-plane conduction and airflow paths may increase temperature non-uniformity and flow resistance. The present study therefore examines a distributed periodic arrangement that divides a long cooling edge into repeated sections.

This study evaluates the thermal-hydraulic performance of a periodic edge-cooling configuration in which equally spaced fans are arranged along the long sides of the stack, as shown in Figure 2b. Because the geometry and airflow pattern repeat along the longitudinal direction, a representative unit containing one complete fan opening and two-half fan openings is extracted to reduce computational cost while preserving the local interaction between adjacent fans. The specific objectives are to: (i) quantify the temperature, pressure, and velocity fields produced by the periodic fan arrangement; (ii) determine how heat flux and fan spacing affect active-area temperature non-uniformity and flow resistance; (iii) evaluate the required cooling airflow, pressure drop, and auxiliary airflow hydraulic power under a common maximum-temperature constraint; and (iv) identify a low-power fan-spacing region within the investigated design space. The identified spacing is interpreted as a configuration-specific design region for the present geometry and operating range, rather than as a universal optimum for PEMFC stacks of different dimensions.

2 Problem definition

2.1 Periodic edge-cooling configuration

Figure 3 shows the drawing of the plate design, as well as the top view and side view of the plate. The configuration represents an interior periodic unit of a stack in which fan openings are distributed at equal intervals along the long cooling edges. As can be seen, heat flux is applied to half of the active zone, which represents the heat generated in the fuel cell during operation. The generated heat conducts to the edge along the graphite plate. Air is supplied through the fan openings into the cooling passage, exchanges heat with the plate edge, and then leaves through the outlet region between adjacent fans, as shown in Figures 3b,c. Because the geometry and fan arrangement repeat in the longitudinal direction, the computational domain contains one complete fan opening at the center and two-half fan openings at the longitudinal boundaries. Together, these openings rep-resent one periodic fan-spacing unit and preserve the local thermal-fluid interaction between adjacent fans while reducing computational cost. The model represents an interior repeated unit of a sufficiently long and regularly arranged stack; end effects near the first and last fans of a finite stack are not included.

2.2 Heat-generation model

The heat-generation rate of a single PEMFC can be approximately estimated as follows (Zhang and Kandlikar, 2012):

q ″ = E t h − E c e l l × i (1)

Here, E_th is the theoretical cell voltage, and Ecell and i are the operating voltage and current density, respectively. Eth depends on the phase of the product water. Under standard conditions (25 °C and 1 atm), Eth is approximately 1.48 V for liquid water and 1.25 V for water vapor (Zaidi and Rauf, 2016). The operating point of 0.65 V and 1.0 A cm⁻² with gaseous product water is used only to obtain a representative baseline thermal load: Equation 1 gives 6000 W m⁻², of which 3000 W m⁻² is assumed to conduct toward the modeled side of the bipolar plate. The five imposed heat fluxes are therefore treated as equivalent thermal loads rather than a complete electrochemical description of all UAV operating conditions.

2.3 Design variables and thermal constraint

In this study, the heat flux (namely, heat generation rate), marked as q ″ (q ″ = 2,400, 2,700, 3,000, 3,300, 3600 W m⁻²), and the clear spacing between adjacent fan openings, marked as d (d = 20, 30, 40, 50, 60 mm), were set as independent variables to show the thermal performance of this edge-cooling strategy. The five heat-flux values are treated as equivalent thermal loads covering changes in PEMFC operating conditions; they do not explicitly resolve altitude-dependent electrochemical or water-management behavior. Figure 3c shows the dimensions of the cooling plate with a plate thickness of 0.5 mm and an active area width of 30 mm. Half of the active area, with a width of 15 mm, was considered the basic simulated unit. Within the simulated unit, the width and height of the cooling fluid field were 20 mm and 0.25 mm respectively with a complete fan in the middle and two-half fans at the edges. The corresponding representative-unit lengths are 160, 180, 200, 220, and 240 mm for d = 20, 30, 40, 50, and 60 mm, respectively. To analyze the fluid flow and the heat transfer processes, CFD software Fluent was used. Accordingly, the dimensional spacing values reported in this study are specific to the stated fan width, cooling-passage height, cooling-edge width, plate thickness, material properties, and thermal constraint; direct scaling to other stack sizes requires reassessment of these geometric and flow ratios.

For each combination of q″ and d, the inlet velocity was successively adjusted in the steady CFD calculations until the maximum active-area temperature satisfied |T _max − 358.15 K| ≤ 0.1 K. This procedure allows the required airflow, pressure loss, and airflow hydraulic power to be compared under the same maximum-temperature constraint. The value 358.15 K (85 °C) is adopted as the design ceiling for the present numerical study and lies within the commonly reported PEMFC operating-temperature range (Ravishankar and Arul Prakash, 2014); it should not be interpreted as a universal degradation threshold or a manufacturer-specific limit for all PEMFC stacks.

3 Numerical simulation

3.1 Governing equations and assumptions

To simplify the calculation, the steady, incompressible, viscous and laminar flow equations were used to study the fluid field, shown from Equations 2-5. To obtain the temperature fields, the energy equation, shown as Equation 6, was employed.

∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z = 0 (2)

u ∂ u ∂ x + v ∂ u ∂ y + w ∂ u ∂ z = − 1 ρ a ∂ p ∂ x + v a ∂ 2 u ∂ x 2 + ∂ 2 u ∂ y 2 + ∂ 2 u ∂ z 2 (3)

u ∂ v ∂ x + v ∂ v ∂ y + w ∂ v ∂ z = − 1 ρ a ∂ p ∂ y + v a ∂ 2 v ∂ x 2 + ∂ 2 v ∂ y 2 + ∂ 2 v ∂ z 2 (4)

u ∂ w ∂ x + v ∂ w ∂ y + w ∂ w ∂ z = − 1 ρ a ∂ p ∂ z + v a ∂ 2 w ∂ x 2 + ∂ 2 w ∂ y 2 + ∂ 2 w ∂ z 2 (5)

u ∂ T ∂ x + v ∂ T ∂ y + w ∂ T ∂ z = 1 ρ a k a C p , a ∂ 2 T ∂ x 2 + ∂ 2 T ∂ y 2 + ∂ 2 T ∂ z 2 (6)

Where u, v, and w are velocities in the x, y, and z directions respectively and T is temperature. To improve durability, the allowed maximum temperature of the stack is 85 °C. The value 358.15 K is used as the common design ceiling for comparing all cases and is close to the upper end of the usual low-temperature PEMFC operating range (Barbir, 2013; Zhang and Kandlikar, 2012). It is not presented as a universal degradation threshold or as a manufacturer-specific limit.

The constant material properties actually adopted in the archived steady Fluent calculations are listed in Table 1. Because the solid-domain equation is steady, only the graphite thermal conductivity affects the solid conduction solution; solid density and heat capacity are therefore not used in the present steady-state results.

TABLE 1. Constant properties adopted in the archived steady CFD calculations.

The baseline model adopts the following assumptions: the airflow is incompressible, viscous, and laminar; gravitational buoyancy is neglected because forced convection is dominant; the graphite plate is homogeneous and isotropic; the active-area heat flux is spatially uniform; thermal contact resistance, radiation, external natural convection, and leakage are not included; and all external solid boundaries outside the prescribed cooling path are treated as adiabatic. The model therefore represents an ideal periodic cooling unit under geometrically symmetric operating conditions.

In the graphite plate, steady heat conduction was represented by ∇·(kp∇Tp) = 0, with the prescribed active-area heat flux and conjugate temperature and heat-flux continuity at the fluid–solid interface.

3.2 Boundary conditions and fan representation

The boundary conditions used for solving the momentum and energy equations are shown below.

Each fan opening was represented by a prescribed uniform velocity inlet normal to the inlet surface. The physical fan blades and rotating flow domain were not explicitly resolved, and no manufacturer-specific pressure-flow curve or efficiency map was imposed. The CFD model therefore determines the airflow and pressure requirement associated with each geometric configuration, which can subsequently be matched to the operating curve of a selected commercial fan.

Inlet: u = ν i n, v = 0, w = 0, T = T_in

Outlet: ∂ u ∂ n = 0, ∂ v ∂ n = 0, ∂ w ∂ n = 0, p = 0

Active area walls: u, v, w = 0, k p · ∂ T ∂ z = q ″

Other walls: u, v, w = 0.

Interface: u, v, w = 0, T p = T a , k p ∂ T p ∂ n = k a ∂ T a ∂ n

A uniform inlet temperature of 300 K was imposed at each fan opening, and the outlet boundaries were specified at the ambient reference pressure. No-slip conditions were applied to the fluid-solid walls, and the prescribed heat flux was imposed on the active-area surface. The right and left boundaries of the representative unit coincide with geometric midplanes that bisect identical half-fan openings. Under equal inlet conditions, these planes were treated as symmetry boundaries to represent an interior periodic unit. This simplification applies to the ideal regularly repeated configuration and does not include fan-to-fan variation, geometric tolerances, or local leakage.

For consistent post-processing, the static-pressure difference, volumetric flow rate, representative-unit airflow hydraulic power, and estimated electrical fan power are defined by Equations 7–10. The same complete inlet and outlet surface sets and the same area-averaging procedure were used for every configuration.

Δ p = p ¯ in − p ¯ out (7)

V ˙ = ∫ A in u · n dA (8)

P hyd = Δ p V ˙ (9)

P fan = P hyd η fan = Δ p V ˙ η fan (10)

Here, p in and p out are the area-averaged static pressures over the complete inlet and outlet surface sets, respectively; V ˙ is the total volumetric flow rate through the complete inlet surfaces of one representative periodic unit; and η fan is the overall fan efficiency. The airflow hydraulic power, P hyd = Δ p V ˙, represents the mechanical power transferred to the cooling air to overcome the pressure loss of the cooling passage. It should not be interpreted as the actual electrical power consumption of a commercial fan. The corresponding electrical fan power may be estimated as P fan , elec = P hyd / η fan. Because no manufacturer-specific pressure–flow–efficiency curve was imposed, P hyd is used only as a relative airflow-power metric for comparing the investigated configurations. Determination of the actual electrical fan power would require a documented commercial-fan performance curve and efficiency data.

ANSYS Fluent was used to solve the steady momentum and energy equations for the coupled fluid and solid domains using the SIMPLEC pressure–velocity coupling algorithm, second-order pressure interpolation, and second order upwind discretization for momentum and energy. Convergence was assessed from the residual histories together with stabilization of the maximum active-area temperature and the inlet and outlet surface-averaged pressures. The same numerical treatment was applied to all heat-flux and fan-spacing cases.

The right and left boundaries of the representative unit coincide with geometric midplanes that bisect identical half-fan openings. These planes were imposed as symmetry boundaries to represent an ideal interior periodic unit. No full non-symmetric or two-fan reference-domain calculation was performed to quantify the sensitivity of the predicted flow and temperature fields to this assumption. Accordingly, the present results apply strictly to a regularly repeated geometry with identical fan inlets, spatially uniform heat loading, and no leakage. Finite-stack end effects, fan-to-fan mismatch, manufacturing deviations, and asymmetric leakage are outside the scope of the model.

3.3 Mesh generation and quality assessment

The solid and fluid domains were discretized using the MultiZone algorithm in ANSYS ICEM, and hexahedral elements were employed. The maximum and average skewness values were below 0.81 and 0.30, respectively. These metrics characterize element quality but do not establish mesh independence.

For fan spacings of 20, 30, 40, 50, and 60 mm, the computational domains contained 8,127,202, 9,191,273, 10,176,000, 12,345,892, and 13,438,286 elements, respectively. The increase in cell count results from the increased domain length while maintaining a common meshing strategy. These geometry-dependent cell counts should not be interpreted as coarse, medium, and fine meshes for a convergence assessment.

A formal three-level mesh-sensitivity assessment of the original three-dimensional conjugate model could not be reconstructed because the archived solution-field and mesh-sensitivity output files are no longer available. Consequently, the spatial-discretization uncertainty of the reported temperature and pressure predictions has not been quantified.

The hydraulic diameter of the rectangular cooling passage and the inlet Reynolds number are defined by Equations 11, 12, where W and H are the passage width and height, Ac is the cross-sectional area, Pw is the wetted perimeter, ūin is the area-averaged inlet velocity, and νa is the kinematic viscosity of air listed in Table 1.

D h = 4 A c P w = 2 WH W + H (11)

Re = u ¯ in D h ν a (12)

For W = 20 mm and H = 0.25 mm, the hydraulic diameter is 0.493827 mm. The maximum inlet velocity of approximately 10 m s⁻¹ gives an inlet Reynolds number of approximately 338. However, the local velocity inside the passage is substantially increased by flow convergence, turning, and local contraction. For the most demanding case, d = 20 mm and q″ = 3600 W m⁻², Figure 4d gives a maximum local velocity of approximately 41.7 m s⁻¹. The corresponding hydraulic-diameter-based local Reynolds number is approximately 1.41 × 10³.

This value remains below the conventional transition criterion for fully developed internal flow, but the maximum velocity occurs in a locally turning and contracting flow region. The calculation therefore supports the order of magnitude of the laminar assumption but does not constitute a turbulence-model sensitivity assessment.

3.4 Verification scope and uncertainty limitations

No directly comparable experimental dataset or independently validated numerical dataset was available for the present periodic fan-and-groove geometry, thin cooling passage, graphite plate, and temperature-constrained inlet-velocity procedure. In addition, a formal mesh-convergence assessment of the original three-dimensional model could not be reconstructed from the archived files. Therefore, the reported temperature and pressure predictions do not have independently quantified validation or spatial-discretization uncertainties.

The numerical results are consequently interpreted as comparative trends obtained using a consistent geometry-dependent meshing strategy, solver setup, convergence procedure, and maximum-temperature constraint. They should not be interpreted as experimentally validated or mesh-independent absolute predictions for a complete PEMFC stack.

4 Results and discussion

A total of 25 baseline cases were simulated using five heat fluxes and five fan spacings. The inlet velocity was independently adjusted for each case until the maximum active-area temperature reached 358.15 ± 0.1 K. Consequently, differences among the cases represent the airflow, pressure loss, and airflow hydraulic power required to satisfy the same maximum-temperature constraint rather than performance under a common inlet velocity.

4.1 Temperature distribution

Figure 5 shows the temperature distributions in the fluid and solid domains for the representative case of q ″ = 3600 W m − 2 and d = 60 mm. As shown in Figure 5a, relatively low air temperatures occur near the fan inlets, whereas the air temperature gradually increases along the flow path as heat is absorbed from the cooling edge. The maximum air temperature is observed near the outlet region between adjacent fans. The maximum temperature of approximately 348 K reported in Figure 5a refers specifically to the fluid domain and should not be confused with the maximum active-area temperature in the solid plate.

As shown in Figure 5b, the solid temperature distribution exhibits a periodic pattern corresponding to the periodic fan arrangement. Lower solid temperatures occur opposite the fan inlets because of the stronger local convective cooling at these locations. In contrast, the highest active-area temperatures occur approximately midway between adjacent fans and away from the cooling edge, where the heat-conduction path is relatively long and the local cooling effect is weaker. The resulting temperature field reflects the combined effects of in-plane heat conduction through the graphite plate and spatially non-uniform forced convection along the cooling edge.

Figure 6a compares the active-area temperature distributions under different heat fluxes while the fan spacing is fixed at d = 40 mm. In all cases, the highest temperature occurs approximately midway between adjacent fan openings and away from the cooling edge, whereas lower temperatures are observed opposite the fan inlets. This spatial pattern results from the non-uniform convective cooling along the edge and the different in-plane heat-conduction distances from the active area to the cooling boundary.

For each heat-flux condition, the inlet velocity was independently adjusted to maintain the maximum active-area temperature at 358.15 ± 0.1 K. Therefore, increasing the heat flux does not produce a substantial increase in T max, because the additional heat load is compensated by a higher cooling-air velocity. Instead, the enhanced local convection reduces the temperature more strongly in the regions opposite the fan inlets. Consequently, the local minimum temperature decreases and the active-area temperature difference increases. The apparent reduction in the high-temperature area should therefore not be interpreted independently as an improvement in overall temperature uniformity.

Figure 6b compares the active-area temperature distributions under different fan spacings at a fixed heat flux of q ″ = 3000 W m − 2. For all investigated spacings, the relatively low-temperature regions are located opposite the fan inlets, while the highest-temperature regions occur approximately midway between adjacent fans. As the fan spacing increases, the distance between the directly cooled regions becomes larger and a higher inlet velocity is required to maintain the prescribed maximum-temperature limit.

The increased airflow strengthens the local cooling near the fan-facing regions and expands the low-temperature area. Meanwhile, the maximum temperature remains close to 358.15 K because it is imposed as the common thermal constraint. Therefore, the increase in the active-area temperature difference with fan spacing is mainly caused by the reduction in the minimum temperature rather than by an increase in the maximum temperature. This result indicates that a larger fan spacing can satisfy the maximum-temperature requirement but may produce a less uniform active-area temperature distribution.

4.2 Fluid pressure distribution

Figure 7 shows the pressure distribution in the cooling passage for the representative case of q ″ = 3000 W m − 2 and d = 40 mm. Relatively high pressure occurs near the velocity-inlet boundaries, whereas lower pressure is observed toward the outlet region between adjacent fans. This pressure variation drives the cooling air from the fan openings toward the discharge boundaries.

The pressure field is governed by the combined effects of flow acceleration, flow turning, viscous friction, and local geometric resistance. In particular, the sealing groove partially obstructs the incoming air and forces the flow to change direction before reaching the outlet, producing local pressure gradients and additional form losses. Therefore, the pressure distribution cannot be interpreted using the Bernoulli relation alone. The periodic pressure pattern along the longitudinal direction is primarily associated with the repeated arrangement of the fan inlets, outlets, and sealing-groove geometry.

Figure 8 compares the pressure distributions at q ″ = 3000 W m⁻² for different fan spacings. At the smallest spacing of d = 20 mm, the available discharge region between adjacent fan openings is relatively restricted, and the incoming streams interact strongly within a short distance. Together with the flow turning around the sealing groove, these effects produce relatively large local pressure gradients.

As the fan spacing increases, the pressure variation does not change monotonically. A moderate increase in spacing enlarges the effective outlet region and weakens the direct interaction between adjacent inlet streams, which tends to reduce the hydraulic resistance. However, a larger spacing also increases the characteristic flow-path length and requires a higher inlet velocity to maintain the prescribed maximum active-area temperature. At sufficiently large spacing, the effects of the increased flow rate, wall friction, and longer flow path become increasingly important. The competition between the enlarged discharge area and the increased flow requirement therefore produces the observed minimum-pressure-drop region.

4.3 Fluid velocity distribution

As shown in Figure 9a, it can clearly be observed that the maximum velocity appears near the outlet between the fans due to the minimum pressure located in the same area. This can be easily explained by the principle of energy conservation that the highest speed occurs where the pressure is the lowest within a horizontal fluid flow. However, the minimum velocity area occurs right in front of the fans and in the middle of the two adjacent fans, where both are near the seal groove. This is attributed to the fact that the edge of the seal groove acts as a block for the cooling air, therefore the air velocity becomes zero and then the air turns its direction near the edge, as shown in Figure 9b. Similarly, the cooling air from the two adjacent fans merges between the adjacent fans in an opposite velocity direction resulting in a zero velocity area.

The fluid velocity contours at q ″ = 3,000 W m⁻² were calculated for all assumptive fan intervals, as shown in Figure 10. It is observed that the area where the fluid velocity is close to zero decreases with the increase in the fan interval. This is due to the fact that the outlet area becomes larger with increasing fan interval and the cooling air flows easier out of the outlet. The velocity gradient of cooling air near the outlet also decreases as the fan interval increases. As the fan spacing increases, the available outlet region becomes larger and the direct interaction between adjacent inlet streams is weakened, which tends to reduce the local velocity gradients near the outlet. However, the longer flow path and the higher inlet airflow required to satisfy the maximum-temperature constraint introduce additional wall-friction and momentum-related losses. Therefore, the velocity and pressure fields are governed by the combined effects of flow redistribution, viscous friction, flow turning, and local geometric resistance, rather than by the Bernoulli relation alone.

4.4 Influence of heat flux

The highest local velocity occurs for d = 20 mm and q″ = 3600 W m⁻² and is approximately 41.7 m s⁻¹. Using the hydraulic diameter and air viscosity reported in Section 3.3 gives a maximum local Reynolds number of approximately 1.41 × 10³. Thus, the maximum local Reynolds number is about four times the inlet-based value but remains of order 10³. Because this maximum occurs near a flow-convergence and turning region, the result should be interpreted as a local flow-regime indicator rather than definitive proof of laminar behavior.

The influence of heat fluxes on pivotal properties was evaluated for all cases and results are shown in Figure 4a . Figure 4a shows that the active-area temperature difference generally increases with heat flux for every investigated fan spacing. The temperature difference is defined as

T a c t i v e = T max − T min

Because the maximum active-area temperature was constrained to approximately 358.15 K in all cases, the increase in Δ T active is not primarily caused by an increase in T max. Instead, a higher heat flux requires a larger cooling-air velocity, which strengthens the local convective cooling opposite the fan inlets and decreases the local minimum temperature.

Consequently, increasing the heat flux enlarges the spatial temperature difference even though the same maximum-temperature limit is satisfied. Under the present constant-T max design criterion, a larger Δ T active therefore indicates reduced temperature uniformity rather than failure to control the maximum temperature. This distinction is important because maximum-temperature control and active-area temperature uniformity represent two different thermal-management requirements.

Figure 4b shows that the pressure drop generally increases with heat flux at a fixed fan spacing. A greater heat load requires a larger inlet velocity and volumetric flow rate to maintain the prescribed maximum-temperature limit. The resulting increase in fluid momentum, wall shear, flow turning, and local geometric losses raises the pressure drop across the cooling passage.

The magnitude of the pressure drop also depends on fan spacing. In particular, the d = 20 mm configuration exhibits relatively high pressure drops over the investigated heat-flux range. This behaviour is associated with the restricted effective discharge region, the strong interaction between adjacent incoming streams, and the local flow turning near the sealing groove. The quantitative pressure drop was obtained from the difference between the averaged inlet and outlet pressures, rather than from the local maximum and minimum values displayed in the pressure contours.

Figure 4c shows that the inlet velocity required to satisfy the maximum-temperature constraint increases monotonically with heat flux for all fan spacings. As the applied heat flux increases, a larger convective heat-removal rate is required. This demand is met by increasing the airflow supplied through the fan openings.

Over the investigated range, the inlet-velocity curves appear approximately linear with heat flux. However, this trend is empirical and is limited to the present geometry, air properties, and thermal-load range. It should not be extrapolated to substantially higher heat fluxes or different stack configurations without further calculations.

Figure 4d shows that the maximum local air velocity also increases overall with heat flux. However, the maximum local velocity is not identical to the imposed inlet velocity. The inlet velocity is a boundary-condition value averaged over the fan opening, whereas the maximum local velocity occurs inside the cooling passage after the airflow has been redistributed by flow convergence, turning, local contraction, and the outlet geometry.

Consequently, the ordering of the maximum local velocities among different fan spacings does not necessarily coincide with the ordering of their inlet velocities. The maximum local velocity should therefore be interpreted as an indicator of local flow acceleration, while the inlet velocity and volumetric flow rate characterize the overall airflow requirement. These quantities should not be used interchangeably when evaluating cooling performance or airflow hydraulic power.

Overall, increasing the heat flux under the constant-T max constraint requires a higher inlet velocity and volumetric flow rate, which in turn increases the pressure drop and local flow velocity. Although the maximum active-area temperature remains controlled, the stronger localized cooling enlarges the active-area temperature difference. Therefore, increasing the cooling airflow alone does not simultaneously improve both maximum-temperature control and temperature uniformity.

4.5 Influence of fan interval

Figure 11a shows that the active-area temperature difference generally increases with fan spacing for all investigated heat fluxes. As the spacing increases, the distance between adjacent fan openings becomes larger, and the regions midway between the fans are located farther from the strongly cooled areas opposite the fan inlets. This produces a more spatially non-uniform cooling pattern along the active area.

Because the inlet velocity was increased separately for each spacing to maintain the same maximum-temperature limit, the stronger airflow mainly reduces the local temperature near the fan-facing regions. Meanwhile, the maximum temperature remains close to 358.15 K. Therefore, the increase in ΔTactive (Equation 13) is mainly caused by a decrease in the local minimum temperature rather than by an increase in the maximum temperature. A larger fan spacing can therefore satisfy the maximum-temperature constraint while producing poorer active-area temperature uniformity.

Δ T active = T max − T min (13)

Figure 11b shows a non-monotonic relationship between pressure drop and fan spacing. At small fan spacings, the effective discharge region between adjacent fan openings is relatively restricted. The incoming air streams interact strongly, and the flow undergoes pronounced turning around the sealing-groove geometry, resulting in relatively large local and overall hydraulic losses.

Increasing the fan spacing initially enlarges the available outlet region and weakens the direct interaction between adjacent inlet streams, thereby reducing the pressure drop. At larger spacings, however, the flow path becomes longer and a higher inlet velocity and volumetric flow rate are required to maintain the prescribed maximum-temperature limit. Wall friction, flow turning, and momentum-related losses then become increasingly important, causing the pressure drop to rise again. The observed minimum therefore results from the competition between the enlarged discharge area and the increased flow-path length and airflow requirement.

As shown in Figure 11c, the inlet velocity required to satisfy the maximum-temperature constraint generally increases with fan spacing. When the distance between adjacent fan openings increases, the cooling air must serve a larger repeated region and the point midway between the fans is located farther from the directly cooled boundaries. A higher airflow rate is therefore required to prevent this thermally disadvantaged region from exceeding the prescribed temperature limit.

The increase in inlet velocity also strengthens the local convection opposite the fan openings. Although this helps maintain the maximum temperature, it lowers the local minimum temperature and consequently contributes to the increase in the active-area temperature difference shown in Figure 11a.

Figure 11d shows that the maximum local air velocity varies non-monotonically with fan spacing. Unlike the imposed inlet velocity, the maximum local velocity is determined by the redistribution of the airflow inside the cooling passage. It is affected by flow convergence, turning around the sealing groove, local contraction, outlet area, and interaction between the streams supplied by adjacent fan openings.

Therefore, the fan spacing corresponding to the minimum maximum local velocity does not necessarily coincide with the spacing corresponding to the minimum pressure drop. The maximum local velocity characterizes local flow acceleration, whereas the pressure drop reflects the integrated hydraulic loss across the cooling passage. For fan selection and auxiliary-power assessment, the pressure drop must be evaluated together with the total volumetric flow rate through the hydraulic-power relation

P hyd = Δ p V ˙

Overall, increasing fan spacing produces both beneficial and adverse effects. A moderate increase in spacing enlarges the effective discharge region and can reduce the pressure drop. However, a larger spacing also increases the distance between the directly cooled regions, requires a higher inlet airflow, and worsens the active-area temperature uniformity. Fan spacing should therefore be selected by jointly considering the temperature difference, pressure drop, volumetric flow rate, and airflow hydraulic power rather than by using any single indicator.

4.6 Airflow hydraulic power and response-surface analysis

4.6.1 Airflow hydraulic-power calculation

The auxiliary power required for the edge-representative periodic unit depends on both the pressure rise supplied by the fan and the volumetric airflow rate. Therefore, the product of inlet velocity and pressure drop used in the original analysis does not directly represent fan power because the inlet area is not included.

In the revised analysis, the volumetric flow rate was obtained directly from the CFD inlet boundaries aswhere A in is the total effective inlet area within the representative computational unit. For a uniform velocity-inlet boundary, this expression can be simplified as

V ˙ = ∫ A in u · n dA

V ˙ = u in A in

The airflow hydraulic power required to overcome the pressure loss in the cooling passage was calculated aswhere Δ p is the difference between the averaged inlet and outlet pressures obtained directly from the CFD results.

P hyd = Δ p V ˙

If a practical fan with an efficiency of η fan is selected, the corresponding electrical input power can be estimated by

P elec = P hyd η fan

Because no manufacturer-specific pressure–flow–efficiency curve was prescribed in the present study, the following analysis reports the airflow hydraulic-power requirement rather than the actual electrical fan power. When the same fan efficiency is assumed for all configurations, the normalized airflow hydraulic-power trend is identical to the normalized electrical-power trend.

Figure 12 shows that the representative-unit airflow hydraulic-power requirement increases with heat flux because a greater thermal load requires a larger volumetric flow rate and generally produces a higher pressure drop. For a fixed heat flux, the airflow hydraulic power exhibits a non-monotonic dependence on fan spacing. Small spacings produce relatively high flow resistance because of the restricted discharge region and strong interaction between adjacent inlet streams. At large spacings, the longer flow path and higher required airflow increase both the volumetric-flow requirement and the associated hydraulic loss.

Among the five discrete fan spacings directly simulated, 40 mm gives the lowest or near-lowest airflow hydraulic power for the investigated heat-flux conditions. However, this discrete comparison alone cannot establish the exact location of the minimum between the simulated design points. A response-surface analysis was therefore performed using all 25 CFD-derived airflow hydraulic-power results.

4.6.2 Response-surface analysis

To characterize the continuous relationship among fan spacing, heat flux, and airflow hydraulic power, a second-order response-surface model was fitted to the 25 CFD-derived airflow hydraulic-power results. The model was expressed aswhere d is the fan spacing, q ″ is the applied heat flux, and β i are the regression coefficients determined using the least-squares method.

P hyd = β 0 + β 1 d + β 2 q ″ + β 11 d 2 + β 22 q ″ 2 + β 12 d q ″

The fitted model achieved a coefficient of determination R 2 of 0.9895, an adjusted R 2 of 0.9868, and a root-mean-square error of 0.0769. The predicted normalized airflow hydraulic-power values were compared with the original CFD-derived values to assess the fitting accuracy. The high coefficients of determination and the relatively low RMSE indicate that the quadratic response-surface model adequately reproduces the overall dependence of airflow hydraulic power on fan spacing and heat flux within the investigated design space.

Within 20 ≤ d ≤ 60 mm and 2,400 ≤ q″ ≤ 3,600 W m⁻², the fitted surface places the predicted minima between 37.7 and 40.0 mm (Table 2). Because the spacing input was sampled at 10 mm intervals and numerical/model uncertainty has finite magnitude, these decimal values are regression estimates rather than manufacturing-level precision. The engineering result is therefore reported as an approximately 38–40 mm candidate low-power region for one representative periodic unit.

TABLE 2. Regression-estimated minimum representative-unit airflow hydraulic-power fan spacings at different heat fluxes.

The response surface refines the trend obtained from the five directly simulated spacing levels, but it does not establish a universal or experimentally validated optimum. The approximately 38–40 mm candidate region depends on fan-opening width, passage height, cooling-edge width, plate properties, ambient conditions, and the imposed temperature constraint. Extrapolation beyond the investigated ranges is not recommended.

As shown in Figure 13a, the fitted surface reproduces the overall variation in the CFD-derived hydraulic power over the investigated design space. The hydraulic power increases substantially with heat flux, whereas its dependence on fan spacing is non-monotonic.

The contour map in Figure 13b shows that the regression-estimated minimum spacing varies from 37.7 to 40.0 mm as heat flux increases from 2,400 to 3,600 W m⁻². These values support the directly simulated 40 mm configuration but are reported as an approximately 38–40 mm candidate region because the underlying CFD spacing increment is 10 mm.

The fitted surface indicates a slight shift of the minimum toward larger spacing as heat flux increases. Nevertheless, 40 mm is the nearest directly simulated spacing at every heat-flux level. The result should therefore be interpreted as a candidate region rather than a single exact optimum.

From an engineering perspective, selecting a low-power region is more appropriate than specifying a single mathematical optimum because commercial fan dimensions, installation tolerances, airflow mismatch, and manufacturing constraints may prevent the exact fitted spacing from being implemented.

5 Conclusion

A three-dimensional steady conjugate heat-transfer model was developed to investigate the thermal and hydraulic performance of a periodic edge-cooling configuration for PEMFC stacks used in unmanned aerial vehicles. The effects of heat flux and fan spacing were systematically investigated under a prescribed maximum-temperature constraint.

The results show that increasing heat flux requires greater airflow to maintain the temperature limit, resulting in increased pressure drop, local air velocity, hydraulic-power demand, and temperature nonuniformity. The temperature distribution is governed by the combined effects of heat conduction through the graphite plate and spatially non-uniform forced convection along the cooling edge.

Fan spacing produces competing thermal and hydraulic effects. Small spacings restrict the discharge region and strengthen interactions between adjacent airflow streams, whereas large spacings increase the flow path and airflow requirement. Consequently, an intermediate fan-spacing range provides a favorable balance between cooling performance and hydraulic-power demand.

The identified favorable spacing region is specific to the present periodic geometry and operating conditions and should not be directly generalized to complete PEMFC stacks. Practical applications should further consider fan characteristics, stack-end effects, leakage, manufacturing tolerances, non-uniform heat generation, and actual UAV operating conditions. The present results should therefore be interpreted as preliminary comparative design guidance rather than validated full-stack predictions. Formal mesh-convergence assessment, independent experimental or numerical validation, and quantitative evaluation of the symmetry assumption were not available and should be addressed in future work.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Footnotes

Edited by: Jinliang Yuan, Ningbo University, China

Reviewed by: Vassilis Kostopoulos, University of Patras, Greece

Senthil S, KCG College of Technology, India