Simulation Research on the Performance of PEMEC with Single Channel Flow Field Established by Orthogonally Sinusoidal Function with Different Amplitudes and Angular Frequencies

. As one of the main hydrogen production devices, the performance of the proton exchange membrane electrolytic cell (PEMEC) is directly related to the geometric design of the ﬂ ow channel. In this paper, a new orthogonally sinusoidal channel with periodic ﬂ uctuations in both horizontal and vertical planes is proposed. The e ﬀ ects of amplitude and angular frequency in two orthogonally sinusoidal channels with the same or opposite frequency on PEMEC performance are studied by numerical simulation. The results show that the spatial amplitude and angular frequency e ﬀ ectively improve the electrochemical performance and mass transfer performance of PEMEC. Compared with Case 1-1, the polarization performance of Case 7-2 is increased by 16.39% because of the largest spatial amplitude. The oxygen content at 40mm away from the inlet in the di ﬀ usion layer is also decreased by 42.63%. The horizontal angular frequency can increase the hydrogen content by 11.07% at most in Case 5-2. However, the increase in vertical angle frequency will inhibit the oxygen removal on the anode side. The sinusoidal ﬂ ow channel in orthogonal spaces will provide a particular reference for the design of PEMEC in the future.


Introduction
In recent years, since fossil energy consumption makes the problem of energy shortage increasingly prominent [1,2], hydrogen energy has begun to receive much attention as a representative of clean energy [3,4].Among many technologies for producing hydrogen as a stable carrier for converting [5][6][7] and storing excess electric energy [8,9] by electrolysis of water, proton exchange membrane electrolysis cell (PEMEC) is considered a promising technology due to its advantages such as small volume of reaction device [10], high purity of hydrogen production [11], and strong adaptability [12].However, the acidic environment in the PEMEC imposes more stringent requirements on the performance of key components such as proton exchange membranes, electrodes, and bipolar plates [13][14][15].Among these components, the flow channel in the bipolar plate is the most critical structural element in the electrolytic cell [16,17].
As a channel for transporting liquid water and discharging product gas, the geometric design of the flow channel determines the spatial distribution of materials [18,19].Gayen et al. [20] placed baffles at the inlet of the microchannel to obtain uniform material spatial distribution and heat transfer.This spatial distribution can be divided into two main kinds: one is in the horizontal plane and the other one is in the vertical plane.In recent decades, researchers have conducted a series of studies on these two main flow behaviors in channels.
The researchers first investigated the effects of periodicities and fluctuations on PEMEC performance in the vertical plane.The vertical fluctuations can not only maintain the electrolytic cell operating environment stable but also make the water transfer along the vertical plane with quick removal of produced bubbles by forced convection to improve the life and mass transfer capacity of PEMEC [21].Since the geometric design of the vertical plane can determine levels of forced convection, Bilgili et al. [22] also set periodic blocks in the vertical plane of the flow channel to cause forced convection.Khatib et al. [23] paid attention to the vertical planes and studied the influence of the fluctuant geometrical parameters in PEMEC with open-cell foam material on forced convection of the fluid.Therefore, the geometrical parameters of the fluctuations should be considered a critical factor in the design of the flow channel.
Besides the geometrical parameters of vertical fluctuations, the same vertically periodic structures were also investigated by Wang et al. [24] with raised periodic structures in an interdigitated flow channel.The presence of periodic structures in the flow channel has been observed to enhance the performance of PEMEC.Wang et al. [25] designed the flow channel with a periodic wave structure in the vertical direction.Compared with the conventional straight flow channel, the maximum power density of the periodic wave structure channel was increased by 4%.Tugirumubano et al. [26] discovered that an increased frequency resulted in more vertical liquid to flush and exhaust bubbles through the sinking periodic structure to improve electrolysis efficiency significantly.Though the necessity of periodic fluctuations in the vertical plane to the performance of PEMEC is proved in the above studies, the attention to these periodic fluctuations is only focused on the vertical plane.Due to the lack of horizontal planes, these studies on periodic fluctuations are insufficient.
Since the design of a horizontal structure also plays a crucial role in guiding fluid fluctuation in the horizontal plane, researchers further investigated the effects of periodicities and fluctuations on PEMEC performance in the horizontal plane.Toghyani et al. [27] found that compared with the parallel flow channel, the more extended serpentine flow channel with a number of corners exhibited horizontal fluctuations to increase current density and hydrogen production efficiency in PEMEC by promoting fluid rotation.In addition, Toghyani et al. [28] also further designed the channel with a spiral structure to enhance horizontal fluctuations.Zhou et al. [29] focused on the influence of structural parameters on PEMFC performance in horizontal fluctuation structure.The results showed that the geometric parameters near the crest of the fluctuation structure greatly influenced the performance of the flow channel, and the performance of the flow channel with the same slope was improved by 19.53%.
At the same time, the influences of the period on PEMEC performance were also found by Majasan et al. [30] in some periodic channels.Based on the electrochemical performance results, the fluid residence time was enhanced in the serpentine channel with periodic structures in the horizontal plane to show a superior electrolytic performance.Similarly, Li et al. [31] also studied the influence of horizontally periodic arrangement on the electrolytic performance of PEMEC.The results showed that a better electrolytic performance was able to be obtained by more periodic structures.But if these periodic structures were too many, the channel would be excessively elongated and the performances of PMEMC would also be reduced by consequent bubble accumulation.Subsequently, with both the fluctuations and periodicities taken together into consideration in the horizontal flow channel, Tijani et al. [32] further found that the fluid slowed down around periodic fluctuations of the channel was more conducive to the interaction between the reactants and the MEA.The importance of the periodic fluctuation structure in the horizontal plane is also shown.At the same time, it should be pointed out that the current research on the fluctuations and periodicities in the horizontal plane is only limited to the horizontal plane, and the lack of periodic fluctuation structure in the vertical plane will also bring the insufficiencies of the research.
The researches on periodic fluctuations are mainly focused on the single horizontal [33,34] or vertical layout [35], while the orthogonal structures in both horizontal and vertical planes at the same time are rarely studied.Therefore, a periodic fluctuation structure can be introduced into the flow channel design of PEMEC.The effect of periodic fluctuations on PEMEC performance is studied by coupling two independent sinusoidal functions on both the horizontal and vertical planes.The sinusoidal function of both planes is shown in Figure 1.
In this study, a three-dimensional and steady-state model is established by coupling the orthogonal structure with both periodic fluctuations in the horizontal and vertical planes.The PEMEC model is constructed and solved utilizing the commercial software COMSOL Multiphysics 5.6.The amplitude and angular frequency of the sinusoidal function are taken into account as geometric variables of the channel in this study.The opposite fluctuations are also considered in the horizontal plane to enhance the degree of fluctuation variation.The effects of horizontally, vertically, and orthogonally sinusoidal channel designs on the electrochemical performance and mass transfer capacity of PEMEC are studied.The four parameters containing the amplitude and angular frequency of the sinusoidal function in both the horizontal and vertical planes are derived from mathematic expressions and discussed with simulation results.Based on these parameters, 24 cases in 8 groups of experiments

