Artificial intelligence-assisted design and CFD re-optimization of graded hybrid metal foams for faster hydrogen desorption in LaNi 5 metal-hydride reactors

1 Introduction

Reversible metal hydrides are attractive for solid-state hydrogen storage because they combine high volumetric density, intrinsic safety, and operation at moderate pressure. Their practical discharge rate, however, is not governed only by the thermodynamic plateau of the alloy. It is controlled by a coupled transient process in which heat must reach the reactive bed while hydrogen leaves through a porous medium. During desorption, insufficient heat supply slows the reaction, leaves part of the bed underused, and reduces the useful response of the reactor. The recent study by (Chibani et al., 2026) provided a validated two-dimensional numerical analysis of L a N i 5 desorption enhanced with additively manufactured AlSi10Mg, CuCr1Zr, and 316L foams. Among the fixed alternatives examined there, CuCr1Zr at 89% porosity gave the strongest response. That result establishes an important reference point, but it also raises a more advanced design question: can the same reactor be improved further if the foam is no longer uniform and is instead treated as a radially graded, materially hybrid structure?

This question is consistent with the broader development of hydride-based hydrogen systems. The relevance of solid-state storage for mobile and stationary applications was emphasized by (Schlapbach and Zuttel, 2001), while (Zuttel, 2003) showed that candidate materials must be assessed through capacity, kinetics, reversibility, and operating conditions rather than by one property alone. The link between material development and reactor-level heat and mass transport was clarified by (Dornheim et al., 2007), and (Sakintuna et al., 2007) identified thermal management as a major obstacle for practical hydride devices. From a gas-reaction perspective (Sandrock, 1999), explained the continued interest in AB 5 alloys, and the review by (Joubert et al., 2021) confirmed the value of L a N i 5 because of its reversibility, fast kinetics, and convenient operating window.

Since thermal resistance remains the main rate-limiting mechanism, hydride-reactor research has gradually shifted from alloy choice toward engineered heat-transfer paths. Aluminium foam was shown by (Laurencelle and Goyette, 2007) to improve heat penetration in hydride beds, while the coupled simulations of (Muthukumar et al., 2009) demonstrated that desorption predictions must include both thermal inertia and gas flow. Thermal-management options were later organized by (Miao et al., 2023) into heat exchangers, conductive inserts, additives, and phase-change routes. Reactor geometry is also decisive, as discussed by (Cui et al., 2022), because the arrangement of the heat-exchange surface affects both charging and discharging. Passive buffering with phase-change materials has been explored by (Ye et al., 2021), who reported a smoother temperature field in a L a N i 5 reactor, by (Arslan et al., 2023), who experimentally observed improved storage behavior in a L a N i 5–H 2 system with PCM, and by (Cheng et al., 2025), who showed that matching the phase-change range to the operating window can enhance absorption and desorption. More recent hybrid concepts are equally relevant (Gao S. et al., 2025): demonstrated the advantage of combining fins with metal foam, and (Gao W. et al., 2025) used numerical optimization for hydride reactors, although within predefined enhancement layouts. Effective-property modelling also remains essential (Mou et al., 2025); highlighted the need for conductivity descriptions that remain valid during both hydrogen uptake and release.

Metal foams are particularly useful because they insert conductive skeletons into the reacting bed while retaining void space for hydrogen transport. The work of (Wang et al., 2021) reported accelerated heat exchange during hydrogen charging and discharging with foam assistance, and (Chibani et al., 2022) confirmed the benefit of foam-assisted latent-heat storage in large hydride reactors. Reduced-order approaches are also important for rapid design studies, as shown by the resistance–capacitance model of (Lamrani et al., 2025), which offers a computationally efficient alternative when repeated full CFD simulations are impractical. At the system level (Alobaid et al., 2025), stressed the need for reliable design criteria for hydride storage in stationary and mobile applications. On the manufacturing side, (Chanraksmey et al., 2026), (Atwater et al., 2018), and (Zhang et al., 2025) collectively show that tailored porous energy materials, high-porosity metallic structures, and laser-additively manufactured three-dimensional lattices are becoming realistic design options. The modelling basis for these architectures is supported by the effective-conductivity model of (Calmidi and Mahajan, 1999), the morphology-sensitive analysis of (Boomsma and Poulikakos, 2001), the thermophysical measurements of (Bhattacharya et al., 2002), the heat-transfer and permeability review of (Liu and Zhao, 2012), and the cellular-solids framework of (Gibson and Ashby, 1997). These studies make graded porous foams both physically meaningful and increasingly feasible.

The literature also shows that graded porosity itself cannot be claimed as a new concept. Bai et al. (2022) optimized a two-layer gradient-porosity metal-foam reactor for hydrogen absorption using a genetic algorithm and reported an 11.5% improvement over a uniform-foam reactor. More recently, Ye et al. (2026) examined radially graded-porosity metal foam in a L a N i 5 absorption reactor and reported about a 20% reduction in the time required to reach 90% absorption relative to a uniform-porosity foam. Zhu et al. (2026) also reported a graded-metal-foam benefit in a metal-hydride/PCM thermal-management configuration. Accordingly, the novelty claimed here is deliberately narrower: the present work concerns desorption and couples a smooth radial porosity profile with a spatial CuCr1Zr/AlSi10 Mg allocation and a transient heat-transfer/pressure schedule in one self-adaptive evolutionary search. It does not claim to introduce graded metal foam itself.

The distinction from the closest graded-foam studies is also process-specific. Absorption is exothermic and is principally limited by removal of reaction heat, whereas desorption is endothermic and requires sustained heat input while the hydrogen pressure is kept sufficiently below equilibrium. The present design problem therefore couples where conductive material is placed with when heat-transfer intensity and outlet pressure are applied. This simultaneous radial material allocation and transient h / P scheduling is the central co-design contribution relative to the absorption-focused graded-porosity studies of (Bai et al., 2022), (Ye et al., 2026), and (Zhu et al., 2026).

The remaining challenge is how to choose such a structure. Comparing a limited number of uniform foams is useful for building a baseline, but it cannot explore the continuous design space created by radial porosity, local alloy fraction, and transient operating conditions. Differential evolution is well suited to this type of nonlinear search. The algorithm introduced by (Storn and Price, 1997) provides a robust global optimization strategy for continuous variables, while (Brest et al., 2006) showed that self-adaptive control parameters can improve reliability by avoiding fixed mutation and crossover choices. Further strategy adaptation was introduced by (Qin et al., 2009), and the review by (Das et al., 2016) confirmed the strength of differential evolution for multimodal engineering problems. Because high-fidelity reactor calculations are expensive, surrogate-assisted optimization is also necessary. The foundations presented by (Queipo et al., 2005; Forrester et al., 2008) support calibrated reduced models for design exploration. In the present workflow, reduced-order screening is not the endpoint: the same decision variables are subsequently re-evaluated and re-optimized with an OpenFOAM finite-volume CFD model.

A clear gap therefore remains. The validated foam-enhanced L a N i 5 study of (Chibani et al., 2026) provides strong evidence for fixed additively manufactured foams, but it does not examine a continuous space in which porosity is graded, alloy composition changes with radius, and early heating and pressure are scheduled together. The present work addresses this missing design level. It replaces the question of which uniform foam is best with a more flexible question: where should conductive material be placed, how much permeability should be preserved, and how should the first minutes of operation be controlled to exceed the best fixed benchmark? In this way, the reactor is treated as an optimized thermal–mass-transfer architecture rather than as a passive insert selected by comparison.

The study uses a three-level evidence hierarchy: (i) the original MATLAB surrogate for inexpensive screening, (ii) a MATLAB finite-volume PDE challenge of the exact surrogate-selected candidate, and (iii) OpenFOAM v2506 re-optimization and fine-grid/full-3D verification. The 7.95-min value is never used as an OpenFOAM calibration, stopping, or optimization target.

