Numerical Simulation of Dynamic Properties of Thin-Film Temperature Sensor Inside PEMFC

Abstract

The dynamic properties of thin-film temperature sensors with different sizes are investigated in detail through numerical simulation and system identification modeling. A one-dimensional transient heat transfer model for the sensor is built based on its location in the proton exchange membrane fuel cell （PEMFC）. The dynamic mathematical model， dynamic performance indicators， and dynamic error are obtained by employing COMSOL simulation and the system identification method. Notably， several significant dynamic parameters including working frequency bands， delay time， rise time as well as dynamic error peak， are determined for insulation layers of 1 μm， 2 μm， 3 μm， 5 μm， and 10 μm thick， and a real thin-film sensor is fabricated and calibrated. The results demonstrate that the sensor dynamic performance reduces with the growth of the insulation layer thickness. This paper reports a novel method to identify whether a thermal probe can capture the internal dynamic temperature variety of PEMFC， thus benefiting the further development of thermal probe on the research for PEMFC dynamic temperature variation under transient conditions， which is likely to inspire the sensor design contained physical parameters selection and structural design.

Keywords

proton exchange membrane fuel cell (PEMFC); temperature sensors; thin-film; dynamic properties; numerical simulation

Proton exchange membrane fuel cell （PEMFC） has become a promising power source due to its high energy conversion efficiency and eco-friendly production^［1］. Since temperature profoundly affects the output power and durability of the fuel cell， it is extremely important to detect the temperature inside the device. The thin-film temperature sensor is widely used for temperature measurement in narrow spaces inside PEMFC thanks to its small size and fast response^［2-8］. Nevertheless， whether the dynamic performance of the sensor meets the requirements for dynamic temperature measurement inside PEMFC is usually neglected. For vehicle-mounted fuel cells， with rapid changes of vehicle operating conditions， the chemical reaction rate of fuel cells constantly varies， thus leading to a violent and fast temperature fluctuation. The poor dynamic properties of the thin-film sensor would cause a great dynamic error compared to the true value when fuel cell temperature changes quickly. Consequently， dynamic characterization of the thin-film temperature sensor in PEMFC is necessary.

Numerical simulations have also been applied in the study of the dynamic properties of temperature sensors^［9-10］. Through the establishment of heat transfer models and the use of simulation methods， complex temperature fields and excitation can be built virtually， thus making it easier to study the dynamic properties of temperature sensors in a relatively feasible way. For dynamic characterization of temperature sensors， the dynamic performance indicators including time and frequency domain are usually presented for easier comparison and evaluation. Besides， dynamic mathematical models of temperature sensors are established， making it easier to acquire performance indexes and calculate the predicted response towards excitation^［11-12］.

Accordingly， numerical simulation is adopted for dynamic characterizations of the thin-film temperature thermal resistor in this paper. Besides， to explore the effects of design sizes on the dynamic behaviors of the sensor， the dynamic performance indexes of various design schemes are investigated. First， the construction and implementation of the heat transfer model for the resistance inside PEMFC are completed. Next， the forms and reciprocal transformations of different dynamic mathematical models are presented. Afterwards， resolutions and discussions are applied to transfer functions， dynamic performance indicators， and especially dynamic errors of the sensor. Finally， the effects of design sizes on the dynamic behaviors are also compared.

1 One-dimensional transient heat transfer model

1.1　Construction of heat transfer model

A one-dimensional transient heat transfer model of the thin-film thermal resistor inside PEMFC is established， whose basic structure is shown in Fig.1a， wherein BP is the bipolar plate， GDB is the diffusion layer backing， MPL is the microporous layer， ACL is the anode catalyst layer， PEM is the proton exchange membrane， and CCL is the cathode catalyst layer. The thin-film resistance is inserted inside the PEMFC， between the ACL and PEM. Fig.1b is a schematic diagram containing a thin-film resistance， where IL and Probe respectively represent the outer insulating layer and the metal probe in the center of the sensor. IL and Probe together form the thin-film thermal resistor. The material of the metal probe is platinum （Pt）， and the insulating layer is polymide （PI） coating. Fig.1c is a simplified one-dimensional model diagram from Fig.1b.

CCL is the most important heat source in PEMFC ^［13］. Therefore， this model ignores the reaction heat of ACL， ohmic heat， and phase change heat， only to consider the heat transfer of CCL. The heat transfer from CCL to PEM occurs at the border on the right. In this paper， this boundary is regarded as the first kind of boundary condition， i.e.， when the time $t>\mathrm{0}$， the boundary temperature $T_w=f(t)$. The case on the left is different， because of the special structure of BP， this boundary condition is the third kind of boundary condition in the gas flow channel， and there occurs convection heat transfer between hydrogen， water， and the wall， i.e.，

