Analysis of Surrogate Models for Vapour Transport and Distribution in a Hollow Fibre Membrane Humidifier

Abstract

To achieve high efficiency and low degradation of a polymer electrolyte fuel cell (PEMFC), it is necessary to maintain an appropriate level of humidification in the fuel cell membrane. Thus, membrane humidifiers are typically used in PEMFC systems. Parameter studies are important to evaluate membrane humidifiers under various operating conditions to reduce the amount of physical tests. However, simulative studies are computationally expensive when using detailed models. To reduce the computational cost, surrogate models are set up. In our study, a 3D computational fluid dynamics (CFD) model of a hollow fibre membrane humidifier is presented and validated using measurement data. Based on the results of the validated CFD model, a surrogate model of the humidifier is constructed using proper orthogonal decomposition (POD) in combination with different interpolation methods. To evaluate the surrogate models, their results are compared against reference solutions from the CFD model. Our results show that a Halton design combined with a thin-plate-spline interpolation results in the most accurate surrogate humidifier model. Its normalised mean absolute error for 18 test points when predicting the water mass fraction in the membrane humidifier is $0.58\%$. Furthermore, it is demonstrated that the solutions of the POD model can be used to initialise CFD calculations and thus accelerate the calculation of steady state CFD solutions.

1. Introduction

Polymer electrolyte fuel cells (PEMFCs) are considered to be a zero emission drive train solution for vehicles that can replace the currently used internal combustion engines [1,2,3]. The current state of research on PEMFCs has been outlined in various recent articles [1,4,5]. One major research area is the water management of PEMFCs [4,6,7]. Our study focuses on a membrane humidifier, which is an important component for the water management of PEMFCs. The purpose of using a membrane humidifier is to ensure sufficient humidification of the proton-exchange membrane (PEM) used in the PEM fuel cell. The humidification of the membrane is indirectly achieved via humidification of the supplied air. In general, the humidification of the PEM is important to reach a high efficiency and service lifetime of the fuel cell, as shown by [7,8,9]. The water used for the humidification is produced by the electrochemical reaction of hydrogen and oxygen in the PEM fuel. The produced water is carried out of the fuel cell by the air flow at the cathode. This wet air flow enters the so-called wet side of the humidifier. Inside the membrane humidifier, the entering water is partially transferred from the wet air flow to the dry supply air flow through the humidifier membrane. This membrane is the key component of the humidifier, which is highly permeable to water vapour but poses a resistance to the transport of other gases. In the literature, mainly two types of membrane humidifiers for PEM fuel cell systems are investigated and discussed: hollow fibre and flat membrane modules. Park et al. [10] and Vu et al. [11] analyse the water transfer of hollow fibre humidifiers experimentally. Furthermore, water transfer in flat membrane humidifiers is studied experimentally by [12,13]. In most cases, experimental work is complemented by modelling and simulation to analyse humidifiers under different operating conditions. Such simulation studies with varying parameter combinations have been carried out by [10,14,15,16], among others. An overview of the application of CFD models and simulations in the field of membrane technology can be found in [17]. CFD models of flat membrane humidifiers are set up and investigated in different studies and are shown to be an accurate tool to investigate the water transfer in membrane humidifiers [15,16,18,19]. Schmitz et al. present a 3D CFD model of a hollow fibre humidifier to analyse the sorption behaviour and the water transport in membrane humidifiers and validate the model with measurement data [20]. Even though only a small part of the flow domain is modelled for the CFD analysis, good agreement with measurement data is achieved [20], which motivates our approach to used a downscaled model. Another 3D CFD model of a hollow fibre membrane module was set up to analyse a gas separation process [21]. This model was successfully employed to analyse local concentrations and to optimise the efficiency of separation [21]. In summary, the potential of using CFD simulations to investigate humidifiers is well documented in the literature. Our study aims to set up a 3D CFD model of a hollow fibre humidifier that allows us to investigate the effect of parameter variations on water transfer and distribution in a membrane humidifier that has not yet been investigated in the literature to the authors’ knowledge. The parameter variation focuses on the operating conditions applied at the inlets of the humidifier, namely the air mass flow rates, the operating temperature and pressure and the amount of water available at the wet stream inlet. The parameters are listed in the vector $\zeta$:

$${\zeta }_{operation}=\left(\begin{array}{c}{\dot{m}}_{air,in}\\ T_{in}\\ Y_{H2O,wet,in}\\ p_{out}\end{array}\right)$$

(1)

Although the developed CFD model is derived from [16] and is built upon physical knowledge, some parameters of the model, especially of the membrane, are unknown and need to be identified using measurements. Therefore, measurements were conducted on a custom tailored test stand and are used for parameter calibration of the membrane model. Another measurement data set is used in the following to validate the calibrated CFD model.

A remaining challenge when making use of excessive simulation studies with 3D CFD models is that many high-fidelity models are time consuming to solve and require large amounts of computational resources [22]. To overcome this issue, surrogate models can be used, yielding a sufficient accuracy while significantly reducing computational cost with respect to time as well as hardware. In the field of aerodynamic research, the use of surrogate models based on proper orthogonal decomposition (POD) is well established, and many publications are available. An overview can be found in [22,23,24]. In addition, POD-based surrogate models have been successfully applied to the analysis of energy converters [25,26]. Bai et al. [25] show that POD in combination with regression splines leads to an accurate model of a PEM fuel cell. Those results motivate the transfer of the method to the humidifier. The current paper is focused on the application of POD as a dimensionality reductions technique in combination with an interpolation model (POD+I) to accelerate parameter studies for steady-state results. For this method, simulation results at different operation conditions are calculated using the CFD model. Thereon, the CFD results are stored in a matrix and POD is applied to it, yielding a set of POD modes and coefficients. The POD coefficients are a function of the parameter combinations, as shown by [23]. This relationship between the parameter vector given in Equation (1) and POD coefficients can be captured with an interpolation or regression model [23]. Using this interpolation or regression model, the POD coefficients for an untried parameter combination, i.e., a combination not included in the training data for the surrogate model, can be obtained [23]. A more detailed explanation of the used methods are given in Section 3. One clear advantage of the investigated POD methods compared to scalar surrogate models is that a 3D solution at untried parameter combinations is acquired, which provides information about the field data and can additionally be used to calculate scalar quantities [22]. The creation of a surrogate model using POD+I requires computational time and resources for calculating solutions of the high-fidelity model, also called full-order model. Those high-fidelity results serve as information upon which the reduced order model is built and should therefore contain as much information about the problem to be solved as possible. To keep this upfront investment to a minimum, the parameter combinations that are investigated by the high-fidelity model are chosen using design of computer experiments (DoCE) methods. In summary, the aims of our study are:

• To build a CFD model representing the hollow fibre membrane humidifier and to validate by experimental data.
• To apply the method of POD+I to the hollow fibre membrane humidifier.
• To examine if POD+I is a suitable tool to accelerate parameter studies for such mass exchangers.
• To identify the most suited DoCE to build the surrogate model upon.
• To find a suited interpolation method to be used by the POD+I model.

To the best knowledge of the authors, POD+I methods have never been applied to membrane mass exchangers. Therefore, the potentials of the aforementioned methods are analysed in the present study.

2. Test Stand and Investigated Hollow Fibre Membrane Humidifier

2.1. Description of the Hollow Fibre Membrane Humidifier

A sketch of a hollow humidifier geometry is schematically depicted in Figure 1. In the sketched geometry, a small number of 12 fibres is depicted to clearly highlight the features and flow situation. A similar geometry is used in the CFD simulations in Section 4. As shown in Figure 1, the wet and dry air streams are in a counterflow arrangement to achieve high vapour transfer rates. As seen in Figure 1, the wet air stream is fed to the fibres, whereas the dry stream flows through the shell. As a result of the manufacturing process, the fibres are placed randomly inside the shell; therefore, the water transfer behaviour is expected to be inhomogeneous. The correlation of fibre placement and water distribution is investigated in the CFD results. A potting technique is used at both ends of the fibres to separate the wet and dry sides of the humidifier.

The humidifier investigated in the experimental studies is commercially available from the company Fumatech [27]. It is designed to be used in combination with a PEM fuel cell, offering an electrical power output of about 700 W rated at an air flow of 50 sLPM. A total of 488 hollow fibres are placed randomly into the housing. The hollow fibre membranes are made of an undisclosed polymer composite consisting of Polyimide and Polysulfone; therefore, the permeation properties of the membrane are unknown. The relevant geometrical data of the fibres and the housing are given in

Table 1. Measured geometric parameters of the hollow fibre membrane humidifier.

Quantity	Symbol	Unit	Value
Number of fibres	$n_F$	1	488
Fibre outer diameter	$d_{F,o}$	mm	1
Fibre inner diameter	$d_{F,i}$	mm	0.9
Fibre length	$l_F$	mm	150.8
Housing inner diameter	$d_H$	mm	39.2
Packing density	$\varphi$	m${}_{\mathrm{mem}}^2$/m${}^3$	1257

.

The packing density of the module is calculated by the following equation:

$$\varphi =\frac{A_{mem}}{V_H}=\frac{n_F\pi \left(d_{F,o}\right)l_F}{\pi {(d_H/2)}^2{l}_F}$$

(2)

This quantity describes how much transfer area is theoretically available for a mass transfer in a given volume and can often be found in the literature.

2.2. Description of the Test Stand

The purpose of the test stand is to investigate the mass transfer capabilities of membrane humidifiers at various operating conditions typical for a PEM fuel cell system. The operating conditions that can be varied and investigated using the test stand are:

• mass flow rates of both, dry and wet, streams;
• temperatures of both, dry and wet, streams;
• pressures of both, dry and wet, streams;
• relative humidity of the wet stream.

The water transfer can be calculated based on sensor values of the dry side:

$${\dot{m}}_{H2O,dry,meas}={\dot{m}}_{air,dry,out}{Y}_{H2O,dry,out}-{\dot{m}}_{air,dry,in}{Y}_{H2O,dry,in}$$

(3)

and on the sensors used on the wet side:

$${\dot{m}}_{H2O,wet,meas}={\dot{m}}_{air,wet,in}{Y}_{H2O,wet,in}-{\dot{m}}_{air,wet,out}{Y}_{H2O,wet,out}$$

(4)

A steady-state operating point is only considered valid when both results agree with a tolerance of 5%. For the comparison with the simulations, both results are averaged.

A P&ID of the humidifier test stand is depicted in Figure 2. All relevant sensors used in the test rig are listed in

Table 2. Sensors used in the humidifier test stand and their specified measurement uncertainties.

Sensor	Measured Quantity	Output Unit	Uncertainty
Vaisala HMT-337	Humidity	%	$\pm (1.5+0.015\phi )$
Omega FMA-1609A	Mass flow rate	g/s	$\pm (0.008\dot{m}+0.00204)$
Omega PXM459	Differential pressure	Pa	$\pm 56$
WIKA P-30	Pressure	bar	$\pm 0.068$
WIKA TR-40	Temperature	K	$\pm 0.15+0.002\ast (T-273.15)$

. Additionally, the uncertainties of the sensors specified by the manufactures are given in Table 2. The presented uncertainties are used to calculate the error propagation used in the validation process in Section 4.1.

The required air is fed to the test stand from a pressurised air storage tank. Behind the air storage, the air path is divided into two lines: one for the wet and one for the dry air flow. The wet path simulates the exhaust gas from the fuel cell and the dry path the supply air. A valve at the inlet of each path is used to control the individual flow rates of the paths. In a next step, both air streams are heated up by electrical heaters to simulate operating temperatures typical for PEM fuel cell systems. The pressure levels of both streams can be set individually by two control valves located near the outlet of the test stand. For additional information, the pressure drop over the humidifier in the test stand is measured for both air streams. In the wet path, vapour is added from a vapour storage tank through a controlled valve to achieve the desired inlet humidity of the air. Typically, the air leaving the PEM fuel cell carries much water produced by the fuel cell reaction. The supply line of vapour is heated to avoid condensation, which is necessary to achieve valid measurement results. The installation of this heater is important because the moisture sensors installed and used to calculate the water transfer are only able to detect gaseous water. Therefore, additional heaters are placed at the outlets of the humidifier to ensure that all present water is in gaseous form. To make sure that all water is evaporated, the temperatures are measured and controlled to be higher than the saturation temperature. Additionally, all tubes as well as the box, wherein the humidifier is positioned, are isolated and heated to avoid condensation of water in the test stand. If condensation of water occurs, this leads to a deviation in the water transfer rates measured on the dry and wet side. Thus, the water transfer rates are continuously checked during the measurement.