The remainder of the paper is arranged as follows. Section 2.1 describes the chemical, physical, and mathematical model used for the hydride reactor. Section 2.2 presents the proposed PG-SADE-GHF optimization procedure and its algorithmic structure. Section 3 provides a figure-by-figure physical interpretation, consolidated numerical evidence, and comparisons with related thermal-management studies. Section 4 discusses the scope of the prediction and the required validation route. Finally, Section 5 summarizes the main conclusions for optimized metal-hydride reactor design. All section, figure, table, equation, citation, URL, and DOI references in the electronic manuscript are active hyperlinks.

2 Materials and methods

2.1 Physical and mathematical model

2.1.1 Hydrogen desorption chemistry and thermodynamics

The reversible hydride reaction is written as shown in Equation 1:where M denotes the storage alloy, M H x the hydride phase, x the stoichiometric coefficient, and Δ H the reaction enthalpy. In desorption mode, external heat is required to release hydrogen from the solid phase. The equilibrium pressure is evaluated with a Van’t Hoff relation,where A and B represent the entropy- and enthalpy-related contributions. Hydrogen density in the gas phase is calculated from the ideal-gas expression, as given in Equation 3:

M + x 2 H 2 ⇌ M H x ± Δ H , (1)

ln p e q = A − B T , (2)

ρ g = p g M g R T . (3)

2.1.2 Mass conservation and Darcy flow

The gas-phase balance in the porous hydride bed is expressed aswhere ε is the bed porosity, ρ g is the hydrogen gas density, u is the superficial velocity, and m ̇ is the desorption source term. The solid-phase mass balance is given by Equation 5:

ε ∂ ρ g ∂ t + ∇ ⋅ ρ g u = − m ̇ , (4)

1 − ε ∂ ρ s ∂ t = m ̇ . (5)

Hydrogen flow through the bed is approximated using Darcy’s law,where K is permeability and μ g is the dynamic viscosity of hydrogen.

u = − K μ g ∇ p g , (6)

2.1.2.1 Scope of the Darcy approximation at high foam porosity

Equation 6 is retained as the low-velocity, bed-scale approximation used by the reference hydride model; it should not be interpreted as proof that all flow through an open-cell metal foam is strictly Darcy flow. For the foam skeleton, a Darcy–Forchheimer form provides the more general screening relationwhere K f and c F are the foam permeability and Forchheimer coefficient. Yang et al. (2021) reviewed experimental pressure-drop data and proposed correlations applicable over approximately 85.7–97.8 % porosity and 10–100 PPI. One recommended permeability form is given by Equation 8:with d f / d p = 2 ( 1 − ε m f ) / ( 3 π ) / [ 1 − exp ( − ( 1 − ε m f ) / 0.04 ) ]. The corresponding inertial coefficient may be screened with c F = 0.00212 ( 1 − ε m f ) − 0.132 ( d f / d p ) − 1.63. Recent high-pore-density measurements further confirm that the Darcy-to-Forchheimer transition depends on pore density and velocity, not on porosity alone (Lv et al., 2024). Because the reduced optimizer does not resolve local pore diameter/PPI or the velocity field of the new graded lattice, these correlations are used as a regime and sensitivity check rather than as a substitute for pore-resolved CFD. The final graded geometry must therefore be re-evaluated with Equation 7 once its manufactured pore scale is fixed.

− ∇ p g = μ g K f u + ρ g c F K f | u | u , (7)

K f d p 2 = 0.00747 1 − ε m f − 0.2329 d f d p 0.00249 , (8)

2.1.3 Reaction kinetics

The desorption source term follows an Arrhenius-type kinetic law driven by the pressure difference between the gas phase and the equilibrium state,where C d is the kinetic constant, E d is the activation energy, and ρ emp is the empty-bed solid density. The normalized hydrogen fraction is defined as shown in Equation 10:

m ̇ = C d p g − p e q p e q exp − E d R T ρ s − ρ emp , (9)

X = ρ s − ρ 0 ρ sat − ρ 0 . (10)

In the reduced reproduction model, desorption is considered complete when X ≤ X crit, with X crit = 0.01 used to maintain numerical robustness.

2.1.3.1 Parameter provenance and thermodynamic consistency

The critical physical constants were not re-identified by the evolutionary optimizer. They were inherited from the validated reference formulation and cross-checked against independent L a N i 5 literature before being used as fixed inputs. In particular, the adopted C d = 9.57 s − 1 and E d = 16.42 kJ mol − 1 coincide with the L a N i 5 values summarized by (Gambini et al., 2024); that source reports a desorption enthalpy magnitude of 30.1 kJ mol − 1 and entropy change of 109.96 J mol − 1 K − 1, close to the 30.8 kJ mol − 1 reaction-enthalpy magnitude used in the present reference model. With the sign convention of Equation 2, A = Δ S / R and B = Δ H / R give A ≈ 13.23 and B ≈ 3620 K for those independent values, yielding an equilibrium pressure of approximately 2.95 bar at 298.15 K under the same pressure-reference convention. This cross-check verifies parameter order and thermodynamic consistency; it is not presented as a new kinetic fit.

The MATLAB reproduction remains a calibrated response model: its fitted time constants reproduce the published CFD trends, whereas the constants above provide the governing physical basis and provenance. Consequently, agreement of these constants with literature should not be confused with an independent validation of the optimized 7.95 min prediction.

2.1.4 Energy equation without and with metal foam

For a hydride bed without foam, the volume-averaged energy balance is given by Equation 11:where the effective heat capacity is defined by Equation 12 and the effective conductivity follows Equation 13:

ρ c eff ∂ T ∂ t + ρ c p g u ⋅ ∇ T = λ eff ∇ 2 T + m ̇ M g Δ H , (11)

ρ c eff = ε ρ c p g + 1 − ε ρ c s , (12)

λ eff = ε λ g + 1 − ε λ s . (13)

When a metal foam skeleton is inserted into the bed, the energy balance becomeswith effective heat capacity and conductivity written as the foam-filled heat capacity is given by Equation 15:

ρ c eff,mf ∂ T ∂ t + ρ c p g u ⋅ ∇ T = λ eff,mf ∇ 2 T + ε m f m ̇ M g Δ H , (14)

ρ c eff,mf = ε m f ρ c eff + 1 − ε m f ρ c m f , (15)

λ eff,mf = ε m f λ eff + 1 − ε m f λ m f . (16)

The external thermal boundary condition is imposed as shown in Equation 17:where h is the convective heat-transfer coefficient and T f is the temperature of the heat-transfer fluid.

− λ ∂ T ∂ n = h T f − T , (17)

2.1.5 Reduced-order calibrated reproduction

The MATLAB implementation follows the same physical assumptions, but it does not use the original Fluent mesh or UDF files. The average temperature response is represented byand the average hydrogen fraction is approximated aswhere τ X is calibrated so that X ̄ ( t d ) = X crit. Spatial contour maps are then reconstructed with smooth core-lag and boundary-heating functions so that the field evolution follows the same ordering and trends as the published Fluent results. The reduced model is used only for rapid optimization screening; its selected design is subsequently challenged and re-optimized with OpenFOAM CFD before high-fidelity performance claims are made.

T ̄ t = T f − T f − T 0 exp − t τ T p T , (18)

X ̄ t = exp − t τ X q X , (19)

A further distinction is essential for interpreting the figures. The average response functions in Equations 18 and 19 are calibrated reduced-order responses, and the two-dimensional temperature and hydrogen maps are smooth reconstructions constrained to reproduce the ordering and qualitative penetration patterns of the published Fluent cases. They are not obtained by re-solving Equations 4, 7 and 14 on a new graded mesh. Accordingly, the Results section uses the terms “surrogate-reconstructed”, “illustrates”, and “is consistent with” for these maps and does not use them as independent CFD evidence for the optimized architecture.

2.2 Developed PG-SADE-GHF optimization method

2.2.1 Graded hybrid foam parameterization

