Document [original]

This version is available at https://doi.org/10.14279/depositonce-11685

This work is licensed under a CC BY-NC-ND 4.0 License (Creative

Commons Attribution-NonCommercial-NoDerivatives 4.0 International). For

more information see https://creativecommons.org/licenses/by-nc-nd/4.0/.

Weigert, J., Hoffmann, C., Esche, E., Fischer, P., & Repke, J.-U. (2021). Towards demand-side

management of the chlor-alkali electrolysis: dynamic modeling and model validation. Computers &

Chemical Engineering, 107287.

https://doi.org/10.1016/j.compchemeng.2021.107287

Joris Weigert, Christian Hoffmann, Erik Esche, Peter Fischer, Jens-Uwe

Repke

Towards demand-side management of the

chlor-alkali electrolysis: dynamic modeling

and model validation

Accepted manuscript (Postprint)Journal article |

Towards demand-side management of the chlor-alkali electrolysis:

dynamic modeling and model validation

Joris Weigerta,∗, Christian Hoffmanna, Erik Eschea, Peter Fischerb, Jens-Uwe Repkea

aTechnische Universit¨at Berlin, Process Dynamics and Operations Group, Sekr. KWT 9, Straße des

17. Juni 135, Berlin 10623, Germany

bVestolit GmbH, Paul-Baumann-Str. 1, 45772 Marl, Germany

Abstract

The chlor-alkali electrolysis promises advantageous application of demand-side management

and several research groups have proposed dynamic models to obtain optimal trajectories

for such applications. However, no dynamic model of those proposed so far has yet been

validated with real, industrial plant data. This is highly important to determine trajectories,

which are indeed feasible. This contribution addresses this issue by proposing a dynamic

model of the chlor-alkali electrolysis that contains those variables relevant for dynamic op-

eration, especially temperature and composition of the electrolytes. This model is validated

with real plant data from an electrolysis located in Marl, Germany, for two scenarios: (1)

several small load changes of around 4 MW and (2) a large load change of more than 50 MW.

Keywords: Chlor-alkali electrolysis, Demand response, Demand-side management,

Flexible operation, Dynamic Modeling

1. Introduction

The increasing share of renewable energy sources and their natural fluctuations pose

significant challenges for grid stability but also offer emerging opportunities as a result of

the likewise fluctuating electricity prices. Demand-side management or demand response

(DR) is therefore a viable possibility to stabilize the power grid (Paterakis et al., 2017). The

application of DR in combination with the resulting highly volatile electricity markets have

encouraged flexible operation of chemical plants and thus their active participation in DR,

e.g., for reverse osmosis plants (Sassi and Mujtaba, 2013), steel plants (Castro et al., 2020),

standard distillation columns (Dowling and Zavala, 2018), and air separation units (Tsay

et al., 2019; Zhao et al., 2019). Electrochemical processes are of particular interest as they

consume a large amount of electrical energy, e.g., in Germany (Klaucke et al., 2017) or in

∗Corresponding author

Email addresses: [email protected] (Joris Weigert), [email protected]

(Christian Hoffmann), [email protected] (Erik Esche), [email protected] (Peter

Fischer), [email protected] (Jens-Uwe Repke)

Preprint submitted to Computers & Chemical Engineering March 17, 2021

Accepted manuscript of: Weigert, J., Hoffmann, C., Esche, E., Fischer, P., & Repke, J.-U. (2021). Towards demand-

side management of the chlor-alkali electrolysis: dynamic modeling and model validation. Computers & Chemical

Engineering, 107287. https://doi.org/10.1016/j.compchemeng.2021.107287

http://creativecommons.org/licenses/by-nc-nd/4.0/

the United States (Worrell et al., 2000), and thus offer an important lever to enable more

flexible operation of chemical plants (Paulus and Borggrefe, 2011; Ausfelder et al., 2018).

The chlor-alkali electrolysis (CAE) plays a decisive role due to its large installed capacity

(Klaucke et al., 2017). For the 1350 MW installed in Germany (Ausfelder et al., 2018) alone,

a flexibility potential of 169 MW for load reduction is estimated (Klaucke et al., 2020).

As a result, the chlor-alkali electrolysis has also been studied extensively concerning DR

applications: Babu and Ashok (2008) solved a mixed-integer nonlinear scheduling problem

to demonstrate the economic benefit of DR wheras Otashu and Baldea (2019) addressed

the nonlinear dynamics during flexible operation and showed that careful control of pro-

cess variables during flexible operation is highly important. They later used their dynamic

model to determine optimal operation under varying electricity prices and also incorporated

frequency-balancing in their framework (Otashu and Baldea, 2020). Recently, Simkoff and

Baldea (2020) used the model by Otashu and Baldea (2019) to identify a Hammerstein-

Wiener model and used this surrogate model to allocate load for the day-ahead and real-

time markets. Richstein and Hosseinioun (2020) studied the impact of various network

tariffs and regulations on the flexibility of the chlor-alkali electrolysis and Baetens et al.

(2020) addressed the implications of exogenous uncertainties on the optimal operation of

the chlor-alkali electrolysis. Another approach to enable flexibility of the CAE was investi-

gated by Br´ee et al. (2019, 2020) and Roh et al. (2019). They used a quasi-stationary model

of the CAE, which described the operation with both a standard cathode and an oxygen

depolarized cathode.

Given the economic benefits of short-term electricity markets (Babu and Ashok, 2008),

the CAE should preferably operate under continuously changing load based on market sig-

nals while maintaining the desired temperature and concentrations of anolyte and catholyte.

Otherwise, especially the membrane between catholyte and anolyte might be damaged

(Wellington, 1992; O’Brien et al., 2007). However, most studies on DR potential and optimal

operation either applied simplified steady-state models (Babu and Ashok, 2008; Richstein

and Hosseinioun, 2020; Baetens et al., 2020) or simplified quasi-stationary (Roh et al., 2019;

Br´ee et al., 2019, 2020). While such simplified models are certainly useful to provide solu-

tions to scheduling problems over larger time horizons due to their lower complexity, they

may not provide more insight into the phenomena within the electrolysis cells. Dynamic

first-principle models have rarely been used yet, but are required to

1. demonstrate that the requirements regarding temperature and cell concentrations can

be met during flexible operation

2. obtain optimal control trajectories for DR scenarios or

3. verify load trajectories that were obtained with a simplified model of the CAE by

solving a scheduling problem

Dynamic models of the CAE have, for example, been proposed by Wang et al. (2014),

Budiarto et al. (2017), and Otashu and Baldea (2019). Wang et al. (2014) developed a dy-

namic model to describe the CAE in combination with a wind farm and fuel cells. However,

only the relationship between cell temperature and power consumption/current density was

shown. Budiarto et al. (2017) developed a dynamic CAE model on laboratory scale, which

was capable of describing the dynamic behavior of the product streams and the power con-

sumption depending on the current density. However, some fundamental effects necessary

to describe the composition of the electrolytes, such as electroosmosis and the evaporation

of water, were not taken into account. Furthermore, no energy balance was included in

the proposed model, which is highly relevant under dynamic operation. A much more de-

tailed description of the CAE and, in particular, of the phenomena in the electrolysis cells

was given by Otashu and Baldea (2019). They modeled the process dynamics of the CAE

at industrial scale in the context of DR. However, the load-dependent description of mass

transport through the membrane is missing as they assumed constant diffusion coefficients

for the components migrating through the membrane, which does not allow for a precise

description of the anolyte and catholyte composition. Furthermore, no validation of the

model on the basis of real plant data was carried out and the accuracy of their detailed

model could thus not be quantified.

To overcome these shortcomings, we present a dynamic model of a real chlor-alkali elec-

trolysis plant, which describes all variables relevant for dynamic operation. In this context,

we discuss the occurring phenomena within such a plant in detail and suggest a modeling

approach for each individual phenomenon. In addition, to the authors’ knowledge, this is

the first time that a dynamic CAE model has been validated with actual, industrial pro-

cess data. The model validation carried out in this contribution is thus a major novelty

and might make the available optimal control and scheduling models more reliable in the

future. The process data is taken from a CAE plant operated by the Vestolit GmbH in Marl

(Germany). Throughout this work, we focus on membrane processes as this process concept

dominates modern chlorine production (Br´ee et al., 2019; Klaucke et al., 2020).

We also note that the CAE cannot be operated flexibly by itself as storing chlorine is

avoided in practice due to safety concerns (Br´ee et al., 2019). In this context, Klaucke et al.

(2020) analyzed the chlorine value chain to determine processes with the largest flexibility

potential. One such process is the production of polyvinyl chloride, in which the intermediate

1,2-dichloroethane can be stored easily. A dynamic model of this subsequent process was

proposed and validated by us in the past (Hoffmann et al., 2020b).

In the next section, the basics of chlor-alkali electrolysis and the modeled plant are

presented. In Section 3, the dynamic model is developed. Section 4 contains the model

validation using real plant data.

2. Process description

The chlor-alkali electrolysis produces chlorine (Cl2), hydrogen (H2), and caustic soda

(NaOH) from sodium chloride (NaCl) brine using electrical power. A schematic represen-

tation of the electrolysis cell and the relevant phenomena modeled in Section 3 is shown in

Figure 1. The cell consists of two compartments: the anode and the cathode, which are sep-

arated by a permselective ion-exchange membrane. At the anode, chloride ions are oxidized

