Non-equilibrium thermodynamic model of water sorption in Nafion membranes

Journal of Membrane Science 540 (2017) 35–49

Contents lists available at ScienceDirect

Journal of Membrane Science journal homepage: www.elsevier.com/locate/memsci

Non-equilibrium thermodynamic model of water sorption in Nafion membranes

T

⁎

Václav Klikaa, , Jan Kubanta,b, Michal Pavelkac,d, Jay B. Benzigere a

Department of Mathematics, FNSPE, Czech Technical University in Prague, Trojanova 13, 120 00 Praha, Czech Republic New Technologies, Research Centre, University of West Bohemia, Univerzitní 8, 306 14 Plzeň, Czech Republic c Mathematical Institute, Faculty of Mathematics and Physics, Charles University in Prague, Sokolovska 83, 186 75 Prague, Czech Republic d Department of Chemical Engineering, University of Chemistry and Technology Prague, Technicka 5, 16628 Prague 6, Czech Republic e Chemical and Biological Engineering, Princeton University, Princeton, NJ 08544, USA b

A R T I C L E I N F O

A B S T R A C T

Keywords: Sorption rates Water sorption Thermodynamic model Nafion membrane Model analysis Swelling

Despite the extensive attention that water transport in Nafion membranes has been subject to, several phenomenal still remain unexplained. It is for instance not clear why desorption rate is much faster than absorption. The purpose of this paper is to suggest an answer to that problem. A simple but thermodynamically consistent model is formulated within the framework of Classical Irreversible Thermodynamics including swelling and interaction of the membrane with the surrounding water vapor. Numerical solutions to the model are in a good match with experimental data of sorption dynamics. By simplifying the model (excluding swelling and assuming instantaneous diffusion within the membrane), we are able to identify the source of the asymmetry between the absorption and desorption rates - the shape of the uptake curve. This natural explanation of the asymmetry is then discussed quantitatively (showing an agreement between both the full and the simplified model but also with the experimental data) and also in the context of active carbon sorption.

1. Introduction Despite the intensive studying of water transport in PEM, its understanding is not yet complete and many questions still remain open. There are many experiments nowadays showing “Non-Fickian” or paradoxical behavior of water in Nafion. For example it is not clear why is water absorption approximately ten times slower than water desorption, why sorption curves for different widths of membranes L colt lapse into one curve when plotted against L (t being time) and not t

against L2 as simple Fickian model would predict or why absorption does not show approximately exponential behavior [1]. From the theoretical point of view, these issues are being addressed by two complementary groups of models. One is trying to capture the physics of membrane sorption in detail to have a rather complete description of involved processes, e.g. [2,3]. The second group consists of simple models allowing for analytical insight and hence revealing the role of the considered physical phenomena on the observed outcomes, e.g. [4,5,1,6]. Hence both groups of models serve its purpose and provide complementary understanding of water sorption. In this work we propose a model from the latter group but which is highly nonlinear and more widely applicable than the typically used

⁎

Corresponding author. E-mail address: vaclav.klika@fjfi.cvut.cz (V. Klika).

http://dx.doi.org/10.1016/j.memsci.2017.06.025 Received 5 August 2016; Received in revised form 22 December 2016; Accepted 9 June 2017 Available online 13 June 2017 0376-7388/ © 2017 Elsevier B.V. All rights reserved.

linear models. Its formulation follows from the non-equilibrium thermodynamics, the choice of a sorption isotherm (which summarizes many physical and chemical phenomena into a single empirical relation), an experimental expression for self-diffusion coefficients, a theoretical relation between self-diffusion and mutual diffusion coefficient and includes swelling. Apart from this more precise formulation of diffusion flux compatible with thermodynamics the proposed model including the boundary conditions is not new [1,6,5] although it has not been formulated in this complete form, at least according to authors best knowledge. Importantly, we successfully analyzed this nonlinear model under certain simplifying yet admissible assumptions revealing the mechanism behind the disparity of the time scales between absorption and desorption. Qualitatively, we show that the cause of disparity can be identified with particular water uptake dependence on water activity. Quantitatively, we show that the presented simple model is able using a single fitted parameter but with one value for all experimental settings to reproduce sorption data. The text is structured as follows. First, a thermodynamically consistent model of water sorption that includes swelling, called the full model, is formulated. Consequently its simplified version for fast diffusion with respect to sorption kinetics is derived, called the lumped

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

Fig. 1. Basic scheme of the sorption experiment. Relative humidity ϕ is kept constant in the chamber surrounding the membrane, L is width of the membrane, λ is water content in the membrane, water transport is illustrated via flux jw and sorption variable s(t) representing normalized total membrane weight (defined in the text). Interfacial layer and the state in chamber is represented via boundary conditions and bulk diffusion via mass balance equation.

diffusivity Dw was used:

model, where the particular choice of the dependency of diffusion coefficient becomes irrelevant and where swelling is neglected. Both these models depend only on one free parameter describing interfacial transport once isotherm (a necessary input for both models), swelling with water content (required only for the full model) and the dependence of diffusivity on water activity (only for the full model) is provided from experiments. Subsequent analysis reveals the source of the disparity of the time scales between absorption and desorption for the lumped model. As a next step, the choice of an isotherm is discussed and the single parameter fitted. Additionally, we consider an experimentally assessed dependency of diffusion coefficient on water activity and swelling on water activity to close the full model and show that it has very similar predictions as the lumped model and compare the results with experimental data. Finally, application of the model to an active carbon material is proposed, where we provide indications that the presented qualitative insight is not limited in applications to PEM but rather to any sorption system with fast equilibration of sorbant in the bulk as is typically the case in thin membranes.

j Dw c w =⎛ w⎞ = L ww . R ⎝ Xw ⎠ jk = 0 ⎜

jw = −Dw (a w )·

c w ∂a w a w ∂x

(4)

which completes the mass balance Eq. (1) once functions c w (a w ) and Dw (a w ) are specified.1 The evolution equation, referred to as the full model in this article, is then described by following differential equation:

∂c w (a w (x , t )) ∂ c (a (x , t )) ∂a w (x , t ) ⎞ = − ⎛Dw (a w (x , t )) w w ∂t ∂x ⎝ a w (x , t ) ∂x ⎠

In this section a model of water sorption is developed. It focuses on modeling of an experimental set-up that is shown in Fig. 1. Water content λ is defined as the ratio of the number of moles of water in the membrane to the number of SO3 groups in the membrane (for details see Appendix A). Width of the membrane L is not constant as it can change with the amount of water in the membrane. General balance of water mass in the membrane reads [7]

⎜

⎟

(5)

accompanied by boundary conditions, following from the balance of mass at the interface,2

jw |0 = k w (ϕ − a w |0 )

(6a)

jw |L = k w (a w |L − ϕ),

(6b)

the phenomenological coefficient k w implicitly contains all the effects of the interface via the boundary condition and hence the diffusion coefficient corresponds to a bulk diffusion. Constant ϕ denotes relative humidity of vapor outside of the membrane. From the thermodynamic perspective these boundary conditions can be constructed from an exponential dissipation potential similarly to the Bulter-Volmer equation, see [10]. The unknown function in the full model is the water activity a w = a w (x , t ) . Alternatively, the evolution equation can be rewritten in terms of the experimentally accessible sorption variable s(t) defined as

(1)

where concentration c w is expressed in mol/ cm3 and water flux jw in mol/ cm2 s . On the level of mechanical equilibrium, total fluxes of protons and water in a membrane are proportional to the corresponding thermodynamic forces (gradients of electrochemical potentials at constant temperature). Assuming isothermal membrane and that proton transport due to coupling beetween proton and water flux is negligible in water sorption dynamics (or no electric current in the membrane [8]), the water flux is given by a constitutive relation provided by classical irreversible thermodynamics [9]

1 D c jw = L ww Xw = − L ww ∇μ w = − w w ∇μ w T RT

(3)

It expresses how the gradient of concentration affects diffusion flux, and is related to the ratio of water flux to its driving force while no other fluxes jk are present. For simplicity we shall restrict the experimental and modeling set-up to (approximately) 1D setting. Next, using this knowledge and the general formula for chemical potential of water in the membrane, water flux can be written (no coupling assumed or in the case of zero electric current) explicitly in terms of water activity a w as

2. Model formulation

∂c w = −div jw ∂t

⎟

1 Note that membrane is treated implicitly via responses to thermodynamic forces that are contained in the relations for fluxes (both in bulk and at the boundary). 2 Expressing boundary fluxes in terms of difference between thermodynamic forces, i.e. k͠ (μ w (ϕ) − μ w |0, L ) , would be more natural. At the same time, however, it can be equivalently expressed using activities (due to monotonicity of ln) and we prefer the latter for simpler analysis.

(2)

where Xw denotes the thermodynamic force being linearly proportional, with a constant L ww , to the flux jw it causes. Chemical potential of water per mole is denoted as μ w . In the equation above, the definition of 36

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

s (t ) =

λ (t ) − λ|t = 0 λ|t =+∞ − λ|t = 0