The proposed framework represents the foam as a radially graded hybrid structure. The normalized coordinate r ∈ [ 0,1 ] is measured from the heated wall toward the reactor core. The metal-foam porosity is described bywhere ε w is the wall-region porosity, ε c is the core-region porosity, r w defines the conductive wall-zone fraction, and δ r controls the transition smoothness. The smoothing function is given by Equation 21:

ε m f r = ε w + ε c − ε w S r ; r w , δ r , (20)

S r ; r w , δ r = 1 2 1 + tanh r − r w δ r . (21)

The local CuCr1Zr fraction follows the same radial transition,

f C u r = f C u , w + f C u , c − f C u , w S r ; r w , δ r . (22)

The resulting local foam conductivity is calculated from

λ m f r = f C u r λ CuCr1Zr + 1 − f C u r λ AlSi10Mg . (23)

This formulation allows the optimizer to create a compact, highly conductive CuCr1Zr-rich region near the wall while keeping the core lighter and more permeable through an AlSi10 Mg-rich composition.

Equation 23 is an interpolation of the ligament/base-material conductivity, not a claim that the effective conductivity of a real mixed open-cell foam is exactly linear in alloy fraction. Porosity enters separately through Equation 16. To quantify the uncertainty introduced by the mixing assumption, the supplementary check brackets the optimized wall composition (f C u , w = 0.65, ε w = 0.861) between arithmetic and harmonic ligament mixing. The arithmetic rule gives λ m f = 274.25 W m − 1 K − 1 and λ eff,mf ≈ 39.26 W m − 1 K − 1, whereas the harmonic rule gives 253.52 and 36.38 W m − 1 K − 1, respectively; the resulting effective-conductivity spread is approximately 7.3%. This bracket is now treated as modelling uncertainty. Contact resistance, anisotropy, strut geometry, and alloy-interface effects require the final CFD/experimental characterization and are not hidden inside the linear rule. The transport closures were also audited over the systematic 44-design CFD dataset. In the actual OpenFOAM property map, replacing arithmetic by harmonic ligament mixing while holding each sampled design coordinate fixed lowers the local effective conductivity by 0.04%–5.68% at the wall and 0.16%–5.24% in the core for the 32-development set; the corresponding held-out ranges are 0.81%–5.54% and 0.18%–5.00%. This complements the earlier 7.3% single-composition analytical bracket and shows that the mixing-rule uncertainty is material but smaller than the original surrogate–CFD completion-time discrepancy. The reduced screening surrogate contains no resolved local Forchheimer velocity correction, whereas OpenFOAM retains the inertial term explicitly. Across the 32 development cases, the maximum Forchheimer/Darcy ratio spans 0.066–0.498; relative to a Darcy-only mobility, the corresponding local Darcy–Forchheimer correction reduces mobility by as much as 33.25%. The held-out set spans 0.086–0.578, corresponding to a maximum local mobility reduction of 36.62%. The uncorrected CFD-minus-surrogate completion-time error is positively associated with the maximum inertial ratio (Spearman ρ = 0.542, p = 0.0014). Thus, omitted inertial resistance is a measurable contributor to the low-fidelity discrepancy, although it is not its sole cause. The closure audit is reported as a momentum/property sensitivity and is not presented as a second full re-optimization under altered constitutive laws.

2.2.2 Time-dependent heating and pressure schedule

In addition to the foam architecture, the optimizer adjusts an early-stage operating schedule. The heat-transfer coefficient and desorption pressure are defined as shown for h(t) in Equation 24:the pressure schedule is then given by Equation 25:with the smooth switching function defined in Equation 26:

h t = h late + h early − h late σ t ; t s , Δ t , (24)

P d t = P late + P early − P late σ t ; t s , Δ t , (25)

σ t ; t s , Δ t = 1 1 + exp t − t s / Δ t . (26)

Thus, the initial period can use stronger heating and a more favorable discharge pressure, after which the conditions gradually return to milder operation.

2.2.3 Decision vector and objective function

The design variables are collected as shown in Equation 27:

x = ε w , ε c , f C u , w , f C u , c , h early , P early , t s , r w . (27)

The objective function iswhere t d is the complete desorption time. The terms Π m, Π u, Π o, and Π f penalize excessive metal loading, thermal non-uniformity, overly severe operating conditions, and manufacturability limitations, respectively. This penalized objective prevents a trivial solution based on maximum copper content and maximum heating throughout the bed.

J x = t d x + w m Π m x + w u Π u x + w o Π o x + w f Π f x , (28)

2.2.3.1 Rationale for bounds and interpretation of active constraints

Several surrogate coordinates were located at search limits; therefore, the range analysis was replaced by an actual two-stage OpenFOAM bound study. The primary CFD search uses exactly the original MATLAB bounds: ε w = 0.86–0.90, ε c = 0.91–0.96, f C u , w = 0.65–1.00, f C u , c = 0–0.50, h early = 1700–2400 W m − 2 K − 1, P early = 0.35–0.75 bar, t s = 2–6 min, and r w = 0.20–0.45. The best fine-mesh OpenFOAM design in this exact envelope is (0.86004, 0.91753, 1.000, 0.500, 2400, 0.36168, 5.6589, 0.44451) with t d = 7.8167 min. Because ε w, f C u , w, f C u , c, and h early remain active, a secondary expanded-envelope sensitivity was then performed using ε w = 0.84–0.92, ε c = 0.89–0.98, f C u , w = 0.55–1.00, f C u , c = 0–0.60, h early = 1500–2800 W m − 2 K − 1, P early = 0.25–0.90 bar, t s = 1.5–7 min, and r w = 0.15–0.55. This gives a separate sensitivity optimum t d = 6.4417 min. We therefore use 7.8167 min as the primary same-problem high-fidelity counterpart to the MATLAB optimization and 6.4417 min only as an expanded-bound extension, avoiding retroactive redefinition of the original problem. Re-scoring the already evaluated OpenFOAM candidates with penalty multipliers 0.5, 1.0, and 2.0 leaves both selected completion times unchanged because the best candidates carry zero penalty.

2.2.4 Self-adaptive differential evolution operators

At each generation g, a current-to-pbest/1 mutation strategy is used to generate the mutation vector,

v i g = x i g + F i x pbest g − x i g + F i x r 1 g − x r 2 g . (29)

The trial vector is produced by binomial crossover,and elitist selection keeps the better candidate,

u i , j g = v i , j g , r a n d j ≤ C R i o r j = j rand , x i , j g , o t h e r w i s e . (30)

x i g + 1 = u i g , J u i g ≤ J x i g , x i g , o t h e r w i s e . (31)

The memory terms F m and C R m are updated from successful trials according to

F m g + 1 = 0.85 F m g + 0.15 ∑ F s 2 ∑ F s , (32)

C R m g + 1 = 0.85 C R m g + 0.15 C R ̄ s . (33)

2.2.4.1 Controlled fixed-DE comparison

The original wording was too broad in stating that self-adaptation “outperforms” fixed-parameter DE without an application-specific control. We now distinguish algorithm selection from proof of physical accuracy. In addition to the literature motivation in (Brest et al., 2006; Qin et al., 2009), a deterministic equal-budget stress test was performed on the same eight-variable bounded decision structure using equal budgets (N P = 28, 30 generations, 30 paired seeds). Fixed DE used F = 0.50 and C R = 0.90, while the adaptive variant used the successful-parameter memories in Equations 32, 33. On this diagnostic normalized objective, self-adaptive DE obtained a lower final value in 29 of 30 paired runs; the median final objective was 0.1063 for SADE and 0.1473 for fixed DE, with standard deviations 0.0073 and 0.0263, respectively. These data support the use of adaptation for the search procedure but are not used to claim that SADE itself validates the reactor physics. The deterministic test and its full run table are supplied in MATLAB_V8_ROBUSTNESS_EVIDENCE.m.