3. Modelling of the Hollow Fibre Humidifier

3.1. CFD Model of the Hollow Fibre Humidifier

The code basis of the developed humidifier model was used in another publication for a flat membrane humidifier [16] and is based on the freely available library OpenFOAM [28,29]. In the CFD model, moist air is described as an ideal mixture consisting of the gas components dry air and water vapour. All necessary thermophysical properties for the fluids are calculated from the JANAF data [30], which are directly available in OpenFOAM. All presented balance equations originate from the solver “reactingFoam” [28]. The mass balance for a compressible fluid in its differential form is:

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \overrightarrow{U}\right)=0$$

(5)

The momentum balance for the compressible fluid reads:

$$\frac{\partial \rho \overrightarrow{U}}{\partial t}+\nabla ·\left(\rho \overrightarrow{U}\overrightarrow{U}\right)=-\nabla p+\nabla ·\left(2\eta \mathit{S}\left(\overrightarrow{U}\right)\right)-\nabla \left[\frac{2}{3}\eta (\nabla ·\overrightarrow{U})\right]$$

(6)

The term $\mathit{S}$ represents the rate of strain tensor. The calculation of the viscosity follows the Sutherland equation [28]:

$$\nu =\frac{A_s\sqrt{T}}{1+T/T_s}$$

(7)

The coefficients $A_s$ and $T_s$ were fitted to match data available in the VDI 4670 [31] for a temperature range from 0 up to $100{}^{\circ }\mathrm{C}$. The viscosity given in the data is reproduced accurately by the Sutherland equation.

In this work, only low Reynolds numbers up to 500 are considered for the flow in the fibres and the shell of humidifier so that the flow should be predominantly laminar. Hence, no turbulence model is applied. The limitation of the Reynolds Number arises from limitations of the mass flow rates that can be investigated by the test rig. The model can be used for higher Reynolds numbers, but no validation data are available from the measurements.

The fluid used in the model is a gas mixture consisting of the components dry air and water vapour. The mass balance of species i in a mixture is:

$$\frac{\partial \rho Y_i}{\partial t}+\nabla ·\left(\rho Y_i\overrightarrow{U}\right)=\nabla ·(\rho D_i\nabla Y_i)$$

(8)

The diffusion coefficient $D_i$ is calculated with Equation (9), taken from [32]:

$$D_i(p,T)=D_{i,ref}{\left(\frac{T}{T_{ref}}\right)}^{3/2}\left(\frac{p_{ref}}{p}\right)$$

(9)

Since only two components are present in the mixture, the binary diffusion coefficient of vapour and air is used. In the reference state, at $T_{ref}=298$ K and $p_{ref}=1$ atm, the diffusion coefficient is $D_{H2O,Air}=0.260$ cm${}^2$/s [32].

The energy balance for the gas mixture is formulated in terms of total enthalpy balance:

$$\frac{\partial \rho h}{\partial t}+\frac{\partial \rho K}{\partial t}+\nabla ·\left(\rho h\overrightarrow{U}\right)+\nabla ·\left(\rho K\overrightarrow{U}\right)=\frac{\partial p}{\partial t}+\nabla ·(\kappa \nabla T)+\nabla ·\sum _i^N\rho D_i{h}_i\nabla Y_i$$

(10)

The second term on the right-hand side represents the heat conduction. In the solution process, this term is reformulated as a function of enthalpy, as shown by Kumar [33]. The last term on the right-hand side is also taken from [33]. It represents the energy transfer due to diffusion of the modelled species.

Model of the Polymer Membrane and Boundary Conditions for the Flow Simulation

The implementation of the polymer membrane boundary condition is based on the class semipermeableBaffle [28]. It is a two-sided boundary condition that connects two fluid regions of the membrane humidifier, typically a fibre with the shell. Inside the membrane model chosen for the hollow fibre humidifier, a mass flow rate is calculated based on a driving concentration difference for each species and for each face of the patch. The calculated mass flow rates are then introduced into the aforementioned equations of the fluid volumes directly next to the membrane.

The permeability of the membrane was not disclosed by the manufacturer; therefore, a calibration using measurement results was conducted to determine the best value for the permeability of vapour. A constant value of $P_{H2O}=4.02\times {10}^{-7}$ m/s resulted in a good calibration and validation result. The calibration and validation processes utilising surrogate models as well as the simulation results are presented in Section 4.1. The membrane area available for mass transfer was calculated as the mean of inner and outer membrane area. The thickness of the membrane is ${\delta }_{mem}=1\times {10}^{-4}$ m.

$${\dot{m}}_{H2O,perm}=P_{H2O}{M}_{H2O}{A}_{mem}/{\delta }_{mem}\Delta c_{H2O}$$

(11)

An accumulation of water in the membrane is not considered.

The permeation of dry air through the membrane is not completely prevented in reality, therefore, the possibility for permeation of dry air is also implemented in the membrane boundary condition. The difference of partial pressure across the membrane is the driving potential for permeation. The entire permeation process is calculated with a constant permeation coefficient $P_{air}=1\times {10}^{-14}$ mol/(m s Pa) obtained form the literature, as suggested by [34]:

$${\dot{m}}_{air,perm}=P_{air}{M}_{air}{A}_{mem}/{\delta }_{mem}\Delta p_{air}$$

(12)

The membrane also connects the wet and dry sides thermally, but, due to the chosen nearly isothermal conditions, the description is omitted here.

The other boundary conditions are shown in