magnitude of the driving force Xw relative to a concentration difference (which is what is measured in sorption s). As we show below, however, one may reveal the precise relation between the isotherm shape and the absorption and desorption rates disparity based on our minimalistic thermodynamically consistent model. Consider an RH interval (ϕini , ϕ) for absorption, i.e. starting with an initial value ϕini (the membrane was initially in equilibrium with water vapor at RH ϕini) and reaching an equilibrium with external vapor at RH ϕ. Analogously, with desorption we consider a RH interval (ϕ, ϕini ) . Sorption curves are denoted as sabs, sdes , respectively. Let us first consider a linear dependence of water uptake on water activity, while assuming that a hysteresis does not occur,

(7)

To complete the description, the relation c w (a w ) can be obtained from the experimentally measured dependence of water uptake λ on water activity a w , see Appendix A and Section 3.2 for more details. The dependence of diffusion coefficient Dw follows from the experimental assessment of self-diffusivity which is subsequently used for calculating Dw, see Appendix B. Note that this constant k w represents the only free parameter of the model which is to be determined by comparing the resulting sorption curves to experiments. Swelling. Nafion membranes undergo swelling during the sorption processes. To include this phenomenon properly, it is convenient to rewrite the evolution equation into Lagrange coordinates corresponding to a dry membrane, see Appendix C. The resulting Eq. (C.13) is somewhat complex although simple in its essence and hence we choose not to include it in the main text.

λ (a w ) = λ1 a w + λ 0 . Thence Δλ = λ1 (ϕ − ϕini ) and the inverse function is equal to λ a−w1 (λ*) = (λ* − λ 0)/ λ1. The lumped model (10) is for both absorption and desorption equivalent to

d α sabs, des (t ) = (1 − sabs, des (t )). dt λ1

3. Lumped model It has been observed that concentration profiles are almost flat during sorption, see [11] and Discussion. This suggests that diffusion is fast with respect to the kinetics of sorption. As a result, we may further simplify the evolution Eq. (5) and get insight into its solution. Integrating the equation with respect to space yields an evolution equation for the total amount of water in the membrane:

∫0

L

∂ cw (aw (x , t )) dx = ∂t

∫0

L

∂jw (x , t ) ∂x

Therefore the sorption curves are also identical

α sabs (t ) = sdes (t ) = 1 − exp ⎛− t ⎞ λ ⎝ 1 ⎠ ⎜

and we may conclude that the model predicts absorption as fast as desorption if the dependence of λ on a w is taken linear and if there is no hysteresis in the considered segment of isotherm. More interesting is the behavior of the model when the dependence is nonlinear. Two criteria stating whether the model predicts faster absorption or desorption can be identified.

dx = (jw |0 − jw |L ) = kw (2ϕ − aw |0 − aw |L ),

(8) where boundary conditions 6 were used. Hence while assuming diffusion is fast enough to equilibrate water activity throughout the membrane (thus we may consider a w (x , t ) = a w (t ) ), Eq. (8) becomes

d 2k EW λ (a w (t )) = w (ϕ − a w (t )), dt ρdry L

• In the case without hysteresis in the isotherm curve the difference between rates can be characterized by the difference between the values of a weighed inverse function h(s) of water uptake on activity λ (a w ) , defined as:

(9)

which we shall refer to as the lumped model and where ρdry is the density of dry Nafion and EW is the equivalent weight, see Appendix A. Note that the relation between concentration and water uptake (A.1) was employed. The description is complete once a closure in the form of experimentally measured isotherm λ (a w ) is provided. Finally, as we are interested in sorption curves it is advantageous to formulate an evolution equation for sorption using the sorption variable s(t) (7). The resulting equation is:

d α s (t ) = (ϕ − λ a−w1 (Δλ·s (t ) + λ|t=0 ) ), dt Δλ where α =

2kw EW , ρdry L

⎟

h (s ) =

λ a−w1 (Δλ·s + λ (ϕini )) − ϕini ϕ − ϕini

with sorption parameter s ∈ [0, 1]. Particularly, if

h (1 − s ) + h (s ) > 1

∀ s ∈ (0, 1),

then absorption is slower than desorption and similarly if

h (1 − s ) + h (s ) < 1

(10)

Δλ = λ|t =+∞ − λ|t = 0 = λ (ϕ) − λ (ϕini ) with ϕini being

the initial water activity in the membrane and λ a−w1 (λ*) denotes the value of the inverse function of λ (a w ) at a point λ* . Note that this simplification completely removes the need of supplying Dw (a w ) . Hence the results of the below analysis of the lumped model are valid for all Dw (a w ) that meet the assumption of fast diffusion with respect to the sorption kinetics. Further we discuss the suitable choices of isotherms as the only necessary closure for the lumped model. Later on in the text we shall numerically compare the full model with the lumped model for a particular choice of Dw (a w ) and check the assumption behind the lumped model.

•

∀ s ∈ (0, 1),

then desorption is slower than absorption. As a special case, absorption is faster than desorption in such regions of isotherm where the dependency of water uptake on water activity is concave (and vice versa for convex isotherm). See Appendix D for the details about the derivation of this criterion. In the case of an isotherm with a hysteresis and when considering an activity interval where both segments of the isotherm are linear, the slopes of the isotherm are the decisive factor. With λd denoting the isotherm curve for desorption and analogically λa for absorption, the criterion reads if

|λd (ϕ) − λd (ϕini )| = |Δλd | < |Δλa | = |λa (ϕ) − λa (ϕini )|

•

3.1. Analysis of the lumped model Before adopting a particular isotherm we shall analyze the effect of the shape of isotherms on the disparity of the sorption rates. It is rather intuitive that the shape of the isotherm affects the absorption and desorption rates differently once we realize that water activity changes with concentration c w (and therefore with λ) and hence modify the 37

then desorption is faster than absorption (and again vice versa for the opposite inequality). See again Appendix D for details. For the general situation we are unable to state an intelligible criterion as the effects described above are competing. However, one can observe that this competition can result in changes of what is the faster process as time evolves. One has ha (s (ϕini )) = hd (s (ϕini )) = 0 , ha (s (ϕ)) = hd (s (ϕ)) = 1, and absorption is faster than desorption at a saturation s if

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

1<

1/|Δλa | 1 − ha (s ) . 1/|Δλd | hd (1 − s )

Hence at the onset of sorption experiments, t = 0 (s = 0), the sorption with flatter average slope, i.e. with larger 1/|Δλa, d |, is faster. On the other hand, the ratio of ad- to de-sorption rates in the ap1 / λ ′ (ϕ) 1 / λ ′ (t → ∞) proach to equilibrium is determined by 1 / λ ′ a(ϕ ) = 1 / λ ′a (t → ∞) d ini d which is the same as the ratio at the onset of sorptions but only for linear segments of isotherms and can be very different otherwise. Note that this is in accordance with the previous findings about the isotherm convexity in the case without hysteresis as the isotherm is a monotonically increasing function for stable processes and as a result the convexity of the isotherm translates into λ′ (ϕini ) < λ′ (ϕ) implying desorption faster than absorption. In summary, if the lumped model described by sorption Eq. (10), where diffusion is fast with respect to kinetics and hence without any effect on the behavior of the lumped model, is sufficient to recover the asymmetry between absorption and desorption, then the asymmetry is caused by the non-linear dependence λ (a w ) . Indeed, if the dependence were linear, the asymmetry would disappear.

Fig. 2. Fits of the considered isotherms (see text for details) to experimental data from [12,14]. They significantly differ only for small and large values of aw but which seem not to be particularly important for the considered data. This corresponding range of experimental data [1] is depicted as the shaded region, see Section 4 for more details. We choose the GAB isotherm for this article, see the text for discussion. (The reader is referred to the online version of this article for figures in full color).

2. a power law fit by Eikerling and Berg [2]

3.2. Sorption isotherm

λpower (a w ) = 3.0a w0.2 + 11.0a w4 ,

The analysis above indicates that the dependence of water uptake λ on activity a w is crucial for the asymmetry in sorption curves and hence the choice of a particular isotherm may affect the sorption rate predictions in a significant way. Further, the role of an isotherm in a phenomenologically based model(the second group of models from the introduction) is that it summarizes many physical and chemical phenomena into a single empirical relation to allow keeping the model simple. Let us, therefore, consider several examples of isotherms. We shall consider a cubic fit [12,6] and power law used for PEM [2], finite layer BET isotherm [13] and a GAB isotherm that we use in the remaining of the article. The choice of the isotherm is also linked to the question of data selection. Two different groups of the experimental results can be identified, one where λ ≈ 14 for saturated water vapor, i.e. ϕ = 100%, see e.g. [12], and the other where λ ≈ 21, see e.g. [14]. Water uptake of a Nafion membrane equilibrated with liquid water is approximately 22, see [15], which is compatible with the latter group, while the former group probably exhibits the so called Schröder's paradox. Moreover, Onishi [16] demonstrated that thermal pretreatment of the membrane affects the maximum water uptake dramatically (λ ≈ 14 for a high temperature pretreatment). In this work we prefer the experiments without Schröder's paradox, i.e. with λ ≈ 21 for saturated water vapor. Particularly, water uptake is adopted from [14] for high RH and from [12] for low RH. As Jeck et al. [14] pointed out, a substantial part of water uptake from vapor appears for very high RH, which is the reason for such combination of sorption curves, see Fig. 2 for the data. Note that according to this paper, Schröder's paradox is just an experimental artefact caused by not reaching sufficiently high RH. Eikerling and Berg [2,3] also believe that the Schröder's paradox does not exist. However, their offered resolution of this paradox is based on a model where they use continuum theory with simple capillarity for pore formation with sizes ∼ 1 nm and the assumption of the presence of capillary pressure between liquid water in pores and saturated vapor in outer region. Finally note that as follows from the above analysis, the particular choice of dataset does not affect the qualitative behavior of the model. We consider the following isotherms:

3. a finite layer BET isotherm fit by Thampan et al [13] with physically meaningful parameters (water loading at monolayer coverage λm = 1.8, total number of water layers in the pores at saturation nL = 13.5 and one fitting parameter c = 150)

λBET (a w ) = λm

ca w [1 − (nL + 1) a wnL + nL a wnL + 1] . 1 ( − a w )[1 + (c − 1) a w − ca wnL + 1]

(13)

4. a GAB isotherm with λG representing the monolayer capacity, cG, kG being energetic constants3 and kG a measure of the difference of free enthalpy of the sobate molecules in pure liquid and sorption state in layers above the monolayer determined via the least squares method

λGAB (a w ) =

λG cG kG a w (1 − kG a w )(1 + (cG − 1) kG a w )

(14)

with fitted values kG = 0.9 ∈ (0.5, 1), λG ≈ 1.93 and cG ≈ 44.3. Fig. 2 shows all four isotherms. All capture the considered experimental data reasonably well. Further, it can be easily shown that all the above considered isotherms satisfy the first criterion from Section 3.1 meaning that it predicts faster desorption than absorption in accordance with the experimental observations. We choose the GAB isotherm because it in addition best fits the experimental sorption rates data as we show below in Section 4. 4. Quantitative modeling of water transport In this section we shall evaluate the proposed full and lumped model together with checking the consistency of the used approximation yielding the lumped model. Results are compared with experimental data from Satterfield et. al. [1] where RH was ramped up as a step function from low humidity to almost fully saturated water vapor.4 Absorption in the calculations started with a membrane equilibrated 3 Measures the difference of chemical potentials of the sorbate molecule in the upper sorption layers and in the monolayer. 4 Membranes were suspended in a controlled atmosphere container while their weight was recorded on a balance. For absorption, ϕ > 0.95 was sustained by hanging membrane just above water level. For desorption, ϕ < 0.05 was sustained in an Erlenmeyer flask partially filled with Drierite. For more details about membrane preparation etc. see [1].

1. a cubic fit by Springer et al. [12,4]

λ cubic (a w ) = 0.043 + 17.81a w − 39.85a w2 + 36.0a w3 ,

(12)

(11)

38

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

Fig. 3. Absorption and desorption in a N115 membrane as predicted by the lumped model for (A) cubic isotherm (11), (B) power law isotherm (12), (C) finite layer BET isotherm (13) and (D) GAB isotherm (14). Although all these fits are capturing the sorption isotherm well, see Fig. 2, their ability to describe the experimental sorption data is very different as can be seen in mol these plots. Fitted values of the single fitted parameter (in 10−7 2 ) are kw = 3.8 for cubic, kw = 3.8 for power law, kw = 4.0 for finite layer BET, and kw = 4.4 for GAB isotherm. We chose cm s

the GAB isotherm (D) for the remainder of this article. Note the one fold difference in sorption rates and see Fig. 6 for full sorption and desorption curves. Legend: fitted sorption curves from the lumped model are displayed as solid smooth curve for absorption, dashed curve for desorption; fuzzy solid lines represent experimental data taken from [1]; a linear curve best approximating desorption data with indicated slope highlights the exponential nature of desorption, contrasting with a clear non-exponential behavior of absorption. (The reader is referred to the online version of this article for figures in full color).

the above fitted value of k w the same, to different membrane thickness and comparing it to experimental data from Satterfield et al. [1]. Figs. 4 and 6b illustrates that model predicts the sorption kinetics well5 .

with water vapor at a low RH ϕini by exposing it at time t = 0 to a controlled environment with RH ϕ > ϕini . We consider the equilibrium to be reached if the amount of water in the membrane changed less than 1% over the period of 4000 s. After the equilibrium was reached the RH of the surrounding water vapor was then suddenly changed to the low value ϕini and the desorption started. Note that the RHs were calculated for each experiment from the reported values of water uptake λ at the beginning and at the end of the experiment using the considered isotherm rather than using the reported values of RH measured during preparation of membranes. These values are different from the values reported in the article [1], which we attribute to a certain delay between the start of sorption and the start of the measurement and also probably to uneven distribution of RH in the chamber. To minimize these unknown influences we propose to use the membrane sample weight as the only measured parameter not only to track sorption but also to set the initial and boundary conditions consistently. RH can be calculated back from a monotonous isotherm if required.

4.2. Full model For the purposes of the full model we consider the measurement of water self-diffusivity in Nafion by Benziger et al [17] and by pulse gradient spin-echo nuclear magnetic resonance (PGSE NMR) [6], where the experimental data were fitted by the curve 0 2 Dwself = D∞ aw,

0 D∞ = 0.265 exp ⎛− ⎝

3343 K cm2 ⎞ . T ⎠ s

(15)

With this knowledge we can use the relation between self-diffusivity and diffusivity, Appendix B, to complete the full model. Note again that the particular functional dependency Dw (a w ) is not important for the dynamics of sorption once diffusion is fast. We shall test this below both by checking the dynamics of concentration profile within membrane and also by the similarity of the solutions of the full and the lumped model. The full model with swelling (C.13), 6 accompanied by relation (B.10) between self-diffusion (15) and diffusion and with the chosen

4.1. Lumped model We have solved the lumped model (10) for the four considered isotherms using the NDSolve solver in Wolfram Mathematica 10. Absorption started with λ 0 = 2.36 and finished at λ1 = 13.12 . Desorption started after the membrane with λ1 = 10.34 was moved into a drying chamber (finished at λ 0 = 1.79). Consequently we chose the GAB isotherm with the best fit to experimental data, see Fig. 3, while assessing the value of interfacial mol transport coefficient to be k w = 4.4·10−7 2 . Finally, we used this value to cm s test our model by applying it, i.e. while keeping the sorption isotherm and

5 As mentioned in the paragraph above, these initial conditions were calculated from experimentally measured λ values via the chosen isotherm for consistency (one measures λ not ϕ to assess initial condition and only consequently via chosen isotherm ϕ is determined). If we use instead (1.79; 13.12) as initial and boundary conditions for λ in both ab- and de-sortpion, ie. desorption directly follows absorption, model yields very similar results (not shown) with only slightly changed value of fitted parameter kw. This observation further supports the above approach of calculating λ0 and λ1 separately.

39

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

Fig. 4. Prediction of absorption and desorption kinetics of the lumped model for N1110 Nafion (B). Note that the single parameter of the model has been assessed from sorption data for N115 (A) and kept fixed. Again note the one fold difference in soprtion rates. (The reader is referred to the online version of this article for figures in full color).

Fig. 5. Absorption and desorption in a N115 membrane as predicted by the full model with swelling for GAB isotherm versus experimental data (A). Concentration profiles within membrane in few indicated time steps all showing rapid flattening out of water activity profiles (B). (The reader is referred to the online version of this article for figures in full color).

GAB isotherm is numerically solved for a thin membrane by finite element method. Initial conditions for absorption and desorption were chosen the same as in the case of the lumped model above. Fig. 5 presents absorption and desorption curves for GAB isotherm from the full model with swelling compared with experimental data from [1]. The fitted value of k w for the desorption data for T = 30° is mol k w = 3.5·10−7 2 . The figure shows qualitative and quantitative agreecm s ment, while the rapid flattening out of the concentration profiles within the membrane is apparent confirming fast diffusion with respect to the sorption kinetics. Therefore the assumption behind the lumped model is plausible and should yield similar predictions as the full model. The fitted values of interfacial transport coefficient k w are slightly different for the lumped and the full model with swelling. We account this difference to explicit inclusion of the diffusion process that slows down the transport within membrane and hence reduces the value of normalized membrane weight s. Hence the fitted value of k w in the lumped model should be considered as an effective value when approximating diffusion as being fast. Results in Fig. 6 strongly suggest that the lumped model is sufficient to recover the disparity of the sorption timescales and shows a qualitative and a quantitative match with the full model. Hence we conclude that our working hypothesis obtained from the analysis of the lumped model, which states that the cause of the disparity can be identified as the non-linear dependence of λ on a w , is a plausible explanation of the origin of the different sorption rates.