The full PG-SADE-GHF pseudocode has been moved to Supplementary Algorithm S1. The main text retains the governing equations, decision variables, adaptive-memory updates, and equal-budget fixed-DE comparison needed to reproduce the method without duplicating procedural detail.

2.3 OpenFOAM v2506 porous-continuum CFD re-optimization and verification

An independent high-fidelity implementation was developed in OpenFOAM v2506, a finite-volume CFD framework based on the FOAM methodology (Weller et al., 1998). An active ANSYS Fluent licence was not available for the present revision, and the original proprietary Fluent project files, mesh, and UDFs associated with the reference calculations were also unavailable for direct reuse. Rather than reconstructing a partially specified proprietary case, the governing equations were implemented directly from their published mathematical definition in OpenFOAM. This choice is not presented as a comparison or criticism of commercial CFD software. The decisive methodological advantage for the revision is reproducibility: OpenFOAM is open source, and the solver source, dictionaries, boundary conditions, numerical settings, execution logic, and post-processing route can be linked directly to editable text-based scripts and archived with the study. The resulting workflow can therefore be inspected, rerun, modified, and extended without relying on hidden project-state information. The LaNi 5 kinetic and thermodynamic constants remain unchanged.

Computational implementation and numerical settings. Calculations were executed on a Windows workstation through Docker Desktop using the fixed container image opencfd/openfoam-dev:2506. Meshes were generated with OpenFOAM blockMesh; the custom solvers pgsadeFoam and pgsadeFoam3D were compiled inside the same v2506 environment and executed in serial mode. ParaView was used for post-processing after VTK/PVD export, with a common clipped cutaway and fixed camera for the full-3D field panels. The software-side route from Windows PowerShell to Docker/OpenFOAM execution and then to ParaView post-processing is documented in Figure 10; the corresponding source files, dictionaries, and execution wrappers are included in the supplementary computational package.

For clarity, the software route used in the revision can be summarized procedurally as follows.

The stepwise PowerShell–Docker–OpenFOAM–ParaView execution pseudocode has been moved to Supplementary Algorithm S1. The main text retains the solver, mesh, discretization, stopping criterion, and post-processing settings that affect the numerical solution.

The CFD-in-the-loop search uses a 5° three-dimensional wedge corresponding to an axisymmetric cylinder of radius 0.06 m and height 0.35 m. The primary wedge mesh is 24 × 1 × 48 with Δ t = 0.5 s; independence checks include 12 × 1 × 24 and 36 × 1 × 72 meshes with time steps from 1.0 to 0.25 s. After optimization, the selected original-bound and expanded-bound configurations are re-solved on a circumferentially complete 360° cylinder using pgsadeFoam3D. The standard full-3D mesh is a structured five-block O-grid-like arrangement with 82,944 hexahedral cells: a 24 × 24 × 48 central block and four 24 × 12 × 48 curved outer blocks. The high-resolution confirmation uses 196,608 hexahedral cells: 32 × 32 × 64 in the central block and 32 × 16 × 64 in each outer block, with Δ t = 0.25 s. OpenFOAM checkMesh reports one connected region, maximum non-orthogonality 38.67° (mean 6.94°), maximum skewness 0.693, maximum aspect ratio 7.67, and “Mesh OK.” The outer cylindrical surface is the heated wall, the top is the pressure outlet, and the bottom is thermally adiabatic and gas-impermeable. The heat-transfer-fluid temperature is 363 K and the initial bed temperature is 290 K.

The temporal derivative is discretized with first-order Euler integration. Spatial gradients use Gauss linear reconstruction; Laplacian terms use Gauss linear corrected discretization; interpolation is linear and the surface-normal gradient uses the corrected scheme. The pressure and temperature linear systems are solved with PCG and DIC preconditioning, with absolute tolerance 1 0 − 10 and relative tolerance 0.01. The Darcy–Forchheimer pressure coupling is updated with three Picard iterations per time step. These settings are identical between the standard and high-resolution full-3D confirmation cases except for mesh density and time step.

The OpenFOAM solver, pgsadeFoam, advances local temperature T, remaining hydride fraction X, hydrogen pressure p, and Darcy–Forchheimer velocity. The equilibrium pressure and local desorption coefficient are given by Equation 34:with C d = 9.57 s − 1, E d = 16.42 kJ mol − 1, | Δ H | = 30.8 kJ mol − 1, and Δ S = 109.96 J mol − 1 K − 1. The hydride fraction is advanced positivity-preservingly as X n + 1 = X n ⁡ exp ( − k d Δ t ). The pressure equation uses the mobility defined in Equation 35:which is iterated three times per time step and therefore retains the Forchheimer correction rather than assuming pure Darcy flow a priori. The energy equation is given by the transient balance in Equation 36:

p e q = exp Δ S R − Δ H R T , k d = C d ⁡ exp − E d R T max p e q − p p e q , 0 , (34)

M = K / μ g 1 + ρ g c F K | u | / μ g , − ∇ ⋅ M ∇ p = m ̇ ρ g , (35)

ρ c p eff ∂ T ∂ t − ∇ ⋅ k eff ∇ T = − m ̇ Δ H M H 2 − ρ g c p , g u ⋅ ∇ T . (36)

At the heated wall, − k eff ∇ T ⋅ n = h ( t ) [ T − T f ] is imposed with the same logistic h ( t ) schedule as the MATLAB design; the outlet pressure is set to the corresponding P d ( t ) schedule.

The radial grading law is identical to Equations 20, 22. The bed conductivity scale k bed,scale = 5.501 and foam-network efficiency η foam = 0.579 are frozen from the two-anchor thermal closure calibration and are never adjusted using 7.95 min. Permeability varies with porosity according to a normalized Kozeny–Carman form, K ∝ ε 3 / ( 1 − ε ) 2, normalized at ε = 0.89 and bounded to 0.05 K bed–50 K bed. The completion criterion is the capacity-weighted volume average X ̄ ≤ 0.01.

For CFD-in-the-loop optimization, each candidate writes the eight design variables to constant/designParameters; the compiled solver is then run without recompilation. The self-adaptive DE outer loop uses a population of 10 and eight generations (up to 90 coarse OpenFOAM evaluations including the initial population) with τ 1 = τ 2 = 0.10 and F ∈ [ 0.1 , 1.0 ]. The five best candidates are independently re-ranked on the 24 × 48 fine mesh. The expanded-bound study repeats the same protocol with a different random seed. A soft penalty is activated only for excessive pressure or Forchheimer contribution; the final best candidates have zero penalty.

The supplementary computational package contains the custom pgsadeFoam and pgsadeFoam3D source files, blockMesh dictionaries, physical-property and boundary-condition files, the Windows PowerShell to Docker execution wrappers, the fixed image identifier opencfd/openfoam-dev:2506, VTK/ParaView export instructions, mesh-quality logs, and machine-readable result summaries. The package therefore preserves both the mathematical model and the executable numerical route. Because the numerical configuration is exposed through source files and scripts, the same route can be used in subsequent work to test alternative graded-foam profiles, materials, operating schedules, kinetic closures, mesh densities, or experimentally calibrated parameters without rebuilding the workflow from an opaque graphical project.

2.4 Systematic surrogate–CFD discrepancy correction and locked challenge

Surrogate accuracy was assessed over the original eight-variable design envelope rather than at the optimized point alone. A 32-case OpenFOAM development set was assembled from the original space-filling designs and the dense follow-on sample. The reported screening surrogate and its continuous pre-floor response were evaluated at the same coordinates. Because the reported surrogate contains a 7.95-min lower safeguard, 29 of 32 development predictions (90.6%) collapse to that value; this saturation is treated explicitly as a limitation rather than as CFD agreement.