and chlorine is produced (Equation (1)). The inner current of the cell is ideally generated

by positive sodium ions migrating through the cation exchange membrane together with

a hydrate shell of approximately four water molecules. At the cathode, hydronium ions

are reduced to hydrogen and water (Equation (4)) and the remaining hydroxide ions from

the dissociation of water (Equation (5)) form caustic soda with the migrating sodium ions

(Equation (6) and Equation (7)). In addition, water leaves the cell along with chlorine or

hydrogen due to evaporation effects.

In the ideal case, only sodium ions and water migrate from anolyte to catholyte. In the

technical process, however, the undesired back-migration of hydroxide ions plays a significant

role regarding the current efficiency of chlorine and caustic soda production. These back-

migrated ions form oxygen instead of chlorine at the anode in an electrochemically preferred

side reaction. This formation of oxygen can be suppressed by acidification of the anolyte.

Nevertheless, a loss of current efficiency for the production of caustic soda occurs (O’Brien

et al., 2007).

Oxidation Reduction

Autoprotolysis

enriched

depleted

Ion exchange membrane

enriched

diluted

Anolyte Catholyte

Figure 1: Schematic representation of the chlor-alkali electrolysis cell. Note that the autoprotolysis of water

also takes place in the anolyte compartment.

Anode

Anode reaction: 2 Cl–Cl2+ 2 e–(1)

Dissociation of sodium chloride: 2 NaCl 2 Na++ 2 Cl–(2)

Overall reaction: 2 NaCl 2 Na++ Cl2+ 2 e–(3)

Cathode

Cathode reaction: 2 H3O++ 2 e–H2+ 2 H2O (4)

Dissociation of water: 4 H2O 2 H3O++ 2 OH–(5)

Overall reaction: 2 H2O + 2 e–H2+ 2 OH–(6)

Overall reaction of chlor-alkali electrolysis

2 NaCl + 2 H2O 2 NaOH + Cl2+ H2(7)

Figure 2 shows a simplified flowsheet of the real-life CAE plant, which is the basis for the

proposed model introduced in Section 3. The electrolysis cell (framed by the dotted line) is

described by two continuously stirred-tank reactors. Reactor 1 (left) represents the anode

and reactor 2 (right) represents the cathode compartment. The membrane only allows for

the components water, sodium ions, and hydroxide ions to flow from one side to the other.

To control the cell temperature, the feed temperature of both compartments is adjusted with

the heat exchangers HE1 and HE2. Downstream to the outlets of both compartments, buffer

tanks (T1 and T2) compensate for strong fluctuations in the outlet streams. In addition, the

catholyte recycle is modeled as it strongly influences the dynamic behavior of the CAE. The

recycle is also required for temperature control of the cell and for control of the catholyte

composition (set points can be found in Table 2). The anolyte recycle is disregarded in the

model since we assume that fresh brine is always supplied with constant composition and

that the feed flow is independent of the outgoing anolyte flow.

Five controllers are employed to control the process and to mimic the behavior of the

real plant: The cell temperature is controlled using the inlet temperature of the catholyte.

Anolyte and catholyte composition are controlled with anolyte feed and water dilution flow

in the catholyte recycle, respectively. The product flow of caustic soda and the level of

the buffer tank T1 are controlled using the outlet flow of the catholyte buffer tank and the

anolyte buffer tank, respectively, as manipulated variables. A detailed description of the

employed controllers can be found in Section 3.5.

3. Process modeling

To model the process displayed in Figure 2, the following assumptions are made:

Assumption 1. The diffusion of chloride ions to the catholyte compartment can be neglected

(O’Brien et al., 2007). In addition, the diffusion of chlorine and hydrogen through the

membrane is negligible.

Assumption 2. Anolyte and catholyte are ideally mixed due to the continuous formation

of chlorine and hydrogen gas bubbles in the cell (Kreysa and Wendt, 2010).

OH-

Na+

H2O

prod

water

Cl2(g), H2O(g) H2(g), H2O(g)

in in

out tank, out

out

tank, out

Anolyte Catholyte

CAE cell

QC QC

HE1

HE2

T1 T2

Figure 2: Flow diagram of the modeled chlor-alkali electrolysis plant.

Assumption 3. The electroylsis cells behave like a bubble column. This indicates that the

gas volume fraction is proportional to the volume flow of the gas (Shah et al., 1982). Hence,

due to Faraday’s law a decreasing current density leads to a reduction of the gas volume as

less gas is formed at the electrodes (Kreysa and Wendt, 2010).

Assumption 4. The membrane does not cause any significant heat transfer resistance given

that it is only a few hundred micrometers thick (O’Brien et al., 2007).

Assumption 5. The internal energy equals the enthalpy as the product of pressure and

temperature is small compared to the enthalpy.

Assumption 6. The product gases chlorine and hydrogen leave the cell saturated with water

(Kreysa and Wendt, 2010).

Assumption 7. Heat loss is caused by free convection and governed by the heat transfer

between cell surface and the surrounding air.

Assumption 8. The correlation between cell voltage and current density can be described

linearly in the relevant operating range of the CAE (O’Brien et al., 2007).

Assumption 9. The influences of the electrolyte compositions and cell temperature on the

cell voltage can be neglected as these vary only within a small operating window.

Assumption 10. Oxygen formation at the anode can be neglected as it is suppressed by

acidification (O’Brien et al., 2007).

Assumption 11. The volume of the buffer tanks is small compared to the inlet flows coming

from the electrolyzer and the temperature of the flows remains constant between the outlet of

the electrolyzer and the inlet of the buffer tanks.

Based on these assumptions, the dynamic mass and energy balances for anolyte and

catholyte are presented. Afterwards, the electrochemical phenomena and the buffer tanks

are described. Finally, the modeling equations for the employed controllers are discussed.

Due to Assumption 1, we do not consider chlorine ions in the model equations of the catholyte

compartment. To describe the occurring components, we define component sets for each cell

compartment. The set for the anolyte compartment is

Can ={Na+,Cl–,OH–,H3O+,H2O}(8)

and the set for the catholyte compartment is

Ccat ={Na+,OH–,H3O+,H2O}(9)

3.1. Mass balances

The mass balances are formulated separately for the anolyte and the catholyte compart-

ment of the electrolysis cell, i.e., the control volumes in Figure 2. As a result of Assumption 2,

the mass fractions in both compartments are equal to the mass fractions in the respective

outlet streams.

3.1.1. Anolyte compartment

The dynamic component mass balances of the anolyte compartment for the considered

components are given by

dHUan

Na+

dt =˙

Van

in ·ρan

in ·wan

in,Na+−˙

Van

out ·ρan ·wan

Na+

| {z }

Convective flow

−˙nmem

Na+·Acell ·MNa+

| {z }

migration Na+

(10)

dHUan

Cl–

dt =˙

Van

in ·ρan

in ·wan

in,Cl–−˙

Van

out ·ρan ·wan

Cl–−2·ranode ·Acell ·MCl–

| {z }

anode reaction

(11)

dHUan

OH–

dt =˙

Van

in ·ρan

in ·wan

in,OH–−˙

Van

out ·ρan ·wan

OH–+ ˙nmem

OH–·Acell ·MOH–

| {z }

back-migration OH–

+Ran

AW ·MOH–

| {z }

autoprotolysis

(12)

dHUan

H3O+

dt =˙

Van

in ·ρan

in ·wan

in,H3O+−˙

Van

out ·ρan ·wan

H3O++Ran

AW ·MH3O+(13)

dHUan

H2O

dt =˙

Van

in ·ρan

in ·wan

in,H2O−˙

Van

out ·ρan ·wan

H2O−˙nmem

H2O·Acell ·MH2O

| {z }

electroosmosis

−Fan,vapor

H2O

| {z }

evaporation

−2·Ran

AW ·MH2O(14)

Therein, ˙

Van

in ,wan

in,c and ˙

Van

out,wan

care the in- and outgoing liquid volume flows and the

respective mass fraction of the considered components. Correlations for calculating the

required electrolyte densities ρan and ρcat and the inlet densities ρan

in and ρcat

in are given

in Appendix A. The computation of both flows of water vapor, Fan,vapor

H2Oand Fcat,vapor

H2O, is

explained in Section 3.2.2. The rates of the electrode reactions, ranode and rcathode, as well

as the material flows of sodium ions and water through the membrane, ˙nmem

Na+and ˙nmem

H2O, are

explained in Section 3.3.1. The extent of the autoprotolysis of water, Ran

AW, is discussed in

Section 3.3.3.

The summation of the component holdups yields the total holdup in the anolyte com-

partment: X

Can

HUan

c=ρan ·Van (15)

As mentioned in Assumption 3, the gas volume in both cell compartments and thus the

liquid volumes Van and Vcat depend on the current density. To mimic the behavior of the

real electrolysis cells, both liquid volumes need to reflect this dependency. We avoid a very

stiff system (would be obtained by assuming an overall mass balance at steady-state) and

a high-index differential-algebraic equation system (would be obtained by just fixing the

liquid volumes of both compartments to their values at the corresponding current density)

by using artificial controllers for this purpose, which determine the liquid outlets of both cell

compartments ˙

Van

out and ˙

Vcat

out . The model equations describing the set points of both liquid

volumes as a function of the current density can be found in Section 3.3.2 and the artificial

controllers are presented in Section 3.5.

The mass fractions are computed from the component holdups by

wan

c=HUan

PCan HUan

∀c∈Can (16)