applicability of the lumped model is not restricted to water transport in PEM. It should apply to any sorption process in a sufficiently small specimen from a gaseous phase (large diffusion coefficient) so that the assumption of fast diffusion yielding the lumped model is a plausible one and hence the actual dependence of diffusion coefficient on state variables is insignificant. The full model with diffusion should be even more general but requires specification of D (a w ) which is not easily experimentally accessible. Let us discuss a reported experiment of water vapor sorption on activated carbon [18]. They measured adsorption and desorption rates for eight water activity ranges (in eight segments of isotherm) instead of one overall step from 0% to 100% RH or vice versa. They used a Ficks’ diffusion model with a constant coefficient and Dirichlet boundary conditions to interpret their findings and to model the performed sorption experiments. Limitations of the considered model are twofold: one, although Dirichlet boundary conditions when accompanied by appropriate nonlinearities in diffusion coefficient can yield distinct adand de-sorption rates [19], we believe that Dirichlet boundary conditions do not apply to sorption experiments (and additionally, the nonconstant diffusion coefficient was not considered [18]); two, Fickian diffusion (constant diffusion coefficient) with plausible boundary conditions was shown being incapable of capturing the observed disparity of time scales of ab- and desorption, see [20] Appendix E. Nevertheless, we may use the measured characteristic times of the absorption and desorption rates, see Table 1, for discussion of the applicability of our lumped model. Further, we discuss below that if one uses the presented lumped model to interpret the experimentally measured isotherm curve, it yields results that qualitatively agree with the observed behavior of sorption rates and provide feedback to experiment.

5. On generality of presented sorption model. Adsorption on active carbon. Isotherm with hysteresis The aim of this section is to show that the plausibility and 40

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

Fig. 6. Comparison of the numerical solutions of the lumped (orange) and full (black) model with swelling. Apart from the fact that the full model explicitly considers diffusion and swelling, they only differ in the fitted value of interfacial transport coefficient kw , see the text for discussion of this observation. The results suggest (A) that the two sets of the solutions are very similar which supports our approach to approximation of the full model by the lumped model and suggests that the qualitative insight from the lumped model is applicable to the full model as well. Fig. B shows sorption plots with experimental data (green for desorption, red for absorption) together with both the lumped and the full model solutions. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

analysis for this section:

Table 1 A table of characteristic times for eight segments of activated carbon isotherm from [18] and comparison with qualitative predictions from the lumped model. The third column lists the reported characteristic times, where the faster sorption is noted by a letter (A stands for adsorption and D for desorption). The fourth column lists a summary of qualitative insight from the lumped model for each segment of isotherm with hysteresis, Fig. 7A, see the summary in the text that lists the conclusions that are being used.

Nr.

pressure range (kPa)

measured characteristic times A (s) ≷D (s)

prediction of the lumped model

1 2 3 4

0–0.3 0.3–0.61 0.61–0.91 0.91–1.22

87 < 303 A 138.9 > 114 D 226.2 < 255.1 A 565 < 621.1 A

5

1.22–1.52

2702.7 > 1818.2D

6

1.52–1.82

1351.4 > 1075.3D

7 8

1.82–2.11 2.11–2.41

505.1 > 381.7 D 370.4 > 318.5 D

A (linearity) A (competition) A (linearity) A significantly faster initially; D eventually (competition) identical rates initially; D eventually (competition) D significantly faster initially; A eventually (competition) D (linearity) D (linearity)

• in linear segments of an isotherm desorption is quicker where the slope of desorption isotherm is flatter; • in nonlinear segments of isotherm there are at least two character-

istic times of sorptions and as a result desorption is faster initially if the average slope of desorption isotherm is smaller while desorption is faster in the approach to equilibrium if the slope of desorption isotherm at equilibrium is flatter.

With these information in mind about the relation of isotherm slope and rate of ad- and de-sorption in the case of hysteresis, we can estimate regions of the isotherm where adsorption is faster than desorption based on the lumped model. We number the studied intervals of water activity in an increasing manner and get the following predictions:

• It seems reasonable to assume linear isotherm in the regions 1, 3, 7, 8 and hence the flatter sorption is the faster one. • In the fifth region where the overall change of water uptake is equal

We shall use the results of the qualitative analysis of the lumped model to predict behavior in the active carbon sorption experiment based on the experimentally measured isotherm displayed in Fig. 7. For convenience we summarize the main implications of the previous

•

for both ad- and de-sorption, we anticipate desorption to be slightly faster due to the faster approach to equilibrium as can be seen from the anticipated slopes at equilibrium (at points 5a, the equilibrium point for adsorption in region 5, and 5d, the equilibrium point for desorption in region 5, in Fig. 7). Similarly in the fourth interval, we expect the adsorption to be faster

Fig. 7. Fig. A shows the experimentally measured water adsorption isotherm of activated carbon at 295 K in eight regions (numbered with shaded numbers) with highlighted estimated slopes close to equilibriums for non-linear regions that determine if ad- or de-sorption occurs at a faster rate. According to the lumped model, in linear regions of isotherm, the sorption with flatter slope is faster. In nonlinear regions (4,5,6th) this holds only initially and at later stages faster sorption is the one with the flatter slope at equilibrium. Consult Table 1 for prediction of convex and concave regions of the isotherm and the text for the relation of the isotherm slope and the rate of sorption. Figure B shows experimental data for desorption in the fifth region and note the presence of two characteristic times in accordance with results from the analysis of the lumped model.

41

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

mentioned sorption phenomena. Similarly, the scaling property of water sorption curves (dependency on t / L ) has already been captured by several authors e.g. [4,5,22]. These studies, however, do not provide a satisfactory explanation of the difference in sorption rates and of the non-exponential behavior of absorption [1]. Further, Satterfield et al. proposed a model that is able to predict both the disparity of the time scales and its magnitude [1] but via an addition of an empirical exponential term (and hence effectively introducing a second characteristic time) to the solution of a Fick's model and attributed its presence to the swelling of the membrane. Hallinan et al. [23] on the other hand proposed two mechanisms (diffusion-reaction and diffusion-relaxation) leading to non-Fickian behavior of sorption, focusing on the initial time lag of the process, which is not discussed here. The model, however, does not explain the disparity of time scales between absorption and desorption. Finally, Goujot and Vitrac have recently studied nonlinear adsorption isotherms via analytical solutions to mass diffusion problems [24]. They considered Fick's diffusion with constant diffusion coefficient for bulk transport accompanied by boundary conditions and involving piecewise linear isotherms. We adopted a different approach where the bulk equations are consistent with nonequilibrium thermodynamics without the assumption of dilute solutions, boundary conditions follow from the balance of mass (where all the boundary and interfacial effects are reduced to a single parameter k w ) and involving nonlinear isotherms. Yet, as we show above, both the full and the lumped model show a good agreement with experimental data and additionally analytical insight in the model behavior can be gained. The simplification of the full model yielding the lumped model is based on the assumption of fast diffusion within the membrane. According to [11], the water concentration profile is flat during the sorption process as determined by SANS. Therefore the sorption is either controlled by the interface properties or by the experimental conditions. In both cases, it does support the admissibility of the lumped model approximation based on fast diffusion yielding flat concentration profiles. Additionally, numerical solutions of the full model confirm rapid transition from the initial non-uniform concentration profile to the resulting almost constant solution, see Fig. 5b. As the lumped model is independent of the particular choice of diffusion coefficient dependency and is a good approximation of the full model once diffusion is a fast process, all the qualitative results hold for any such functional dependency. We adopted the conclusion of [6] where D ∝ a w2 6 but QENS study or conductivity measurements through the Nernst-Einstein equation suggest almost linear dependency with respect to a w . Such distinctions should not be important in our study and for the obtained model due to the observed fast diffusion. Finally, this highlights the rather intuitive fact that for fast bulk diffusion the interfacial transport (represented by k w and BCs) in sorption kinetics is indeed the rate limiting process as observed before [6,8]. Both presented models have only a single fitting parameter k w . However, with a single value it is able to capture sorption kinetics for various membrane thicknesses. The fitted value for the lumped model is mol mol k w = 4.4·10−7 2 and k w = 3.5·10−7 2 for the full model with swelling cm s cm s when using the data for T = 30°. Monroe et al. [5] estimated the value mol of a ‘vaporization exchange rate’ as k w = 3.2·10−6 2 for T = 50° and cm s discussed that it can be judged accurate within an order of magnitude. Note that their model and experimental set up included flow in chamber surrounding the membrane with gas flows which may be

initially but with slightly faster desorption eventually.

• In the second interval, desorption can be considered linear but likely

not the adsorption isotherm where significant concavity is expected. However, adsorption is expected to be faster both initially and near the equilibrium based on the model prediction. Finally, in the sixth interval lumped model provides a prediction that desorption is significantly faster initially but slower as the equilibrium is approached.