$$-\lambda {(\frac{\mathrm{d}T}{\mathrm{d}x})}_w=h(T_w-T_f)$$

（1）

wherein： $T_w$ is the wall temperature； $T_f$ and $h$ are the temperature and convection heat transfer coefficient of hydrogen. However， at the ridge of the bipolar plate， it is not the third type of boundary condition， and it is more difficult to obtain the temperature data here. To simplify the calculation， the established model is located at the flow channel， i.e.， the left boundary adopts the third type of boundary conditions. Assuming that the hydrogen temperature on the left boundary remains constant， the specific heat capacity， thermal conductivity， convective heat transfer coefficient， and other parameters used in the model do not change with temperature. The structural and physical parameters used in the model are shown in Tab.1， and the parameters are from the literature ［14-22］.

  Tab.1　Parameters of the model

Component	Thermal conductivity/ （W·m-1·K-1）	Specific heat capacity/（J·kg-1·K-1）	Density/（kg·m-3）	Thickness/μm	CHTC/ （W·m-2·K-1）
PEM	0.177	1050	2 076.2	25.4	—
ACL	0.061±0.006	3300	473.3	7.5（0.2 mg·cm-2）	—
GDB	0.300	568	2 786.4	140.0	—
MPL	0.150	568	3 485.9	60.0	—
IL	0.084	712	1 289.0	1， 2， 3， 5， 10.0
Pt	70.000	130	21 460.0	0.1	—
Hydrogen	—	—	—	—	2×103
CHTC： Convective heat transfer coefficient.

1.2　Implementation of heat transfer model

After the construction of the one-dimensional transient heat transfer model， COMSOL software is used to simulate the model. A temperature excitation is applied to the bounder on the right， and the insulation layer thickness of the resistor is changed to simultaneously achieve the temperature response of the probe at the same temperature excitation and with different sizes. The average temperature of all grid points in the metal probe is used as the probe temperature. The thermal resistor of Pt is relatively small， and the metal probe is very thin， therefore the temperature difference within the probe is small at the same time. It is reasonable to take the average temperature of all grid points as the probe temperature.

To get an accurate dynamic mathematical model of the sensor， the appropriate excitation signal must be available， which owns a wide frequency range and enables the inspiration for all modes of the sensor^［23］. The ideal step signal has a zero rise time and rich frequency components. The ideal step signal is used as the excitation signal here.

2 Results and discussion

2.1　Identification of sensor transfer function

Fig.2a shows the response of different insulating layer thicknesses to the same step signal. It can be found that as the thickness of the insulating layer increases， the sensor responds to the same step signal more slowly and the final stable value is also smaller. The reason for this is that， as the insulating layer thickness increases， the distance and hindrance of heat transfer increase， and the temperature rise decreases. With the growth of thermal resistance， a larger temperature gradient is caused when the stable stage is reached， making the final stability value decreases as the thickness of the insulating layer increases.

  

  

Fig.2  Response of different insulating layer thicknesses to the same step signal， and transfer functions and model adaptation rates at different insulation thicknesses

a Step response of sensors at different insulation thicknesses b Amplitude frequency characteristics of sensors at

According to the input/output data， the System Identification Toolbox is used to achieve the continuous transfer function of the sensor. Tab.2 depicts the transfer functions and model adaptation rates at different insulation thicknesses. The amplitude frequency characteristic curves are drawn through the transfer function， as demonstrated in Fig.2b.

Tab.2  Transfer functions of sensors ai different insulation thicknesses

Insulation thickness/μm	Continuous Transfer function	Model adaptation rate/%
1	$H(s)=\frac{1.46\times {10}^4{s}^4+6.371\times {10}^7{s}^3+3.228\times {10}^{10}{s}^2+2.73\times {10}^{12}s+2.709\times {10}^{13}}{s^5+2.849\times {10}^4{s}^4+8.399\times {10}^7{s}^3+3.61\times {10}^{10}{s}^2+2.841\times {10}^{12}s+2.739\times {10}^{13}}$	99.39
2	$H(s)=\frac{7775s^4+2.932\times {10}^7{s}^3+1.322\times {10}^{10}{s}^2+1.027\times {10}^{12}s+9.858\times {10}^{12}}{s^5+2.052\times {10}^4{s}^4+4.63\times {10}^7{s}^3+1.615\times {10}^{10}{s}^2+1.108\times {10}^{12}s+1.008\times {10}^{13}}$	99.48
3	$H(s)=\frac{3897s^4+1.077\times {10}^7{s}^3+4.098\times {10}^9{s}^2+2.864\times {10}^{11}s+2.651\times {10}^{12}}{s^5+1.183\times {10}^4{s}^4+1.911\times {10}^7{s}^3+5.365\times {10}^9{s}^2+3.198\times {10}^{11}s+2.739\times {10}^{12}}$	99.72
5	$H(s)=\frac{1301s^3+7.383\times {10}^5{s}^2+6.224\times {10}^{7}s+6.015\times {10}^8}{s^4+3413s^3+1.17\times {10}^6{s}^2+7.465\times {10}^{7}s+6.344\times {10}^8}$	99.02
10	$H(s)=\frac{346.6s^3+9.399\times {10}^4{s}^2+4.934\times {10}^{6}s+4.063\times {10}^7}{s^4+1100s^3+1.731\times {10}^5{s}^2+6.69\times {10}^{6}s+4.498\times {10}^7}$	98.96