3.1.2. Catholyte compartment

The mass balances for the catholyte compartment read:

dHUcat

Na+

dt =˙

Vcat

in ·ρcat

in ·wcat

in,Na+−˙

Vcat

out ·ρcat ·wcat

Na++ ˙nmem

Na+·Acell ·MNa+(17)

dHUcat

OH–

dt =˙

Vcat

in ·ρcat

in ·wcat

in,OH–−˙

Vcat

out ·ρcat ·wcat

OH–+rcat ·Acell ·MOH–

−˙nmem

OH–·Acell ·MOH–+Rcat

AW ·MOH–(18)

dHUcat

H3O+

dt =˙

Vcat

in ·ρcat

in ·wcat

in,H3O+−˙

Vcat

out ·ρcat ·wcat

H3O+−rcat ·Acell ·MH3O+

+Rcat

AW ·MH3O+(19)

dHUcat

H2O

dt =˙

Vcat

in ·ρcat

in ·wcat

in,H2O−˙

Vcat

out ·ρcat ·wcat

H2O+ ˙nmem

H2O·Acell ·MH2O

−Fcat,vapor

H2O−2·Rcat

AW ·MH2O(20)

No component mass balance is required for chloride ions in the catholyte compartment due to

Assumption 1. The variables ˙

Vcat

in ,wcat

in,c,˙

Vcat

out , and wcat

care the in- and outgoing liquid flows

and the respective mass fractions of the considered components. Additionally, expressions

analogue to Equation (15) and Equation (16) are included. Here, only the components in

Ccat have to be taken into account.

3.2. Energy balance

Due to Assumption 2, Assumption 4, and Assumption 5, an identical temperature pre-

vails in both cell compartments and the internal energy of the cell can be represented by

the cell’s enthalpy Hcell. Thus, the cell’s dynamic energy balance reads:

dHcell

dt =˙

Hcell

in −˙

Hcell

out −˙

Qreac −˙

Hevap −˙

Qloss +Pcell,(21)

which is analogous to the formulation given by Otashu and Baldea (2019). Here, ˙

Hcell

in and

Hcell

out are the in- and outgoing overall enthalpy streams. The variable ˙

Qreac describes the heat

of reaction caused by Equation (7). Furthermore, ˙

Hevap considers the required enthalpy to

evaporate water in the cell compartments, ˙

Qloss describes the heat loss, and Pcell is the input

of electrical energy. All these variables are defined in the following sections.

3.2.1. Enthalpy flows and heat of reaction

The total enthalpy of the cell is given by

Hcell =Van ·ρan ·Can,liquid

p,NaCl ·(Tcell −Tref) + Vcat ·ρcat ·Ccat,liquid

p,NaOH ·(Tcell −Tref) (22)

while the enthalpy flows of the inlet and outlet are determined by

Hcell

in =˙

Van

in ·ρan

in ·Can,liquid

p,in,NaCl ·(Tan

in −Tref) + ˙

Vcat

in ·ρcat

in ·Ccat,liquid

p,in,NaOH ·(Tcat

in −Tref) (23)

and

Hcell

out =˙

Van

out ·ρan ·Can,liquid

p,NaCl ·(Tcell −Tref) + ˙

Vcat

out ·ρcat ·Ccat,liquid

p,NaOH ·(Tcell −Tref)

+˙

Nan

Cl2·MCl2·Can

p,Cl2·(Tcell −Tref) + ˙

Ncat

H2·MH2·Ccat

p,H2·(Tcell −Tref) (24)

Therein, Tcell represents the cell temperature and Tref is the reference temperature of the

enthalpies of formation (25 °C). The product streams of hydrogen ( ˙

Ncat

H2) and chlorine ( ˙

Nan

Cl2)

are defined in Equation (37) and Equation (38), respectively. Since the operating range with

respect to temperatures of the modeled plant is strongly limited (Wellington, 1992; O’Brien

et al., 2007), temperature-independent heat capacities are assumed. Their composition-

dependent expressions for all streams are given in Appendix B. Moreover, no reference

enthalpies need to be considered in Equation (22), (23), and (24) as the heat of reaction is

explicitly considered in the energy balance.

The heat flow generated by the overall cell reaction is defined as (Kreysa and Wendt,

2010): ˙

Qreac =˙

Nan

Cl2·∆hreac,ref (25)

The enthalpy of reaction ∆hreac,ref is given by the enthalpies of formation at reference state

of the involved components (Gmehling et al., 2019):

∆hreac,ref =νNaCl ·∆hf,ref,NaCl(aq) +νH2O·∆hf,ref,H2O(l) +νNaOH ·∆hf,ref,NaOH(aq)

+νCl2·∆hf,ref,Cl2(g) +νH2·∆hf,ref,H2(g) (26)

Using the enthalpies of formation given by Masterton et al. (1983) and Chase (1998),

∆hreac,ref is set to 446.8 kJ mol−1.

3.2.2. Evaporation of water

Due to the high operating temperature (ca. 90 °C) of the cell, a significant amount of

water from the anolyte and catholyte is evaporated and leaves the cell with the product gases

chlorine and hydrogen, respectively. Due to Assumption 6, the partial pressure of water in

the product gases is equal to the vapor pressure in the respective electrolyte solution. Since

both electrolytes are highly concentrated, the vapor pressure is lower than that of pure

water. (O’Brien et al., 2007; Kreysa and Wendt, 2010). The associated enthalpy stream is

Hevap =Fan,vapor

H2O+Fcat,vapor

H2O·∆hevap,H2O(27)

The calculation of the enthalpy of evaporation of water for both compartments is also de-

scribed in Appendix B. The mole fractions of water in the product gases of both cell com-

partments are calculated by assuming that the product gases chlorine and hydrogen leave

the cell saturated with water (Kreysa and Wendt, 2010):

yan,vapor

H2O·Pan ·HUan

H2O

MH2O

+HUan

Na+

MNa+=HUan

H2O

MH2O

·γan

H2O·PLV

H2O(28)

Therein, Pan is the pressure of the anolyte compartment. Based on plant characteristics, its

value is fixed to 1.28 bar. PLV

H2Oand γan

H2Oare the vapor pressure and the activity coefficient

of water, respectively. The vapor pressure of water in anolyte and catholyte is described as

PLV

H2O=aLV ·Tcell +bLV (29)

in which the parameters aLV and bLV were fitted to data from Green and Perry (2008) in the

temperature range 80 to 95 °C. They are given in Table 1.

0.16 0.17 0.18 0.19 0.2

0.9

0.905

0.91

0.915

0.92

0.925

0.93

0.935 eNRTL 80 °C

eNRTL 85 °C

eNRTL 90 °C

simplified model

0.16 0.18 0.2 0.22 0.24 0.26 0.28

3.1

3.2

3.3

3.4

3.5 Tanner and Lamb (1978) 70 °C

Tanner and Lamb (1978) 80 °C

Tanner and Lamb (1978) 90 °C

simplified model

Figure 3: Temperature influence on two properties of the NaCl solution in the relevant composition range:

(a) comparison of computed activity coefficient of water for eNRTL at three different temperatures and the

simplified expression from Equation (30); (b) comparison of heat capacity of aqueous NaCl for correlation

and simplified expression from Equation (B.1).

Since the electrolytes are highly concentrated solutions of both aqueous brine and caus-

tic soda, they cannot be assumed ideal. Instead, the electrolyte non-random two-liquid

(eNRTL) model (Chen and Evans, 1986) is employed to calculate the activity coefficient

of water. The required eNRTL parameters were retrieved from the database APV 100

ENRTL-RK in Aspen Properties V10. We confirmed that eNRTL in combination with the

retrieved parameters can accurately predict vapor-liquid equilibria of aqueous solutions of

both sodium chloride and sodium hydroxide by comparing vapor pressure measurements

at various temperatures and compositions of both solutions, taken from Green and Perry

(2008), with predictions from eNRTL. It was found that the temperature dependence of

the activity coefficients of water in NaCl and NaOH solutions is negligible in the relevant

temperature range of the CAE. This is shown in Figure 3 (a) for the activity coefficient of

water in aqueous NaCl. Therefore, a simplified expression of the activity coefficient of water

in a NaCl solution was set up to increase the robustness of the process model by removing

the highly nonlinear eNRTL equations from the model. This simplified expression is given

γan

H2O=aan

act ·

HUan

Na+

MNa+

HUan

Na+

MNa++HUan

Cl–

MCl–+HUan

H2O

MH2O

+ban

act (30)

The parameters aan

act and ban

act were estimated using data generated by evaluating eNRTL at

90 °C in the relevant composition range from 16 to 20 wt.% of NaCl and 29 to 35 wt.% of

NaOH. The parameters are given in Table 1. Equation (30) is valid as long as NaCl or

NaOH are completely dissociated. Based on these thermodynamic considerations, the mass

flow of evaporated water from the anolyte compartment is

Fan,vapor

H2O= ˙

Nan

Cl2

1−yan,vapor

H2O

−˙

Nan

Cl2!·MH2O(31)

The mole fraction of water in the catholyte gas phase ycat,vapor

H2O, the activity coefficient γcat

H2O,

and the mass flow of evaporated water Fcat,vapor

H2Oare described analogously to Equation (28),

Equation (30), and Equation (31). The pressure in the gas phase of the catholyte com-

partment Pcat is fixed to 1.30 bar. The parameters for calculating γcat

H2Ocan be found in