Model Description
2.1.Design of Orthogonally Sinusoidal Flow Channel.As the flow channel exhibits periodic fluctuations along its length, the minimum repetitive units of the orthogonally sinusoidal flow channel are selected as the research object in Figures 2(a) and 2(b).The width direction of the flow channel is along the x-axis direction, while the height direction of the flow channel is along the y-axis one.All the repetition units will be periodically repeated along the z-axis direction from z = 0 to 50 mm.With the left and top views of the orthogonally sinusoidal flow channel projected to the YOZ and the XOZ planes, respectively, the curves of the walls projected can be seen as sinusoidal (S) or linear (L) curves with independent variable z.Since both the fluctuations and periodicities are governed by the crucial geometric characteristics of these walls, V upper and V lower are used to represent the curves of upper and lower walls in the left view, while H left and H right stand for the ones of left and right walls in the top view.
The flow channel design is described based on four curves.Then, the sinusoidal curve controlled by amplitude (A) and angular frequency (ω) can be expressed in a more distinguishable way as follows: where A z 1 and ω j are the amplitude and the angular frequency of the left wall curve function, respectively.The height of the channel bottom is defined as zero.Limited by the width and height of the channel, the equations of the lower, upper, and right wall curve functions are as follows: where A z 1 and A z 2 are used to represent the amplitudes of the right and upper wall curve functions, respectively.ω i and ω k represent the angular frequencies of the right wall and upper wall curve functions, respectively.Considering that the number of equations is too complicated to describe the characteristic curve of the minimum repetitive unit, the difference between the curve equations of the right and left walls is H, and the difference between the curve equations of the upper and lower is V, these two differences can be explained by After substituting Eqs.(1) and (2) into Eqs.( 5) and (6), Eq. ( 7) is received as follows: Since the angular frequencies of the left and right walls curve in the horizontal plane can be the same or the opposite of each other, the same angular frequencies Front view Left view Front view Left view  In order to explore the influence of the amplitude (A) and angular frequency (ω) on the performance of singlechannel PEMEC, the channels with different geometrical parameters are divided into 8 groups and 24 cases.The characteristics of the channel are further described by some symbols.S and L are used to represent the sinusoidal and linear curves of the walls projected, respectively.In order to avoid ambiguity in the following discussion, the descriptions of the same or opposite angular frequency of the left and right wall curves are called the same or opposite frequency.Therefore, when the frequency of the left wall curve is opposite to that of the right wall curve, the left wall curve can be further named Ŝ.
The naming rule equation for groups is as follows: Group no of group = ω with different positions and sizes no of case where no of group can be 1, 2, 3, 4, 5, 6, 7, or 8; no of case can be 1, 2, 3, 4, or 5; V upper can be S or L; V lower are always L; H left can be Ŝ, S, or L; H right can be S or L; and ω with different positions and sizes can be 2ω i or 2ω k .The amplitude and angular frequency are considered as variables to study the fluctuation and period, respectively.The grouping of the main groups in the above naming is shown in Table 1, which can be combined with Eq. ( 7) for a detailed explanation.
(1) When A z 1 = 0 and A z 2 = 0, the horizontal and vertical planes are linear geometric features in Group 1 (2) When A z 1 = 0, A z 2 ≠ 0, and ω k ≠ 0, the flow channels have sinusoidal geometric features in the vertical plane, while in the horizontal plane, there are still linear geometric features in Group 2 (3) When A z 1 ≠ 0, A z 2 = 0, and ω i = ω j , the vertical plane of the flow channels is a linear geometric feature, and the horizontal plane is sinusoidally geometric features of the left and right walls with the same frequency in Group 3 (4) When A z 1 ≠ 0, A z 2 ≠ 0, and ω i = ω j = ω k , the angular frequencies of the flow channels in both vertical and horizontal planes are equal, and the left and right wall curve frequencies in the horizontal planes are the same in Group 4 (5) When A z 1 = A z 2 ≠ 0 and ω i = ω j ≠ ω k , the amplitudes in both vertical and horizontal planes are equal, and the left and right wall curve frequencies in the horizontal plane are the same.However, the angular fre-quency in the vertical plane is different from that in the horizontal ones in Group 5 (6) When A z 1 ≠ 0, A z 2 = 0, and ω i = −ω j , the vertical plane of the flow channels is a linear geometric feature, and the left and right wall curve frequencies in the horizontal plane are opposite to each other in Group 6 (7) When A z 1 ≠ 0, A z 2 ≠ 0, and ω i = −ω j = ω k , the angular frequencies of the flow channels in both the vertical and horizontal planes are the same, too.However, the left and right wall curve frequencies in the horizontal plane are opposite to each other in Group 7 (8) When A z 1 = A z 2 ≠ 0 and ω i = −ω j ≠ ω k , the amplitudes in both vertical and horizontal planes are the same, while the frequencies of the left and right wall curves are opposite in the horizontal plane.The angular frequency in the vertical plane is different from those in the horizontal ones in Group 8 Based on the groups above, the groups are further subdivided into 24 cases with amplitudes and angular frequencies in orthogonal space.The abbreviation for each case is named as follows: Case abbreviation = case no of group − no of case 10 The detailed naming and geometric parameters of the single channel in each case are shown in Table 2.