Table 3. Boundary conditions for the simulation of the hollow fibre membrane humidifier.

	p	$\overrightarrow{\mathit{U}}$	T	${\mathit{Y}}_{\mathit{H}2\mathit{O}}$	${\mathit{Y}}_{\mathbf{Air}}$
inlets	zero gradient	fixed value	fixed value	fixed value	fixed value
outlets	fixed value	zero gradient	zero gradient	zero gradient	zero gradient
walls	zero gradient	no slip	zero gradient	zero gradient	zero gradient

and their specification is given in the OpenFOAM documentation [28,29].

3.2. Down-Scaled Geometry Used for the Simulation of the CFD Model

The discretisation of the hollow fibre humidifier geometry given in Table 1 would result in many volume elements, which, as a consequence, would make the solution of the full humidifier very time consuming and expensive. To overcome this issue, only a small fraction of the humidifier is simulated. This is done under the assumption that the flow inside the shell and the fibres is distributed homogeneously. In the scaling process, the number of fibres to be simulated is reduced from 488 to 12, which results in a scaling factor of 40.6. Furthermore, the shell diameter is adjusted in order to obtain the same packing density in the scaled and the full-size humidifier. The scaling process does not affect the lengths and diameters of the individual fibres. The parameters of the scaled geometry are listed in

Table 4. Geometric parameters of the down-scaled hollow fibre membrane humidifier.

Quantity	Symbol	Unit	Value
Number of fibres	$n_F$	1	12
Fibre outer diameter	$d_{F,o}$	mm	1
Fibre inner diameter	$d_{F,i}$	mm	0.9
Fibre length	$l_F$	mm	150.8
Housing inner diameter	$d_H$	mm	6.16
Packing density	$\varphi$	m${}_{\mathrm{mem}}^2$/m${}^3$	1257
Scaling factor	-	1	40.6

.

Additionally, the mass flow rates applied to the model are scaled in the same manner so that the mass fluxes in the shell and in the fibres are the same in measurement and simulation. The results achieved with the down-scaled model are then scaled up to the full size of the humidifier.

3.3. Surrogate Model of the Hollow Fibre Humidifier

The process to build a POD+I surrogate model follows the steps described by Franz [23] for an aerodynamic problem and is presented in Figure 3.

In the first step, the parameter limits have to be defined. Next, a DoCE is calculated to place the parameter samples inside the defined parameter limits, and the high-fidelity model is evaluated at those parameter combinations. Both aforementioned steps are discussed in the next section. To obtain the maximum amount of information out of the computational expensive full-order-model (FOM) calculations, the parameter sets for those simulations are selected using a special design for computer experiment introduced in Section 3.4. The surrogate model is then created based on these steady-state results, so-called snapshots. These snapshots can also be called training data for the surrogate model. A snapshot contains the result variables of the steady-state solution in all the cell centres of the finite volumes in the fluid domain at a given parameter combination ${\zeta }_i$. Following the method of snapshots [35], snapshots are grouped into the matrices W of size ${\mathbb{R}}^{m\times n}$ for each result variable of interest (e.g., temperature, water mass fraction or pressure), where n is the number of collected snapshots and m is the number of volumes in the fluid domain:

$$W=\left(\begin{array}{cccccc}{w}_1^1& w_1^2& \cdots & w_1^i& \cdots & w_1^n\\ w_2^1& w_2^2& \cdots & w_2^i& \cdots & w_2^n\\ ⋮& ⋮& ⋮& ⋮& ⋮& ⋮\\ w_m^1& w_m^2& \cdots & w_m^i& \cdots & w_m^n\end{array}\right)$$

(13)

Each snapshot depends on the operating conditions $W^i\left({\zeta }^i\right)$ and represents a column of the given matrix W. Often, the overall mean of the snapshot matrix is subtracted from each snapshot without loss of generality.

To obtain the POD-base, the singular value decomposition of the snapshot matrix W is computed:

$$W=\Phi \Sigma V^T$$

(14)

where $\Phi \in {\mathbb{R}}^{m\times m}$ and $V\in {\mathbb{R}}^{n\times n}$ are orthogonal matrices. The columns of $\Phi$ are called POD modes. The resulting diagonal matrix $\Sigma \in {\mathbb{R}}^{m\times n}$ contains the singular values of W on the diagonal, ordered by their magnitude in decreasing order.

A certain snapshot $W^i\left({\zeta }^i\right)$ results from a linear combination of the POD modes and the POD coefficient $b^i\left({\zeta }^i\right)$, which, itself, depends on the vector of operating conditions:

$$W^i=\sum _{j=1}^m{b}_j^i{\Phi }^j$$

(15)

Assuming that the magnitude of the singular values decreases rapidly, the number of modes, k, to be used can be determined by looking at the variance captured by the first k singular values compared to all singular values. This metric is called the Relative Information Criterion and is calculated as follows:

$$RIC\left(k\right)=\frac{{\sum }_{i=1}^k{\sigma }_i^2}{{\sum }_{i=1}^p{\sigma }_i^2}.$$

(16)

In order to obtain an approximate solution $W^{*}$ of a specific parameter combination ${\zeta }^{*}$, that is not included in the snapshot matrix, an arbitrary interpolation scheme can be used to calculate the POD coefficient $B^{*}=B\left({\zeta }^{*}\right)$. The approximate solution of the surrogate model at the specified parameter combination is determined using the POD modes and the interpolated POD coefficients:

$$\widehat{W}=\sum _{j=1}^m{B}_j^{*}{\Phi }^j$$

(17)

The vector of parameters depends on the problem to be solved. For the analysis at hand, the parameters to be analysed are the mass flow rates, the operating temperature, the operating pressure and the mass fraction of water at the wet inlet.

To evaluate the surrogate model at untried parameter combinations, regression methods for the POD coefficients can be used [23,36]. Franz [23] uses radial basis functions as a surrogate model based on a thin plate spline (TPS) due to its good approximation capability. The linear, thin-plate-spline and the cubic method employed in this study are taken from the library SMARTy [37].

3.4. Design of Computer Experiments