A discrepancy model was trained on the difference between the reduced prediction and the OpenFOAM completion time. Model selection used repeated cross-validation on the 32 development cases, and the locked predictor is a weighted ensemble of the three best kernel-ridge discrepancy models. Development-set nested cross-validation gives R 2 = 0.955, RMSE = 0.271 min, MAE = 0.198 min, and Spearman ρ = 0.968. A second set of 12 Sobol/maximin designs was then generated inside the same original bounds. Their corrected predictions and 80%/95% prediction intervals were written and SHA-256 locked before the H01–H12 OpenFOAM outcomes were computed. The held-out designs, locked predictions, lock receipt, and CFD result tables are supplied as machine-readable supplementary data. The supplementary data also include casewise arithmetic–harmonic effective-conductivity diagnostics and the Darcy–Forchheimer mobility-reduction audit used for the closure-sensitivity rows of Table 1. This correction quantifies and reduces the low-/high-fidelity discrepancy; the final optimized reactor values continue to come from direct OpenFOAM re-optimization.

TABLE 1. Systematic surrogate–CFD validation, transport-closure sensitivity, and held-out uncertainty. Panel (a) reports development and locked-.challenge predictive metrics; Panel (b) reports arithmetic–harmonic conductivity and Darcy–Forchheimer diagnostics; Panel (c) reports bootstrap and interval-coverage uncertainty.

3 Results

3.1 Benchmark reproduction, exact-candidate challenge, and OpenFOAM re-optimization

The analysis separates design-point validation from optimization-problem validation. The original MATLAB surrogate predicted t d = 7.95 min at

ε w , ε c , f C u , w , f C u , c = 0.861 , 0.946 , 0.65 , 0 , h early , P early , t s , r w = 1700,0.75 , 2.27 , 0.30 .

Direct OpenFOAM evaluation of those exact coordinates gives t d = 12.1917 min; therefore the exact surrogate coordinates are not presented as validated. This negative result is retained transparently because it demonstrates that the CFD model was not tuned to recover the desired value.

The crucial additional test is then performed: the same eight variables are re-optimized with actual OpenFOAM evaluations inside the exact original MATLAB bounds. The uniform CuCr1Zr-89% reference is independently reproduced at 9.975 min. The fine-mesh original-bound OpenFOAM optimum is t d = 7.8167 min, with X ̄ = 0.009936, p max = 0.9993 bar and a maximum Forchheimer/Darcy ratio of 0.1587. This value is 1.677% below the original 7.95-min MATLAB magnitude and 21.833% below the 10-min benchmark. The top five independently re-ranked fine-mesh candidates span only 7.8167–7.9583 min, showing that the result is not an isolated single point.

Because several original-bound coordinates remain active, an expanded-bound study is reported separately as sensitivity rather than as retroactive validation. The expanded optimum is t d = 6.4417 min with X ̄ = 0.009997, p max = 0.9353 bar and maximum Forchheimer/Darcy ratio 0.0955. It is 35.583% faster than the 10-min reference and 18.973% faster than 7.95 min. The primary high-fidelity CFD counterpart to the original MATLAB problem remains 7.8167 min; the 6.4417-min value is an extension enabled by the wider engineering envelope.

3.2 Reference reproduction and material-response comparison

To conserve main-text artwork, the numerical-independence evidence is now reported directly in the text rather than as a standalone figure. For the expanded-bound OpenFOAM optimum, the wedge sequence gives 6.4500 min on 12 × 24 with Δ t = 1.0 s, 6.4500 min on 24 × 48 with Δ t = 1.0 s, and 6.4417 min on the 24 × 48/Δ t = 0.5 s and 36 × 72/Δ t = 0.25 s checks; the maximum deviation from the finest result is 0.129%. The subsequent full-360° verification is even more restrictive: the expanded-bound full-3D calculation and its 196,608-cell high-resolution repeat both give 6.4417 min. Thus, the performance differences discussed below are not attributable to the chosen spatial or temporal resolution.

The comparison in Figure 1 addresses whether the reduced response has the characteristic time scale and monotonic behavior of the aluminium-foam hydride benchmark before optimization; it does not validate the optimized graded reactor. The reconstructed curve follows the reference numerical/experimental trend used for calibration. This behavior is consistent with Laurencelle and Goyette (Laurencelle and Goyette, 2007), who showed that a conductive porous skeleton shortens the thermal path inside a hydride bed, and with the coupled heat-and-mass-transfer analysis of (Muthukumar et al., 2009). We therefore use this figure as a trend/calibration check. High-fidelity optimization and verification of the final design space is provided separately by the OpenFOAM results in Figures 8, 9, 11.

Figure 2 is retained as a surrogate-screening comparison. The plotted responses are the bed-averaged temperature and hydrogen-fraction histories for the fixed 316L, AlSi10Mg, and CuCr1Zr foams together with the PG-SADE-GHF candidate. It preserves the reproduced fixed-material ordering associated with the adopted bulk conductivities (approximately 19, 180, and 325 W m − 1 K − 1 for 316L, AlSi10Mg, and CuCr1Zr) and shows why the reduced optimizer favored spatial redistribution of conductivity. The 7.95 min graded trajectory in this figure must now be read strictly as the surrogate-selected trajectory. The independent frozen PDE reassessment summarized in the lower two rows of Figure 4 and in Figure 7 challenges those exact coordinates; the subsequent same-bound OpenFOAM re-optimization in Figure 8 recovers t d = 7.817 min within the original design envelope.

The porosity trends reinforce this interpretation. Reducing porosity from 93% to 89% lowers the completion time by 33.3% for 316L, 33.3% for AlSi10Mg, and 37.5% for CuCr1Zr. A larger solid fraction therefore improves conduction, but a uniform reduction in porosity also increases metal mass and can restrict the connected gas volume. This conductivity–permeability conflict is well established in metal-foam studies (Calmidi and Mahajan, 1999; Boomsma and Poulikakos, 2001; Bhattacharya et al., 2002; Liu and Zhao, 2012) and in the cellular-solids framework of (Gibson and Ashby, 1997). The proposed wall–core grading resolves the conflict more effectively: it accepts a dense wall layer where conduction is valuable and preserves high core porosity where transport and thermal mass become more important. Table 2 reports the complete fixed-case matrix used to produce these comparisons, rather than only the best value from each material family.

TABLE 2. Combined numerical basis: (a) traceability from the original MATLAB screening result to the OpenFOAM v2506 high-fidelity calculations; (b) complete fixed-foam surrogate-screening matrix and within-material porosity sensitivity. This panelled table retains the full contents of the former Table 3 while reducing the number of main-text floats.

3.3 Thermal and hydrogen-field behavior

Average curves alone do not resolve the local wall-to-core progression. In Figure 3, panel (a) represents reconstructed temperature evolution and panel (b) reconstructed hydrogen-fraction evolution for the fixed cases and the optimized reduced-model trajectory. Both panels are therefore retained only as surrogate-reconstructed visualizations of the calibrated wall-to-core thermal and depletion patterns, rather than as local CFD solutions of the graded geometry. Within that reconstruction, the foam-free and weakly conducting reference cases retain a larger cold-core region, whereas the higher-conductivity CuCr1Zr reference exhibits faster apparent penetration. The optimized reconstructed field is consistent with the intended wall-focused conductive architecture. However, these panels are not obtained by solving the local energy equation on a new graded CFD mesh and therefore do not verify the local temperature distribution of the optimized reactor. Their role is explanatory: they illustrate the mechanism assumed by the reduced model and identify a hypothesis that must be checked in the requested high-fidelity graded-geometry CFD calculation.

The distinction from other thermal-management routes is also important. PCM-assisted systems primarily buffer or recycle reaction heat (Ye et al., 2021; Arslan et al., 2023; Cheng et al., 2025); fin–foam systems increase external and internal exchange area (Gao S. et al., 2025); and improved heat-exchanger arrangements shorten the distance between the heat-transfer fluid and the hydride (Cui et al., 2022; Miao et al., 2023). The present approach is complementary but different: it redistributes the internal solid conduction pathway itself. In the reduced screening model, that spatial redistribution is associated with a surrogate-predicted 20.5% improvement beyond a fully uniform CuCr1Zr-89% reference. The independent frozen PDE reassessment does not reproduce this gain; it predicts 12.233 min for the same candidate. The design hypothesis therefore remains interesting, but the quantitative superiority claim is withdrawn pending PDE-in-the-loop re-optimization and subsequent CFD/experimental confirmation.