Table 1.

3.2.3. Heat loss

Due to Assumption 7, we can describe the heat loss of the cell with the following equation:

Qloss =αloss ·Asurf ·(Tcell −Tamb) (32)

Herein, αloss is the overall heat transfer coefficient between the cell surface and the surround-

ing air for the case of free convection, Asurf is the outward-facing surface area of the cell and

Tamb is the ambient temperature.

As VDI e.V. (2013) offers different correlations for the heat transfer coefficient for the

various surfaces at the side, the top, and the bottom, an average value is used. This average

value is obtained by assuming constant surface and ambient temperatures of 90 and 20 °C, re-

spectively. The applied characteristic lengths are based on the geometry of a real electrolyzer

operated by Vestolit GmbH in Marl (personal communication). A selection of typical values

for cell geometries of different electrolyzer types is provided in O’Brien et al. (2007). The

used value for the product of heat transfer coefficient and area of the cell surface is given in

Table 1.

3.2.4. Electrical power

Finally, the input of electric power must be described (Kreysa and Wendt, 2010). For

one cell, it is given by

Pcell =Ucell ·j·Acell ·10−3(33)

Therein, Ucell is the voltage drop over one electrolysis cell due to the electric current, i.e.,

the product of current density jand cell area Acell. The factor 10−3is used for conversion

to kW. Based on Assumption 8 and Assumption 9, the cell voltage is only a function of the

current density:

Ucell =aU·j+bU(34)

in which the parameters aUand bUare regressed using measurement data from the real

plant (personal communication) in a wide load range (the used measurements are shown in

Section 4). These parameters are given in Table 1. Other approaches to determine the cell

voltage can, for example, be found in Otashu and Baldea (2019) and Wang et al. (2014).

The overall electric power consumption of the CAE is given by

Pel =Pcell ·ncell (35)

where ncell is the number of electrolysis cells used in the plant.

3.3. Electrochemical phenomena

The anode and cathode reaction rates considered in the mass balances in Equation (11),

(18), and (19) are calculated using Faraday’s law (Kreysa and Wendt, 2010):

ran =rcat =j

z·F(36)

where zis the number of transferred electrons per molecule of product and Fis Faraday’s

constant. The parameter zhas been set to 2 for both reactions.

In principle, current inefficiencies in the mass balances could be neglected if we assume

that all transferred electrons involved in the anolyte and catholyte reactions participate in

the production of chlorine or hydrogen. In the anolyte compartment, this is true in case no

significant aniodic oxygen formation takes place. This is given via Assumption 10, which

holds in case a proper anode material is selected and the anolyte is operated at low pH

value (O’Brien et al., 2007). In the catholyte compartment, the aforementioned assumption

is often fulfilled as the produced hydrogen is typically of high purity and only contains

evaporated water. (O’Brien et al., 2007).

3.3.1. Current efficiencies and electroosmosis

Since highly accurate results with respect to the produced amount of hydrogen are not

required here, the current efficiency for hydrogen, ξH2, is set to 100 % so that the mole flow

of produced hydrogen is ˙

Ncat

H2=rcat ·Acell (37)

At the anode, we need to consider current efficiencies to describe the usable amount of

produced chlorine because of the finite solubility of chlorine in the anolyte solution and

occurring reactions involving the dissolved chlorine. The mole flow of chlorine gas that

leaves the cell is defined as ˙

Nan

Cl2=ξCl2·ran ·Acell (38)

where ξCl2is the current efficiency with respect to chlorine.

To correctly predict the amount of produced caustic soda, the migration of sodium ions

through the membrane must be described. This migration is governed by the applied current

and could be computed using only Faraday’s law in case of an ideal membrane. To increase

the model accuracy regarding the produced caustic soda, we additionally need to consider

the back-migration of hydroxide ions. Thus, the migration of sodium ions is given by

˙nmem

Na+=ξNaOH ·j

zNa+·F(39)

where ξNaOH is the current efficiency with respect to caustic soda. Because electroneutrality

must be ensured in both cell compartments, the back-migration of hydroxide ions is given

˙nmem

OH–= (1 −ξNaOH)·j

zNa+·F(40)

Another important phenomenon is the mass transfer of water through the membrane.

Here, we only consider water migration due to electroosmosis since the influence of diffusion

can be neglected at industrially relevant current densities (O’Brien et al., 2007). Depending

on the membrane condition and the operating mode of the CAE, three to four molecules of

water are dragged through the membrane with each migrating sodium ion (O’Brien et al.,

2007). As the electroosmotic flow of water ˙nmem

H2Ois only coupled to the flow of sodium ions

˙nmem

Na+, this is described as

˙nmem

H2O=αelosm ·˙nmem

Na+(41)

For a detailed description of how to determine the parameters ξCl2,ξNaOH, and αelosm, the

reader is referred to O’Brien et al. (2007). In this contribution, they are determined by

parameter estimation using concentration and flow measurements at steady-state from a real

plant, operated by Vestolit GmbH in Marl (personal communication), at different operating

conditions. The estimated parameters are listed in Table 1 and are within the range predicted

by O’Brien et al. (2007). In accordance with O’Brien et al. (2007), these parameters are

assumed load-independent.

3.3.2. Liquid volume of the cells

As mentioned in Section 3.1, we use artificial controllers (see Section 3.5) to mimic

the dynamic behavior of the liquid volumes of both cell compartments and their reaction

to changes in the current density. Since the total volume of the gas and liquid phases is

kept constant by a weir, the theoretical liquid volumes used as set points for the artificial

controllers are

Van,SP =Vcat,SP =Vcell ·(1 −Φgas) (42)

Here, Vcell is the volume of each cell compartment and Φgas is the volume fraction taken

by chlorine or hydrogen in the respective compartment. Since the mass flow of the product

gases is linearly correlated to the current density due to Faraday’s law and the pressure is

held constant in both cell compartments, the linear relationship to the current density also

applies to the volume flow of the product gases (Assumption 3). Additionally, the behavior

of the product gases can be described as being similar to that in a bubble column, where

the gas fraction is proportional to the volume flow of the gas (Shah et al., 1982). Based

on these assumptions, the following linear correlation describing the relationship between

volume fraction and current density is used:

Φgas =agas ·j(43)

The parameter agas is determined based on measurement data from the real plant and its

value can be found in Table 1. For a detailed description of the design and the product gas

treatment of different membrane cell types, the reader is referred to O’Brien et al. (2007).

3.3.3. Autoprotolysis of water

To describe the extent of reaction of the autoprotolysis of water Ran

AW from Equation (12),

(13) and (14) without formulating a high-index differential-algebraic equation system, we

use the following expression for the equilibrium constant in the anolyte compartment:

KAW =˙

Van

in ·ρan

in ·wan

in,OH–

MOH–+Ran

AW + ˙nmem

OH–·Acell·˙

Van

in ·ρan

in ·wan

in,H3O+

MH3O++Ran

AW

˙

Van

in ·ρan

in ·wan

in,H2O

MH2O−2·Ran

AW −˙nmem

H2O·Acell −Fan,vapor

H2O

MH2O2(44)

Here, we use the outgoing molar flows of OH–, H3O+and water instead of the respective

molar holdups. The autoprotolysis of water in the catholyte compartment can be formulated

analogously to Equation (44) by replacing the outgoing molar flows with those from the

catholyte mass balances from Equation (18), (19), and (20).

The expression for the equilibrium constant is

ln (KAW) = aAW +bAW

Tcell + 273.15 +cAW ·ln Tcell + 273.15(45)

The included parameters were regressed with measurement data taken from K¨uster and

Thiel (2009) and are listed in Table 1.

3.4. Buffer tanks and catholyte recycle

Following Assumption 11, we use the same liquid density in the buffer tanks and the

respective cell compartments, which leads to the following expression:

dHUan,tank

dt =˙

Van

out ·ρan ·wan

c−˙

Van,tank

out ·ρan ·wan,tank

c∀c∈Can (46)

The total holdup and the mass fractions in the anolyte buffer tank are defined analogously

to Equation (15) and Equation (16), respectively. Additionally, the level (in %) of the buffer

tank is defined as

Lan,tank =Van,tank

Van,tank

ges

·100 % (47)

where Van,tank

ges is the overall usable volume of the buffer tank.

In addition, corresponding equations are formulated for the component mass balances,

the total holdup, the mass fractions and the level for the buffer tank in the catholyte recycle.

Again, only c∈Ccat has to be taken into account for the catholyte compartment.

In Section 2, we pointed out the relevance of the recycle section of the catholyte compart-

ment. As illustrated in Figure 2, the recycle can be described with the following steady-state

component balance:

Vcat

in ·ρcat

in ·wcat

in,c =˙

Vcat,tank

out ·ρcat ·wcat,tank

c−˙

Vcat

prod ·ρcat ·wcat,tank

+˙

Vcat

water ·ρcat

water ·wcat

water,c ∀c∈Ccat (48)

Therein, ˙

Vcat,tank

out and wcat,tank

care the outgoing liquid flow and the liquid mass fraction of

all components of the buffer tank (see Section 3.4). We assume that the density of caustic

Table 1: Regressed model parameters. Data source for parameter regression are: Eq. (29): Green and Perry

(2008); Eq. (30): Chen and Evans (1986) and APV 100 ENRTL-RK in Aspen Properties V10; Eq. (32):

VDI e.V. (2013); Eq. (34), (38), (39), (41), and (43): Vestolit GmbH (personal communication); Eq. (45):