2.2　Model test

The mathematical models need to be tested to determine the applicability. By changing the step excitation， different excitation signals can be obtained， and the response towards excitations can be calculated in COMSOL as the true value. Next， the Tustin transform and Z transform are applied to the continuous transfer function for achieving Difference Equations. The output response for the excitation signal can be computed utilizing Difference Equation. Comparing the simulation results from COMSOL with model results from the difference equation， the model reliability can be verified. Fig. 3 manifests the excitation signal and two responses of the model test， as well as the error （absolute value of the difference between two responses）. It reveals that the response calculated by Difference Equation is quite consistent with that by COMSOL.

a  1 μm

b  2 μm

c  3 μm

d  5 μm

e  10 μm

Fig.3  Simulation output, Difference Equation output, and error(absolute value of the difference between the two outputs) for unit step signal of sensors at different insulation thicknesses

Notably， in an extremely short time （0.1‒0.2 ms） after the start of excitation， the error is large， even when the insulation layer thickness is 1 μm， the maximum error reaches 0.33 （the simulation result is greater than Difference Equation result）， and the maximum error in the initial stage decreases with the rise of the insulation layer thickness. Since transient heat transfer is divided into the initial irregular regime and the subsequent regular regime. The heat transfer properties of the two stages are different. The former is mainly affected by the initial temperature， while the latter is mainly affected by the boundary conditions， which leads to different structures of the temperature data in the two stages. Consequently， the one mathematical model cannot completely match the data of the two stages. As can be seen from Fig. 3， the results from Difference Equation and COMSOL agree well in the regular regime， but differ in the initial irregular regime， which leads to a large initial error. In general， the output from the Difference Equation model has a good agreement with the simulation results， which is applicable to study the dynamic properties.

2.3　Calculation of dynamic performance indexes

Tab.3 lists the frequency and time domain indicators containing working frequency bands whose amplitudes errors are less than ±5% and ±10%， the delay time， and the rise time. The delay time refers to the time at which the output value reaches 50% of the steady-state value from zero， and the rise time refers to the time interval at which the output value reaches 90% of the steady-state value from 10%. It can be found that the working frequency band of the sensor reduces gradually with the increase of the thickness of the insulating layer. When the thickness of the insulating layer is 1 μm， the working frequency bands whose amplitudes are less than ±5% and ±10% are 249.04 Hz and 85.91 Hz respectively. The response time has a positive correlation with the insulation thickness. When the insulating layer thickness is 1 μm， the delay time is 48 μs and the rise time is 3.04 ms.

Tab.3  Frequency and time domain indicators of sensors at different insulation thicknesses

Insulation thickness/μm	Frequency domain/Hz		Time domain/ms
	±10%	±5%	delay time	rise time
1	249.04	85.91	0.048	3.040
2	88.42	31.32	0.150	8.270
3	49.80	17.04	0.390	14.760
5	23.42	9.89	1.080	30.178
10	10.09	5.68	3.880	71.550

Comparing the frequency domain and time domain dynamic performance indicators with the designed indexes， the design scheme can be evaluated and optimized. But sometimes there are no dynamic performance indexes. Only the size of the dynamic error is listed. Therefore， the dynamic error is discussed in Section 2.4.

2.4　Dynamic error calculation

