In this study, the multi-parameter constrained optimization procedure integrating the design of experiments (DOE), full factorial experimental design (FFED), genetic algorithm (GA) and computational fluid dynamics (CFD) are proposed to design a gas diffusion layer (GDL) linked to a wavy channel. The Forchheimer–Brinkman extended Darcy model and two-equation energy model are adopted to describe the fluid flow and heat transfer characteristics in the porous media. Numerical computations are performed with a wavy channel for the parameters studied. That is, Reynolds number (Re), wave amplitude (α), wave number (β) and gas diffusion layer porosity (ε) in the range of 200≤Re≤1000, 0.1≤α≤0.3, 2≤β≤10, 0.4≤ε≤0.5, respectively. The interactive effects of Re, α, β, and ε on the local Nusselt number (Nu) and the average Nusselt number (Nu¯) are discussed in this study. The results show that computations of local and average Nusselt numbers within the wavy channel are enhanced compared to the straight channel.In addition, the optimization of this problem is also presented by using a full factorial experimental design and the genetic algorithm (GA) method. The objective function defined as thermal performance factor (TPF) has developed a correlation function with three design parameters. That is, wave amplitude α, wave number β and gas diffusion layer porosity ε. The results show that the thermal performance factor (TPF) prediction can be approached under the maximum error of 6.2% from the present regression model, which is compared with CFD results.