K¨uster and Thiel (2009).

Parameter Eq. Parameter Eq.

aLV 0.0247 bar/°C (29) bU2.48 V (34)

bLV −1.5182 bar (29) ξCl20.984 −(38)

aan

act −1.9301 −(30) ξNaOH 0.983 −(39)

ban

act 1.0345 −(30) αelosm 3.812 −(41)

acat

act −3.8857 −(30) agas 5·10−5m2/A (43)

bcat

act 1.1772 −(30) aAW 156.48 −(45)

αloss ·Asurf 2.41 ·10−3kW/K (32) bAW −14 537 K (45)

aU1.46 ·10−4V m2/A (34) cAW −25.971 −(45)

soda solution leaving the system via ˙

Vcat

prod is roughly equal to the density in the buffer tank.

The density of water ρcat

water used to dilute the recycled caustic soda via ˙

Vcat

water is set to its

value at the approximate process temperature of 20 °C, i.e., 998.21 g L−1.

As shown in Table 2, the catholyte’s inlet temperature is used in the model to control

the temperature in the cell. However, the inlet temperature is not measured in the real

plant (see Section 4). Instead, the temperature behind the heat exchanger HE1 is measured

(see Figure 2). To allow a validation of the temperature control, the following steady-state

energy balance is used to determine the temperature in the catholyte recycle (Tcat

rec ):

Vcat

in ·ρcat

in ·Ccat,liquid

p,in,NaOH ·(Tcat

in −Tref) = ρcat(˙

Vcat,tank

out −˙

Vcat

prod)·Ccat,liquid

p,NaOH ·(Tcat

rec −Tref)

+˙

Vcat

water ·ρcat

water ·Ccat,liquid

p,water ·(Tamb −Tref ) (49)

Herein, the mass balance for the product splitter is considered. Since the temperature de-

pendence of the heat capacity of aqueous NaOH was found to be negligible (see Appendix B),

the same value as for the cell is used. The water used to dilute the catholyte is assumed to

be at ambient temperature (20 °C) with a heat capacity of 4.185 ·10−3kJ g−1K−1(VDI e.V.,

2013). Moreover, no reference enthalpies are considered in Equation (49) as no reaction is

expected to take place during the mixing process.

3.5. Control

According to Figure 2, five PI controllers are considered in the model to ensure that the

real plant behavior can be simulated. We use the following generic equation to implement

the controllers in the model:

u(t) = uOP +KP·v(t)−vSP+1

·Zt

0v(τ)−vSPdτ(50)

Table 2: Control pairing and base parameterization of the implemented controllers. The asterisk marks

artificial controllers for the liquid volume in both cell compartments.

v(t)u(t)vSP uOP KP1/TI

Tcell Tcat

in 90.0 °C 75.0 °C 7.0 °C/°C 0.005 1/s

wan

Na+˙

Van

in 0.073 g/g 0.057 L/s 100.0 L/s 0.001 1/s

wcat

Na+˙

Vcat

water 0.184 g/g 0.005 L/s 6.0 L/s 0.001 1/s

Lan,tank ˙

Van,tank

out varies 0.043 L/s 1.0 L/s 0.1 1/s

Vcat

prod ˙

Vcat,tank

out varies 0.085 L/s 1.0 L/s 0.1 1/s

Van(*) ˙

Van

out Van,SP 0.043 L/s 0.05 1/s 0.01 1/s

Vcat(*) ˙

Vcat

out Vcat,SP 0.085 L/s 0.5 1/s 0.01 1/s

Therein, u(t) is the manipulated variable, v(t) is the controlled variable, vSP is the respective

set point, and KPand TIare tuning parameters. An additional parameter uOP represents the

value of u(t) at the desired operating point at steady-state. It is used to stabilize the system

in the desired domain. The pairing and parameterization of the implemented controllers are

given in Table 2. The tuning parameters listed in Table 2 have been determined by solving

a dynamic optimization problem. In this optimization problem, the power consumption was

increased by following a heaviside step function and the control deviations of all controllers

were minimized simultaneously. Afterwards, a manual tuning procedure was carried out

based on the validation scenarios from Section 4.1 and Section 4.2 to improve the results.

The given values can be used as a base parameterization and a good control performance

should be achieved in the majority of possible simulation scenarios. Nevertheless, in some

cases it may be necessary to further tune the parameters.

To solve this equation with standard algorithms for differential-algebraic equation sys-

tems, which do not accept integrals in the equations, the controller equations are reformu-

lated by replacing the integral term with a new variable I, given by

I=Zt

0v(τ)−vSPdτ(51)

and by adding another differential equation for this new variable:

dt =v(t)−vSP (52)

As mentioned in Section 3.1, the liquid volumes Van and Vcat depend on the current density.

To mimic the real dynamic behavior of the liquid volumes in our model, we use artificial PI

controllers. The liquid volume set points Van,SP and Vcat,SP are given by Equation (42). The

controllers are defined analogously to Equation (50). The parameterization of these artificial

controllers can be found at the bottom of Table 2. Due to abrupt changes in the set points

of the artificial controllers, the manipulated variables ˙

Van

out and ˙

Vcat

out can take negative values

causing a back-flow from the buffer tanks into the cell compartments. To avoid this non-

physical behavior, we apply a smooth reformulation of a max operator (Hoffmann et al.,

2020a) to the manipulated variable determined by the PI controller. The smooth max

operator is defined as

umodel =uPI +pu2

PI + 0.01

2(53)

where umodel is the value of the manipulated variable used in the model equations and uPI is

the value determined by Equation (50). Using this procedure, positive values of uPI will be

directly applied to the model whereas negative values in uPI will result in umodel switching

to zero.

3.6. Implementation

As part of this work, the dynamic model of the CAE has been implemented in MO-

SAICmodeling, a web-based modeling, simulation, and optimization environment (Merchan

et al., 2015; Esche et al., 2017). Using its code generator for various programming languages

and modeling environments, the dynamic system is exported to the gPROMS model builder

(Process Systems Enterprise, 1997-2020).

4. Model validation

In the following section, the model is validated using two dynamic data sets from an

industrial CAE plant operated by Vestolit GmbH in Marl (Germany). The considered

plant produces ca. 730 t of chlorine per day in more than 1000 CAE cells and consumes

around 75 MW of electrical power. The first validation set contains multiple load changes

of approximately 4 MW (ca. 5.5 % of the maximum load) with load ramps of 0.13 MW s−1

over a time horizon of 175 min. These specifications meet the requirements of the primary

balancing market in Germany for which a load change must be carried out within 30 s. In the

second validation, one large load change of approximately 50 MW (ca. 67 % of the maximum

load) with a load ramp of 0.08 MW s−1over a time horizon of 160 min is performed. These

specifications are rather representative of an application in the day-ahead market for which

larger load changes occur. While such large load changes are unlikely to be realised on

a regular basis in a real plant, the scenario still allows for an assessment of the model’s

performance in regions far away from the typical operating point. It may thus serve as a

worst-case scenario and model-plant mismatch is expected to be even smaller for less drastic

load changes.

The validation focuses on the evaluation of the model performance regarding the cell

temperature and the concentrations in both cell compartments. Since the reaction rates

of chlorine and hydrogen only depend on the current and the accuracy of the respective

measurement sensors of both streams in the real plant is low, results for these variables

are not shown in the following. Instead, we focus on the composition and temperature of

both electrolytes as they are much more relevant for membrane stability and quality of the

produced caustic soda. Their measurements are provided with a sampling rate of 15 s. To

increase the clarity of the results, the total flow measurements are divided by the number

of installed cells to obtain a representative average cell.

Table 3: Design specifications used in both validation simulations. The aqueous NaCl feed and the incoming

water are specified with a pH value of 3 and 7, respectively.

Parameter Parameter

Acell 2.72 m2wan

in,Cl–0.157 g/g

Tan

in 63.0 °Cwan

in,OH–1.706 ·10−13 g/g

Tamb 20.0 °Cwan

in,H3O+1.908 ·10−5g/g

Vcell 100.0 L wcat

water,Na+0.0 g/g

Van,tank

ges 27.8 L wcat

water,OH–1.706 ·10−9g/g

Vcat,tank

ges 27.8 L wcat

water,H3O+1.908 ·10−9g/g

wan

in,Na+0.102 g/g

The design and feed specifications for both dynamic simulations are shown in Table 3.

The aqueous NaCl feed ˙

Van

in and the incoming water ˙

Vcat

water are specified with a pH value of 3

and 7, respectively. The other geometric (e.g., cell area and cell and buffer tank volumes) or

operating specifications (e.g., anolyte feed temperature Tan

in ) are taken from real plant data.

Both simulations are initialized at the steady-state defined by the initial set of measurements.

4.1. Validation 1: small load changes

To simulate the first load change scenario, the catholyte feed flow ˙

Vcat

in and the current

density jare set to their measured values during operation at every time step during the

simulation. In addition, the set points of the implemented controllers (Section 3.5) are set

to the measured values of the respective variables. In order to fully retain the information

content regarding the highly dynamic load changes, the unfiltered values of the measured

data are used to specify the inputs and the set points.

Both dynamic inputs are shown in Figure 4 (a and b). Figure 4 (c and d) shows the result-

ing cell voltage and the power consumption of all cells. The comparison with measurement