The operating sweeps in the upper two rows of Figure 4 remain useful for interpreting the surrogate objective. Panels (a)–(d) combine the heat-transfer-coefficient and desorption-pressure sweeps with the candidate-response trends and the adaptive heating/low-pressure schedule selected during surrogate screening. Increasing constant h shows diminishing returns, whereas lower fixed pressure increases the modeled desorption driving force. However, the previous statement that the adaptive graded design “wins” over these fixed cases has been removed. The 7.95 min value shown in this screening composite belongs to the reduced model; the independent frozen PDE evaluation of the same candidate is 12.233 min. The operating sweep therefore motivates the coupled design variables but does not establish a verified performance advantage.

This interpretation is consistent with the multi-step optimization direction reported by (Gao W. et al., 2025), but the present framework extends the idea by optimizing material grading and operating schedule simultaneously. It also responds to the effective-property concern raised by (Mou et al., 2025): the benefit of a foam cannot be represented by a single nominal conductivity when the local solid fraction, reaction state, and transport role differ across the bed. The complete deterministic sweep values and the release milestones are summarized in the text and supplementary computational package; a separate main-text table was removed to respect the journal artwork limit. The comparative design context and the PG-SADE-GHF decision variables and search mechanisms are summarized in Table 3

TABLE 3. Combined design context: (a) positioning relative to the closest graded-metal-foam reactor studies; (b) decision variables, ranges, and search mechanisms used by PG-SADE-GHF. The panel format consolidates two related tables without removing their original information.

The hydrogen-fraction maps in Figure 3 are surrogate-reconstructed fields used to illustrate how the calibrated reduced model represents the expected wall-to-core depletion sequence. In the fixed benchmark cases, the reconstruction retains a high-X core after the wall region has begun to deplete; for the optimized surrogate trajectory, the reconstructed front moves inward earlier. This pattern is qualitatively consistent with the coupled behavior described by (Muthukumar et al., 2009), but it is not an independently solved hydrogen field for the new graded structure. We therefore use these maps only as mechanistic visualization and do not treat them as evidence beyond the scalar surrogate prediction.

3.3.1 Comparison with the closest graded-foam studies

The improvement predicted here should be interpreted relative to, rather than independently of, prior graded-foam reactor work. (Bai et al., 2022) demonstrated during absorption that a two-layer porosity gradient optimized by a genetic algorithm can outperform a uniform foam, while (Ye et al., 2026) showed that radial grading can reduce the 90%-absorption time by about 20% and that very high foam conductivity can reduce the incremental benefit of grading. (Zhu et al., 2026) further reported a graded-foam advantage in a metal-hydride/PCM configuration. These studies support the physical premise that conductivity should be placed non-uniformly, but they also mean that “graded foam” alone is not the novelty of the present paper. The contribution examined here is the desorption-side co-design of continuous porosity, hybrid alloy allocation, and a transient h / P schedule. The 20.5% improvement over the reproduced CuCr1Zr-89% benchmark is therefore identified explicitly as a reduced-model screening result. The exact-candidate PDE/OpenFOAM challenge does not confirm the original coordinates; however, same-bound OpenFOAM re-optimization independently recovers a 7.817-min optimum and the expanded-bound study reaches 6.442 min. The novelty claim is therefore the coupled structure–material–schedule optimization plus a transparent surrogate-to-CFD re-optimization workflow, not the assertion that one surrogate coordinate vector is universally optimal.

3.4 Optimization performance and staged verification

The lower two rows of Figure 4 are retained deliberately as the intermediate challenge of the exact surrogate-selected coordinates. Panels (e)–(h) progress from the calibration anchors and held-out optimized prediction, through the surrogate-screening reductions relative to the foam-free baseline, to the surrogate and frozen-model PDE hydrogen-depletion histories and the final parity comparison between target/surrogate times and PDE completion times. This sequence shows why simply reporting the 7.95-min reduced-model value was insufficient: the frozen MATLAB PDE verifier gives approximately 12.233 min for those coordinates. This mismatch is retained explicitly. Instead, Figure 8 adds the decisive next stage: actual OpenFOAM re-optimization of the same eight variables. The same-bound CFD optimum is 7.8167 min, whereas the separately reported expanded-bound optimum is 6.4417 min. Therefore, the final claim is not that OpenFOAM reproduces the exact original candidate, but that OpenFOAM independently recovers the original performance scale when the original optimization problem is solved at higher fidelity.

Figure 5 documents the optimization diagnostics rather than physical validation. Panel (a) gives the equal-budget fixed-DE versus SADE comparison; the adaptive method attains the lower final diagnostic objective in 29 of 30 paired runs, with median 0.1063 versus 0.1473 for fixed DE. Panel (b) now visibly plots the selected radial porosity and CuCr1Zr-fraction profiles, including the wall/core transition. Panel (c) shows the early thermal-pressure schedule, and panel (d) retains the surrogate release comparison as an optimization-stage output. These panels demonstrate search behavior and design interpretability, not correctness of the 7.95 min physical prediction.

Figure 6 documents the MATLAB-stage diagnostics, while Figures 8, 9 replace the former diagnostic-only bound and penalty argument with actual OpenFOAM evidence. Within the exact original bounds, the fine optimum is 7.8167 min; expanded bounds reduce it to 6.4417 min. Penalty multipliers 0.5, 1.0, and 2.0 do not change either selected completion time because the best evaluated candidates are feasible and have zero penalty. In the OpenFOAM transport solution the Forchheimer correction is retained directly; the maximum Forchheimer/Darcy ratio is 0.1587 for the primary original-bound optimum and 0.0958 on the finest expanded-optimum calculation. Thus, the high-fidelity analysis does not rely on a Darcy-only assumption at high porosity.

Finally, Figure 7 retains the independently solved MATLAB PDE fields for the exact original surrogate-selected candidate (The frozen-model prediction is 12.233 min versus 7.95 min from the surrogate, a 53.88% increase). Their 12.233-min completion time is now interpreted as the trigger for CFD-in-the-loop redesign rather than as the final performance conclusion. The OpenFOAM workflow then re-solves the optimization problem itself and produces the same-bound and expanded-bound results summarized in Figures 8, 9.

3.5 OpenFOAM traceability and same-bound recovery of the MATLAB performance scale

Figure 8 summarizes the full traceability chain. The uniform CuCr1Zr-89% benchmark is reproduced at 9.975 min. The exact MATLAB-selected coordinates give 12.192 min in OpenFOAM, confirming that the CFD model is not being forced toward 7.95 min. Most importantly, OpenFOAM re-optimization within the exact original MATLAB bounds yields 7.8167 min. The difference from 7.95 min is only − 1.677 %, while the reduction relative to the 10-min uniform benchmark is 21.833%. This is the primary quantitative result used to connect the original MATLAB optimization problem to the high-fidelity solver.

3.6 Architecture and operation fairness decomposition

A four-case matched-condition study was added to prevent the total speed-up from being attributed incorrectly to grading alone. Within the original MATLAB bounds, the graded optimum under its optimized schedule completes at 7.8167 min, while uniform CuCr1Zr-89% under the same schedule completes at 7.9417 min. Thus, the architecture-only advantage under identical optimized operation is 1.574%. On the same graded architecture, changing from the reference schedule to the optimized schedule improves completion by 20.911%. Under the expanded optimum, the graded and uniform designs under the same optimized schedule complete at 6.4417 and 7.4417 min, respectively, giving a 13.438% architecture effect, while the schedule effect on the same graded architecture is 26.938%. The physically correct interpretation is therefore a coupled architecture–operation benefit whose relative contributions depend on the admissible design envelope.