The Design of Computer Experiments is a special branch of the classical Design of Experiments (DoE) introduced by Fisher [38]. A major difference between computer and physical experiments is that many computer experiments are fully deterministic. This means that a given set of input parameters, often called factors, always produces exactly the same output in a computer experiment, whereas, on a physical test stand, the results will vary due to random factors influencing the experiment [38]. The CFD model presented above is such a deterministic model. Therefore, no repetitions of certain input sets are necessary nor specifically wanted in computer experiments [39,40]. Due to the fact that, often, a small subset of input parameters causes the major effects in the response variables, a DoCE design should avoid repetitions in all dimensions [40]. Garud et al. [41] studied different designs and analysed different space-filling criteria to judge the performance of the individual DoCEs. They conclude that Quasi-Monte-Carlo methods such as the Halton and Sobol sequence are well suited for simulations with a number of input parameters below eight [41]. In their analysis, they analyse sample counts greater than 100 [41], which are not suitable for the evaluation of the membrane humidifier CFD model. The goal of our study is, therefore, to reduce the sample points to a minimum due to the high amount of resources necessary for the evaluation of the full-order model. In order to achieve a specific design of experiments for the investigation of the humidifier, the parameters are given in the introduction. For each of the parameters, an upper and a lower limit must be defined.

Each simulation is done with one parameter vector containing a set of these variables:

$${\zeta }_{operation}^i=\left(\begin{array}{c}{\dot{m}}_{i,air,in}\\ T_{i,in}\\ Y_{i,H2O,wet,in}\\ p_{i,out}\end{array}\right)$$

(18)

A matrix of parameter combinations, $Z\in {\mathbb{R}}^{d\times n}$, is generated when applying a DoCE method.

Three different experimental designs are used and analysed in this study: a latin hypercube design (LHD), a Halton sequence and a fractional factorial design from classical DoE. The LHD is often used in engineering literature [41] for DoCE purposes. As mentioned before, the Halton and Sobol sequences showed good space-filling abilities for a number of four input parameters and small sample counts [41]. Both the Halton sequence and Latin-Hypercube designs can be generated using Scipy [42]. The fractional factorial design poses the advantage that the minima and maxima of the parameter domain are occupied. Due to its screening capabilities, the fractional factorial design is often used in early stages of an experimental investigation [38]. The number of samples used for the Halton and LH designs are 8 and 17. Eight points were chosen, because this number of runs is considered to be the lower limit for a statistical analysis [38]. On the other hand, 17 points are arguably the minimum number of runs that are necessary for building a nonlinear model in classical DoE when four factors are taken into account [38]. The fractional factorial design consists of eight points at the outer bounds of the parameter domain and is well suited to set up a model containing the main effects. In order to achieve better space filling, a centre point is added to the fractional factorial design. This extension was carried out because it involves little computational effort. The low number of model evaluations is intended in this work to minimise the calculation time for the creation of the POD+I model. In addition, it should be shown which results can already be achieved with a small data set. A graphical overview is given in Figure 4.

The values of temperature, pressure, air mass flow rate and water mass fraction are determined by each of the DoCEs within the limits specified in

Table 5. Parameter ranges for simulation of the membrane humidifier. Water mass fraction and relative humidity are applied to the wet inlet only. The air flow rates are scaled to match the model.

Parameter	Minimum	Maximum
Temperature	52 °C	93 °C
Pressure	1.4 bar	2.1 bar
Air flow rate	5 mg/s	15 mg/s
Water mass fraction	$0.03$	$0.25$

. Water mass fraction is first set by the DoCE and, in subsequent steps, is limited to the saturation mass fraction valid for the pressure and temperature at that point to prevent a relative humidity above 100%. This is necessary because the model does not take the phase change from vapour to liquid into account. This bias is not considered to be a drawback, because it reflects the operating conditions in the fuel cell system, where the relative humidity of the wet inlet is also often close to saturation humidity.

4. Results and Discussion

4.1. Calibration and Validation of the CFD Model

The aim of the calibration is to identify the undisclosed membrane permeability coefficient. In the calibration process, the model was simulated with different values for the membrane permeability $P_{H2O}$ at four selected input parameter combinations that are characterised by having different operating pressure, temperature, water mass fraction and mass flow rates. The coefficient of determination was used as metric for calibration and is defined as:

$$R^2({\dot{m}}_{H2O,diff,sim},{\dot{m}}_{H2O,diff,meas})=1-\frac{{\sum }_{i=1}^n{({\dot{m}}_{H2O,diff,sim}-{\dot{m}}_{H2O,diff,meas})}^2}{{\sum }_{i=1}^n{({\dot{m}}_{H2O,diff,sim}-{\overline{\dot{m}}}_{H2O,diff,meas})}^2}$$

(19)

It measures the amount of variance present in the data that can be explained by the model. The fitted permeability coefficient was then used in validation simulations at further operating points. In

Table 6. Results of the calibration process of the membrane permeability for vapour.

${\mathit{P}}_{\mathit{H}2\mathit{O}}$ in m/s	${\mathit{R}}^2$ for Calibration	${\mathit{R}}^2$ for Validation
$4.02\times {10}^{-7}$	0.9990	0.995

, the coefficient for calibration and validation is given.

In Figure 5, the diffusion flow rates measured in the experiments are plotted against the simulated ones. The measurement uncertainty was calculated by the GUM method [43], and the sensor data are given in Table 2. The diffusion flow rate is one of the most important metrics when evaluating the humidifier performance. In order to compare model and simulation, the results of the simulation model were scaled to the size of the experimentally investigated humidifier. The operating temperature was varied between 60 and 80 °C, the pressure took values between 1.5 and 2.0 bar and the water mass fraction at the wet inlet was controlled between 0.0278 and 0.1728 (kg H${}_2$O)/(kg Dry Air). Overall, the simulation and experimental results agree well over the whole investigated parameter range. In

Table A1. Comparison of measured and simulated results, upscaled to the full size humidifier. The points used for the calibration of the diffusion coefficient of the membrane are marked with a *.