data shows that the model describes the cell voltage (and therefore the power consumption

as well) with high precision. In Figure 5 (a and b) and Figure 6 (a), the results for quality

and temperature control of the cell are shown. The dynamic set points for both quality

and temperature controllers are continuously set to the measured values of the controlled

variables and are met with high precision. The results for the manipulated variables of the

aforementioned quality controllers and the catholyte recycle temperature are displayed in

Figure 5 (c and d) and Figure 6 (b). Concerning quality control, it is evident that the

simulated control actions show the same dynamics as the real plant. Possible explanations

for the larger deviations between the flow measurements and the simulation results of ˙

Van

in Figure 5 (c) could be the assumed constant anolyte inlet temperature Tan

in or small de-

viations in the electroosmotic water flow through the membrane. The dynamic behavior of

the catholyte’s recycle temperature is not completely consistent with the measurement data.

Instead, less aggressive control actions are taken, which indicates that the simulated system

0 50 100 150

0.075

0.0755

0.076

0 50 100 150

5200

5300

5400

5500

5600

5700

0 50 100 150

3.22

3.24

3.26

3.28

3.3

0 50 100 150

Figure 4: Dynamic inputs and results for validation 1: (a) anolyte feed flow, set to the unfiltered measure-

ments; (b) current density, set to the unfiltered measurements; (c) simulated and measured cell voltage; (d)

overall power consumption.

responds faster to temperature changes. Additionally, a small offset between simulated and

measured temperature is present. A reasonable explanation for this offset is the deviation in

the anolyte inlet flow (Figure 5 (c)), which leads to a stronger cooling effect in the anolyte

compartment. This cooling must be compensated for by a higher catholyte recycle tempera-

ture. To test this hypothesis, we removed the quality controller for the anolyte composition

from the model and fixed the anolyte inlet flow according to its measured values (see Fig-

ure 5 (c)). The results of this procedure are shown in Figure 6 (b) (blue line), where the

deviations between simulation and measurements are reduced to less than 1 % on average.

However, removing the quality controller from the model leads to larger deviations for the

anolyte composition (around 1 % with a maximum offset of 2 %).

Finally, Figure 7 (a) demonstrates that the load-dependent liquid volumes in both cell

compartments are accurately obtained by the artificial controllers. We do not show the

associated manipulated variables here as they are not recorded in the real plant. Figure 7 (b)

shows the results for the product flow of caustic soda. Again, the manipulated variable of

this controller – the liquid outlet of the buffer tank – is not measured. The controlled flow

shows fluctuations compared to its measurement. This behavior can be attributed to the

specified input variables ˙

Vcat

in and j. Since this flow is not a manipulated variable of any

controller, it is used in the model to fulfill the mass balances. Consequently, the fluctuations

in the unfiltered dynamic inputs are reflected in this variable.

In summary, the dynamic trajectories from scenario 1 can be reproduced with high

accuracy by the formulated model.

0 50 100 150

0.071

0.072

0.073

0.074

0.075

0 50 100 150

0.18

0.182

0.184

0.186

0 50 100 150

0.052

0.054

0.056

0.058

0.06

0.062

0.064

0 50 100 150

4.6

4.8

5.2

10-3

Figure 5: Product quality results for validation 1: (a) mass fraction of Na+in the anolyte, unfiltered

measurements used as set points for the quality controller; (b) mass fraction of Na+in the catholyte,

unfiltered measurements used as set points for the quality controller; (c) inlet volume flow for anolyte,

manipulated variable of the quality controller in (a) compared to its measured values; (d) inlet volume flow

of water for catholyte, manipulated variable of the quality controller in (b) compared to its measured values.

4.2. Validation 2: large load change

In the following, the results of the second validation scenario are presented. As described

in Section 4.1, the catholyte feed flow ˙

Vcat

in , the current density j(Figure 8 (a and b)), and

the set points of the implemented controllers are set to their measured values. Again, no

filtering of the measurements was carried out to retain the full information content of the

highly dynamic load changes. The cell voltage from Figure 8 (c) is again accurately described

given the small deviations to the measurement data. In Figure 8 (d), the load change of

more than 50 MW (67 % of the maximum load) is displayed.

The set points of the implemented quality controllers (Figure 9 (a and b)) and of the

temperature controller (Figure 10 (a)) are tracked just as well as in scenario 1. The minor

deviations can be attributed to the controller parameterization. The manipulated volume

flows for quality control shown in Figure 9 (c and d) and the catholyte’s recycle temperature

(Figure 10 (b)) coincide well with their corresponding measurements.The strong fluctuations

in the measured water flow (Figure 9 (d)) between minute 120 and 130 are probably caused

by a malfunction of the measuring device, since no effects of the strong fluctuations on the

catholyte composition in the cell (see Figure 9 (b)) can be observed. Again, less aggressive

control actions are taken by the temperature controller, which is particularly well visible

after 40 minutes. At this point, the temperature in the real plant changes at a much

0 50 100 150

87.2

87.4

87.6

87.8

88.2

0 50 100 150

Figure 6: Temperature results for validation 1: (a) cell temperature, unfiltered measurements used as set

points for the temperature controller; (b) catholyte recycle temperature (normal simulation and test case

with anolyte inlet flow fixed to measurements) compared to its measured values.

higher rate compared to the simulation results. To illustrate the influence of these different

recycle temperatures on the cell temperature, another test was performed by removing the

temperature controller from the model and setting the recycle temperature to its measured

values in Figure 10 (b). The results of this test case are shown in Figure 10 (a) (blue line)

and show the relatively small impact of this mismatch for the recycle temperature. The

deviation between simulation and measurements does not exceed 1.1 % at any time.

Figure 11 (a and b) also shows a good agreement between set point and actual value of

the load-dependent liquid volumes in both cell compartments and the caustic soda product

flow. However, the high dynamic changes in liquid volume cannot be fully reproduced by the

implemented controllers without producing non-physical behavior. Again, we do not show

the associated manipulated variables here as they are not recorded in the real plant. The

fluctuations in the unfiltered dynamic input specifications are reflected in the product flow

of caustic soda. However, due to the scaling of Figure 11 (b) these fluctuations are invisible.

Nevertheless, the deviations between the manipulated variables and the catholyte recycle

temperature are somewhat larger compared to scenario 1. This can be explained by the fact

that the large and highly dynamic changes in liquid volume in both cells cannot be fully

described by the developed model. This may have a large influence on the energy and the

mass balance, leading to inaccuracies in volume flows and temperatures, especially for large

load changes.

As already mentioned in Section 1, the cell temperature and the concentrations need

to be maintained at a constant value when implementing load changes in the context of

demand response to ensure membrane stability. While this shows certain limitations of the

model’s accuracy for large load changes, we note that such a scenario would be irrelevant for

DR as such drastic load changes would rarely appear. In summary, the dynamic trajectories

from scenario 2 can be reproduced with high accuracy by the proposed model.

0 50 100 150

71.5

72.5

73.5

0 50 100 150

0.013

0.0132

0.0134

0.0136

0.0138

Figure 7: Control results with unmeasured manipulated variables for validation 1: (a) liquid volume of

anolyte and catholyte compartment, set points from Equation (42) used for the volume controller; (b)

product flow of caustic soda, unfiltered measurements used as set points for the flow controller.

0 50 100 150

0.06

0.07

0.08

0 50 100 150

2000

4000

6000

0 50 100 150

2.4

2.6

2.8

3.2

0 50 100 150

Figure 8: Dynamic inputs and results for validation 2: (a) anolyte feed flow, set to the unfiltered measure-

ments; (b) current density, set to the unfiltered measurements; (c) simulated and measured cell voltage; (d)

overall power consumption.

0 50 100 150

0.071

0.072

0.073

0.074

0.075

0 50 100 150

0.18

0.182

0.184

0.186

0 50 100 150

0.02

0.04

0.06

0 50 100 150

810-3

Figure 9: Product quality results for validation 2: (a) mass fraction of Na+in the anolyte, unfiltered

measurements used as set points for the quality controller; (b) mass fraction of Na+in the catholyte,

unfiltered measurements used as set points for the quality controller; (c) inlet volume flow for anolyte,

manipulated variable of the quality controller in (a) compared to its measured values; (d) inlet volume flow

of water for catholyte, manipulated variable of the quality controller in (b) compared to its measured values.

0 50 100 150

Figure 10: Temperature results for validation 2: (a) cell temperature (normal simulation and test case

with recycle temperature fixed to its measured values), unfiltered measurements used as set points for the

temperature controller; (b) catholyte recycle temperature compared to its measured values.

4.3. Evaluation of model performance

Table 4 lists the average and maximum deviations between simulation and measurements

for the measured variables that were not used as controller set points. The average deviation

is calculated as

|yi−y∗

y∗

·1

·100 % (54)

where yiand y∗

iare the simulation result and the associated measurement at time point i,

respectively, and nmis the total number of measurement points.

In scenario 1, the largest average deviations occur for the inlet flow of the anolyte and the

inlet temperature of the cathode and amount to 2.44 and 1.59 %, respectively. Compared to

scenario 1, the deviation between simulation results and measurements increases for every

considered variable in scenario 2 with the anolyte inlet flow showing the largest average

deviation. Although the large load change cannot be fully described by the developed

model, the average deviation of all variables is below 6.5 %. However, significantly larger

deviations occur during drastic load reduction, which is shown by the maximum deviations

in Table 4 for the anolyte inlet flow and especially for the inlet flow of water.