These observations are summarized in Table 1 where we also list the measured characteristic times of ad- and de-sorption from [18]. Note that in the nonlinear regions this “competition”s of convexity and overall slope makes it difficult to predict the overall relation of sorption rates just from the qualitative insights provided by the lumped model analysis. Additionally, as was discussed above, it is not sufficient to characterize the sorption process by a single characteristic time and more datapoints are necessary - either in time or to split non-linear regions into smaller ones where linearity to a certain degree can be assumed. Nevertheless, a reasonable agreement of the sorption characteristic times and the presented model has been reached demonstrating that the application of the presented model is not restricted to PEM. 6. Discussion and conclusions Our aim was to develop a macroscopic approach and to identify a minimalistic but physically sound model that has a clear interpretation and may help to identify the causes of the observed phenomena. The toll that is paid for this approach is that only some physical and chemical phenomena are captured by the model. Many relevant processes are omitted or only implicitly included in the used empirical relations (e.g. the sorption isotherm). As mentioned in the introduction, there are models taking a complementary approach to the studied phenomena by formulating actual physical and chemical phenomena taking place on the microscale and inferring the macroscopic behavior of sorption kinetics and water transport [2,3]. This bottom-up approach faces different challenges like the plausibility of the usage of continuum description for pore formation with sizes ∼ 1 nm or what physical phenomena should be considered at this lenghtscale (capillarity is probably too crude) but offers a promising way of addressing the discussed problems. The macroscopic framework that we chose allows analytical insight and hence indicating the effect of the considered physical phenomena on the experimentally observed data. In this sense we showed that neither swelling nor diffusion (and its functional form) plays an important role in the sorption characteristics and that the nonlinearity of sorption isotherm is the likely reason behind the disparity of the sorption time scales. We briefly summarize properties of the proposed models in this article in relation to sorption in Nafion: 1. sorption curves for different t widths of membranes collapse into one when plotted against L (a direct consequence of the lumped model and visible from the agreement of the model for experiments with different widths of membranes); 2. absorption is approximately ten times slower than desorption; 3. the likely cause of this disparity in sorption time scales is linked to nonlinearity in sorption isotherm; 4. for description of absorption at least two characteristic times are necessary. Several other models have already been proposed that capture (water) sorption kinetics but, at least to our best knowledge, none is possessing all these four properties nor providing an explanation for disparity of time scales. We shall briefly discuss the most relevant ones. Some authors propose a hypothesis that the disparity of time scales can be explained via the choice of boundary conditions (which demonstrates the importance of the choice of appropriate boundary conditions in a model formulation) by setting unequal rates of flow over boundary in the case of absorption and desorption [4,21]. This approach has been discussed in detail by us in [20], summary of which is given in Appendix E. The conclusion is that even though this hypothesis could explain the disparity of the time scales, it is unable to explain the other

6 The diffusion coefficient was measured by means of PGSE NMR. The experimental value differs with chosen diffusion delay time. In the limit of long diffusion delay times, water molecules have enough time to pass through the porous network, and the heterogeneous morphology of the membrane is thus taken into account. Estimating the average distance a molecule travels during the PGSE NMR measurement (2 = Dt ) in the long-time limit (t = 1 s ), the distance is 23μm , which is lower than dry thickness of N115 membrane (121μm ). Therefore, the measured diffusion coefficient should reflect bulk of the membrane rather than the interfaces and the heterogeneous morphology of the membrane is thus taken into account. Therefore the way we implement this diffusion coefficient into the full model is consistent with the experimental technique used for its assessment.

42

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

responsible for the difference in the fitted values. Similarly, Zhao et al. [6] reported the value of a ‘lumped mass transport coefficient’ for water evaporation rate with liquid water on the other side of the membrane mol k w = 4.48·10−5 2 at T = 50°. Given that k w is highly temperature decm s pendent (value in Zhao et al. study nearly doubled when raising the temperature from T = 30° to T = 50°) and that the here considered experimental data were from static (not dynamic) vapor sorption, we believe that the here estimated value of k w is plausible. The lumped model analysis and results matching both experimental observation and data are promising but the model may be easily considered as a too crude description of sorption process in a membrane. For this reason we presented also a model that we refer to as the full model which additionally includes both swelling and diffusion coefficient dependency on water activity. As we have shown in this article, the agreement between the full and the simplified lumped model is very good (at least within the studied experimental ranges of parameter values) which suggests that the detailed and uneasily accessible functional dependency Dw (a w ) is redundant for description of sorption in thin Nafion membranes. Similarly, the role of swelling seems to be of secondary importance7 although it has been previously suggested as a possible reason for the difference in rates of absorption and desorption [1]. An article supporting our observation is by Gebel et al. [11] where they suggested that water transport in Nafion can be governed by a flow over the boundary. They found out using a SANS method that one can identify a single characteristic time-scale of water sorption in bulk of the membrane. Consequently, water transport is suggested to be Fickian inside the membrane and strongly influenced by water transport through the membrane/gas interface as the strong disparity of time scales of absorption and desorption is well known. This is not in accordance with the hypothesis that water sorption is governed by Fickian diffusion and polymer relaxation since another time-scale corresponding to the polymer relaxation would be present. Nevertheless, such hypothesis cannot be ruled out yet as it can be dependent on the region in physical parameter space. Note, however, that the role of swelling could become substantial as it could, in certain situations, result into a built up of an

inhomogeneous internal pressure that in turn would act as an additional force for water transport. Although contribution of such force can be in principle included in the proposed model via CIT, we have not pursued this idea due to the complexity of the resulting model but will be subjected to further research. As both qualitative and quantitative results of the lumped and the full model are very similar to each other and match the experimental sorption curves rather well, we suggest that the reason for disparity of the time scales in ab- and de-sorption is due to the nonlinearity of isotherms as was shown to be the case for the lumped model. Detailed experimental verification of this theoretical result will be a subject of our future research. Finally, we demonstrated in the activated carbon example that the applicability of the presented qualitative insights may not be restricted to the studied sorption in Nafion but rather to any sorption process in a system with fast equilibration of sorbant in the bulk as is typically the case in thin membranes.8

Acknowledgements We are grateful to Frantiek Mark for revealing the world of thermodynamics to V.K., M.P. and J.K. and for generously supporting this research. The work was partially developed within the POLYMEM project, reg. no CZ.1.07/2.3.00/20.0107, that is co-funded from the European Social Fund (ESF) in the Czech Republic: “Education for Competitiveness Operational Programme”, from the CENTEM project, reg. no. CZ.1.05/2.1.00/03.0088, co-funded by the ERDF as part of the Ministry of Education, Youth and Sports OP RDI programme and, in the follow-up sustainability stage, supported through CENTEM PLUS (LO1402) by financial means from the Ministry of Education, Youth and Sports under the “National Sustainability Programme I”. Further, the study was partly supported by the Charles University in Prague, project GA UK No 70515. The work was supported by Czech Science Foundation (project no. 14-18938S)

Appendix A. Equilibrium water sorption in Nafion membranes The amount of water in a Nafion membrane is usually expressed in terms of water uptake λ as defined in Section 2. Water uptake can be converted to concentration of water cw (with units mol cm−3 ) via

c w = λ·

1 ·ρ EW N

(A.1)

where EW is the equivalent weight, which expresses weight of the Nafion polymer per 1 mol of SO3 groups, see [29], and ρN is density of the Nafion polymer in the swollen membrane, i.e. mass of the polymer divided by volume of the membrane. This equation can be used to relate cw to a w as λ = λ (a w ) is experimentally accessible (isotherm). Global relation between water concentration and total water content (which is being measured in sorption experiments) then follows

1 L

∫0

L

c w dx = λ·

1 ·ρ . EW N

(A.2)

Assuming isotropic swelling, the experimental dependence of strain εon relative humidity (or activity of water in the membrane) was published in [6]. The experimental data can be fitted as

ε (a w ) = 0.168426a w − 0.314683a w2 + 0.287563a w3 .

(A.3)

Volume of the membrane is related to volume of the dry membrane through

V = Vdry (1 + ε )3,

(A.4)

and density of the polymer itself, i.e. without water, can be calculated from the dry density as

ρN =

ρdry (1 + ε )3

.

(A.5)

7 However, to soundly compare this hypothesis with experiments it would be necessary to formulate a mathematical model coupling elasticity or viscoelasticity of the polymer with diffusion of water. Such a model lies beyond classical irreversible thermodynamics, presented for example in [9,25], but also beyond the usual modeling of transport processes in chemistry, for example [26], although several models describing diffusion within viscoelastic porous media could be adopted, see for example [27,28]. 8 In fact Goujot and Vitrac in their similar analytical study of piecewise linear isotherms were motivated by assessment of the risk of food contamination by substances originating from polymer materials in contact [24].

43

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

Having the swelling and water uptake as functions of water activity (or relative humidity) it is now easy to calculate thickness and density of dry Nafion membranes. Such calculation is needed because thicknesses and densities of the membranes are usually measured at RH = 50%. For example in [30] thickness and density of Nafion 115 were determined at RH = 50% and 23°C to be ρN115 = 1.969 gcm−3 and tN115 = 127μm . Nafion N1110 has the same density as N115 and its thickness is tN1110 = 254μm . Using formula (A.3), thickness of dry Nafion N115 is given by

tN 115 = 121.6μm , 1 + ε (0.5)

tN 115, dry =

(A.6)

and thickness of dry Nafion N1110 is

tN 1110, dry =

tN 1110 = 243.3μm . 1 + ε (0.5)

(A.7)