Physical Model.
Figure 3 shows a single channel in the PEMEC.The main structure from top to bottom includes an anode bipolar plate (ABP), anode flow channel (ACH), anode gas diffusion layer (AGDL), anode catalytic layer (ACL), proton exchange membrane (PEM), cathode catalytic layer (CCL), cathode gas diffusion layer (CGDL), cathode flow channel (CCH), and cathode bipolar plate (CBP).The geometric parameters and operating conditions of PEMEC are shown in Tables 3 and 4, respectively.
2.3.Numerical Model.Due to the complex chemical reactions and physical changes in PEMEC, the simulation study cannot be carried out individually in a short time.Therefore, with the ignoring of the nonessential factors, some assumptions of the simplified model are applied as follows [36,[41][42][43]: (1) The water in PEMEC is always the liquid state without phase transition during the reaction where i v,a and i v,c represent the local electrochemical reaction rates.a v,a and a v,c denote the electrochemically specific active surface areas.i 0,a and i 0,c represent the reference where φ s is the solid overpotential, φ m represents the electrolyte overpotential, and E eq denotes the equilibrium potential.The equilibrium potential is calculated by the Nernst equation [45] as follows: where E eq,ref T can be obtained from the equation below [46].
where T represents the operating temperature of the electrolytic cell, R is the universal gas constant, F denotes the Faraday constant, and p i (i stands for O 2 , H 2 , and H 2 O) represents the respective equilibrium pressure of each substance.

Proton and Electron Transport
Equations.The proton and charge transport equations are as follows [37]: where σ s represents the electron conductivity of the solid phase electrode, σ m is the proton conductivity of the proton exchange membrane, and S s and S m refer to the source terms of solid overpotential and electrolyte overpotential, respectively.
The proton conductivity of a proton exchange membrane is a function of water content and temperature [47].
Conservation of Mass and Momentum.The continuity equation is as follows [48]: where ρ l is the density of liquid water, u l represents the velocity vector of liquid water, and S l denotes the source term of liquid-phase mass conservation.
∇ ρ g u g = S g , 18  where ρ g is the density of the gas, u g represents the velocity vector of the gas, and S g denotes the mass source term of a gas produced by a chemical reaction.
The momentum conservation equation is as follows [49]: where κ is the permeability of porous media, ε p denotes the porosity of porous media, I denotes the transposition symbol, p l is the liquid water pressure, and μ l represents the dynamic viscosity of liquid water.
where κ is the permeability of porous media, ε p denotes the porosity of porous media, I denotes the transposition sym-bol, p g is the gas pressure, and μ g represents the dynamic viscosity of the gas.

Transport of Energy.
The energy equation in all computational domains is defined as follows [50]: where ρ eff , C p,eff , and k eff stand for the effective density, the effective heat capacity, and the effective thermal conductivity, respectively, which can be calculated as follows [51]: where ρ s , C p,s , and k s are the density, heat capacity, and thermal conductivity of the solid matrix, respectively.ρ f , C p,f , and k f are the density, heat capacity, and thermal conductivity of the fluid, respectively.S T denotes the source term of the energy equation, which includes the irreversible reaction heat and entropic heat generated during the electrolytic reaction, as well as ohmic heat generated in the electrode and membrane.Table 5 shows the source terms in the governing equation.

Boundary Conditions.
In this study, the PEMEC channel is operated at a temperature of 353.15 K, and the working pressure of the outlet is set to 1 atm.Only liquid water is introduced into the inlets of both the anode and cathode channels, Cathode reference specific surface area, a v,c (1/m) 1 × 10 6 Thermal conductivity of the membrane, Heat capacity of the membrane, C p,m (J•kg -1 •K -1 ) 1090 8 International Journal of Energy Research which are located on the same side with an inlet velocity of 0.125 m/s.The liquid water in the anode channel serves as the reactant and is used to promote the removal of oxygen generated, while the liquid water in the cathode channel is only used to remove hydrogen.The bipolar plate of the anode side is set to 2.2 V, and that of the cathode side is set to 0 V.

Numerical Computation
Procedure.In this study, the numerical model is built in COMSOL Multiphysics 5.6 with free and porous medium flow modules and electrolyzer modules to solve the problems of charge transfer and material transfer.The iterative process is completed when the residuals of free and porous media flow modules are less than 10 -3 , and those of the electrolyzer coupling with multiphysical fields are all less than 10 -4 .
2.6.Model Verification.In order to verify the accuracy of the numerical model, the polarization curve obtained by numerical simulation is compared with the experimental results received by Majasan et al. [30] in Figure 4.The maximum error between the experimental results and the simulated ones is 3.1%, which proves that the numerical model established is accurate and reliable.Due to the application of the finite element method, the quality and number of grids become necessary factors affecting calculation time and accuracy [52].When the voltage is 1.8 V, different grid numbers are selected for the simulation study, and the curves obtained are shown in Figure 5.When the number of grids is lower, the current density decreases from 446.8 mA/cm 2 to 446.67 mA/cm 2 with the grids increasing from 48500 to 106700.Unlike the lower number of grids, the current density will remain almost constant when the number is increased from 155200 to 242500.So the number of grids is selected as 155200.

Results and Discussions
3.1.Effect of Sinusoidal Curve in Vertical Plane on PEMEC Performance.The polarization curves and oxygen content curves of different Groups 1 and 2 are shown in Figures 6(a) and 6(b).When 2.0 V is selected as the standard for comparison, the polarization curves of Group 2 are better compared to Group 1 with the maximum 5.1% increment.Also, from Case 2-1, Case 2-2, and Case 2-3, the polarization performances are increased with the increasing amplitude in the vertical plane.Also, since the amplitude is too small in Case 2-1, the polarization performance of Case 2-1 in Group 2 is a little improved compared with Case 1-1.The oxygen content of Case 2-3 with the highest amplitude is reduced by about 22.3% at 40 mm away from the inlet compared with Case 1-1.Furthermore, oxygen in the diffusion layer  Figure 7(a) shows the flow channel cross-sections of Groups 1 and 2 in the YOZ plane.The overall effect of the local current density distribution in Figure 7(b) is the same as that of the polarization curve.Interestingly, it can be seen from Figure 7(b) that the current density of Group 2 is significantly higher than that of Group 1 near the inlet.The increase of outlet current density is not apparent, which can be attributed to the fact that the design of the vertically sinusoidal curve improves the efficiency at the initial stage of the electrochemical reaction.However, due to the continuous consumption of liquid water in the process of passing through the flow channel, the effect of amplitude on improving current density at the outlet is weakened.In the vertical plane, appropriately increasing the amplitude of the sinusoidal curve is conducive to enhancing the polarization performance of PEMEC.
In order to explore the effect of the vertically sinusoidal curve on PEMEC's mass transfer capacity, Figures 7(c  Figure 7(e) shows the influence of amplitude in the vertical plane on the vertical velocity.In Group 1, after the vertical velocity with 0.04 m/s is distributed at a small distance from the inlet, the vertical velocity with 0 m/s will be almost uniformly distributed along the rest of the channel.It can be proved that the flow in the noninlet region of the straight channel is parallel flow.High-speed regions positively correlated with amplitudes will appear near the trough of the vertical sinusoid in Group 2. The same finding can be obtained in the research proposed by Yan et al. [53].In Group 2, by increasing these vertical velocities up to 0.079 m/s towards GDL at most, the mass transfer process of liquid water involved in the electrochemical reaction is accelerated to produce more oxygen and hydrogen.Therefore, it is proved that the electrolytic result can be changed by changing the sinusoidal function parameters of the periodic fluctuation structure in the vertical plane.

Effect of Sinusoidal Curve in Horizontal
Plane on PEMEC Performance.To investigate the effect of the horizontally sinusoidal amplitude on the performance of PEMEC, the polarization curves of horizontally sinusoidal channels with the same and opposite frequencies are compared in Figure 8(a).Overall, the polarization performances of Groups 3 and 6 are better than those of Group 1 because of the horizontally sinusoidal curve.At 2.0 V, the polarization performance of Case 6-3 is the best in Group 3 and Group 6, which is improved by 14.1% compared with Case 1-1.Interestingly, the polarization curve of PEMEC in horizontally sinusoidal channels is also positively correlated to amplitude.Though the frequencies of Groups 3 and 6 are different, the polarization performances with the same amplitude are almost the same.
It can be seen from Figure 8(b) that the oxygen content in the AGDL is related to the structure of the horizontally sinusoidal curve.With increasing amplitude, the oxygen content in the horizontal structure with opposite frequency decreases more significantly than that with the same frequency.Among all of the cases with opposite frequency, Case 6-3 has the lowest oxygen content.This phenomenon may be due to the varying width of the flow channel.This channel facilitates the transformation of large bubbles into smaller ones to some extent, and the oxygen can be removed more efficiently to reduce its content.11 International Journal of Energy Research same and opposite frequencies at the ACL-PEM interface, respectively.In all the cases, with the increasing amplitude, the local current density at the inlet is increased to enhance the distributed area of high current density across the interface.It should be pointed out that the degree of changes caused by the horizontally sinusoidal curve is much higher than that caused by the vertically sinusoidal one.When the amplitude of the sinusoidal channel with the same frequency is 0.8 mm, the maximum local current density occurs at 0.1183 A/cm 2 .The amplitude changes the inlet local current density of the same frequency sinusoidal curve more significantly than that of the opposite frequency one.In general, Case 6-3 shows the best polarization performance among the horizontally sinusoidal PEMEC structures.Moreover, Case 3-3 exhibits higher local current density at the inlet.
Figures 10(a)-10(d) show the influence of the horizontally sinusoidal amplitude on the internal mass transfer of PEMEC.Since the transportation laws of both oxygen and water are similar, flow channels with opposite frequencies also have a stronger ability to supply water in Figures 10(a) and 10(b).In Figures 10(c) and 10(d), the variation of hydrogen content is closely related to the variation of polarization performance.No matter whether the frequencies of horizontally sinusoidal curves are the same or opposite, the hydrogen production capacities are improved at the same level.Among them, the hydrogen content of Case 6-3 is increased by 16.13% compared with Case 1-1.In general, compared with Case 1-1, the horizontally sinusoidal curve has a great influence on the mass transfer performance of PEMEC.
The horizontal velocity distribution employed by Amit et al. [54] is used to investigate the states of fluid flow.Since the structure is designed on the horizontal plane, the effect of the horizontally sinusoidal amplitude on the performance of PEMEC can also be studied by the horizontal velocity at the ACH-AGDL interface, as shown in Figures 11(a Comparing the distribution of horizontal velocities in Group 3 and Group 6, the high-speed regions positively correlated to amplitude will be formed near both the trough of the horizontal sinusoid.This conclusion is the same as that obtained by Zhou et al. [35].However, this region is more prominent in the same frequency channel than the opposite one.Among Group 3 with the same frequency, Case 3-3 has the largest amplitude and high-speed region, with a maximum velocity of 0.0368 m/s, which enhances the horizontal mass transfer capacity to accelerate reactant supply and product gas removal.

Effect of Orthogonally Sinusoidal
Curve in Space on PEMEC Performance.The polarization curves of orthogonally sinusoidal channels with the same and opposite frequencies are shown in Figures 12(a) and 12(b), respectively.In Figure 12(a), the larger amplitude is positively related to the polarization performance of PEMEC in both horizontal and vertical planes.Different from the amplitude, the influence of angular frequencies on the polarization performance is not specific.This conclusion can be drawn from the same polarization performances in Case 5-1 and Case 4-5.In addition, by comparing both Figures 12(a) and 12(b), it can be seen that similar laws of the impact of the amplitude and angular frequency on the polarization performance of PEMEC can be obtained in orthogonally sinusoidal channels with the same frequency or opposite frequency.Overall, Case 4-2 and Case 7-2 have     As one of the products, the oxygen content distribution will affect the mass transfer capacity of PEMEC.Figures 12(c) and 12(d) show the oxygen content curve of orthogonally sinusoidal structures with the same and opposite frequency at the center of GDL along the length, respec-tively.The laws that the oxygen content distributions in Group 4 and Group 7 are more deeply affected by the vertical amplitudes than the horizontal ones with the same frequency, while the horizontal ones with the opposite frequency seem more significant to the oxygen content distribution obtained.Since there are the largest vertical and horizontal amplitudes in Case 4-4 with the same frequency and Case 7-2 with the opposite frequency, respectively, the  In order to further investigate the effects of amplitude and angular frequency on local current density in PEMEC with orthogonally sinusoidal channels, both Figures 13(a Similar to the variation of oxygen content, the distribution of liquid water is also shown in Figures 14(a) and 14(b).Whether the frequencies are the same or opposite, the optimal distributions of liquid water content in horizontal planes will occur in the cases with the largest amplitudes, such as Case 4-2, Case 4-4, Case 7-2, and Case 7-4.The liquid water content of Case 7-2 surpasses that of the other three channels.Since the spatial amplitudes are the largest, the volumes of the channel are also the biggest.Then, both the transportation of the liquid water to GDL and the dispersion of the oxygen bubbles are enhanced.
The mole fraction distribution of hydrogen at the CCL-PEM interface is shown in Figures 14(c) and 14(d).From these figures, the amplitudes and angular frequencies in both horizontal and vertical planes are positively related to the mole fraction of hydrogen.In general, there will be a higher hydrogen content in the orthogonally sinusoidal structure with the same frequency and opposite frequency.Compared with Case 1-1, the hydrogen contents of Case 4-2 and Case 7-2 are increased by 10.94% and 8.76%, respectively.In an orthogonally sinusoidal flow channel with the same frequency, the hydrogen content of Case 5-2 with high horizontal angular frequency is increased by 11.07%.The results show that the horizontal angular frequency greatly influences the hydrogen content in the same frequency orthogonally sinusoidal flow channel.
The vertical velocity in the middle cross-section vertical to the XOZ plane along the channel is illustrated in Figures 15(a) and 15(b) with contours.Firstly, the area of the high-speed region at the trough also increases with the increase of the amplitude in both horizontal and vertical planes.Because of the forced convection, the highest vertical velocities in these high-speed regions are 0.18 m/s and 0.2 m/ s in Case 4-4 with the same frequency and in Case 7-4 with the opposite frequency, respectively.Secondly, with the increase of angular frequency in the vertical plane, the stronger periodicity in the flow channel makes the distributions of velocity in these high-speed regions more uniform, such as Case 5-1 and Case 8-1.However, the change of the angular frequency in the horizontal plane is not able to significantly improve the area and value of the high-speed region.
The horizontal velocity at the ACH-AGDL interface is illustrated in Figures 15(c) and 15(d) with contours.Similar to vertical velocity, the high-speed region area of the horizontal velocity also varies with the amplitude.Moreover, the maximum area of the high-speed region appears in Case 4-2 and Case 7-2, but the difference is that the area of the horizontal high-speed region with the same frequency is large enough to influence the adjacent high-speed regions [55].With the adjacent high-speed regions overlayed with each other, the horizontal transmission of materials is extremely strengthened.The maximum horizontal velocity can be almost 0.0767 m/s in Case 4-2.In Case 5-2, more high-speed regions will be shown in the whole region, and the maximum value of 0.09 m/s will be also found near the inlet.It is shown that the horizontal angular frequency   International Journal of Energy Research impacts the horizontal mass transfer of PEMEC with the same frequency.On the contrary, the vertical angular frequency only affects the number of high-speed regions but does not obtain the maximum horizontal velocity.
The lifetime and performance of the proton exchange membrane are greatly affected by the internal heat transfer performance of PEMEC.Therefore, in order to study the accuracy of heat transfer performance, the temperature distribution curves of the single channel and orthogonally sinusoidal channel along the y-axis direction of the PEMEC are shown in Figure 16.In Case 1-1, the temperature of the anode side flow field drops to a minimum of 353.70 K due to the rapid flow rate of the fluid center.The MEA is the main region of the electrolysis process, where the temperature reaches the highest value, 355.62 K.Although water is supplied to the cathode for cooling, the temperature on the cathode side is still higher than that on the anode side.The temperature distribution law of this model is the same as that obtained by Zhou et al. [48].In addition, the overall temperature is lower than the maximum working temperature of PEMEC [56].The results indicate that the internal heat transfer process of the model is accurate.
Case 4-5 is selected to study the influence of the orthogonally sinusoidal structure flow field on the internal heat transfer process of a single channel.Compared with Case 1-1, the maximum temperature at the MEA decreases by 0.57 K.The average temperature recorded at ACL is 0.08 K lower than the highest average temperature reported by Wang et al. [36].The minimum temperature of the anode side is almost the same as in Case 1-1.Besides, the introduction of the orthogonally sinusoidal structure on the cathode side significantly decreases the temperature of the cathode side flow field.In analysis of both mass and heat transfer, the orthogonally sinusoidal structure enhances fluid velocity to promote material transport and heat transfer.Sui et al. also found and proved these conclusions [57].Since some other practical engineering problems can be solved in this novel orthogonally sinusoidal structure by effective heat dissipation performance [58,59], the structure can be also regarded as an advantage in PEMEC.

Conclusions
In this study, a three-dimensional orthogonally sinusoidal flow channel model with periodicity and fluctuation is proposed.The different amplitudes and angular frequencies affect the electrochemical and mass transfer performance of PEMEC with vertically, horizontally, and orthogonally sinusoidal flow channels.The main conclusions are as follows: (1) The amplitude has a positive effect on the performance of the vertically sinusoidal flow channel PEMEC.Case 2-3 with the largest amplitude improves the polarization performance by 5.1% at 2.0 V, reduces the oxygen content by 22.3% at 40 mm away from the inlet, increases the maximum hydrogen content by 8.3%, and reaches the maximum vertical velocity of 0.079 m/s.The vertically sinusoidal flow channel accelerates the material transport process to supply the reactant and remove the products (2) The amplitude greatly influences the hydrogen production performance of the horizontally sinusoidal flow channel.Compared with the same frequency structure, the hydrogen content of Case 6-3 is higher than that of 17 International Journal of Energy Research 16.39%, respectively.Meanwhile, the oxygen content at 40 mm away from the inlet also decreases by 33.87% and 42.63%, respectively.The hydrogen content of Case 5-2 with the maximum horizontal angular frequency increases by 11.07%compared with Case 1-1.However, Case 5-1 and Case 8-1, with the maximum vertical angular frequency, inhibit oxygen removal compared with Case 4-5 and Case 7-5

BP:
Bipolar plate CH: Channel GDL: Gas diffusion layer CL: Catalytic layer PEM: Proton exchange membrane V upper : Curve of the upper wall in the left view V lower : Curve of the lower wall in the left view H left : Curve of the left wall in the top view H right : Curve of the right wall in the top view H: The difference between the curve equations of the right and left walls V: The difference between the curve equations of the upper and lower S: The sinusoidal curve of the walls projected Ŝ: The sinusoidal curve of opposite frequency L: The Electrolyte overpotential (V) E eq : Equilibrium potential (V) E eq,ref T : Temperature-dependent value of reversible cell voltage (V) T: Temperature (K) σ s : Electron conductivity of the solid phase electrode (S/m) σ m : Proton conductivity of the proton exchange membrane (S/m) ρ l : Liquid water density (kg/m 3 ) ρ g : Gas density (kg/m 3 ) u l : Liquid water velocity vector (m/s) u g : Gas velocity vector (m/s) R: Universal gas constant (C/mol) F: Faraday constant (J/(mol•K)) P: Pressure (Pa)

Figure 1 :
Figure 1: Sinusoidal functions in the horizontal and vertical planes.

2
International Journal of Energy Research are conducted to compare these effects in the horizontal and vertical planes, respectively.Finally, an orthogonally sinusoidal flow channel is formed by coupling the horizontally and vertically sinusoidal functions.The impact of these four parameters on all the performances of PEMEC, such as polarization curves, current density, molar fractions of hydrogen and oxygen, liquid water content, and velocity distribution, are assessed through simulations.Through the comprehensive study of spatially periodic fluctuations in this paper, a new design idea for the flow channel in PEMEC is provided to potentially help in further research.

Figure 2 :
Figure 2: Schematic diagram of the minimum repeating unit of the flow channel in space.(a) Fluctuation with the same frequency.(b) Fluctuation with the opposite frequency.

3
Figure 2(a) and the opposite ones ω i = −ω j in Figure 2(b) are both shown together.In order to explore the influence of the amplitude (A) and angular frequency (ω) on the performance of singlechannel PEMEC, the channels with different geometrical parameters are divided into 8 groups and 24 cases.The characteristics of the channel are further described by some symbols.S and L are used to represent the sinusoidal and linear curves of the walls projected, respectively.In order to avoid ambiguity in the following discussion, the descriptions of the same or opposite angular frequency of the left and right wall curves are called the same or opposite frequency.Therefore, when the frequency of the left wall curve is opposite to that of the right wall curve, the left wall curve can be further named Ŝ.The naming rule equation for groups is as follows:

( 2 )
All fluids in PEMEC are laminar because of the small Reynolds number (3) The diffusion layer, catalytic layer, and proton exchange membrane are homogeneous and isotropic mediums 4 International Journal of Energy Research (4) The gaseous products obtained by the electrochemical reaction are all incompressible fluids (5) The cross-permeation of hydrogen and oxygen in the proton exchange membrane is ignored (6) The contact resistances between each part are ignored 2.3.1.Transfer and Conservation of Electric Charges.The Butler-Volmer equation is as follows [44]:

Table 1 :
Symbolic description and minimum unit diagram of groups.
0.5 4π/5 4π/5 0.5 8π/5 4π/56 International Journal of Energy Research exchange current densities at the anode and cathode catalytic layers, respectively.α a and α c denote the charge transfer coefficient.η a and η c are the activation overpotential of the anode and cathode, respectively.

Figure 4 :
Figure 4: Comparisons of the polarization curve obtained from the present study and reference data.

Figures 9 (
a) and 9(b) are schematic diagrams of horizontally sinusoidal flow channels in the XOZ plane.Figures 9(c) and 9(d) show the local current density of PEMEC with the

Figure 7 :
Figure 7: (a) Cross-section of the flow channel in the YOZ plane.(b) Contour of local current density distribution at the ACL-PEM interface (A/cm 2 ).(c) Contour of H 2 O molar fraction at the AGDL-ACL interface.(d) Contour of H 2 molar fraction at the CCL-PEM interface.(e) Contour of vertical velocity in the middle cross-section vertical to the XOZ plane along the channel (m/s).

Figure 8 :Figure 9 :
Figure 8: PEMEC simulation results with sinusoidal curves in the horizontal plane.(a) Polarization curve of horizontally sinusoidal curve flow channel.(b) Curve of oxygen content at the center of GDL along the length.

Figure 10 :
Figure 10: (a) Contour of H 2 O molar fraction at the AGDL-ACL interface for Group 3. (b) Contour of H 2 O molar fraction at the AGDL-ACL interface for Group 6. (c) Contour of H 2 molar fraction at the CCL-PEM interface for Group 3. (d) Contour of H 2 molar fraction at the CCL-PEM interface for Group 6.

Figure 11 :
Figure 11: (a) Contour of horizontal velocity distribution in the XOZ plane below the flow channel for Group 3 (m/s).(b) Contour of horizontal velocity distribution in the XOZ plane below the flow channel for Group 6 (m/s).

13
International Journal of Energy Research the largest spatial amplitudes, the performances are lower compared to Case 4-2 and Case 7-2 due to the relatively weaker promotion from vertical plane amplitudes.

Figure 12 :
Figure 12: Simulation results of PEMEC with orthogonally sinusoidal structures of the same and opposite frequency.(a) Polarization curve of orthogonally sinusoidal structures with the same frequency.(b) Polarization curve of orthogonally sinusoidal structures with the opposite frequency.(c) Oxygen content curve of orthogonally sinusoidal structures with the same frequency at the center of AGDL along the length.(d) Oxygen content curve of orthogonally sinusoidal structures with the opposite frequency at the center of AGDL along the length.
) and 13(b) are shown.In Group 4 and Group 7, similar with the polarization performance, since the local current density at the inlet of the flow channel is positively correlated with the amplitude, Case 4-2 and Case 7-2 with the largest amplitudes in the horizontal plane have the best uniform distribution of the global current density.Secondly, the effects of horizontal and vertical angular frequencies in Group 5 and Group 8 are also not specific.Compared with Case 4-5 and Case 7-5, although the angular frequencies in the vertical plane of Case 5-1 and Case 8-1 are both increased, the local current density at the inlet will be decreased.On the contrary, due to the increase of the angular frequency in the horizontal plane, the local current density at the inlet of Case 5-2 is greatly improved, while there is almost no change in Case 8-2.

Figure 13 :Figure 14 :
Figure 13: (a) Local current density distribution contour of orthogonally sinusoidal structures with the same frequency at the ACL-PEM interface (A/cm 2 ).(b) Local current density distribution contour of orthogonally sinusoidal structures with the opposite frequency at the ACL-PEM interface (A/cm 2 ).

Figure 15 :
Figure 15: (a) Contour of vertical velocity in the middle cross-section vertical to the XOZ plane along the channel with the same frequency (m/s).(b) Contour of vertical velocity in the middle cross-section vertical to the XOZ plane along the channel with the opposite frequency (m/s).(c) Contour of horizontal velocity distribution at the ACH-AGDL interface with the same frequency (m/s).(d) Contour of horizontal velocity distribution at the ACH-AGDL interface with the opposite frequency (m/s).

Figure 16 :
Figure 16: Temperature distribution curves of Case 4-5 and Case 1-1 in the y-axis direction of the PEMEC.
Liquid water pressure (Pa) p g : Gas pressure (Pa) μ g : Dynamic viscosity of gas (Pa•s) μ l : Dynamic viscosity of liquid water (Pa•s) k: Thermal conductivity (W/(m•K)) C p : Heat capacity (J/(kg•K)) S l : Source term of liquid-phase mass conservation (kg/(m 3 •s)) S g : Source term of gas-phase mass conservation (kg/ (m 3 •s)) S s : Source term of solid overpotential (A/m 3 ) S m : Source term of electrolyte overpotential (A/m 3 ) S T : Source term of energy equation (W/m 3

Table 2 :
Dimension of a single channel in different cases.

Table 5 :
Source terms in the governing equations.
v η c + σ s ∇φ s 2 + σ m ∇φ m 2 , CCL i v η a + σ s ∇φ s 2 + σ m ∇φ m 2 − T dE eq,ref /dT , ACL ) and 7(d) are selected to measure the material transport in PEMEC's hydrogen production process.In Group 2, the linear curve of the walls projected i v,a : ).VariablesA z 1 : Amplitudes of the left and right wall curve functions (mm) A z 2 : Amplitude of the upper wall curve function (mm) ω i : Angular frequency of the left wall curve function (rad/s) ω j : Angular frequency of the right wall curve function (rad/s) ω k : Angular frequency of the upper wall curve function (rad/s).