Nevertheless, we deem the model successfully validated with real plant data for load

changes of typical magnitude in demand response and it may thus be used for determining

optimal trajectories in demand response scenarios. For larger load changes, qualitatively

correct results can still be expected.

5. Conclusion and outlook

In this contribution, a dynamic model of an industrial chlor-alkali plant was presented,

which describes all sensors and actuators relevant for dynamic operation. The CAE is

described by two continuously stirred-tank reactors representing both cell compartments.

The membrane is modeled by describing the migration of sodium ions, the back-migration

of hydroxide ions, and the electroosmotic flow of water between both compartments. The

reactions at anode and cathode are described based on Faraday’s law. Additionally, autopro-

tolysis and evaporation of water are considered in the cell. Furthermore, the load-dependent

liquid volumes in both cell compartments, the catholyte recycle, and buffer tanks down-

stream of both cell compartments are taken into account. To increase the robustness of the

model and to enable online applications (e.g. dynamic real-time optimization), we simpli-

fied as many of the nonlinear expressions as possible to obtain linear correlations for process

characteristics and fluid properties. For this purpose, we assumed that the model is used in

a small operation window regarding the cell temperature and the concentrations in the cell

as this procedure is suggested for an economic plant operation.

The model was successfully validated with dynamic plant data from an industrial CAE

plant. Under load change conditions that meet the requirements for DR operation in the

primary balancing market in Germany, the industrial plant’s behavior can be reproduced

with a maximum average deviation of 2.44 %. Under differing operating conditions, for

which a drastic load reduction of around 67 % occurs, the average deviation is larger but

0 50 100 150

10-3

Figure 11: Control results with unmeasured manipulated variables for validation 2: (a) liquid volume of

anolyte and catholyte compartment, set points from Equation (42) used for the volume controller; (b)

product flow of caustic soda, unfiltered measurements used as set points for the flow controller.

Table 4: Average deviation D(calculated using Equation (54)) and maximum deviation between simulation

and measurements. In validation 2 values between minute 120 and 130 are not taken into account for the

inlet water flow, due to a measurement malfunction in this region..

Variable Average dev. in % Maximum dev. in %

Val. 1 Val. 2 Val. 1 Val. 2

Cell voltage Ucell 0.04 0.32 0.57 2.93

Inlet flow of anolyte ˙

Van

in 2.44 5.11 6.28 22.24

Inlet flow of water ˙

Vcat

water 0.64 5.06 3.65 87.94

Recycle temperature Tcat

rec 1.59 1.97 2.71 6.78

a value of ca. 6.5 % is never exceeded. However, in the areas of the drastic load reduction

the real plant behavior cannot be fully described with the developed model. Such drastic

load changes are, however, unlikely to occur in DR applications and the applicability of the

developed model in DR is not limited.

In the future, the model will be used to determine optimal load change trajectories for

the CAE in demand response scenarios. For this purpose, an optimal control approach

will be applied where the implemented controllers will be removed from the model and

the manipulated variables will be determined by dynamic optimization. We will focus on

demand response scenarios based on day-ahead and balancing markets. Additionally, the

results of the optimal control can be used to evaluate the feasibility of the given load changes

by showing whether the specified set points can be fulfilled or whether purity constraints

may be violated.

Acknowledgements

The authors acknowledge the financial support by the Federal Ministry of Economic

Affairs and Energy of Germany in the project ChemEFlex (project number 0350013A).

References

Alexandrov, A.A., 2005. The equations for thermophysical properties of aqueous solutions of sodium hy-

droxide, in: Proceedings of the 14th International Conference on the Properties of Water and Steam, pp.

86–90.

Ausfelder, F., Seitz, A., von Roon, S. (Eds.), 2018. Flexibilita¨atsoptionen in der Grundstoffindustrie

Methodik, Potenziale, Hemmnisse. DECHEMA Gesellschaft f¨ur Chemische Technik und Biotechnolo-

gie e.V., Frankfurt am Main.

Babu, C.A., Ashok, S., 2008. Peak Load Management in Electrolytic Process Industries. IEEE Transactions

on Power Systems 23, 399–405. doi:10.1109/tpwrs.2008.920732.

Baetens, J., Kooning, J.D.M.D., Eetvelde, G.V., Vandevelde, L., 2020. A two-stage stochastic optimisation

methodology for the operation of a chlor-alkali electrolyser under variable DAM and FCR market prices.

Energies 13, 5675. doi:10.3390/en13215675.

Br´ee, L.C., Bulan, A., Herding, R., Kuhlmann, J., Mitsos, A., Perrey, K., Roh, K., 2020. Techno-economic

comparison of flexibility options in chlorine production. Industrial & Engineering Chemistry Research

59, 12186–12196. doi:10.1021/acs.iecr.0c01775.

Br´ee, L.C., Perrey, K., Bulan, A., Mitsos, A., 2019. Demand side management and operational mode

switching in chlorine production. AIChE Journal 65, e16352. doi:10.1002/aic.16352.

Budiarto, T., Esche, E., Repke, J.U., Leksono, E., 2017. Dynamic model of chloralkali membrane process.

Procedia Engineering 170, 473–481. doi:10.1016/j.proeng.2017.03.076.

Castro, P.M., Ave, G.D., Engell, S., Grossmann, I.E., Harjunkoski, I., 2020. Industrial demand side man-

agement of a steel plant considering alternative power modes and electrode replacement. Industrial &

Engineering Chemistry Research 59, 13642–13656. doi:10.1021/acs.iecr.0c01714.

Chase, M.W., 1998. NIST-JANAF Themochemical Tables, Fourth Edition. Journal of Physical and Chemical

Reference Data, Monograph 9 .

Chen, C.C., Evans, L.B., 1986. A local composition model for the excess gibbs energy of aqueous electrolyte

systems. AIChE Journal 32, 444–454. doi:10.1002/aic.690320311.

Dowling, A.W., Zavala, V.M., 2018. Economic opportunities for industrial systems from frequency regulation

markets. Computers & Chemical Engineering 114, 254–264. doi:10.1016/j.compchemeng.2017.09.018.

Esche, E., Hoffmann, C., Illner, M., M¨uller, D., Fillinger, S., Tolksdorf, G., Bonart, H., Wozny, G., Repke,

J.U., 2017. MOSAIC - enabling large-scale equation-based flow sheet optimization. Chemie Ingenieur

Technik 89, 620–635. doi:10.1002/cite.201600114.

Gmehling, J., Kleiber, M., Kolbe, B., Rarey, J., 2019. Chemical thermodynamics for process simulation.

Wiley-VCH, Weinheim.

Green, D.W., Perry, R.H., 2008. Perry’s chemical engineers’ handbook. McGraw Hill professional, McGraw-

Hill, New York.

Hoffmann, C., Weigert, J., Esche, E., Repke, J.U., 2020a. A pressure-driven, dynamic model for distillation

columns with smooth reformulations for flexible operation. Computers & Chemical Engineering 142,

107062. doi:10.1016/j.compchemeng.2020.107062.

Hoffmann, C., Weigert, J., Esche, E., Repke, J.U., 2020b. Towards demand-side management of the chlor-

alkali electrolysis: Dynamic, pressure-driven modeling and model validation of the 1,2-dichloroethane

synthesis. Chemical Engineering Science 214, 115358. doi:10.1016/j.ces.2019.115358.

Klaucke, F., Hoffmann, C., Hofmann, M., Tsatsaronis, G., 2020. Impact of the chlorine value chain on

the demand response potential of the chloralkali process. Applied Energy 276, 115366. doi:10.1016/j.

apenergy.2020.115366.

Klaucke, F., Karsten, T., Holtrup, F., Esche, E., Morosuk, T., Tsatsaronis, G., Repke, J.U., 2017. Demand

response potenziale in der chemischen industrie. Chemie Ingenieur Technik 89, 1133–1141. doi:10.1002/

cite.201600073.

Kreysa, G., Wendt, H., 2010. Electrochemical Engineering. Springer Science & Business Media.

K¨uster, F.W., Thiel, A., 2009. Rechentafeln f¨ur die Chemische Analytik. Walter de Gruyter.

Masterton, W.L., Slowinski, E.J., Stanitski, C.L., 1983. Chemical principles: alternate edition with a

qualitative analysis supplement. CBS College Publishing, Philadelphia.

Merchan, V.A., Esche, E., Fillinger, S., Tolksdorf, G., Wozny, G., 2015. Computer-aided process and plant

development. a review of common software tools and methods and comparison against an integrated

collaborative approach. Chemie Ingenieur Technik 88, 50–69. doi:10.1002/cite.201500099.

O’Brien, T.F., Bommaraju, T.V., Hine, F., 2007. Handbook of Chlor-Alkali Technology. volume 1. Springer-

Verlag GmbH.

Otashu, J.I., Baldea, M., 2019. Demand response-oriented dynamic modeling and operational optimization

of membrane-based chlor-alkali plants. Computers & Chemical Engineering 121, 396–408. doi:10.1016/

j.compchemeng.2018.08.030.

Otashu, J.I., Baldea, M., 2020. Scheduling chemical processes for frequency regulation. Applied Energy 260,

114125. doi:10.1016/j.apenergy.2019.114125.

Paterakis, N.G., Erdin¸c, O., Catal˜ao, J.P.S., 2017. An overview of demand response: Key-elements and

international experience. Renewable and Sustainable Energy Reviews 69, 871–891. doi:10.1016/j.rser.