To assess density of dry Nafion, one can start with realizing that mass of water in the swollen Nafion can be expressed as, see (A.1),

m w = Mw

λ λ mN , dry = Mw ρ Vdry EW EW dry

(A.8)

where ρdry and Vdry are density and volume of the dry Nafion, respectively, and Mw is the molar mass of water. Vdry is related to volume of the membrane at RH = 50% by

VRH = 50% = (1 + ε (0.5))3Vdry.

(A.9)

Mass of the swollen Nafion is then equal to the sum of mass of the polymer and mass of water contained in the swollen in the polymer

ρRH = 50% VRH = 50% = ρdry

VRH = 50% M λ (0.5) VRH = 50% + w ρdry (1 + ε (0.5))3 EW (1 + ε (0.5))3

(A.10)

so that density of the dry Nafion is

ρdry =

(1 + ε (0.5))3ρRH = 50% 1+

Mw λ (0.5) EW

. (A.11)

Density of dry Nafion N115 is equal to

ρN 115, dry = ρN 1110, dry = 2.12

g , cm3

(A.12)

which is also equal to density of dry Nafion N1110 as these two types of Nafion only differ by thickness, see [29] for the naming conventions. Appendix B. Diffusivity vs self-diffusivity Several theories have been proposed to convert self-diffusivity into the thermodynamic diffusivity, see [31–33]. A mutual diffusion coefficient defined as

jw = −Dmutual

∂c w ∂x

(B.1)

was introduced in [31] (Eq. (5c)). The mutual diffusion coefficient for a polymer-solvent mixture was estimated by Vrentas and Duda in [32], Eqs. (25) and (26) therein, as

Dmutual =

Dwself cN vN c w ⎛ ∂μ w ⎞ RT ⎝ ∂c w ⎠T , p ⎜

⎟

(B.2)

where cN and vN are concentration and partial molar volume of the polymer. By combining this last equation with Eq. (B.1) and assuming only isothermal and isobaric conditions, we obtain that

jw = −

Dwself cN vN c w ∂μ w ∂x RT

(B.3)

and the thermodynamic diffusion coefficient from Eq. (2) thus becomes

Dw = Dwself cN vN .

(B.4)

Let us now calculate the quantity cN vN . It follows from the assumption of local thermodynamic equilibrium, see [9,25] or [34], that (B.5)

1 = c w vw + cN vN ,

where concentrations are expressed with respect to the total volume of the membrane (i.e. water and Nafion together) and partial molar volumes are defined as

∂V ⎞ vα = ⎛ , ∂ ⎝ nα ⎠T , p, nβ ≠ α ⎜

⎟

α, β ∈ {w, N }, (B.6)

where nw and nN being number of moles of water and Nafion in the membrane, respectively. Partial molar volumes express how the total volume (of the whole membrane regarded as a mixture of water and Nafion) changes when a small amount of a constituent (water or Nafion) is added while keeping temperature, pressure and amount of the other constituent constant. From Eq. (B.5) it follows that the sought quantity cN vN can be calculated once we know c w vw . The latter can be extracted from experimental data on membrane swelling. Indeed, volume of the whole membrane is given by Eq. (A.4). Taking derivative of the whole volume with respect to mass of water leads to 44

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

∂ε ⎞ vw = Vdry·3(1 + ε )2 ⎛ . ⎝ ∂n w ⎠T , p, nN ⎜

⎟

(B.7)

From the relation between water content and concentration, (A.1), it follows that

nw =

λ ρ Vdry EW dry

and c w = λ

ρdry

1 . EW (1 + ε )3

(B.8)

Therefore,

∂ε ∂λ ∂ε EW ⎛ ∂ε ⎞ = = . ∂ n ∂ λ ∂ n ∂ λ ρdry Vdry w ⎝ w ⎠T , p, nN ⎜

⎟

(B.9)

Plugging these relations back into Eq. (B.7) leads to

c w vw =

3 ∂ε 3 ε′ (a w ) λ = λ , 1 + ε ∂λ 1 + ε λ′ (a w )

(B.10)

where the experimental dependencies ε (a w ) and λ (a w ) are employed. The right hand side of this last equation is between 0 and 0.15 for all activities, and thus the term cN vN can perhaps be approximated by unity, i.e. cN vN ≈ 1. In such situation, the thermodynamic diffusion coefficient is approximately equal to self-diffusivity,

Dw ≈ Dwself .

(B.11)

We shall, however, use the full form of the revealed relation (B.10) for calculations based on the full model. Appendix C. Transport of water with swelling C.1. Lagrangian coordinates The membrane is supposed to be a very thin rectangular parallelepiped, which is, however, always the case for Nafion membranes, and it is homogeneous in the longer directions, let us say y and z, so that a w or any other physical quantity only vary in the x direction (in which the membrane is thin). Therefore, the sorption can be considered as a one-dimensional problem. The membrane is swelling when absorbing water. For the sake of simplicity of numerical simulation, it is advantageous to transform the spatial (Eulerian) coordinates to Lagrangian coordinates attached to the dry membrane. The physical spatial (Eulerian) coordinates have already been denoted by x, y and z, and let us denote by xm(t), ym (x , t ) and z m (x , t ) boundaries of the membrane. In other words,

x ∈ [0, x m (t )], y ∈ [0, ym (x , t )]andz ∈ [0, zm (x, t)].

(C.1)

The Lagrangian coordinates (C.2)

X ∈ [0, Xm ], Y ∈ [0, Ym]andZ ∈ [0, Zm] which are attached to the membrane in the dry state are related to the Eulerian coordinates as follows

x (X , Y , Z , t ) =

y (X , Y , Z , t ) =

z (X , Y , Z , t ) =

∫0

X

∫0

Y

∫0

Z

(1 + ε (Aw (X ′, t ))) dX ′

(C.3)

(1 + ε (Aw (X , t ))) dY ′ = (1 + ε (Aw (X , t ))) Y

(C.4)

(1 + ε (Aw (X , t ))) dZ ′ = (1 + ε (Aw (X , t ))) Z

(C.5)

where Lagrangian activity of water is defined as (C.6)

Aw (X , t ) = a w (x (X , t ), t ). C.2. Lagrangian concentration Lagrangian concentration of water in the slice of the membrane with Lagrangian coordinate X at time t is equal to

C (X , t ) =

∫0

ym (x (X , t ), t )

∫0

zm (x (X , t ), t )

= c w (x (X , t ), t )(1 + ε (Aw (X ,

c w (x (X , t ), y, z , t ) dy dz

(C.7)

t )))2 Y ·Z

ΔA ⏟

(C.8)

where area of the dry membrane ΔA was introduced. Lagrangian area concentration of water is then defined as

Cs (X , t ) =

C (X , t ) = c w (x (X , t ), t )(1 + ε (Aw (X , t )))2 . ΔA

(C.9)

Concentration c w (x , y, z , t ) can be expressed in terms of water activity a w . Indeed, as c w (x , y, z , t ) gives the number of moles of water per cm3 of the swollen membrane, it can be related to water uptake λ by Eq. (A.1), and water uptake λ can be expressed in terms of activity according to a chosen isotherm, e.g. (14) Converting also density of the polymer into density of the dry Nafion by using Eqs. (A.1) and (A.5) then leads to the Lagrangian expression for concentration

45

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

ρdry λ (Aw (X , t )) . (1 + ε (Aw (X , t )))3 EW

c w (x (X , t ), t ) =

(C.10)

Comparing this formula with Eq. (C.9) yields a relation between the Lagrangian area concentration and the Lagrangian activity

Cs (X , t ) =

λ (Aw (X , t )) ρdry . 1 + ε (Aw (X , t )) EW

(C.11)

Finally, the total amount (in moles) of water in the membrane can be calculated as follows

nw = = =

∫0

xm (t )

∫0

Xm

ρdry EW

y (x , t )

∫0 m

∫0

zm (x , t )

c w (x , y, z , t ) dx dy dz = ∫0

xm (t )

C (X (x , t ), t ) dx =

∂x X dX = ∫0 m C (X , t )(1 + ε (Aw (X , t ))) dX = C (X , t ) ∂X

∫0

Xm

λ (Aw (X , t )) dX ·ΔA.

(C.12)

C.3. Lagrangian evolution equation Let us now rewrite the evolution Eq. (1) in the Lagrangian coordinates. Plugging relations (4) and (C.10) into Eq. (1), the evolution equation becomes

ρdry ∂⎛ λ (a w (X , t )) 1 ∂ ⎛ λ (Aw (X , t )) Aw (X , t ) ρdry ∂Aw (X , t ) ⎞ ⎞= Dw ∂t ⎝ (1 + ε (Aw (X , t )))3 EW ⎠ 1 + ε (Aw (X , t )) ∂X ⎝ (1 + ε (Aw (X , t ))) 4 EW ∂X ⎠ ⎜

⎟

⎜

⎟

(C.13)

where derivatives with respect to x were replaced with derivatives with respect to X by formula

∂ ∂X ∂ 1 ∂ = = . ∂x ∂x ∂X 1 + ε ∂X

(C.14)

It is, moreover, convenient to introduce a dimensionless coordinate