3.7 OpenFOAM numerical robustness, penalty sensitivity, and final design profiles

The expanded-bound optimum was subjected to an independent mesh/time-step sequence. The completion times are 6.4500 min on 12 × 24 with Δ t = 1.0 s, 6.4500 min on 24 × 48 with Δ t = 1.0 s, 6.4417 min on 24 × 48 with Δ t = 0.5 s, 6.4417 min on 36 × 72 with Δ t = 0.5 s, and 6.4417 min on 36 × 72 with Δ t = 0.25 s. The maximum difference from the finest result is 0.129%. The finest maximum Forchheimer/Darcy ratio is 0.0958. Re-scoring the available CFD candidates with penalty multipliers of 0.5, 1.0, and 2.0 leaves the original-bound and expanded-bound solutions unchanged at 7.8167 and 6.45 min, respectively, because the selected candidates have zero penalty.

The optimized coordinates clarify why the high-fidelity solution differs from the original surrogate point. OpenFOAM places more CuCr1Zr in the core-side transition, extends the conductive wall zone, uses stronger early heat transfer, and maintains a lower early outlet pressure for longer. The same-bound optimum remains within the original engineering envelope; the expanded solution moves further in the same physical direction. These trends are consistent with a stronger need for heat delivery and desorption driving force than estimated by the reduced model. These numerical and design trends are consolidated in Figure 9. Panel (a) restates the benchmark-to-CFD traceability, panels (b) and (c) separate architecture and schedule contributions for the original and expanded envelopes, panels (d) and (e) summarize mesh/time-step and penalty robustness, and panels (f) to (i) compare the radial porosity, CuCr1Zr-fraction, heat-transfer, and outlet-pressure profiles selected by the surrogate, the same-bound CFD optimum, and the expanded-bound CFD optimum.

3.8 Full-360° three-dimensional verification and field structure

The final optimized configurations were additionally solved on the complete cylindrical domain to test whether the axisymmetric wedge used during CFD-in-the-loop search biases the result. For the original-bound optimum, the wedge prediction is 7.8167 min and the full-360° result is 7.8083 min, a difference of only − 0.1066 %. For the expanded-bound optimum, both the wedge and standard full-3D calculations give 6.4417 min. The high-resolution 196,608-cell full-3D case with Δ t = 0.25 s again gives 6.4417 min, with final X ̄ = 0.009948, T ̄ = 346.913 K, p max = 0.9361 bar and maximum Forchheimer/Darcy ratio 0.09575. These results provide a direct circumferential-domain check in addition to the earlier wedge mesh/time-step study.

The execution and post-processing route is shown first in Figure 10, providing the software-level provenance of the calculation before the resolved fields are interpreted.

Panel (a) records the actual Windows PowerShell-to-Docker/OpenFOAM execution of the high-resolution full-360° case in the fixed opencfd/openfoam-dev:2506 environment and reports the completion metrics directly at run time. Panel (b) documents the ParaView post-processing stage used to interrogate the exported three-dimensional solution with a clipped cutaway and fixed camera; the representative view is colored by Darcy-velocity magnitude. One terminal diagnostic in panel (a) retains the legacy text “wedge domain”; this is a print-label artifact only, and the corrected pgsadeFoam3D source changes the wording to “full 360-degree domain” without altering the calculation. Figure 11 then visualizes the final three-dimensional field structure. The temperature field shows the wall-to-core thermal penetration that motivates the conductive near-wall region; the remaining-hydride field identifies the spatially nonuniform completion process underlying the capacity-weighted X ̄ ≤ 0.01 stopping criterion; and the pressure field shows the driving gradient for hydrogen transport toward the outlet. The porosity and CuCr1Zr-fraction maps show the material grading imposed by the optimized design. Taken together with the quantitative completion-time agreement, these fields support the paper’s central mechanism: performance is improved by coordinating where thermal conductivity is placed with when heat-transfer intensity and outlet pressure are applied. This is high-fidelity numerical confirmation of the design concept and its three-dimensional consistency, not experimental validation.

The software route summarized in Figure 10 closes the methodological chain between solver execution and scientific interpretation. The PowerShell/Docker panel confirms the reported completion metrics directly at run time, while the ParaView panel shows the clipping-based interrogation environment used to recover the same wall-heated, core-delayed desorption structure seen in Figure 11. This behavior is consistent with the heat-delivery limitations and graded-material advantages discussed for metal-foam hydride reactors in the related studies of (Bai et al., 2022), (Ye et al., 2026), and (Zhu et al., 2026).

3.9 Locked held-out multifidelity validation and uncertainty

The systematic comparison explains why the 53.35% mismatch of the exact MATLAB candidate is not an isolated anomaly. For the 32 development cases, the original reported surrogate gives R 2 = − 4.894, RMSE = 3.116 min, MAE = 2.831 min, and maximum absolute error = 5.990 min; 29/32 predictions are pinned to the 7.95-min safeguard. The discrepancy arises from three interacting simplifications: response compression by the lower safeguard, calibration on fixed reference cases rather than the complete graded design space, and omission of the locally resolved energy/pressure/Forchheimer coupling that is present in OpenFOAM. Accordingly, the original surrogate is retained as a rapid screening model, not as a globally accurate CFD emulator.

The stronger test is the locked 12-case challenge in Figure 12. Predictions were generated before the corresponding OpenFOAM outcomes and archived with a SHA-256 lock receipt. The original screening surrogate gives R 2 = − 6.910, RMSE = 3.354 min, and MAE = 3.143 min. The pre-locked CFD-corrected multifidelity model improves all 12 cases and achieves R 2 = 0.982, RMSE = 0.159 min, MAE = 0.121 min, maximum absolute error = 0.389 min, Pearson r = 0.995, and Spearman ρ = 0.986. This corresponds to a 95.27% RMSE reduction and a 96.14% MAE reduction. The nominal 80% and 95% prediction intervals cover 91.7% and 100% of the held-out CFD outcomes, respectively. Bootstrap 95% confidence intervals are R 2 = [ 0.951 , 0.994 ], RMSE = [ 0.079 , 0.228 ] min, and MAE = [ 0.071 , 0.184 ] min. The correction removes the dominant systematic bias while preserving the rule that final optimum values are established by direct CFD.

3.10 Experimental realization pathway

A concrete validation path is now defined rather than leaving future experiment unspecified. A first prototype should manufacture the conductive outer and more porous inner lattice regions as separately controlled additively manufactured inserts, avoiding an assumption of reliable direct CuCr1Zr/AlSi10 Mg co-melting before interface compatibility is established. X-ray micro-CT and gravimetric measurements should verify porosity, strut dimensions, and pore connectivity; effective-conductivity and pressure-drop tests should then identify k eff, permeability, and inertial coefficients before the lattice is packed with LaNi 5. The assembled reactor should be instrumented with radial/axial thermocouples, inlet/outlet pressure sensing, hydrogen mass-flow measurement, and controlled HTF temperature/flow. The primary comparison is measured t d, T ( r , z , t ), pressure history, and released hydrogen against the same-bound OpenFOAM prediction, followed by repeated absorption/desorption cycling. The critical uncertainties are foam–wall and foam–hydride contact resistance, manufactured pore size/PPI, inter-zone thermal resistance, alloy oxidation and dimensional tolerance, hydride settling/decrepitation, and kinetic drift with cycling. These quantities should be measured and propagated back into the open computational workflow.

4 Scientific position and validation route

The final scientific position distinguishes coordinate-level validation from validation of the optimization problem. The exact MATLAB surrogate-selected coordinates are not validated: they give 12.192 min in OpenFOAM rather than 7.95 min. However, this does not invalidate the optimization concept or the predicted time scale. When the same eight variables and exact original MATLAB bounds are optimized directly with OpenFOAM evaluations, the high-fidelity optimum is 7.8167 min, only 1.68% from the original 7.95-min magnitude and 21.83% faster than the reproduced 9.975-min CuCr1Zr-89% reference. We therefore distinguish “validation of the exact candidate” from “high-fidelity confirmation of the optimization problem and performance scale.” Only the latter is claimed.

