Generic Dynamical Model of PEM Electrolyser under Intermittent Sources

: Proton Exchange Membrane (PEM) water electrolysis system is one of the promising technologies to produce green hydrogen from renewable energy sources (wind and solar). However, performance and dynamic analysis of PEM water electrolysis systems are challenging due to the intermittent nature of such sources and involved multi-physical behaviour of the components and subsystems. This study proposes a generic dynamical model of the PEM electrolysis system represented in a modular fashion using Bond Graph (BG) as a uniﬁed modelling approach. Causal and functional properties of the BG facilitate the formal PEM electrolyser model to adapt and to ﬁt the different conﬁgurations of the electrolyser ranging from laboratory scale to industrial scale. The system-speciﬁc key parameter values are identiﬁed optimally for a laboratory-scale electrolyser system running on a multi-source energy platform using experimental data. The mean absolute percentage error between simulation and experimental data is found to be less than 5%. The performance characteristic curves of the electrolyser are predicted at different operating temperatures using the identiﬁed key parameters. The predicted performance is in good agreement with the expected behaviour of the electrolyser found in the literature. The model also estimates the different energy losses and the real-time efﬁciency of the system under dynamic inputs. With these capabilities, the developed model provides an economical mean for design, control, and diagnosis development of such systems.


Introduction
Using hydrogen as storage for surplus electricity from Renewable Energy Sources (RES) has gained popularity in recent years due to its numerous advantages.Hydrogen can be produced from water, which is available in abundance.It has the highest energy content while being the lightest and clean fuel [1].Hydrogen production through electrolysis of water is a well established technique and can be easily integrated with suitable electric sources [2].
In electrolysis, water molecules are dissociated into hydrogen and oxygen gases with the application of electrical energy [3].Alkaline electrolyser and Proton Exchange Membrane (PEM) electrolyser are the two well-matured electrolysis technologies for hydrogen generation which have been utilised extensively on a commercial scale [4].Anion Exchange Membrane (AEM) electrolysis is also an emerging technique for water electrolysis and it tries to combine the advantages of both alkaline and PEM electrolysis techniques [5].PEM based water electrolysis has gained popularity in recent years due to its high performance.Various installations all over the world have been done to integrate the RES and the electrolyser as a means of storage [6].PEM electrolyser has superiority over other electrolysers for load following when running on intermittent energy sources [7].Due to the intermittent nature of the renewable power sources, the dynamics of the PEM electrolyser is affected and it becomes necessary to study the effect of the intermittency of the power sources on the overall performance of the electrolyser.
Figure 1 demonstrates the basic working principle of the PEM electrolyser.PEM electrolyser cell consists of two half cells separated by a thin PEM.Half cell reaction for each electrode (anode and cathode) and the overall reaction are written as  The water is fed to the anode side of the cell where water is reduced to oxygen, positively charged hydrogen atoms (protons) and electrons.The oxygen produced during this half-cell reaction is removed with the unconsumed water.The protons move through the electrolyte membrane towards the cathode where they combine with the electrons from the Direct Current (DC) power source and form the hydrogen gas.The water feed to the cathode side is optional as it is only there to facilitate the efficient removal of the hydrogen.
These reactions take place at the catalyst layer coated on each electrode.The diffusion layer on each side helps in efficient current distribution and also connects the Membrane Electrode Assembly (MEA) to the distribution plates.Distribution (bipolar) plates aid toward structural integrity of the cell.It also separates one cell from the other when cells are assembled together in the form of a stack to deliver the required hydrogen flow rate.In addition to the cell/stack, an electrolyser has auxiliary components to ensure the proper functioning of the stack.This includes power supply/voltage controller, input water conditioning system, water circulation system, water-gas separators for hydrogen side as well as oxygen side (optional), heat exchanger, controlled valves and safety devices.These components also consume energy and directly affect the performance and efficiency of the electrolyser.The general schematic of a PEM electrolyser is shown in Figure 2a.Depending on the size and specification of the electrolyser, one or more components may or may not be present in the configuration.Figure 2b shows various multi-physics phenomena involved in a PEM electrolysis cell.PEM electrolyser is a complex system that involves the non-linear relations to describe its dynamics due to coupling between various physical phenomena such as electro-chemical, thermo-electrochemical and thermo-fluidic.Multi-physics phenomena in PEM electrolyser cell/stack (inspired by Reference [2]).
High cost of the PEM electrolyser has always been a challenge for the researchers.New materials and cell design are being developed to reduce the overall cost of the electrolyser [8,9].There is also a need for developing the tools in order to test the electrolysers under different conditions to understand various phenomena taking place inside the electrolyser to predict the performance of the electrolyser in the cost effective way.Mathematical modelling plays a very crucial role towards this objective and acts as a dynamic connection between the electrolyser and the intermittent power source [10].It also enables the researchers for the design optimization and developing the control for the system [2].
Exhaustive reviews about PEM electrolyser modelling show that there is not much work done on the modelling of PEM electrolyser as compared to the PEM fuel cell [2,10].However, the basic concept related to fuel cell can be adapted for modelling, control and optimisation of electrolyser due to structural and process similarity [11].An important contribution for fuel cell modeling is found in Reference [12], where three-dimensional multi phase model is developed in computational fluid dynamics to study the influence of various parameters on I-V characteristic curve of the system.In Reference [13], optimal parameters are found for MEA using numerical model to enhance the performance of the of PEM fuel cell.Some other works as in References [14,15] can be found which are dedicated for fuel cell control analysis and its performance optimisation.Contrary to fuel cell modelling and analysis, some works for PEM electrolyser based on analytical models have been developed for studying the effect of different key variables on the polarization curve which can simulate it with fair accuracy [2,10].As per literature review in Reference [16], mostly the empirical and semi-empirical models for PEM electrolyser exist that can predict the behaviour of the electrolyser under varying operating conditions (for example, different temperatures and pressures).These models are mainly focused on phenomenon understanding and for developing control algorithms.Some recent works related to modelling and performance analysis of the PEM electrolyser are presented in References [9,[17][18][19][20][21].The authors of Reference [17] developed the power-hardware-in-loop simulator to study the effect of variable solar power input to the PEM electrolyser by considering it as a part of smart grid.The characteristics of several power supply electronics are mainly analysed based on simple electro-chemical model of the PEM electrolyser.In Reference [18], a mathematical model for the PEM electrolysis cell has been developed that incorporates electrochemical reactions, gas transport mechanism and various physical phenomena.Authors have also defined the thermal efficiencies to analyse the cell performance.In Reference [19], a simplified mathematical model based on zero-dimensional dynamics and multi-physics approach has been proposed to avoid complex model with too may parameters.The model incorporates different multi-physics phenomena such as electrochemical model, thermal model and fluid-dynamic model.The power consumption by the auxiliaries of the electrolyser has also been incorporated to study the expenditure of thermal and electrical energy.Recent work in Reference [20], equivalent electrical circuit model for a PEM electrolyser has been proposed for studying the influence of different operating conditions and power electronics ripples effect on the cell voltage.The work is limited only for study of cell voltage modelling under static and dynamic operating condition.In Reference [21], an analytical dimensionless model for a PEM electrolyser (single cell) has been proposed.A mathematical model in closed form for overpotential variation, current density distribution and water content distribution in membrane are obtained in non-dimensional form.The developed method offers a tool for study of water management through the PEM electrolyser.
Some authors have also used graphical modelling techniques to develop the models of the PEM electrolysers to couple the multiple phenomena into one model.In Reference [22], a model based on causal ordering graph has been proposed for control analysis of multi-source hybrid system, where electrolyser with its auxiliaries are modeled for controlling and managing the flow of power for desired rate of hydrogen production.In Reference [23], other graphical approach has been presented for modelling of PEM electrolyser system using energetic macroscopic representation.This technique incorporates the different interaction of thermal, fluidic, electro-chemical phenomena in the model that occur simultaneously.However, authors have neglected the fluidic aspects in the modelling.In Reference [24], a graphical approach based on Bond Graph (BG) has been proposed for global modelling of PEM electrolyser system with its auxiliaries.It also incorporates the interaction of different phenomena such as fluidic, mass transport, electro-chemical and thermal for behavioural and performance analysis of the electrolyser.However, the model is complex and it cannot be directly sized and used as a generic tool to represent electrolysers of different scales.BG is a very good tool for modelling multi-physics dynamical systems [25,26].BG technique is not only suitable for multi-physics systems modelling, but also for developing diagnosis, prognosis, actuator sizing and control analysis algorithms [14,[27][28][29][30].In Reference [27], BG technique is proposed for modelling and fault diagnosis for chemical reactor; while in Reference [28], fault diagnosis for hybrid system is proposed.Then, the BG approach is further extended for hybrid system prognosis in References [29,30].There are some models for the study of the durability of the electrolyser for long term diagnosis and prognosis [2].The main disadvantage of these models is that they are limited to specific design of the electrolyser and work for known operating conditions only [10].
The modelling of complex multi-physics systems like PEM electrolyser is one of the most crucial tasks for studying the real behaviour of components, subsystems and overall system behaviour under dynamic intermittent operations.Moreover, system performance degrades at different levels (components and subsystems levels) due to ageing and dynamic operational behaviour.Such studies require a modular design approach, breaking the complete model into different sub-model levels, so that variation of different component parameters on the monitoring variables of the system can be analysed and tested easily and effectively.From the industrial perspective, the model should be adaptable with real monitoring of the system behaviour on suitable supervision platform.Current study considers such aspect of modelling for the PEM electrolyser system running on the intermittent sources (solar and wind energy) using BG as a unified modelling approach [25,26].The proposed approach is not focused only on stack/cell level but also considers the auxiliaries of the electrolyser system.Also, the developed model incorporates the real-time efficiency estimation of the electrolyser.Causal and functional properties of the BG facilitate the formal PEM electrolyser model to adapt and to fit the different configurations of the electrolyser ranging from laboratory scale to industrial scale.Thus, in this study, developed model is adapted for laboratory-scale electrolyser system running on a multi-source energy platform by identifying the system-specific key model parameter values using experimental data.The predicted behaviour of the electrolyser using developed model provides a good agreement with the expected behaviour of the electrolyser found in the literature.
Apart from the introduction, a brief description about the BG modelling approach for multi-physics systems is presented in Section 2. Based on this approach, the modelling of the PEM electrolyser in a modular fashion is discussed in Section 3. Details of the platform of the electrolyser running on intermittent sources used for experimental validation and results are explained in Section 4. In Section 5, conclusion and perspective of the current study is derived.

BG Technique for Model Building
BG is a well-adapted multi-disciplinary and unified graphical modelling approach to describe complex process having multiple energy exchanges [25,31,32].A brief introduction to BG methodology for multi-physics systems is presented in this section.

BG Elements and BG Variables
BG can be denoted as G(N, B), where nodes N represent the BG elements that correspond to energetic physical elements (inertia, resistance and capacitance), source elements (battery, pump, etc.), power/energy constraint BG connecting elements and technological elements (subsystems) and bonds B represent the set of oriented edges that correspond to power/energy exchange among nodal elements.Bond is labeled by two power variables: called flow ( f ) and effort (e) variables whose product provides the value of physical power flow in that bond.Note that for some energy domains like chemical, the product of these variables on a bond may not be the representation of physical dimension of a power, however it is used for systematically model such energy domains, called pseudo BG model [24,26,33].All the BG theory obviously remains valid in such cases also.
In Figure 3a, a physical link between subsystem A and subsystem B using the half arrow power bond is presented that represents the direction of energy/power flow among them.In the represented configuration, energy flow is considered towards the destination node B, only when the product of power variables is positive; otherwise, energy flow is in the reverse direction, that is, away from node B. In BG, causality (represented by cross-stroke) is an important structural property that determines the relationship between power/energy variables based on cause and effect analysis.The node which receives the effort information (or gives the flow information), cross-stroke is marked near to that node.For example, according to causality in Figure 3a and corresponding block diagram as in Figure 3b, B receives effort information from A as cross-stroke is marked near B, while A receives flow information from B. In Table 1, various flow and effort variables belonging to different energy domains are shown.Different systems can be modelled in BG by employing a handful of elements (R, C, I, S f , D f , Se, De, TF, GY and J) [34].Here, R is resistive element (represents dissipation phenomena); I is the inertia element (stores kinetic energy); C is capacitance element (stores potential energy); S f and Se are, respectively, the source of flow and source of effort; D f and De are, respectively, detectors for flow and effort (virtual sensors); TF is the transformer and GY is the gyrator.These elements represent the energy transformation between different domains.J represents the junction that accounts for energy conservation laws.It is used to connect the elements based on their interaction with other elements or the sub-systems.Its value is 0 (to connect elements having same effort) and 1 (to connect elements having same flow).R element may also be used as a multi-port resistive-source element (noted RS) to represent the active resistance that generates entropy.For example, an electric coil that is used for heating acts as a resistance from electricity point of view, but also acts as a source of thermal energy from thermal viewpoint.The vectorial representation is used to represent the coupling variables of a complex system [26].The coupled power variables in vectorial form are depicted as: where F and E, respectively, are the flow and effort vectors.The suffices f l, th, ch, respectively, denote the fluidic/hydraulic, thermal and chemical energy domain.Figure 4 shows different ways to illustrate a vector bond.

Different Levels of Modelling Abstraction
To grasp complex systems, BG technique integrates the different levels of modelling, that is, technological, physical, mathematical and algorithmic levels using the common tool.In technological level, different subsystems of a dynamical system are first identified and then interconnected through energy or power variables.This level is represented by using word BG that shows the architecture of the system in a modularized form.In physical level, the system is modelled in the form of lumped parameters (showing different physical phenomena such as energy storage, dissipation and transformation) using generalized BG elements with power/energy bonds.In mathematical level, dynamical behaviour of the system is represented in the form of mathematical equations such as differential and algebraic equations or in state-space form.The constitutive relations of the constraints and the components of the BG model provide these equations in mathematical form.In algorithmic level, the causal property of BG technique decides how these mathematical equations are algorithmically derived from the graphical BG model.Thus, the structural and causal properties of BG technique enable the dedicated software to systematically and algorithmically obtain the system dynamical equation either in the form of state-space equations or in the form of differential equations for the simulation and analysis purpose.Based on the BG approach, a number of software have been developed and are available for use such as Symbols-Shakti [35] and 20-Sim [36].In this work, model builder of Symbol-Shakti is used for developing the generic PEM electrolyser model where the structural integrity of different components and sub-systems for the global system modelling is checked.Once the models of different sub-systems (capsules) are built, their corresponding Matlab-Simulink models are systematically derived from implementation point of view.

Modular Building (Capsules)
In Symbols-Shakti, the capsules are properly modelled subsystems that have single or multiple input and output ports.They already have their partially derived equations.When the capsules are assembled to form the complete model, the equations are linked together to form behavioural equations for the complete model.For example, the capsule of PEM stack is developed by coupling the BG sub-models with different energy interactions (electro-chemical, thermal and fluidic).The internal model of the capsule in Symbols-Shakti for a PEM stack is shown in Figure 5.

Inputs
Outputs Capsule (PEM stack)  Depending upon the nature of interaction of the capsules with its environment, the ports can be defined as effort or flow input port and effort or flow output port.Based on the requirements these ports can be attributed as optional and hence are not needed to be connected in order to create the global model.Coupled energy interaction between the sub-models can be represented by a vector bond as shown by bond number 26 in Figure 5.A number of capsules are already available in the database of the software for the commonly used components such as tanks, valves, pipes, heaters and sensors.User defined capsules can be stored in the database and can be used and modified as per the model requirement.Figure 6 shows the user defined capsules, developed for the PEM electrolyser system.

Grammar and Connectivity Rules
Using the capsules, a global model of a complex system can be assembled in the model builder of the Symbols Shakti software.Figure 6 shows the graphical user interface of the model builder.The global model is assembled in the form of piping and instrumentation diagrams using the capsules.For example, the architectural model of a PEM electrolyser is shown in Figure 6.The model builder also generates the behavioural equations for the whole system automatically.To successfully connect the capsules with each other there are certain rules that have to be respected.The connectivity of the capsules is automatically checked by the software for validity.Only the like ports can be connected.The flow (resp.effort) output of a capsule can only be connected to the flow (resp.effort) input of the other capsule and vice versa.For example, Figure 7 demonstrates the requirement of the coherence of causality for the connection of the capsules of two tanks.From Figure 7a, it can be seen that it is not possible to directly connect two tank capsules together as the input for the tank capsule is flow and the output is effort.In order to connect the two tank capsules the resistance of the pipe needs to be considered which takes effort as input and flow as output (as shown in Figure 7b).References [37][38][39] and the reference manual of the software can be referred to for further details.

PEM Electrolyser Modelling
In this section, a generic dynamical model of the PEM electrolysis system represented in a modular fashion (in the form of subsystems/capsules) using BG approach is presented.The assumptions taken for the current study are as follow: • The cells constituting the stack are identical in nature and connected in series.Thus, the stack with N cells can be modelled as an equivalent single cell that has the same dynamics of the stack.

•
Uniform fluid flows and current distribution are considered between cells.• Overpotential due to mass transport or diffusion is negligible with the assumption that PEM system usually operates at low current density.• Electrolysis reaction kinetics is assumed firmly as a Faradic process and considers that there is no mass limitation problem in the system.• Gases produced are assumed to have similar properties as that of an ideal gas and the partial pressures of these gases are governed by Dalton's law.

•
Temperature is homogenous throughout the stack.

•
Cell is operated below the boiling temperature of the water.

•
The system parameters are considered as lumped parameters.Pumps and fans are assumed as perfect mass flow sources.

Technological Representation
The schematic architecture of the PEM electrolyser and the physical variables exchanged between its subsystems are shown by the word BG presented in Figure 8.In ,

Modular Representation
BG capsules are used, from one part, to describe the inner components of the electrolyser, and from another part, to explicit the coupling between the different energies hightlighted by the word BG.

Stack Model
The stack is the heart of the electrolyser.The stack model is the integration of various coupled phenomena which have been modelled as different sub-models to provide more insight.

Electrochemical Sub-Model of the Stack
The electrochemical sub-model is one of the crucial sub-models of the electrolysis system.It forms the basis to provide a relation between the cell current and cell voltage at different operating conditions, that is, at different operating pressures and temperatures.The sum of all the electrochemical phenomena leads to a sub-model linking cell voltage and the current density, which is usually represented by a polarization curve.This sub-model explains the real kinetics of the reaction that occurs and gives the information about the amount of product flows according to the water consumed in the electrolysis system and also the production of heat in the stack system.This sub-model considers that electrochemical phenomenon occurs at steady state and thus the sub-model is responded with no time delay which means that there is instantaneous response with respect to change in any input to this sub-model.It is assumed that transient response is very fast and dies out quickly, so it is neglected in model which is also well experimentally demonstrated by the work presented in References [40,41].The electrochemical BG sub-model in the form of a capsule is shown in Figure 9.
This sub-model can also be used to find the required cell voltage, E cell , at any operating condition with respect to reversible potential, E rev , and different overpotentials which occur due to current flows in the cell results in irreversible heat dissipation.The reason of irreversible heat dissipation in the cell is due to the following: (i) the activation (Faradaic losses) voltages E act,a and E act,c , respectively, occur at anode-electrolyte and cathode-electrolyte interfaces due to disturbance from the equilibrium chemical reaction and involved activation energies barriers in the preferred reaction.(ii) the electrical ohmic overpotential (non-Faradaic losses) due to internal cell resistance, that is, it is proportional to the current which flows through electrodes, current collectors, bipolar plates and corresponding interconnections between them.Also, this overpotential occurs due to resistance offered by both electrolyte and membrane to the ions flow that separates both the electrodes.Generally, the electrode material in PEM electrolysers has high conductivity, so the flow of electrons is much faster as compared to ionic flow, and thus, usually the ohmic resistance due to ionic transport is considered.(iii) Overpotential due to mass transport or diffusion (non-Faradaic losses).It should be included in the model for a system having higher current density.Nernst equation is used to define the overpotential due to diffusion and it reveals that the diffusion limitation increases with increase in concentration of the product species at the reaction interface [42].
In this sub-model, E cell is modelled as a modulated source of effort BG element, that is, Mse : E cell and all the overpotentials are modelled by generalized RS coupled resistive elements as RS : R act,a , RS : R act,c and RS : R ohm corresponding to anodic activation resistance, cathodic activation resistance and ohmic resistance, respectively.Here, the generalized RS coupled resistance element shows the coupled energy dynamics and links the electrical domain with the thermal domain.It acts as a source and transmits the generated heat or entropy to the thermal domain due to the current flow in electrical domain [26,43].Also, the transformation of electrical energy into chemical energy and vice a versa is modelled by transformer T f generalized BG element, which describes the Faraday's law by linking the rate of reaction flow ( ξ) with the current flow through cell (I cell ), number of moles of electron transferred n (here n = 2) and Faraday's constant F; and, also by linking reversible cell voltage(E rev ) with free energy of water dissociation (∆G R ), n and F, which are represented in Equations ( 3) and (4).
Here E rev is the reversible potential which is the minimum potential required for the electrolysis reaction.As per Reference [2], the reversible potential E rev is expressed as in Equation ( 5) at any operating temperature and pressure condition.
where E 0 rev is the standard reversible cell potential at standard operating conditions (at standard temperature and pressure), R is the ideal gas constant, p i is the partial pressure of the i th species and a H 2 O is the chemical activity of water.The E 0 rev is usually temperature dependent and is empirically defined as [2,23,44,45] Activation overpotential: There is a relation between the kinetic part of the reaction, that is, the current density, with the thermodynamic part, that is, the overpotential.At equilibrium current exchange density (J 0 ) is defined and it is related with the exponential function of negative free Gibbs energy at equilibrium (∆G 0 ).The value of the parameter J 0 or ∆G 0 should be obtained from the experiment.For the electrolysis process, that is, possibility of occurring the dissociation of water, the cell voltage is usually more than the reversible standard voltage which results in overpotential because of current flow, and thus, it contributes to polarization.This current flow is usually the difference between anodic and cathodic currents at the non-equilibrium condition.The relation between the actual current and the overpotential is usually obtained by using Buttler-Volmer equation as given in Equation (7) and used in much of the literature [14,17,44,46,47].
where J cell = I cell /A M is current density of cell, A M is cross-sectional area of membrane.Also, J 0,k and α k represent the current exchange density and charge transfer or symmetry factor coefficient, respectively, at anode or cathode.Usually, symmetry factor is taken as 0.5; however it lies between 0 and 1.Thus, at α k = 0.5, the overpotential E act,k at anode or cathode can be expressed by Equation (8).
In thermochemical model (Figure 9), Equation ( 8) is the constitutive relation for nonlinear resistive element RS : R act,k .The parameter J 0,k for each electrode is expressed as where J re f 0,k and ∆G 0,k are obtained from the experimental data using model fitting technique and the temperature T is obtained from the thermal model of the stack.
Ohmic overpotential: This overpotential is modelled by the dissipative element, that is, generalized BG resistive element R : R ohm .According to the resistive causality in the BG sub-model (Figure 9), the constitutive relation for R : R ohm is given as where E ohm is the ohmic overpotential.The ohmic resistance R ohm mainly includes the resistance offered by the membrane to ions flow and also includes the other ohmic resistance R other offered by the internal cell components except the membrane.Thus, R ohm is obtained from Equation ( 11) that depends upon the properties (σ M ) and parameters of the membrane (L M , A M ).σ M represents the membrane conductivity, L M and A M , respectively, represent the length and cross-sectional area of the membrane.Moreover, the resistance R other can be obtained from the experiments using model fitting technique.
where d M is the ratio of L M to A M and σ M can be empirically defined as in Equation (12).Membrane conductivity depends on cell temperature T and the membrane water content parameter γ.
The parameter γ value varies from 14 to 25, whose value depends on the nature of membrane hydration.For a poorly hydrated membrane, it is taken as 14 while for fully hydrated it is taken as 25 [42].
Based on BG causal properties and constitutive relations for different elements, different relations can be systematically obtained from the sub-model of electrochemical phenomenon.At junction 1 1 (refer Figure 9), the conservative phenomenon can be written as in Equation ( 13): Since, at junction 1 1 all the flows are equal, that is, Thus, as per causality in sub-model (Figure 9), e 3 = E ohm can be obtained from Equation ( 13) as where value of e 1 = E cell , e 2 = E act,a , e 4 = E act,c , e 5 = E rev .So, Equation ( 14) becomes where E cell is derived from BG source element Mse : E cell as input causality.E act,a , E act,c are derived from constitutive relations for BG elements RS : R act,a , RS : R act,c as in Equation ( 8) and E rev is obtained from T f : 1/2F as in Equation ( 4) or (5).Thus, using BG submodel, E ohm can be systematically derived as Thus, using Equation ( 10), the cell current I cell can be derived from the electrochemical BG submodel by using constitutive relation for R : R ohm , which is used for controlling the hydrogen production by the cell.The required input cell voltage can also be derived by rearranging the Equation ( 16) Likewise, using the property of junction 0 1 (refer Figure 9), the conservative phenomenon can be written as Since, at junction 0 1 all the efforts are equal, that is, e 6 = e 7 = e 8 = e 9 .Thus, as per causality in sub-model (Figure 9), f 9 = Qirr can be obtained from Equation (18) as where value of where Qirr denotes the irreversible heat flow which is systematically derived from BG submodel and obtained as the summation of activation losses at anode, cathode and ohmic loss.

Chemical-Fluidic Sub-Model of the Stack
In chemical-fluidic sub-model, the amount of hydrogen and oxygen production rates can be predicted with respect to consumed water.The BG chemical-fluidic sub-model in the form of a capsule is shown in Figure 10.

Electrochemical
Sub-model In the sub-model, transformers T f : ν i and T f : M i generalized BG elements are used to find the rate of production of the products with respect to consumed reactant.Also, C : C i represent the storage of the matter for the i th species.The rate of produced mass for the i th species, that is, ṁi , is obtained from the rate of reaction flow ξ from Equation (21)

Fluidic and Mass Transfer
where ν i and M i , respectively, denote the coefficient of stoichiometry and molar mass of i th species.
Thermal Sub-Model of the Stack Thermal sub-model predicts the dynamic behaviour of the temperature evolution inside the stack which ultimately affects the relation between cell current and voltage.Thus, durability and efficiency of the electrolysis system is also affected.The complete thermal BG sub-model in the form of a capsule is represented in Figure 11.The sub-model considers the major contribution of the heat from the different sources and phenomena such as heat input to stack by the water inflow towards anode and cathode side, heat taken away by the water outflow from stack towards hydrogen and oxygen separators, thermodynamics of chemical components during reaction (dissociation of water and production of gases) or endothermic nature of chemical reaction, Joule effect phenomenon due to circulation of charge/current, heat due to entropy change in reaction, thermal activated phenomenon due to mass transfer and diffusion and thermal effects due to system enclosure temperature.This thermal sub-model is considered as a first order nonlinear model having one dynamic generalized BG capacitance C element.Here generalized C : C stack element, considered as a constant lumped parameter, models the thermal capacitance of the stack, and thus, it is assumed that temperature is homogenous throughout the stack.This lumped parameter can be estimated from the experiment using model fitting technique.From the BG thermal sub-model, the temperature of stack using constitutive relation of C : C stack and junction 0 2 is obtained as In Equation ( 22), Ḣrec,a , Ḣrec,c and Ḣa,Osep , Ḣc,Hsep respectively, represent the enthalpy rate due to water inflow from recirculation circuit towards anode side, cathode side and water outflow from the stack towards oxygen and hydrogen separator side.These enthalpy flows are related with the molar enthalpy (H i ) and molar flow rate ( ṅi ) of the involved species in the reaction.Qirr is irreversible heat flow due to activations and ohmic losses, QH that is, QH = (I cell × ∆H R )/(2 F) denotes heat flow rate due to chemical components during reaction or endothermic nature of reaction (dissociation of water and production of gases), QS that is, QS = (−I cell × ∆S R )/(2 F) denotes the heat flow rate due to entropy change and Qst,enc denotes heat transfer from stack to system enclosure due to temperature difference.

Fluidic and Mass Transfer Sub-Model
The global behaviour and the efficiency of cell can be affected due to the phenomena of fluidic motion and mass transfer and these cannot be ignored while developing the model of the electrolysis system.Thus, the phenomenon of fluidic motion and mass transfer must be integrated within thermal and electrochemical sub-models.It is assumed that there is no accumulation of fluid inside the stack and the product produced is continuously evacuated from the stack.This sub-model includes the transfer of water from anode to cathode by considering the electro-osmosis phenomenon and also considers the transfer in reverse direction, that is, from cathode to anode, by the process of diffusion.This sub-model also includes the crossover of produced gases, that is, oxygen and hydrogen gas flows towards the opposite electrode with respect to the electrode where they produce.The BG sub-model of fluidic and mass transfer phenomenon in the form of a capsule is shown in Figure 12.

Anode Side Recirculation
Sub-model  In this sub-model, resistive BG field element R : and transformer element T f : n eo M H 2 O are used, respectively, to model the crossover or diffusion flows due to water diffusion, oxygen diffusion, hydrogen diffusion and electro-osmosis diffusion.R Hyst,a and R Hyst,c are the hydraulic resistances at anode and cathode side, respectively.Thus, the model provides the output mass flow of the mixture from the stack at the anode and cathode sides using constitutive relation of junction 0 3 and 0 4 as given by Equations ( 23) and (24), respectively.
where ṁrec,a , ṁrec,c , respectively, denote rate of inflow masses from recirculation circuit to anode and cathode side.ṁprod,O 2 , ṁprod,H 2 and ṁcons,H 2 O , respectively, denote production of oxygen, hydrogen and water consumed which are obtained from electrochemical model.Also, ṁeo,H 2 O denotes the rate of water mass flow due to electro-osmosis, however, this phenomenon does not depend on thickness of the membrane.The electro-osmosis flow, ṁeo,H 2 O , mainly depends on the magnitude of electric current flow during electrolysis process and the number of dragged water molecules by a proton ion as given in Equation ( 25).This number is denoted by n eo and called as electro-osmosis coefficient [44,48].This flow, ṁeo,H 2 O , is obtained by using the generalized transformer element T f : n eo M H 2 O in the fluidic model which is coupled with the rate of reaction ξ of the thermochemical model.
Also, ṁdi f f ,i denotes the flow of diffusion for the i th species through membrane and this flow for i th species is obtained from the constitutive relation of the generalized resistive element R : R di f f ,i .This sub-model shows that the diffusion flow for the i th species is proportional to the change in partial pressure at the both side of the membrane ∆p i and also the resistance R di f f ,i depends on the length of membrane L M , diffusion parameter D i and the Henry's parameter H i [49,50] as represented in Equations ( 26) and ( 27), respectively.
Thus, the global stack model, as shown in Figure 13, is obtained by integrating all the sub-models (capsules) such as electrochemical, chemical-fluidic, thermal, fluidic and mass transfer sub-models.
Oxygen Separator Sub-model  The state-space equation for the stack can now be derived systematically and automatically from the BG sub-model of the stack (shown in Equation ( 28)), where x O 2 , x H 2 , x H 2 O and x st are the states of the stack and function f shows the non linear relation for the state x st .Similarly, the state-space equations for the complete model and each sub-model can also be generated for developing control algorithms.Further, BG modelling of other subsystems and components as prescribed in word BG (see Figure 8) such as: electrical converter, hydrogen and oxygen separator vessels, heat exchanger, cooling and recirculation systems, purification system, the pneumatically controlled hydraulic valves and the stack enclosure, is required for the complete modelling and dynamical study of the electrolysis system.

Converter Sub-Model
Converter provides the control power input to the stack according to the operational point of the electrolyser.Usually, converter has a very fast response time (less than 0.1 s); thus, the dynamic model of the converter has instantaneous power output.This sub-model is shown in Figure 14.Here, physical phenomena of the converter is modelled by using a modulated transformer (MT f : β) and resistive (R : R conv ) generalized BG elements.In the model, Mse : V in is the input voltage to the converter from the electrical source (solar and/or wind source) and V stack is the output of the converter which is fed in to the electrolyser for water dissociation.Thus, using dynamic BG sub-model (Figure 14), constitutive relations can be derived as P in = P stack + P conv (30) where β is the coefficient of transformer, P in = V in × I in denotes input power, P stack = V stack × I stack denotes output power from the converter which is input to the stack, P conv denotes the dissipation power of the converter due to its internal resistance and ε conv denotes converter efficiency.

Separator Sub-Models
Usually gas is collected over water in a vessel and due to the light weight of the gas it goes up and can be separated from the water.When the gas pressure reaches a certain fixed preset value in the separator it comes out from the separator vessel to the storage tank.In the electrolysis system, two separators are used for gas liquid separation: (i) hydrogen separation and (ii) oxygen separation.The separator subsystem includes the interaction of fluidic, thermal and chemical phenomena.In the modelling of each separator, it is assumed that the vessel is in cylindrical shape having a cross-sectional area A sep,i .The BG sub-models of the hydrogen and oxygen separators are shown in Figures 15a,b, respectively.
Separator fluidic Phenomena: The hydraulic pressures or the water levels in the Hydrogen Separator Vessel (HSV) and Oxygen Separator Vessel (OSV) are determined from the continuity equations and these phenomenons are modelled in the sub-models as in Figures 15a,b, respectively. ..
..  The hydraulic capacity as a lumped parameter for i th separator (i = H 2 , O 2 ) is denoted by where g is gravity.First the entire phenomena are discussed for HSV, and then likewise these are presented for OSV.Thus, for HSV from Figure 15a, using constitutive relation of sep,H 2 and junction 0 6 , HSV pressure can be obtained as where P sep,H 2 = ρ w × g × L sep,H 2 is water pressure in HSV, L sep,H 2 is water level in HSV, ρ w is water density, ṁcw,Hsep is the water mass flow which is the product of known mass fraction of water x c,w and the mass flow of mixture, ṁc,Hsep (obtained in Equation ( 24)), ṁsepv is the water mass flow from the HSV to OSV through the separator valve (modelled using modulated resistive element MR : R sepv ).The sub-model of the separator valve [39] is highlighted in the Figure 15a.ṁsep,H 2 is the water mass flow out from the HSV to cathode side recirculation circuit.
Separator thermal Phenomena: The thermal capacity of the vessel is modelled by using lumped parameter, C : C th sep,i , whose value depends on cross-sectional area A sep,i and water level L sep,i in the vessel.Also, it is assumed that the coefficient of heat transfer of the vessel is constant.Thus, for HSV from Figure 15a, using constitutive relation of C : C th sep,H 2 and junction 0 7 , HSV temperature can be obtained as where T sep,H 2 temperature inside HSV, Ḣc,Osep is the enthalpy flow of the water and gas mixture from the cathode side of the stack to HSV and Ḣsep,H 2 is the enthalpy flow out from the HSV to cathode side recirculation circuit and QHenc is the rate of heat loss due to temperature gradient between HSV and system enclosure.Note that all the enthalpy flow is correlated with their respective mass flow of water and mixture.The value of the enthalpy flows Ḣgas,i and Ḣj coupled with mass flows ṁgas,i and ṁj (i = O 2 , H 2 ; j = tank) can be obtained from the expression as presented in Equations ( 34) and (35), respectively.This coupling between thermal and fluidic phenomena is represented by generalized where C pj and ∑ i C pi , respectively, denote specific heat constant of water and gas.The thermal loss, QHenc , is modelled by generalized BG resistive element R : R Hsep,enc .Also, the system enclosure temperature is modelled as BG modulated source of effort element as Mse : T enc , which is obtained from the thermal sub-model of the enclosure.Moreover, generalized detector element De : L sep,H 2 is used which shows the measurement of separator water level that is further utilized by the controller for controlling the percentage opening of separator valve as per level set value.Separator chemical Phenomena: For calculation of partial pressure, P ch Hsep,i , of the i th species, the chemical capacitance is modelled by using lumped parameter, C : C ch Hsep,i at cathode side separator, whose value depends on molar mass M i , volume V Hsep,i of the i th species, temperature T Hsep,i and the gas constant R which is obtained from Equation (36).
Thus, for HSV from Figure 15a, using constitutive relation of C : C ch Hsep,i and junction 0 8 /0 9 / 0 10 , gases partial pressures can be obtained as where ṅci denotes gas mass flow rate at cathode side, x c,i and x c,mix , respectively, denote mass fraction of i th species and gas mixture (hydrogen, oxygen and water vapor) leaving from cathode side of stack to HSV, ṁgas,H 2 is the output gas flow towards purification system.The total pressure P sep,H2 of the oxygen separator can be obtained from the summation of all the partial pressures of the gases according to Dalton's law.The valve and the pipe resistances are modelled using R elements.Likewise, the fluidic, thermal and chemical phenomena for the OSV can be represented by Equations ( 38)- (40), respectively.For OSV from Figure 15b, using constitutive relation of C : C f l sep,O 2 and junction 0 12 , OSV pressure can be obtained as can be written as where P sep,O 2 = ρ w .g.L sep,O 2 is water pressure in OSV, L sep,O 2 is water level in OSV, ṁtank is the water mass flow from water tank to OSV. ṁaw,Osep is the water mass flow from the anode side of the stack to OSV and ṁsep,O 2 is the water mass flow out from the OSV to anode side recirculation circuit.The water mass flow ṁaw,Osep can be obtained using the product of known mass fraction of water x a,w and the mass flow of mixture, ṁa,Osep , obtained in Equation ( 23).
For OSV from Figure 15b, using constitutive relation of C : C th sep,O 2 and junction 0 13 , OSV temperature can be obtained as where T sep,O 2 temperature inside OSV, Ḣtank is the enthalpy flow from water tank to OSV, Ḣsepv is the enthalpy flow from the HSV to OSV through the separator valve, Ḣa,Osep is the enthalpy flow of the water and gas mixture from the anode side of the stack to OSV and Ḣsep,O 2 is the enthalpy flow out from the OSV to anode side recirculation circuit, Ḣgas,O 2 is the heat flow output from OSV towards oxygen circuit and QOenc is the rate of heat loss due to temperature gradient between OSV and system enclosure.For OSV from Figure 15b, using constitutive relation of C : C ch Hsep,i and junction 0 14 /0 15 / 0 16 , gases partial pressures can be obtained as

.4. Cooling and Recirculation Circuits
The desired temperature of the stack is automatically maintained by the cooling and the recirculation system by getting the information of hydrogen and oxygen production along with the information of heat flow in this system.The recirculation system has its own process and instrumentation diagram which is used to feed the water at desired rate to the stack.In the current model the recirculation of water is done in each side of the cell, that is, anode and cathode side using a recirculating controlled pump.However at anode side, a heat exchanger along with its own cooling system is also incorporated which mainly maintains the stack temperature by controlling the input water temperature for the stack.The BG sub-model of this system for anode and cathode sides in capsule forms are shown in Figures 16a,b, respectively.The modelling of the cooling unit of the heat exchanger is also shown in Figure 16a.In the sub-model presented in Figure 16a, different pumps are model as the modulated source of flows Ms f : ṁrec,a and Ms f : ṁcool , respectively, represent the anode side recirculating pump and the pump used for cooling of heat exchanger itself.The thermal capacity is modelled by the lump parameter C : C th rec,a .Thus, from the Figure 16a, using constitutive relation of C : C th rec,a and junction 0 19 , anode side recirculation temperature can be obtained as where the enthalpy flows in the anode side recirculation system Ḣsep,O 2 and Ḣrec,a are calculated from the Equation (35).The heat dissipation, QOrec,enc , due to temperature gradient between recirculation system and enclosure is modelled by BG generalized resistive element R : R Orec,enc and Qcool is the extracted rate of heat flow by the heat exchanger which is modelled by using resistive element R : R htex .In order to calculate the Qcool , the Number of Transfer Unit (NTU) technique is used for heat exchanger [51].Likewise, for the sub-model of cooling unit of the heat exchanger in Figure 16a, using constitutive relations of C : C th cool , C : C th cold and junctions 0 20 , 0 21 temperatures T cool , T cold can be obtained as where T cool , T cold and C th cool , C th cold are the temperatures and thermal heat capacities of cooling system of heat exchanger and cooling tank for the coolant, respectively.In sub-model of cooling unit in Figure 16a, Ms f : ṁcool represents the coolant mass flow which is used for the calculation of Ḣcool and Ḣcold using Equation ( 35) with temperatures T cool and T cold , respectively.
Likewise, using the sub-model of recirculation system in cathode side as shown in Figure 16b, using constitutive relation of C : C th rec,c and junction 0 23 , cathode side recirculation temperature can be obtained as where the enthalpy flows in the cathode side recirculation system Ḣsep,H2 and Ḣrec,c are calculated from the Equation (35).The lump parameter C : C th rec,c is used to model its thermal capacity towards the cathode side.

Hydrogen Purification Subsystem
This unit purifies the gas according to the preset hydrogen purity level required for the particular application.A dryer is used which removes the unwanted moisture from the gas by adsorbing the fraction of water content in it.The modelling of this unit requires the coupling of thermal, chemical and fluidic phenomena.The BG sub-model of the purification subsystem is shown in Figure 17.The lump parameters C : C th dry , R : R dry,enc , R : R dry , R : R exhaust , C : C ch dry , C : C ch ads and RS : R ads are, respectively, used to model the thermal capacity, thermal resistance, internal pneumatic resistance, exhaust resistance, dryer chemical capacity, water adsorption capacity and the adsorption resistance of the unit.Here, element RS : R ads is used to couple the chemical and thermal phenomena of the unit.According to Figure 17, using constitutive relation of C : C th dry and junction 0 28 , dryer temperature T dry can be obtained as where Ḣhpcv is the enthalpy flow coming through hydrogen pressure control valve, Ḣads = ṅH2O ads × ∆H ads is the enthalpy flow due to fraction of water adsorption with molar flow ṅH2O ads in the reaction, Ḣpuri enthalpy flow out from the purification system toward hydrogen production, Qdry,enc is the rate of heat loss due to temperature gradient between dryer unit and system enclosure and Ḣexha is the enthalpy flow when the exhaust valve is on (a e = 1).
Likewise, according to Figure 17, using constitutive relations of C : C t dry,i and junction 0 24 /0 26 for partial pressures p dry,i of oxygen and hydrogen gases, respectively, and C : C ch dry,j and junction 0 25 for partial pressure p dry,j of water vapour can be obtained as where ṁhpcv is the mass inflow coming through hydrogen pressure control valve, ṁpuri mass outflow from the purification system towards hydrogen production, x i and x j are the respective mass fractions.The capacities C ch dry,i and C ch dry,j can be similarly obtained from Equation (36) for the different species.The total pressure P dry of the dryer is the summation of partial pressure of the gases.

System Enclosure
The major subsystems of the electrolysis system are kept inside an enclosure which is designed according to the environmental requirements and the different functions of the subsystems.The enclosure temperature is maintained at desired temperature by designing the proper venting system and providing the fan/blower for its cooling.The BG sub-model of the system enclosure is shown in Figure 18.In this sub-model, the fan is modelled by the source of flow generalized BG element (MS f : ṁ f an ) to represent the mass flow of air and the atmospheric temperature is modelled as the source of effort generalized BG element (MSe : T atm ).The thermal capacity and the thermal resistance are modelled by the lumped parameters using C : C th enc and R : R enc .Thus, from Figure 18, using constitutive relation of C : C th enc and junction 0 29 , enclosure temperature T enc can be obtained as where Qsbs is the summation of all the rate of heat losses from the different subsystems to the enclosure, Ḣin, f an and Ḣout, f an are the rate of heat enthalpy in to the enclosure and out from it which depend on the specific heat (C P,air ) and the mass flow of air ( ṁ f an ).Also, Qenc,atm is the rate of heat loss to the atmosphere.

Efficiency of the PEM Electrolysis System
For analyzing the performance of the electrolysis system, efficiency is also the one of the most important parameters whose definition mainly depends on the different operating conditions and the system designs.Here, the efficiency is defined in the two levels: (i) cell/stack level (ii) system level including the auxiliaries.Here, it is also assumed that the cell operating voltage is always greater than the thermal-neutral voltage.Also, operating temperature is below the boiling point of water and water supplied in the liquid form.The real output of the system is only the useful hydrogen produced.Oxygen is not considered as a real output of the system, however it is also produced along with hydrogen.

Efficiency of Cell/Stack
For the dissociation of water in electrolysis process, a fixed amount of energy is required, which is equal to the summation of Gibbs free energy, ∆G R , and energy due to entropy change, T∆S R , and this summed energies is called the enthalpy, ∆H R , of the electrolysis reaction.Here ∆G R is the minimum amount of energy that must be supplied by the electrical input by assuming the rest of energy is thermally contributed by T∆S R term.The change in enthalpy (∆H R ) and Gibbs free energy (∆G R ) can be obtained from the chemical kinetics of the reaction as 286 kJ/mol and 237 kJ/mol, respectively, at standard temperature and pressure conditions, that equivalently represented in the form of standard electrode potential as 1.23 V and thermo-neutral potential as 1.48 V, respectively [52].In the calculation of standard electrode potential in reversible condition, it is assumed that all the thermal energy along with electrical energy are contributed in the electrolysis reaction based on Lower Heating Value (LHV) that does not include enthalpy of evaporation of water.However, in the calculation of thermo-neutral potential it is assumed that only electrical energy is contributed in the electrolysis reaction based on Higher Heating Value (HHV) that includes enthalpy of evaporation of water.Thus, the efficiency calculation of an electrolysis system is based on either the HHV or the LHV and it should be quoted while mentioning the efficiency of the system.Thus, the efficiency can be defined as the ratio of energy or power content of the hydrogen produced based on HHV or LHV and the input electrical energy or power to the system.Usually, HHV is preferred over LHV to calculate the efficiency of an electrolysis system supplied with liquid water, as the enthalpy of evaporation has to be provided by the process and it is represented as in Equation ( 49) [42].
where ηH 2 is molar flow of hydrogen and P elec is electrical DC input power to the electrolysis cell.Efficiency (ε diss ) for water dissociation can also be defined as the ratio of energy requirement in reversible condition (En rev ) to the energy requirement in irreversible condition (En irrev ).
Thus, based on Equation (50), other important efficiencies such as: voltage efficiency and current or Faradaic efficiency are also defined.Voltage efficiency can be obtained as the ratio of thermal-neutral voltage to the actual cell voltage at any operating condition, assuming the cell is always operating at a voltage above the thermal-neutral voltage.Thus, voltage efficiency is given as Voltage efficiency as defined in Equation ( 51) can be analytically represented in the form of any operating temperature T < 393.15K as [53,54] The voltage efficiency as presented in Equation ( 51) is valid with the assumption that the current supplied to the cell is fully converted into electro-chemical reaction of water dissociation.However, this is not the real situation in the cell operation due to the presence of some stray current in the cell and the unintended side reactions.Furthermore, gases permeation through PEM and successive recombination to water and leakage of gas leads to a loss of actual production of hydrogen.To consider this, the current or Faradaic efficiency ε curr is defined, which is the ratio of actual hydrogen produced ( ηH 2 ,act ) to the theoretically hydrogen produced ( ηH 2 ,th ) based on Faraday's law.
Thus, the overall cell efficiency can be obtained by multiplying voltage and current efficiency as

Efficiency of System Including the Auxiliaries
Electrolysis system includes the various other supporting subsystems and auxiliaries for the production of green hydrogen.Every auxiliaries and subsystems have their own efficiency which can be estimated or supplied by the manufacturer.It is assumed that an imaginary boundary is considered that includes the electrolysis cell and the supporting auxiliaries for defining the system efficiency.Here, solar and/or wind subsystems are excluded for defining the system efficiency.Thus, the efficiency of the complete system is ratio of the energy content of hydrogen produced and the total amount of energy consumed and it is given as where P syst is the total amount of energy or power consumed in the considered system, P elec denotes electrical DC input power and ε conv denotes converter efficiency, P pump denotes pump input power, P htex denotes input power to heat exchanger and P other denotes input power to other auxiliaries.However, there is always confusion and misunderstanding which heating value, that is, LHV or HHV, should be used in efficiency calculation.To eliminate this confusion, efficiency can be represented in terms of power consumed by the electrolysis system, which represents the amount of electrical energy consumed by the system to produce one kilogram (kWh/kg of H 2 ) or one normal cubic meter (kWh/Nm −3 of H 2 ) of hydrogen.It is noted that the efficiency calculation is constant when the electrical power input is constant.But, due to use of intermittent and time varying input sources, the value of efficiency changes with time.Thus, this needs the instantaneous calculation of efficiency [55].

Experimental Validation and Results
The developed model is adapted for laboratory-scale electrolyser system running on a multi-source energy platform.For implementation viewpoint, Matlab-Simulink environment is used for simulation and sensor data acquisition through a programmable logic controller.The experimental platform of PEM electrolyser system running on intermittent sources is shown in Figure 19.The platform consists of two photo-voltaic panels (200 W power per panel), permanent magnet type wind turbine (350-400 W power) and two batteries (act as buffer with 55 Ah capacity per battery) running a commercial 300 W single cell PEM electrolyser.The electrolysis cell is fed with water at the anode side only and the produced hydrogen (pressure range 1.4-10.7 bar) is stored in the metal hydride canister (H 2 bottle) of 760 standard liter capacity.
In order to reduce the complexity, only the model of the electrolyser is considered and the intermittent sources are considered as an electrical source with signal noise (to represent intermittent nature).For model validation, the key model parameters needs to be tuned with the experimental data.For example, the key parameters for the cell are identified using non-linear least square error optimization technique, which are shown in Figure 20.The bounds for the parameters values are taken from the literature for better estimation [8,24,42,54].Other parameter values such as recirculation pump flow, height and cross-sectional area of hydrogen separator are taken based on the system specification provided by the manufacturer.The polarisation (characteristic) curve from the actual measurements and the one estimated using the developed model is shown in Figure 21a.The mean absolute percentage error between simulation and experimental data is found to be 4.8% which is reasonable and well within the acceptable limit.Figure 21b shows the predicted characteristic curve at different temperatures.It can be seen that with the increase in the cell temperature the required cell voltage for electrolysis decreases.The predicted performance is in good agreement with the expected behaviour of the electrolyser found in the literature [18,24,42].Figure 22a shows the contribution of different overvoltages in the overall cell voltage.It can be seen that with increase in the current the ohmic losses increases linearly.The contribution of the activation losses in the overall cell voltage remains almost constant.The losses due to mass transportation is negligible for small electrolysis cells and hence it is neglected.Figure 22b shows the variation of the cell efficiency with respect to the cell voltage.The efficiency of the cell is constant till the cell voltage is less than thermo-neutral voltage.Once the cell voltage exceeds the thermo-neutral voltage, the efficiency decreases continuously.However, the efficiency at the system level can be calculated as discussed in Section 3.3.2,but, due to the limited number of the sensors placed in the experimental platform, the power consumed by different auxiliaries cannot be measured.Hence, the system level efficiency is not shown.Figure 22c shows the prediction of the temperature evolution of the cell from the thermal sub-model.The temperature starts rising from the ambient temperature and attains a constant value when the steady state of the operation is reached.

Number of Iterations
To predict the performance of the electrolyser powered by intermittent sources, a simulation of the electrolyser directly powered by the combination of solar and wind energy is performed.Figure 23a shows the total power fed to the cell by the converter and this figure also shows the individual contribution of each intermittent source.Figure 23b shows the input voltage to the electrolyser model with the objective of achieving high rate of hydrogen production.The current drawn by the cell, the hydrogen production and the cell efficiency are shown in the Figures 24a-c respectively.It can be seen that the efficiency of the cell is about 40%.This is due to the fact that the current electrolyser is designed to run at particular set point with the objective of delivering continuous hydrogen flow at higher rates.This is achieved by running the electrolyser at higher voltage (around 4 to 5 V) which results in reduced efficiency.In order to improve both efficiency and hydrogen production rate, the number of electrolysis cells should be increased.
In the current study, intermittent power is directly taken as an input for the model to predict the electrolyser's performance.However, in real electrolysis system running on RES, the controller with buffer power source (battery) with operating mode management is always present [43].This is important to maintain the constant input power to the cell in reference to the set point for stable operation and to avoid faster degradation of the cell membrane.

Conclusions
A generic multi-physics dynamical model of PEM electrolyser (including the auxiliaries) is developed using the BG approach.The developed model consists of sub-models (capsules) that can be easily modified to fit any PEM electrolyser.Because of the generic nature, it can be easily sized and scaled as per the requirement.The model incorporates the efficiency of stack and auxiliaries in order to estimate the performance of the electrolyser.Also, the system dynamics and the different thermal losses can be predicted which cannot be done in the real process.With the increase in current density, efficiency of the electrolyser decreases due to the reduction in the performance and corresponding heat dissipation.Final implementation and simulations are performed on MATLAB-Simulink.The validation the model has been done against the experimental data single electrolyser running on (solar and wind).accuracy of the validation corresponds to the industrial requirements as the error is found to be 4.83%.However, the model is capable of capturing the effect of system dynamics with the parameters variation; but, due to the industrial constraints on the parameters identification, the model validation in different dynamic configurations is not performed.
The developed model is economical in nature.Due to the use of specific formalism (causal and structural properties of the BG), the developed model can not be used only for dynamic analysis but can also be used for control, diagnosis and prognosis studies.Therefore, this work is a foundation stone for the future work for developing model-based diagnosis and prognosis for PEM electrolysers.Since, the approach is modular and generic, the developed model can be extended to other types of electrolysers.

Figure 4 .
Figure 4. Representation of vector bond (a) with a small ring around the power bond, (b) with multi-bonds and (c) with separate power bond for different energy domains.

Figure 6 .
Figure 6.Graphical user interface of the model builder.

Figure 8 ,
line with half arrow represents the power bond, where yellow bonds show the electrical energy (current and voltage variables), red bonds show the thermal energy (heat flow and temperature variables), blue bonds show the fluidic energy (mass flow and pressure variables), green bonds show the thermal-fluidic energy coupling (heat flow, temperature and mass flow, pressure variables) towards hydrogen side.Likewise, orange bonds show the thermal-fluidic energy coupling towards the oxygen side.

Figure 8 .
Figure 8.The word BG model of the electrolysis system.

Figure 10 .
Figure 10.The BG chemical-fluidic sub-model in capsule form.

Figure 11 .
Figure 11.The BG thermal sub-model in capsule form.

Figure 12 .
Figure 12.BG sub-model of fluidic and mass transfer phenomenon in capsule form.

Figure 13 .
Figure 13.The global model of the stack of the PEM electrolysis system.

Figure 16 .
Figure 16.BG recirculation sub-models (a) anode side including cooling circuit and (b) cathode side.

Figure 17 .
Figure 17.BG sub-model purification subsystem in capsule form.

Figure 18 .
Figure18.BG sub-model of system enclosure in capsule form.

Figure 19 .
Figure 19.Experimental platform of the PEM electrolyser running on intermittent sources.

Figure 20 .
Figure 20.Parameter estimation for the sub-model of the stack.

Figure 22 .
Figure 22.(a) Contribution of overvoltages, (b) efficiency of electrolysis cell and (c) temperature evolution of cell.

Figure 23 .
Figure 23.(a) Powers consumed by the cell for 12 h run and (b) Input voltage for the cell running on intermittent sources.

Figure
Figure Current drawn by the during of operation, (b) and (c) corresponding cell efficiency the cell.

1 C f l−th sep,H 2 Field− 2 J 1 R 1 R 1 R 1 T 3 V
Dryer's chemical capacitance of the purification unit for the i th species, mol•Pa −1 C th dry Dryer's thermal capacitance of the purification unit, J•K −1 Chemical capacitance of hydrogen separator for the i th species, mol•Pa −1 Chemical capacitance of oxygen separator for the i th species, mol•Pa −1 C th rec,k Thermal capacitance of recirculation circuit (anode/cathode side), J•K −capacitance element representing fluidic capacitance and thermal capacitance of hydrogen separator C f l−th sep,O 2 Field capacitance element representing fluidic capacitance and thermal capacitance of oxygen separator C stack Thermal capacitance of the stack, J•K −1 D i Parameter for diffusion, m 2 •s −1 d M Ratio of length to the cross-sectional area of the membrane, m −1 E act,k Activation overvoltages for k th electrode, V 0,k Standard current exchange density for k th electrode, A•m −2 L M Length of the membrane, m L sep,i Water level in Separators (HSV and OSV), m L sep,i Water level of the i th separator, m M i Molar mass for i th species, kg mol −1 Partial pressure of i th species, Pa R c Coupling element for fluidic flow to thermal flow R act,k Non linear activation resistance for k th electrode R ads Coupling resistance of adsorbed water molar flow and the enthalpy flow towards dryer R di f f ,i Diffusion resistance of the i th species, Pa•s•kg −1 R dry,enc Thermal resistance between purification unit and enclosure, K•s•J −1 R dry Internal pneumatic resistance of the dryer, Pa•s•kg −1 R enc Thermal resistance between the enclosure and the atmosphere, K•s•J −1 R exhaust Pneumatic resistance of the exhaust valve, Pa•s•kg −Hrec,enc resistance between hydrogen recirculation circuit and enclosure, K•s•J −Hsep,enc Thermal resistance between hydrogen separator and enclosure, K•s•J −1 R htex Thermal resistance of heat exchanger, K•s•J −1 R hyst,k Internal fluidic resistance of the stack at th k th electrode side, Pa•s•kg −Lrec,k Hydraulic resistance representing leakage in recirculation circuit (anode/cathode side), Pa•s•kg −1 R ohm Total ohmic resistance of the cell, Ω R Orec,enc Thermal resistance between oxygen recirculation circuit and enclosure, K•s•J −1 R Osep,enc Thermal resistance between oxygen separator and enclosure, K•s•J −1 R others Ohmic resistance of the cell except membrane, Ω R sep,hc Pneumatic resistance between hydrogen separator and hydrogen circuit, Pa•s•kg −1 R sep,oc Pneumatic resistance between oxygen separator and oxygen circuit, Pa•s•kg −1 R sepv Hydraulic resistance of the separator valve, Pa•s•kg −1 R stack Thermal resistance of the stack, K•s•J −1 R tank Hydraulic resistance between tank and oxygen Separator, Pa•s•kg −Temperature, K V Hsep,i Volume of the i th species in HSV, m Osep,i Volume of the i th species in OSV, m 3 x i Mass fraction of the i th species µ i Chemical potential of the i th species, J•kg −1 A i Chemical affinity of the i th species, J•mol −1 C ch ads Water adsorption capacity of the purification unit, mol•Pa −1 F Faraday's constant, C•mol −1 R Ideal gas constant, J•mol −1 •K −1 H 2 O → 2H + + 1 2 O 2 + 2e − ; Cathode : 2H + + 2e − → H 2 ; Overall : H 2 O(l) → H 2 (g) + Anode :

Table 1 .
Flow and effort for various energy domains. 3 Charge transfer or symmetry factor coefficients for k th electrode β Transformer coefficient of the converter ∆G R Gibb's free energy of water dissociation reaction, J•mol −1 ∆H R Enthalpy Change of water dissociation reaction, J•mol −1 ∆S R Entropy change of water dissociation reaction, J•mol −1 •K −1 ξ Rate of reaction flow, mol•s −1 Ḣ Enthalpy rate, J•s −1 ṁi Mass flow rate of the i th species, kg•s −1 ṅi Gas mass flow rate for i th species, kg•s −1 Q Rate of heat J•s −1 Matter storage capacity of the i th species (i = H 2 , O 2 , H 2 O), kg 2 •J −1