ID	${\dot{\mathit{m}}}_{\mathit{wet},\mathit{in}}$	${\dot{\mathit{m}}}_{\mathit{dry},\mathit{in}}$	${\mathit{T}}_{\mathit{dry},\mathit{in}}$	${\mathit{T}}_{\mathit{dry},\mathit{in}}$	p	${\mathit{\xi }}_{\mathit{H}2\mathit{O},\mathit{wet},\mathit{in}}$	${\dot{\mathit{m}}}_{\mathit{diff},\mathit{exp}}$	${\dot{\mathit{m}}}_{\mathit{diff},\mathit{sim}}$
$[-]$	$\mathbf{g}/\mathbf{s}$	$\mathbf{g}/\mathbf{s}$	$\mathbf{K}$	$\mathbf{K}$	$\mathbf{bar}$	$\frac{\mathit{kg}{\mathit{H}}_{\mathbf{2}}\mathbf{O}}{\mathit{kg}\mathit{Dry}\mathit{Air}}$	$\mathbf{g}/\mathbf{s}$	$\mathbf{g}/\mathbf{s}$
0	0.44	0.40	346.65	343.15	1.75	0.0928	0.0279	0.0275
1	0.41	0.40	336.21	333.13	1.75	0.0278	0.0071	0.0077
2 *	0.41	0.40	333.36	333.30	1.75	0.0296	0.0077	0.0083
3	0.44	0.40	345.41	343.13	1.50	0.0977	0.0281	0.0283
4 *	0.43	0.41	344.56	343.14	2.00	0.0778	0.0232	0.0232
5	0.42	0.40	344.50	343.15	2.00	0.0437	0.0118	0.0128
6	0.32	0.30	347.66	343.14	1.75	0.0486	0.0098	0.0111
7	0.32	0.30	343.81	343.13	1.75	0.0669	0.0141	0.0155
8	0.42	0.40	356.24	353.15	1.75	0.0584	0.0152	0.0168
9 *	0.48	0.40	357.71	353.11	1.75	0.1698	0.0535	0.0535
10	0.48	0.40	356.53	353.14	1.75	0.1728	0.0543	0.0551
11	0.55	0.50	343.45	353.16	1.75	0.0835	0.0279	0.0289
12	0.42	0.40	333.77	333.48	1.75	0.0573	0.0161	0.0164
13	0.45	0.40	353.64	343.32	2.00	0.1092	0.0338	0.0337
14	0.42	0.40	343.91	343.28	2.00	0.0542	0.0148	0.0159
15	0.43	0.40	347.00	343.15	1.75	0.0751	0.0212	0.0219
16 *	0.53	0.50	344.90	343.15	2.00	0.0602	0.0203	0.0209
17	0.32	0.30	342.92	343.13	2.00	0.0611	0.0138	0.0142

, all results are displayed. The maximum of the mean absolute error between measurement and simulation is 0.0016 g/s at point 8 in Table A1. Furthermore, the maximum relative deviation is 13.1% for point 6 in Table A1. At points with high vapour transfer rates, the measurement uncertainty is the highest due to the fact that the uncertainty of vapour measurements increases when the relative humidity increases. All simulation results of the vapour diffusion rate are within the range of the measurement uncertainty. An additional metric to assess humidifier performance is the mass transfer efficiency, as described by Brandau and Köhler [8]:

$${\eta }_1=\frac{c_{H2O,1}^{″}-c_{H2O,1}^{\prime }}{c_{H2O,2}^{\prime }-c_{H2O,1}^{\prime }}$$

(20)

This quantity describes the ratio of the actual to the maximum possible water transfer. The mass transfer efficiency predicted by the simulation model is plotted against the measured one in Figure 5, and is solely used for validation purposes. Hence, it is not accounted for during the fitting process. It can be seen that the efficiency for all measured points is in the range from 0.62 to 0.72. The mean absolute error in the efficiency is 0.023. It it clearly visible in Figure 5 that the measurement uncertainty is quite high for this quantity. This can be explained by the fact that three moisture measurements, which are subject to significant uncertainty, are included in the calculation of the mass transfer efficiency, cf. Equation (20). The highest deviations that can be observed in Figure 5 occur at operating points with relatively low vapour transfer rates. Due to this low value for the absolute vapour transfer rate, small absolute deviations cause a high relative deviation.

4.2. Set-Up and Evaluation of the Surrogate Models

The first step in the procedure is to calculate the sampling plan, where the FOM is evaluated. Based on these results, the POD-basis is calculated. In Figure 4, the samplings of the different DoCEs are shown together with the measurement points used for calibration and validation. The mass flow rates of the CFD model are scaled to match the experiment. It is clearly visible that the fractional factorial design, which was extended by a centre point, places points in the corners of the design space. Only five points can be seen in the pair plots, except for the ones including the water mass fraction, due to the fact that always two points are sampled with the same temperature, pressure and mass flow rate. In terms of absolute humidity, seven points can be seen, which is because of the limitation posed to relative humidity, as explained above. As a consequence, the points are shifted from the right to the left part of the pair plots using the water mass fraction as the X-axis. In contrast to the fractional factorial design, the sampling points generated by the Halton sequence and the Latin-Hypercube are more evenly distributed in the design space. In the pair plots of temperature and water mass fraction, unfilled areas can be observed that result from the limitation to 100% relative humidity at the wet inlet. Both the Halton sequence and Latin-Hypercube do not guarantee that the corners of the design space are filled with points. In fact, without adjustments, no points will be placed in any corner of the design space. Hence, the choice of parameter limits should be adjusted, or additional corner points should be added manually in order to avoid extrapolation with the surrogate model.

After the solutions of the FOM are calculated at the given sampling points, the POD-basis can be calculated. To evaluate the surrogate models using the different DoCEs in combination with the interpolation methods, the surrogate model is evaluated at the points that were used above for validation, but this is not included in the DoCEs, cf. Table A1. To assess and compare the set-up surrogate models, the following performance metrics are used:

The maximum absolute error that occurs in the whole fluid domain:

$$maxAE=max\left(\right|y_{surr,i}-y_{hf,i}\left|\right)$$

(21)