The expanded-bound optimum of 6.4417 min is reported as a sensitivity extension, not as retroactive validation. Its mesh/time-step independence and zero-penalty status make it numerically robust, but several coordinates remain at expanded limits. The subsequent full-360° calculations strengthen the numerical verification: the original-bound full-3D result differs from its wedge counterpart by only − 0.1066 %, while the expanded-bound wedge, full-3D, and high-resolution full-3D calculations all give 6.4417 min. This confirms that the circumferential domain reduction used during optimization does not materially alter the reported completion time under the present axisymmetric boundary conditions. Experimental realization must still establish the actual foam pore size/PPI, contact resistance, alloy-interface behavior, additive-manufacturing tolerances, decrepitation and cycling durability. The OpenFOAM model is a porous-continuum CFD representation rather than a pore-resolved strut simulation. These limitations are stated explicitly. The systematic validation further changes the interpretation of the reduced model. Its role is now explicitly limited to low-cost screening, while the 32-case discrepancy model and the pre-locked 12-case challenge quantify how sparse OpenFOAM information can correct systematic low-fidelity bias. This distinction removes the previous implication that similarity of one completion-time magnitude establishes surrogate accuracy. It also provides an uncertainty-aware route for future optimization: screen cheaply, correct against selected high-fidelity cases, and confirm final candidates directly with CFD and, ultimately, experiment.

5 Conclusion

This study develops PG-SADE-GHF as a coupled structure–material–operation optimization framework for LaNi 5 hydrogen desorption and closes the high-fidelity verification gap through OpenFOAM v2506 CFD-in-the-loop re-optimization. The original MATLAB reduced model predicted t d = 7.95 min. Direct OpenFOAM evaluation of those exact coordinates gives 12.192 min, so the manuscript no longer equates one surrogate design vector with a validated CFD optimum.

The decisive result is the same-bound re-optimization. Using the exact original MATLAB decision ranges, OpenFOAM gives t d = 7.8167 min, 21.833% faster than the 10-min uniform CuCr1Zr-89% benchmark, while a wider sensitivity envelope reaches 6.4417 min. The systematic design-space study shows that the original surrogate should not be characterized as a quantitatively accurate CFD predictor: 90.6% of the 32 development predictions reach its 7.95-min lower safeguard. Instead, the reduced model supplies inexpensive screening, and the CFD-discrepancy layer supplies high-fidelity correction. The pre-locked 12-case challenge independently verifies this separation, with the corrected model reaching R 2 = 0.982, RMSE = 0.159 min and Spearman ρ = 0.986 while the final optimum remains a direct OpenFOAM result.

The numerical evidence is internally consistent: the expanded optimum varies by at most 0.129% across the tested wedge meshes/time steps, its finest maximum Forchheimer/Darcy ratio is 0.0958, and penalty re-scoring from 0.5 to 2.0 does not alter the selected solution. The full-360° verification provides an additional and more direct geometric check: 7.8083 min versus 7.8167 min for the original-bound full-3D and wedge calculations, respectively, and identical 6.4417-min completion for the expanded-bound wedge, standard full-3D, and 196,608-cell high-resolution full-3D cases. The matched-condition study also prevents over-attribution: in the original-bound problem most of the gain comes from operating-schedule optimization, while the expanded design exhibits a larger independent grading benefit. The strongest conclusion is therefore that the best desorption performance emerges from co-design of the radial conductive architecture and the transient pressure/heating schedule.

The resulting workflow provides a transparent hierarchy from MATLAB screening, through an independent PDE challenge, to OpenFOAM CFD re-optimization and full-360° numerical verification. The OpenFOAM solver source, case dictionaries, execution scripts, full-3D checks and visualization/export workflow are included with the supplementary computational package. The remaining requirement is experimental validation of the manufacturable graded foam and its long-term hydrogen-cycling behavior.

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.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fenrg.2026.1956543/full#supplementary-material

Supplementary material Supplementary material at the source

Nomenclature

A Van’t Hoff entropy-related coefficient B Van’t Hoff enthalpy-related coefficient (K) C d Desorption kinetic rate constant (1/s) c p Specific heat capacity at constant pressure (J/(kgK)) E d Desorption activation energy (J/mol) F i Mutation factor of the ith differential-evolution candidate F m Adaptive memory value of the mutation factor f C u Local CuCr1Zr fraction in the hybrid foam f C u , w CuCr1Zr fraction near the heated wall f C u , c CuCr1Zr fraction in the reactor core G max Maximum number of optimization generations h Convective heat-transfer coefficient (W/(m²K)) h early Early-stage heat-transfer coefficient (W/(m²K)) h late Late-stage heat-transfer coefficient (W/(m²K)) J Penalized optimization objective function K Hydride-bed permeability (m²) M Hydrogen-storage metal alloy M g Molecular mass of hydrogen gas (kg/mol) m ̇ Hydrogen desorption source term (kg/(m³s)) N P Population size in the differential-evolution algorithm p g Hydrogen gas pressure (Pa) p e q Equilibrium pressure (Pa) P d Desorption pressure (bar) P early Early-stage desorption pressure (bar) P late Late-stage desorption pressure (bar) p T Temperature-response shape exponent in the reduced model q X Hydrogen-fraction shape exponent in the reduced model r Normalized radial coordinate from wall to core r w Conductive wall-zone fraction R Universal gas constant (J/(molK)) S Smooth radial transition function t Time (s) t d Complete desorption time (min) t s Schedule switch time (min) T Local or average temperature (K T 0 Initial temperature (K) T f Heat-transfer-fluid temperature (K) u Superficial hydrogen velocity vector (m/s) w f Manufacturability penalty weight w m Metal-mass penalty weight w o Operating-intensity penalty weight w u Thermal-nonuniformity penalty weight x Optimization decision vector X Normalized hydrogen concentration or fraction X crit Desorption completion criterion x Stoichiometric coefficient in the metal-hydride reaction Δ H Metal-hydride reaction enthalpy (J/mol) Δ t Smooth duration of the operating schedule transition (min) δ r Smoothness parameter of the radial grading function ε Hydride-bed porosity ε m f Metal-foam porosity ε w Wall-region foam porosity ε c Core-region foam porosity λ Thermal conductivity (W/(mK)) λ eff Effective thermal conductivity of the hydride bed (W/(mK)) λ eff,mf Effective thermal conductivity of the foam-enhanced bed (W/(mK)) λ m f Local metal-foam thermal conductivity (W/(mK)) μ g Hydrogen dynamic viscosity (Pas) Π f Manufacturability penalty term Π m Metal-mass penalty term Π o Operating-intensity penalty term Π u Thermal-nonuniformity penalty term ρ 0 Initial solid density (kg/m³) ρ emp Empty-bed solid density (kg/m³) ρ g Hydrogen gas density (kg/m³) ρ s Solid hydride density (kg/m³) ρ sat Saturated solid density (kg/m³) σ Logistic operating-schedule function τ T Temperature time constant in the reduced model (min) τ X Hydrogen-fraction time constant in the reduced model (min)

Footnotes

Edited by: Rahul R. Bhosale, University of Tennessee at Chattanooga, United States

Reviewed by: Bahador Abolpour, Sirjan University of Technology, Iran

Bogu Liu, North China Electric Power University, China

Footnotes

AM, Additive manufacturing; CFD, Computational fluid dynamics; DE, Differential evolution; GHF, Graded hybrid foam; HTF, Heat-transfer fluid; MF, Metal foam; PCM, Phase-change material; PG-SADE-GHF, Physics-guided self-adaptive differential evolution for graded hybrid foam design; SADE, Self-adaptive differential evolution; UDF, User-defined function.