= X . X Xm

(C.15)

Activity can be then expressed in the dimensionless coordinate

w (X  , t ) = Aw (X (X  ), t ), A

(C.16)

 are given by and derivatives with respect to X ∂ 1 ∂ = .  ∂X Xm ∂X

(C.17)

This way Eq. (C.13) becomes

w (X w (X w (X w (X ρdry ⎞ ρdry λ (A  , t ))  , t )) A  , t ) ∂A , t) ⎞ ∂⎛ λ (A 1 ∂ ⎛ ⎜ ⎟ = ⎜D w ⎟. 3 4         ∂t ⎝ (1 + ε (Aw (X , t ))) EW ⎠ ∂X Xm (1 + ε (Aw (X , t ))) ∂X ⎝ EW Xm (1 + ε (Aw (X , t ))) ⎠

(C.18)

Note that this equation can be simplified by neglecting the swelling, i.e. putting ε = 0 , so that it becomes

  ∂ w (X w (X w (X  , t ))) = 1 ∂ ⎜⎛Dw λ (A  , t )) A  , t ) ∂Aw (X , t ) ⎟⎞. (λ (A   ∂t Xm2 ∂X ∂X ⎝ ⎠

(C.19)

Numerical calculations then demonstrate that solutions to this simplified equation are qualitatively the same as solutions to the full Eq. (C.13), e.g. desorption is still faster than absorption (see the results). On the other hand, the solutions are not completely the same, and quantitative agreement between experimental data and theoretical prediction requires using the full Eq. (C.13). In Section 3 Eq. (C.19) is further approximated by an ordinary differential equation providing further qualitative insight and giving also qualitatively the same results as the full Eq. (C.13). In summary, the swelling effect of the membrane is not essential for the observed significant difference in sorption rates and hence this disparity of time scales is not swelling-driven. It is, however, providing more precise predictions in the sense that a physical process necessarily present during sorption, swelling, is considered. See [35] for more details. Appendix D. Detailed analysis of the lumped model In this section the lumped model from this article is analyzed in detail and the conclusions presented in the article are derived. The assumptions and notation are the same as in Section 3.1: We analyze the behavior of the lumped model for different shapes of isotherm on an interval of RH (ϕini , ϕ ) for absorption and on (ϕ, ϕini ) for desorption. The lumped model enables us to study the reason why the rate difference between desorption and absorption appears in the full model. Let us compare the differential equations for sabs and sdes from (10)

d α −1 sabs (t ) = (ϕ − λabs (Δλabs sabs (t ) + λabs (ϕini ))), dt Δλabs d α −1 sdes (t ) = (ϕ − λdes (Δλdes sdes (t ) + λdes (ϕini ))). dt Δλdes

46

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

The introduction of an auxiliary function h: (0, 1) → (0, 1) with λ a−w1; abs, des denoting the inverse of functions λabs (a w ) and λdes (a w ) , respectively

habs, des (y ) ≔

λ a−w1; abs, des (y·Δλabs, des + λabs, des (ϕini )) − ϕini ϕ − ϕini

.

gives simpler version of the differential equations:

ϕ − ϕini d sabs (t ) = α (1 − habs (sabs )) ≔ fabs (sabs ), Δλabs dt

(D.1)

ϕ − ϕini d (hdes (1 − sdes )) ≔ fdes (sdes ), sdes (t ) = α dt Δλdes

(D.2)

The initial conditions are exactly the same for both differential equations (D.1) and (D.2):

sabs, des (0) = 0. Therefore, it can be deduced what causes the inequality of sorption rates. It is the difference between the functions fabs and fdes. Based on generalization of Groenwall-Bellman inequality from [36] it is possible to conclude that if

fdes (x ) ≥ fabs (x )

∀ x ∈ [0, 1],

(D.3)

then

sdes (t ) ≥ sabs (t )

∀ t ≥ 0,

which means that desorption is faster than absorption. The inequality (D.3) can be rewritten to

1 1 (hdes (1 − x )) ≥ (1 − habs (x )) Δλdes Δλabs

∀ x ∈ [0, 1].

(D.4)

The opposite inequality would imply faster absorption. All the criteria mentioned in Section 3.1 can be deduced from this inequality (D.4). Firstly, if hysteresis is neglected (hdes = habs ), the inequality becomes:

h (1 − x ) + h (x ) ≥ 1

∀ x ∈ [0, 1],

(D.5)

which gives the first criterion. The special case for convex/concave isotherm can be derived as follows. For the inverse function h

h−1 (x )

=

d2 h (x ) dx 2

= =

−1

it holds that:

λ ((ϕ − ϕini ) x + ϕini ) − λ (ϕini ) Δλ d 1 (h−1) ″ (h (x )) = − −1 dx (h−1) ′ (h (x )) ((h ) ′ (h (x ))3 λ″ (h (x )) , −C (λ′ (h (x )))3

where C is a positive constant. Using the fact that an isotherm is growing on the whole interval, it follows that

sgn

d2 h (x ) = −sgn(λ″ (h (x ))), dx 2

which implies that if the isotherm is convex, h is concave, and vice versa. From this we get:

(λ convex) ⟹ (h concave) ⟹

d h decreasing dx

⎛ dh decreasing ∧ (h (0) = 0, h (1) = 1) ⎞ ⟹ ( ∀ x ∈ (0, 1))(h (x ) ≥ x ) ⎝ dx ⎠ The fact that last implication holds is intuitive and can be shown rigorously by Rolle theorem. In conclusion we know that if an isotherm is convex then ∀ x , h (x ) ≤ x and similarly if it is concave than ∀ x , h (x ) ≥ x . Using inequality (D.5), it is proven that if the isotherm is concave, absorption is faster than desorption while if it is convex, desorption is faster (as in the case of Nafion). The second criterion in Section 3.1 directly follows from (D.4) if linearity of h is assumed. Appendix E. Fickian diffusion A Fick's model (FM) as a generalization of models from articles [5,1,4] is considered below and in detail in [20]. This model was created with the aim to decide whether such models are capable of explaining the experimental observations that are discussed in this article.9 It is also a phenomenological model with the following additional assumptions:

9 A general comment is in place although it may sound trivial: the experimentally measured parameters (as diffusion coefficient) are assessed only via a theoretical model (typically Fick's law of diffusion) and hence the obtained experimental results have to be interpreted in accordance with and within the limits of the used model. Hence if one measures a diffusion coefficient of Nafion membranes using dynamic water sorption DVS (implicitly assuming Fick's law of diffusion with certain boundary conditions) and observes that diffusion coefficient is a function of membrane thickness or humidity one needs to check if the model itself permits such dependencies. For example these dependencies are not admissible in the discussed Fick's law neither with Dirichlet nor with evaporation (also known as flux) boundary conditions and the only implication is that the model is not well suited for the carried out experiment. The carried out analysis in this paper offers a partial resolution to this problem as it turns out that the particular form of the uneasily accessible functional dependency D w (aw ) is insignificant once the diffusion is fast when compared to the sorption kinetics. This is likely to be the case for thin objects like membranes.

47

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

• Transport of water in the membrane follows Fick's law diffusion which states that flow is proportional to the gradient of concentration: ∂c w (x , t ) ∂ ⎛ ∂c (x , t ) ⎞ = D (c w ) w , ∂t ∂x ⎝ ∂x ⎠

• •

where D (c w ) is a diffusion coefficient. the membrane is separating two chambres with constant relative humidities ϕ0 and ϕ L , the following boundary conditions10 corresponding to flows of water ( jw |0, L ) are considered:

∂c (x , t ) ∂x

x=0

∂c (x , t ) jw |L = −D ∂x

x=L

jw |0 = −D

where 0, L c∞

̂0 k int

and

≔ c (a w =

̂L k int

0

0 ̂ (c∞ = k int − c (0, t )) ∀ t ∈ [0, + ∞ ) L

̂ (c∞L − c (L, t )) ∀ t ∈ [0, + ∞ ), = −k int

(E.1)

are interfacial mass transfer coefficients (having different units than k int ) and where

ϕ0, L).

• initial water concentration in the membrane is given by: c (0, x ) = c0

∀ x ∈ [0, L].

• coefficients k ̂

0 int

(E.2)

L

̂ can depend on whether the water is being absorbed or desorbed on the corresponding side of the membrane (as is and k int assumed by Ge et al. in [4]), but are otherwise constant.

Note that Satterfield et al. in [1] added an empirical term to such a model to explain some of their experimental observation. Monroe and Ge consider diffusion coefficient to be concentration dependent, we, however, assume that D (c w ) is a constant to be able to analyze the model rigorously and gain some insight in model behavior and its limitations. The analytical solution of the model can be found. It is in the form of an infinite sum: ∞

s (t ) = 1 −

∑ A (βi) e

β 2L − i2 t D ,

(E.3)

i=1

where the exact relations for A (βi ) and βi can be found in [20] and where the only the first few terms of the sum play significant role. By studying the solution it was shown under which conditions are members of the sum (E.3) negligible and conditions under which the leading term behaves exponentially [20]. a d < k int We have hypothesized that the assumption of k int is going to explain experimental observations of [1]. This is what can be concluded from the analysis [20]:

• The model predicts a collapse of sorption curves s(t) when plotted against t /L in certain ranges of parameters. This agrees well with claims and observations from [1]. • The model is able to explain the difference in rates of absorption and desorption by setting k < k appropriately. However, it is unable a int

•

d int

(because of the independence of s(t) on boundary conditions) explain that the rate ratio of absorption and desorption changes significantly depending on interval of RH as was observed by [37] and in experiments [20]. The model is unable to predict the observed non-exponential behavior of absorption. There are two reasons for that: – For a range of parameters that are experimentally relevant (e.g. [1]) the model predicts exponential behavior for both absorption and desorption. – Although the experimental assessment of parameter values vary greatly (as e.g. the reported values of diffusion coefficients [38,1]), it is possible to deduce that the model is unable to predict the magnitude of non-linearity measured in [1]. Details can be found again in [20].

These conclusions were in some manner already found in [1], but they were not derived rigorously and into such detail. Additionally, our conclusion provides more information than the ones from the mentioned article. In conclusion, a very general model based on Fick's diffusion was considered. However, from the detailed analysis we are forced to conclude that the model is unable to explain many fundamental experimental observations about dynamics of water transport in Nafion. Therefore, different models have to be searched that may have the capacity to explain the behavior and we propose such a model in the main text. For some special cases FM can fit well experimental data but not under variety of conditions.

and swelling in polymer electrolyte membranes: Diagnostic applications, J. Phys. Chem. B, 2015. [4] S. Ge, X. Li, B. Yi, I.-M. Hsing, Absorption, desorption, and transport of water in polymer electrolyte membranes for fuel cells, J. Electrochem. Soc. 152 (6) (2005) A1149–A1157. [5] C.W. Monroe, T. Romero, W. Mérida, M. Eikerling, A vaporization-exchange model for water sorption and flux in nafion, J. Membr. Sci. 324 (1) (2008) 1–6. [6] Q. Zhao, P. Majsztrik, J. Benziger, Diffusion and interfacial transport of water in Nafion, J. Phys. Chem. B 115 (12) (2011) 2717–2727.

References [1] M.B. Satterfield, J.B. Benziger, Non-fickian water vapor sorption dynamics by Nafion membranes, J. Phys. Chem. B 112 (12) (2008) 3693–3704 (PMID: 18314970). [2] M.H. Eikerling, P. Berg, Poroelectroelastic theory of water sorption and swelling in polymer electrolyte membranes, Soft Matter 7 (13) (2011) 5976–5990. [3] M. Safiollah, P.-É. A. Melchy, P. Berg, and M. H. Eikerling, Model of water sorption

10

(generalization of the considered boundary conditions in the main text of this article)

48

Journal of Membrane Science 540 (2017) 35–49

V. Klika et al.

Hydrogen Energy, 2015. [22] L. Karpenko-Jereb, P. Innerwinkler, A.-M. Kelterer, C. Sternig, C. Fink, P. Prenninger, R. Tatschl, A novel membrane transport model for polymer electrolyte fuel cell simulations, Int. J. Hydrog. Energy 39 (13) (2014) 7077–7088. [23] D.T. Hallinan Jr., M.G. De Angelis, M.G. Baschetti, G.C. Sarti, Y.A. Elabd, Nonfickian diffusion of water in nafion, Macromolecules 43 (May (25)) (2010) 4667–4678. [24] D. Goujot, O. Vitrac, Extension to nonlinear adsorption isotherms of exact analytical solutions to mass diffusion problems, Chem. Eng. Sci. 99 (2013) 2–22. [25] S. Kjelstrup, D. Bedeaux, Non-Equilibrium Thermodynamics of Heterogeneous Systems. Series on Advancesin Statistical Mechanics, World Scientific, 2008. [26] J. Newman, K. Thomas-Alyea, Electrochemical Systems. Electrochemical Society series, Wiley, 2004. [27] A. El Afif, M. Grmela, Non-fickian mass transport in polymers, J. Rheol. (1978Present) 46 (3) (2002) 591–628. [28] Q. Liu, X. Wang, D.D. Kee, Mass transport through swelling membranes, Int. J. Eng. Sci. 43 (19–20) (2005) 1464–1470. [29] K.A. Mauritz, R.B. Moore, State of understanding of Nafion, Chem. Rev. 104 (10) (2004) 4535–4586. [30] Nafion® membranes, extrusion-cast. http://www2.dupont.com/FuelCells/en_US/ assets/downloads/dfc101.pdf. (accessed 27 January 2015). [31] R. Bearman, On molecular basis of some current theories of diffusion, J. Phys. Chem. 65 (11) (1961) 1961–1968. [32] J. Vrentas, J. Duda, Diffusion in polymer - Solvent Systems.1. Re-examination of free-volume theory, J. Polym. Sci. Part B: Polym. Phys. 15 (3) (1977) 403–416. [33] J. Vrentas, J. Duda, Diffusion in Polymer-Solvent Systems.2. Predictive Theory for Dependence of Diffusion-Coefficients on Temperature, Concentration, and Molecular-Weight, J. Polym. Sci. Part B: Polym. Phys. 15 (3) (1977) 417–439. [34] M. Pavelka, F. Maršík, V. Klika, Consistent theory of mixtures on different levels of description, Int J. Eng. Sci. 78 (0) (2014) 192–217. [35] M. Pavelka, Thermodynamic analysis of processes in Hydrogen fuel cells. dissertation, Charles University in Prague, 2015. [36] I. Bihari, A generalization of a lemma of bellman and its application to uniqueness problems of differential equations, Acta Math. Hung. 7 (1) (1956) 81–94. [37] P. Krtil, A. Trojánek, Z. Samec, Kinetics of water sorption in Nafion thin films-quartz crystal microbalance study, J. Phys. Chem. B 105 (33) (2001) 7979–7983. [38] D. Rivin, C.E. Kendrick, P.W. Gibson, N.S. Schneider, Solubility and transport behavior of water and alcohols in Nafion, Polymer 42 (2) (2001) 623–635.

[7] V. Klika, A guide through available mixture theories for applications, Crit. Rev. Solid State Mater. Sci. 39 (2) (2014) 154–174. [8] J.B. Benziger, M.J. Cheah, V. Klika, M. Pavelka, Interfacial constraints on water and proton transport across nafion membranes, J. Polym. Sci. Part B: Polym. Phys. 53 (22) (2015) 1580–1589. [9] S.R. de Groot, P. Mazur, Non-equilibrium Thermodynamics, Dover Publications, New York, 1984. [10] M. Pavelka, V. Klika, P. Vágner, F. Maršík, Generalization of exergy analysis, Appl. Energy 137 (0) (2015) 158–172. [11] G. Gebel, S. Lyonnard, H. Mendil-Jakani, A. Morin, The kinetics of water sorption in Nafion membranes: a small-angle neutron scattering study, J. Phys.:Condens. Matter 23 (234107) (2011). [12] T. Springer, T. Zawodzinski, and S. Gottesfeld, Polymer Electrolyte Fuel-Cell Model, J Electrochem Soc, vol. 138, pp. 2334–2342, AUG 1991. [13] T. Thampan, S. Malhotra, H. Tang, R. Datta, Modeling of conductive transport in proton-exchange membranes for fuel cells, J. Electrochem. Soc. 147 (9) (2000) 3242–3250. [14] S. Jeck, P. Scharfer, M. Kind, Absence of Schroeder's paradox: experimental evidence for water-swollen Nafion®membranes, J. Membr. Sci. 373 (1–2) (2011) 74–79. [15] A.Z. Weber, J. Newman, Transport in polymer-electrolyte membranes i. physical model, J. Electrochem. Soc. 150 (7) (2003) A1008–A1015. [16] L.M. Onishi, J.M. Prausnitz, J. Newman, Water-Nafion equilibria. absence of Schroeder's paradox, J. Phys. Chem. B 111 (34) (2007) 10166–10173 (PMID: 17685645). [17] J. Benziger, A. Bocarsly, M.J. Cheah, P. Majsztrik, B. Satterfield, Q. Zhao, Mechanical and transport properties of nafion: effects of temperature and water activity, Struct. Bond 141 (2011) 85–113. [18] N. Foley, K. Thomas, P. Forshaw, D. Stanton, P. Norman, Kinetics of water vapor adsorption on activated carbon, Langmuir 13 (7) (1997) 2083–2089. [19] J. Crank, M. Henry, Diffusion in media with variable properties. part i. the effect of a variable diffusion coefficient on the rates of absorption and desorption, T Faraday Soc. 45 (1949) 636–650. [20] J. Kubant, Modelling of water transport in a proton exchange membrane together with its experimental verification, Master’s thesis, Faculty of Nuclear Sciences and Scientific Engineering, Czech Technical University in Prague, Czech Republic, 2015. [21] V. Liso, S. S. Araya, A. C. Olesen, M. P. Nielsen, and S. K. Kær, Modeling and experimental validation of water mass balance in a pem fuel cell stack, International J.

49