This metric is used to identify if there are outliers where the surrogate model predictions are very inaccurate. The mean absolute error is given by:

$$MAE=\frac{1}{n}\sum _{i=1}^n|y_{surr,i}-y_{hf,i}|$$

(22)

where n is the number of mesh elements. The most important field variable in the membrane humidifier is the water mass fraction driving the water transport from the wet to the dry side. The water mass fraction is therefore used to calculate the MAE and maxAE. Additionally, the normalised mean absolute error (nMAE) is used, where the MAE is divided by the maximum driving potential for the water transfer:

$$nMAE=\frac{MAE}{max\left(y_{hf}\right)-min\left(y_{hf}\right)}$$

(23)

Using this quantity, it is easier to compare simulations with different boundary conditions for the water mass fraction.

An overview of the MAE, maxAE and nMAE for the different surrogate models combining an interpolation method and a DoCE is given in Figure 6. The results given in Figure 6 are averaged over the 18 validation points to assess the overall performance of the surrogate models.

In Figure 6, it can be seen that the surrogate models based on the fractional design produce the largest deviations in all metrics. Comparing the averaged metrics in Figure 6 to each other, one can see that the results are similar to each other. Only for the fractional design using a linear interpolation an outlier in terms of maxAE can be identified. The MAEs, maxErr and nMAEs of the surrogate models using the fractional factorial design are significantly higher compared to the MAEs and nMAEs produced by surrogate models using LHD and Halton designs. Even though the deviations of the surrogate models using the fractional factorial sampling as basis are much higher than the other ones, the nMAE is below 2% for the water mass fraction in all cases. Furthermore, the results of the surrogate models using the LHD are more sensitive to the choice of the interpolation method than the Halton sequence-based ones. Overall, the Halton sequence-based models are quite robust in terms of the chosen interpolation method, since average results do not deviate much. The best combinations of an interpolation method with a given DoCE for the surrogate model of the humidifier are:

• the thin-plate spline with the fractional factorial design,
• the cubic spline with LHD,
• the thin-plate spline with the Halton sequence,

as depicted in Figure 6. On average, the model using the combination of the Halton sequence and thin-plate spline gives in the best surrogate model. Choosing this combination yields a MAE for a surrogate model below 1%. Moreover, Halton yields a deterministic sampling plan, which can be easily extended if more sample points are required, which is generally not the case for LHD-based sampling plans.

Figure 7 shows the results given in Figure 6 broken down by the operating points given in Table A1. However, only the best surrogate models for a given DoCE are considered in further analysis. For each operating point in Table A1, three values expressing the MAE, the maxErr and the nMAE of the water mass fraction predictions of the best surrogate models based upon the three DoCEs are shown. From the results given in Figure 7, a correlation of the MAE and the maximum error can be identified. The correlation coefficient between the MAE and the maxErr is higher than 0.96 for all surrogate models. This points out that those metrics are highly correlated and no extreme outliers are to be expected when the MAE is low. On the other hand, the nMAE behaves differently due to the scaling.

In Figure 8, a box plot for each error metric is displayed. The box covers the range from the first to the third quartile. A one and a half interquartile range was used for the whiskers. In each of the subplots, a box plot is drawn for the best surrogate model based upon the previously introduced DoCEs. The plotted mean values correspond to the ones given in Figure 6. As additional information, the median values are shown. When the median of the prediction errors is considered, the same order of surrogate models, as already shown for the mean values, is obtained. Again, the Halton+TPS model performs best overall. For all metrics but the MAE of the LHD+Cubic model, the median is lower compared to the mean. For the LHD+Cubic and the Halton+TPS models, an outlier in terms of MAE and maxAE can be identified. Moreover, the box plot of the Halton+TPS model displays an outlier in the nMAE. The highest range in all metrics can be observed for the fractional design, whereas the lowest range is produced by the LHD+Cubic model in all metrics.

Overall, the fractional factorial-based surrogate model achieves the minimum MAE in only two operating points. The LHD- and Halton-based models perform best at 10 and 6 operating points, respectively, when considering MAE. The same number of best results is obtained when considering the maximum error and the nMAE. Moreover, it can be seen in Figure 7 that the model based on the Halton design produces the highest deviation in the 18th operating point, but the nMAE stays below 2%. Even the surrogate model using the fractional factorial design achieves an nMAE lower than 3% in all validation points. Therefore, a fractional factorial-based surrogate model may be a good choice in the early stages of an investigation or if a sampling plan has to obey additional constraints, as it might be the case if a physical experiment is set up.

In Figure 9, the results of the surrogate model are displayed when only eight points are available for the POD-basis generation. The purpose of this investigation is to analyse the effect of varying available results, which is crucial when the proposed method should be employed in early stages of a parameter study. A clear trend can be seen in Figure 9: the more results are available, the more accurate the surrogate solution becomes. The nMAE averaged for the 18 validation cases, and only 8 design points available, is ${\mathrm{nMAE}}_{\mathrm{Halton},\mathrm{n}8}=0.013$. This result is similar to the results produced by the fractional factorial design-based surrogate models.

Figure 10 shows a comparison of the water mass fractions predicted by the FOM and the surrogate model on a cutting plane in the middle of the hollow fibre module for cases 16 and 18 of the validation data. The shown results correspond to the validations points in Table A1, where the mean absolute errors between FOM and surrogate models reached the minimum and maximum of the 18 points investigated. Overall, a high similarity between the results of the FOM and surrogate model can be observed from the contours. For both cases, the highest water mass fraction occurs in the fibres. Therefore, the fibres are well visible in Figure 10 in all four simulation results. The random placement of the fibres results in higher fibre density in the upper centre part of the humidifier, whereas the lower part is less densely packed. Due to the laminar flow, only a small amount of mixing due to convection occurs in the shell. If the packing density of the fibres within a certain region is high and the convective mixing is low, a high wet volume flow is present, opposed to a small dry volume flow, in this region. As a result, the water mass fraction in the shell along the flow length increases faster in regions with high packing density than in regions with low packing density. This behaviour can be observed in Figure 10 and is also present in FOM and surrogate results.