The dynamic error refers to the difference between the output and the input caused by the thermal inertia during dynamic measurements^［24］. To ensure that the sensor can meet the requirements of dynamic measurement， and the dynamic error is within the allowable range， it is necessary to determine the dynamic error. To solve the dynamic error， the input signal is required. The excitation temperature signal used in this paper is derived from Ref.［25］. The signal is the temperature change of the cathode catalyst layer in PEMFC， as exhibited in Fig. 4a. Taking the temperature signal as the input， the output of the sensor is calculated in COMSOL， and the difference（absolute value） between the input and output is the dynamic error. Fig. 4b and 4c respectively shows the excitation and COMSOL simulation output， as well as the evolution of dynamic error with time， i.e.， the difference （absolute value） between the temperature excitation and the simulation output. It can be found from Fig. 4b that， for the sensor with the same insulating layer thickness， the dynamic error first grows to a peak value， and then gradually reduces to a steady value. When the insulating layer thickness is 1 μm， the maximum dynamic error is 0.094 K， and the final stable value is 0.08 K. The dynamic error grows with the insulation layer thickness rising. The evolution of dynamic error is caused by the variety of temperature signal change rates. The change rate of temperature is shown in Fig. 4c， and the calculation method is

$$\Delta =\frac{T_k-T_{k-\mathrm{1}}}{T_{k-\mathrm{1}}}$$

（2）

a  Dynamic temperature variety of cathode catalyst layer in PEMFC

  

  

Fig.4  Dynamic temperature variety, temperature excitation-simulation, and the differences between them

b Temperature excitation and COMSOL simulation output of c Difference (absolute value) between excitation and simulation output of sensors at different insulation thicknesses, temperature change rate of the excitation, and an insulating layer of 10 μm thick sensors at different insulation thicknesses

wherein k is the sampling point， k=2，3，…，N （N is the total number of sampling points）. Fig.4 shows the change rate of temperature excitation and 10 μm （thickness of insulating layer） sensor simulation output. Due to the thermal inertia， the temperature change rate of the response is delayed compared with that of the excitation. It can be seen from Fig.4c that the change rate of the excitation before the peak time is greater than that of the response of the 10 μm sensor， i.e.， the growth of the excitation is faster than that of the response. Consequently， before the peak time， the dynamic error must be increasing； after the peak time， the response temperature change rate is greater than the excitation， so the dynamic error starts to decrease from the maximum. Finally， the temperature change rate of excitation and response tends to zero， and the dynamic error also tends to a fixed value. When it tends to be stable， the thicker the insulating layer， the greater the thermal resistance. According to the Fourier law of heat conduction， the temperature gradient between the right boundary and the sensor is also greater， so the error after stabilization is also larger.

Moreover， the peak time of dynamic error for sensors with different thicknesses is the same in general， but there are some differences. The peak time of 1， 2， 3， 5 10 μm is 0.794， 0.797， 0.801， 0.806， 0.817 s respectively. With the growth of thickness， the peak time has a few milliseconds delay. It is because the thicker the insulation layer is， the greater the thermal inertia is， so the later the response change rate reaches the same time as the excitation change rate， that is， the later the peak time is. The dynamic error in Fig.4c clearly illustrates the dynamic error of the sensor with different designed sizes during virtual measurement. By comparing with the design target， it is very convenient to evaluate and optimize the design scheme.

Except for numerical simulation， detailed experimental plans and preliminary preparation has been completed. Fig.5a-5c shows the design layout， partial details， and a photograph of the self-designed sensor made by magnetron sputtering. The calibration curve of the sensor can be seen in Fig.5d. And R² value is 0.995 69， which demonstrates good linearity. Further experimental dynamic calibration and validation will be attempted in the later work.

Fig.5  Design layout， partial details， and a photograph of the self-designed sensor made by magnetron sputtering

3 Conclusions

In this paper， a numerical simulation is conducted to study the dynamic properties of sensors with different design sizes. The following conclusions are reached：

（1） The transfer functions of sensors with different sizes are obtained and tested. It is found that the model is in good agreement with the simulation results， and the model is applicable.

（2） With the growth of the insulating layer thickness， the dynamic performance of the sensor worsens. When the thickness of the insulating layer is 1 μm， the working frequency bands with an amplitude error of less than ± 5% and ± 10% are 249.04 Hz and 85.91 Hz respectively， the delay time is 48 μs， and the rise time is 3.04 ms.

（3） As the thickness of the insulation layer increases， the sensor response to the excitation slows down and the dynamic error enlarges. Due to the influence of the change rate of the measured temperature signal， for the same insulation thickness， the dynamic error first increases to a peak and subsequently decreases to a steady value. When the insulation thickness is 1 μm， the dynamic error peak value is 0.094 K， which finally stabilizes to 0.08 K.

Different from the experimental calibration， this paper proposes an approach to study the sensor dynamic properties through simulation computing. Noticeably， dynamic properties with temperature sensors of different sizes aimed at detecting thermal variety inside PEMFC are evaluated. Additionally， it proposes a practicable way for sensor design from the perspective of dynamic property research.