2016.11.167.

Paulus, M., Borggrefe, F., 2011. The potential of demand-side management in energy-intensive industries

for electricity markets in germany. Applied Energy 88, 432–441. doi:10.1016/j.apenergy.2010.03.017.

Process Systems Enterprise, 1997-2020. gPROMS. URL: https://www.psenterprise.com/gproms.

Richstein, J.C., Hosseinioun, S.S., 2020. Industrial demand response: How network tariffs and regulation

(do not) impact flexibility provision in electricity markets and reserves. Applied Energy 278, 115431.

doi:10.1016/j.apenergy.2020.115431.

Roh, K., Br´ee, L.C., Perrey, K., Bulan, A., Mitsos, A., 2019. Flexible operation of switchable chlor-alkali

electrolysis for demand side management. Applied Energy 255, 113880. doi:10.1016/j.apenergy.2019.

113880.

Sassi, K.M., Mujtaba, I.M., 2013. Optimal operation of RO system with daily variation of freshwater

demand and seawater temperature. Computers & Chemical Engineering 59, 101–110. doi:10.1016/j.

compchemeng.2013.03.020.

Shah, Y., Kelkar, B.G., Godbole, S., Deckwer, W.D., 1982. Design parameters estimations for bubble column

reactors. AIChE Journal 28, 353–379.

Simkoff, J.M., Baldea, M., 2020. Stochastic scheduling and control using data-driven nonlinear dynamic

models: Application to demand response operation of a chlor-alkali plant. Industrial & Engineering

Chemistry Research 59, 10031–10042. doi:10.1021/acs.iecr.9b06866.

Tanner, J.E., Lamb, F.W., 1978. Specific heats of aqueous solutions of NaCl, NaBr, and KCl: Comparisons

with related thermal properties. Journal of Solution Chemistry 7, 303–316. doi:10.1007/bf00644277.

Tsay, C., Kumar, A., Flores-Cerrillo, J., Baldea, M., 2019. Optimal demand response scheduling of an

industrial air separation unit using data-driven dynamic models. Computers & Chemical Engineering

126, 22–34. doi:10.1016/j.compchemeng.2019.03.022.

VDI e.V. (Ed.), 2013. VDI-W¨armeatlas. 11, Springer Vieweg.

Wang, X., Teichgraeber, H., Palazoglu, A., El-Farra, N.H., 2014. An economic receding horizon optimization

approach for energy management in the chlor-alkali process with hybrid renewable energy generation.

Journal of Process Control 24, 1318–1327. doi:10.1016/j.jprocont.2014.04.017.

Wellington, T.C., 1992. Modern Chlor-Alkali Technology : Volume 5. Springer Netherlands, Dordrecht.

Worrell, E., Phylipsen, D., Einstein, D., Martin, N., 2000. Energy use and energy intensity of the U.S.

chemical industry. Technical Report. Lawrence Berkeley National Lab., CA (US). doi:10.2172/773773.

Zhao, S., Ochoa, M.P., Tang, L., Lotero, I., Gopalakrishnan, A., Grossmann, I.E., 2019. Novel formulation

for optimal schedule with demand side management in multiproduct air separation processes. Industrial

& Engineering Chemistry Research 58, 3104–3117. doi:10.1021/acs.iecr.8b04964.

Nomenclature

Abbreviations

CAE Chlor-alkali electrolysis

eNRTL Electrolyte Non-Random

Two-Liquid model

Greek Symbols

Φ Volume fraction

αelosm Electroosmosis ratio for water

γActivity coefficient

ξCurrent efficiency

ρDensity, g L−1

Latin Symbols

AArea, m2

CComponent set

CpHeat capacity, J g−1K−1

DAverage deviation, %

FMass flow, g s−1

FFaraday’s constant,

96 485.3 A s mol−1

HEnthalpy, kJ

HEnthalpy flow, kJ s−1

HU Holdup, g

IIntegral of controller offset,

unit varies

KAW Equilibrium constant of water

autoprotolysis

KPProportional parameter of PI

controller, unit varies

MMolecular weight, g mol−1

NMole flow of product gases,

mol s−1

PElectrical power, kW

QHeat flow, kJ s−1

PPressure, bar

RAW Extent of reaction for

autoprotolysis, mol s−1

TTemperature, °C

TIIntegration time of PI

controller, s

UVoltage, V

VVolume, L

VVolume flow, L s−1

aParameter in various

correlations

bParameter in various

correlations

cParameter in various

correlations

jCurrent density, A m−2

nNumber of

˙nMolar flux through membrane,

mol s−1m−2

rReaction rate, mol s−1m−2

tTime, s

uManipulated variable of PI

controller, unit varies

vControlled variable of PI

controller, unit varies

wMass fraction, g g−1

yVapor mole fraction, mol mol−1

yvariable placeholder

zNumber of transferred electrons

Indices

cComponent index

Subscripts

act Activity coefficient

amb Ambient

evap Evaporation

in Inlet

loss heat loss

m measurements

model in the model eqeuations

out Outlet

PI PI controller

prod Product flow

reac Reaction

rec Recycle

ref Reference state

water Water feed

Superscripts

an Anolyte

cat Catholyte

cell Electrolysis cell

elosm Electroosmosis

gas gas (chlorine or hydrogen) in

the cells

liquid Liquid phase

LV Liquid-vapor

mem Membrane

OP Feed-forward value of PI

controller

SP Set point

U Voltage

surf Surface

tank In buffer tank

vapor Vapor

* measurement

Appendix A. Density correlations

The densities of aqueous NaCl and aqueous NaOH are modeled as functions of mass

fraction of sodium ions and temperature. To describe the respective density in the model,

the corresponding temperature and mass fractions must be entered into Equation (A.1). For

example, the density of the anolyte is:

ρan =aNaCl +bNaCl ·



wan

Na+

MNa+

wan

Na+

MNa++wan

Cl–

MCl–+wan

H2O

MH2O



+cNaCl ·Tcell + 273.15

+dNaCl ·



wan

Na+

MNa+

wan

Na+

MNa++wan

Cl–

MCl–+wan

H2O

MH2O





+eNaCl ·Tcell + 273.152

+fNaCl ·



wan

Na+

MNa+

wan

Na+

MNa++wan

Cl–

MCl–+wan

H2O

MH2O



·Tcell + 273.15(A.1)

The parameters for the brine correlation (aNaCl to fNaCl) and the caustic soda correlation

(aNaOH to fNaOH) were regressed using density data from Green and Perry (2008) and are

given in Table A.1.

Table A.1: Parameters of simplified density correlation for aqueous NaCl and aqueous NaOH, valid in the

temperature range 0 to 100 °C, 1 to 26 wt.% of NaCl, 1 to 50 wt.% of NaOH. The last column refers to both

entries in a row.

Parameter Parameter Unit

aNaCl 908.249 aNaOH 1030.296 g/L

bNaCl 2786.685 bNaOH 2883.396 g/L

cNaCl 9.170 ·10−1cNaOH 1.782 ·10−1g/(L K)

dNaCl −312.883 dNaOH −1290.705 g/L

eNaCl −2.088 ·10−3eNaOH −9.989 ·10−4g/(L K2)

fNaCl −1.740 fNaOH −1.140 g/(L K)

Appendix B. Heat capacities and enthalpy of evaporation

The heat capacities of the molecular components Can

p,Cl2and Ccat

p,H2, were calculated with

correlations from Chase (1998). Since the variations in the relevant temperature range were

found to be negligible, the heat capacities have been fixed to their respective values of

0.4923 ·10−3kJ g−1K−1for Cl2and 14.56 ·10−3kJ g−1K−1for H2(reference state: 90 °C).

The heat capacity of aqueous NaCl was calculated using the correlation by Tanner and

Lamb (1978). An analysis in the relevant temperature and composition range showed that

Table B.2: Parameters of simplified correlation for heat capacities for aqueous NaCl, valid in the temperature

range 70 to 95 °C and 17 to 26 wt.% of NaCl.

Parameter

aan

Cp −9.0937 ·10−3kJ/(g K)

ban

Cp 4.0248 ·10−3kJ/(g K)

the influence of the temperature is negligible and the heat capacity depends linearly on the

composition (see Figure 3 (b)). Therefore, the following simplified expression is used:

Can,liquid

p,NaCl =aan

Cp ·



wan

Na+

MNa+

wan

Na+

MNa++wan

H2O

MH2O



+ban

Cp (B.1)

The heat capacity of the anolyte inlet stream Can,liquid

p,in,NaCl can also be described using Equa-

tion (B.1). The estimated parameters are given in Table B.2.

The heat capacities of NaOH solution Ccat,liquid

p,in,NaOH and Ccat,liquid

p,NaOH are described based on

correlations from Alexandrov (2005). Since the variations in the relevant temperature and

composition range were found to be negligible, the heat capacities of aqueous NaOH have

been fixed to 4.155 ·10−3kJ g−1K−1for Ccat,liquid

p,in,NaOH (reference state: 80 °C and 30 wt.%) and

4.109 ·10−3kJ g−1K−1for Ccat,liquid

p,NaOH (reference state: 90 °C and 32 wt.%).

The enthalpy of evaporation of water ∆hevap,H2Ois described based on the correlation by

VDI e.V. (2013). Again, the variations in the relevant temperature range are negligible and

the parameter has been fixed to 2.3148 kJ g−1(reference state: 88 °C).