The lowest values of water mass fraction occur in the lower section of the contour shown. This can again be seen in both predictions. In the lower region of the humidifier, the distance between individual fibres is higher than in the upper section and, therefore, the water mass fraction in the shell is lower than in the densely packed regions. Both models predict the highest values of water mass fraction in the region where the fibre density is the highest. Deviations in the contours produced by the FOM and the surrogate model can be seen in places where a transition between contour levels in the high-fidelity model takes place. This applies to both fibres and the shell, with the greatest deviations occurring in the shell. For case 18, the highest deviations occur in the bottom region of the shell. The water mass fraction in fibres is well reproduced by the surrogate model.

Overall, it can be concluded that the surrogate model reproduces the results of the FOM with high accuracy, even at the parameter combination with the highest deviations. Our results suggest that the best combination of interpolation method and DoCE is the Halton design and TPS interpolation. However, the differences between the Halton and Latin-Hypercube DoCEs are relatively small. The same applies to the interpolation methods investigated in this study.

4.3. Using POD Solutions to Initialise the CFD Model

Even though the predictions of the surrogate models are as accurate as shown above, there might be situations where a full-order solution is desired. For this case, a POD+I solution can be used to accelerate the convergence of the CFD simulation. To show the benefit of the POD+I initialisation in terms of computational time saving, we use the Halton+TPS model to calculate the initial solutions and compare the computational time necessary for reaching convergence to a normal solution process.

For the comparison of computational times to reach converged simulation, the validation data set was used again. In the first step, the simulations are done using a uniform initialisation of the model. The water mass fractions in the whole field are set equal to the dry inlet, the velocity is set to zero and the pressure equals the outlet conditions.

To demonstrate the advantage of the POD+I initialisation compared to the uniform one, the two cases with the lowest, Case 16, and highest, Case 18, nMAE from the validation data set are considered. For both cases, each initialisation method is applied and the convergence behaviour with respect to the overall water transfer from wet to dry side is analysed. The residuals were checked to be below $1\times {10}^{-4}$ for each quantity. Additionally, the water transfer is considered as the most important integrated quantity and, therefore, is analysed in more detail. The water transfer rate is considered converged if the relative change between two iterations is below $1\times {10}^{-6}$ for 1000 iterations. In Figure 11, the overall water transfer is plotted against the iterations. A zoom is provided to show the convergence behaviour of the POD+I-initialised simulation. The black dotted line marks the steady-state result and the grey area a relative deviation of 1 % from the steady state. It can be seen that the POD+I-initialised cases converge after fewer iterations for both cases. For case 16, the very low deviations between the POD+I solution and the full-order model result in a fast convergence after just 437 iterations and 1770 s of computing time on 8 cores. The normally initialised solution takes about 14,002 iterations and 22,602 s on 8 cores to converge. In the same figure, in the bottom plot, again, one curve for each of the initialisation methods is shown for case 18. In the zoom box, it can be seen that the POD+I-initialised solution converges after 5522 iterations. The time needed to reach convergence on 8 cores is 9015 s. In contrast, the uniformly initialised solution needs more than twice the iterations, about 20,728, and 35,489 s of computational time to converge.

The mean simulation time of the FOM averaged over the 18 validation points is 12.6 h per simulation to solve for 30,000 iterations using 8 CPU cores of an AMD EPYC 7302. This high amount of iterations was chosen to definitely achieve convergence. The necessary time to calculate the POD+I solution is much lower than the time required to solve the FOM. The calculation of the 18 validation points is done in 107 s. The average prediction time for a single operating point is 5.95 s. This time already includes the writing of the OpenFOAM files on the disk. In fact, the isolated evaluation time for POD+I lies below one second. The average solution time for a POD+I-initialised solution is 1.8 h. In summary, using POD+I solutions for initialisation can save significant computing time in case of the membrane humidifier.

5. Conclusions and Outlook

In this paper, it was demonstrated that accurate and fast surrogate models of hollow fibre membrane humidifiers can be set up by using proper orthogonal decomposition in combination with an interpolation scheme. The presented surrogate models are able to accurately predict the water mass fraction at previously untried parameter combinations. It has been shown that parameter variations of the humidifier operating conditions can be significantly accelerated by using the proposed methods.

In order to obtain a database for the construction of the surrogate models, a high-fidelity 3D model of a hollow fibre humidifier was set up using OpenFOAM. In a first step, we presented a new set of measurement data of a commercially available hollow fibre humidifier. The measurement data were used to calibrate the vapour permeability of the membrane used in our high-fidelity model and to validate the same model in the following step. The developed high-fidelity model accurately predicts the water transfer from the wet to the dry side of the hollow fibre membrane humidifier. In the following, a POD-based surrogate modelling approach was proposed and applied to accelerate the parameter studies while preserving the three-dimensional character of the high-fidelity model. The best-suited DoCE for the presented problem set-up is the Halton design to generate the samples in combination with the cubic TPS interpolator as the approximation scheme for the water mass fraction distribution. However, it can be concluded that all tested DoCEs in combination with the interpolators performed well in terms of prediction accuracy. Our research shows that even with this small data set, a precise 3D surrogate model of the membrane humidifier can be generated. The maximum deviation of the best surrogate model configuration, a Halton DoCE combined with TPS interpolation, compared to the detailed CFD model is 1.1 g Water/kg Air. With the same configuration, the mean absolute error was 0.38 g Water/kg Air. The mean simulation time of this surrogate model is 6.0 s compared to 12.6 h necessary to calculate a full-order CFD solution for the same operating point. In a further investigation, it was demonstrated that the POD+I-based surrogate results can be used effectively as an initial solution for the full-order model, when the most accurate results of the high-fidelity model are required. In summary, our research demonstrates that fast and accurate surrogate models of membrane humidifiers can be obtained using POD+I. Further investigations should focus on the prediction of different humidifier geometries using surrogate models in order to support faster design processes of membrane humidifiers. Another field of application of the presented combination of DoCE and POD+I methods can be a PEM fuel cell that involves additional equations, making the solution process harder than for the humidifier.