Generalized model for the diffusion of solvents in glassy polymers: From Fickian to super Case II

TL;DR

This paper introduces a unified kinetic model for solvent diffusion in glassy polymers that can describe Fickian, anomalous, Case II, and super Case II diffusion behaviors based on two key parameters.

Contribution

The authors develop a single generalized model capable of capturing all known diffusion types in glassy polymers using just two parameters, simplifying previous specialized models.

Findings

The model accurately reproduces various diffusion behaviors.

The dimensionless parameter γ determines the diffusion type.

Time-delayed swelling has minimal impact on diffusion behavior.

Abstract

The diffusion of small solvent molecules in glassy polymers may take on a variety of different forms. Fickian, anomalous, Case II, and super Case II diffusion have all been observed, and theoretical models exist that describe each specific type of behavior. Here we present a single generalized kinetic model capable of yielding all these different types of diffusion on the basis of just two parameters. The principal determinant of the type of diffusion is observed to be a dimensionless parameter, $\gamma$, that describes the influence of solvent-induced swelling in lowering the potential barriers separating available solvent sites. A second parameter, $\eta$, which characterizes the effect of the solvent in reducing the potential energy of a solvent molecule when at rest at an available site, only influences the type of diffusion to a lesser extent. The theoretical analysis does not include any effects that are explicitly non-local in time, an example of which is the inclusion of polymer viscosity in the Thomas-Windle model; it thus represents a variant of Fick's second law utilizing a concentration-dependent diffusivity when $\eta$ is small. To check the significance of time-delayed swelling, a simulation was performed of a modified model that contained a history-dependent term. The results were found to be very similar to those in the time-local model.

Tables

Table 1. Table 1: The fitting parameters found by fitting the function f ​ ( t ) = k ​ t n 𝑓 𝑡 𝑘 superscript 𝑡 𝑛 f(t)=kt^{n} to the mass-time curves in Figs.  13 (d) (second row), 14 (d) (third row), and 15 (d) (fourth row).

$\eta$	0	2.5	5.0
$n$ for $\gamma =0$	0.502	0.507	0.531
$n$ for $\gamma =5.0$	0.599	0.610	0.610
$n$ for $\gamma =9.6$	1.024	0.974	1.000

Equations

$$\frac{\partial \phi }{\partial t}=D{\nabla }^2\phi ,$$

$$\frac{\partial \phi }{\partial t}=\nabla \cdot \left[D(\phi )\nabla \phi \right]$$

$$\phi (x,t)={\phi }_0{\left(e^{q(x-vt)}+1\right)}^{-1},$$

$$V_2\left(\mathbf{l}+\frac{\mathbf{a}}{2}\right)-\frac{1}{2}\left[V_1(\mathbf{l})+V_1(\mathbf{l}+\mathbf{a})\right]\equiv V_0\left(\mathbf{l}+\frac{\mathbf{a}}{2}\right)$$

$$V_1(\mathbf{l}+\mathbf{a})-V_1(\mathbf{l})\equiv W\left(\mathbf{l}+\frac{\mathbf{a}}{2}\right).$$

$$\begin{split}\mbox{\boldmath$J$}_{{\boldsymbol{l}},{\boldsymbol{l}}+{\boldsymbol{a}}}=&\mbox{\boldmath$a$}\nu e^{-\beta V_{0}}[\phi_{l}(1-\phi_{l+a})e^{-\frac{1}{2}\beta W}\\ &-\phi_{l+a}(1-\phi_{l})e^{\frac{1}{2}\beta W}],\end{split}$$

$$W=-K(\phi)\mbox{\boldmath$a$}\cdot\nabla\phi,$$

$$K(\varphi )=-\frac{d{V}_1}{d\varphi },$$

$$V_0(\varphi )=V_b-(V_b-V_a){\varphi }^4,$$

$$\gamma \equiv \beta (V_b-V_a).$$

$$\eta \equiv \frac{1}{2}\beta K=-\frac{1}{2}\beta \frac{d{V}_1}{d\varphi }.$$

$$\begin{split}\mbox{\boldmath$J$}_{{\boldsymbol{l}},{\boldsymbol{l}}+{\boldsymbol{a}}}=&\mbox{\boldmath$a$}\nu e^{-\beta V_{b}}e^{\gamma\phi_{l+\frac{a}{2}}^{4}}[\phi_{l}(1-\phi_{l+a})e^{\eta(\phi_{l+a}-\phi_{l})}\\ &-\phi_{l+a}(1-\phi_{l})e^{-\eta(\phi_{l+a}-\phi_{l})}]\\ =&\mbox{\boldmath$a$}\nu e^{-\beta V_{b}}e^{\gamma\phi_{l+\frac{a}{2}}^{4}}\{(\phi_{l}-\phi_{l+a})\cosh[\eta(\phi_{l+a}-\phi_{l})]\\ &+(\phi_{l}+\phi_{l+a}-2\phi_{l}\phi_{l+a})\sinh[\eta(\phi_{l+a}-\phi_{l})]\}.\end{split}$$

$$\frac{\partial\phi}{\partial t}=-\nabla\cdot\mbox{\boldmath$J$},$$

$${\varphi }_l(t+\Delta t)-{\varphi }_l(t)=\frac{\Delta t}{a}\left(J_{l-a,l}-J_{l,l+a}\right).$$

$$\varphi (x,t){|}_{t=0}=0,(x>0)$$

$$\varphi (x,t){|}_{x=0}=1,$$

$$\varphi (x,t){|}_{x\to \infty }=0.$$

$$\varphi (x,t){|}_{x=L}=1,$$

$$\delta {\varphi }^{(1)}={\varphi }_{l+a}-{\varphi }_l$$

$$\begin{array}{rl}\delta {\varphi }^{(5)}=& \frac{1}{16}[({\varphi }_{l+3a}-{\varphi }_{l+2a})+4({\varphi }_{l+2a}-{\varphi }_{l+a})+6({\varphi }_{l+a}\\ & -{\varphi }_l)+4({\varphi }_l-{\varphi }_{l-a})+({\varphi }_{l-a}-{\varphi }_{l-2a})]\\ =& \frac{1}{16}[{\varphi }_{l+3a}+3{\varphi }_{l+2a}+2{\varphi }_{l+a}-2{\varphi }_l-3{\varphi }_{l-a}-{\varphi }_{l-2a}]\end{array}$$

$$\begin{array}{rl}\delta {\varphi }^{(9)}=& \frac{1}{256}[{\varphi }_{l+5a}+7{\varphi }_{l+4a}+20{\varphi }_{l+3a}+28{\varphi }_{l+2a}+14{\varphi }_{l+a}\\ & -14{\varphi }_l-28{\varphi }_{l-a}-20{\varphi }_{l-2a}-7{\varphi }_{l-3a}-{\varphi }_{l-4a}]\end{array}$$

$$\begin{array}{rl}\delta {\varphi }^{(13)}=& \frac{1}{4096}[{\varphi }_{l+7a}+11{\varphi }_{l+6a}+54{\varphi }_{l+5a}+154{\varphi }_{l+4a}\\ & +275{\varphi }_{l+3a}+297{\varphi }_{l+2a}+132{\varphi }_{l+a}-132{\varphi }_l\\ & -297{\varphi }_{l-a}-275{\varphi }_{l-2a}-154{\varphi }_{l-3a}-54{\varphi }_{l-4a}\\ & -11{\varphi }_{l-5a}-{\varphi }_{l-6a}].\end{array}$$

$${\varphi }_{l+\frac{a}{2}}^{(1)}=\frac{1}{2}({\varphi }_l+{\varphi }_{l+a})$$

$$\begin{array}{rl}{\varphi }_{l+\frac{a}{2}}^{(5)}=& \frac{1}{2}\times \frac{1}{16}[({\varphi }_{l+3a}+4{\varphi }_{l+2a}+6{\varphi }_{l+a}+4{\varphi }_l+{\varphi }_{l-a})\\ & +({\varphi }_{l+2a}+4{\varphi }_{l+a}+6{\varphi }_l+4{\varphi }_{l-a}+{\varphi }_{l-2a})]\\ =& \frac{1}{32}({\varphi }_{l+3a}+5{\varphi }_{l+2a}+10{\varphi }_{l+a}+10{\varphi }_l+5{\varphi }_{l-a}\\ & +{\varphi }_{l-2a})\end{array}$$

$$\begin{array}{rl}{\varphi }_{l+\frac{a}{2}}^{(9)}=& \frac{1}{512}({\varphi }_{l+5a}+9{\varphi }_{l+4a}+36{\varphi }_{l+3a}+84{\varphi }_{l+2a}+126{\varphi }_{l+a}\\ & +126{\varphi }_l+84{\varphi }_{l-a}+36{\varphi }_{l-2a}+9{\varphi }_{l-3a}+{\varphi }_{l-4a})\end{array}$$

$$\begin{array}{rl}{\varphi }_{l+\frac{a}{2}}^{(13)}=& \frac{1}{8192}({\varphi }_{l+7a}+13{\varphi }_{l+6a}+78{\varphi }_{l+5a}+286{\varphi }_{l+4a}\\ & +715{\varphi }_{l+3a}+1287{\varphi }_{l+2a}+1716{\varphi }_{l+a}+1716{\varphi }_l\\ & +1287{\varphi }_{l-a}+715{\varphi }_{l-2a}+286{\varphi }_{l-3a}+78{\varphi }_{l-4a}\\ & +13{\varphi }_{l-5a}+{\varphi }_{l-6a}).\end{array}$$

$$\begin{array}{r}{\varphi }_{l+\frac{a}{2}}=\frac{1}{2}({\varphi }_{l+a}+{\varphi }_l)\\ \nabla \varphi \cdot \mathbf{a}={\varphi }_{l+a}-{\varphi }_l\end{array}$$

$$\begin{split}&\mbox{\boldmath$J$}_{{\boldsymbol{l}},{\boldsymbol{l}}+{\boldsymbol{a}}}=\mbox{\boldmath$a$}\nu e^{-\beta V_{b}}e^{\gamma\phi_{l+\frac{a}{2}}^{4}}\Bigl{\{}-{\boldsymbol{\nabla}}\phi\cdot{\boldsymbol{a}}\cosh(\eta{\boldsymbol{\nabla}}\phi\cdot{\boldsymbol{a}})\\ &+\left[-2\phi_{l+\frac{a}{2}}^{2}+2\phi_{l+\frac{a}{2}}+\frac{1}{2}({\boldsymbol{\nabla}}\phi\cdot{\boldsymbol{a}})^{2}\right]\sinh(\eta{\boldsymbol{\nabla}}\phi\cdot{\boldsymbol{a}})\Bigr{\}}.\end{split}$$

$$\begin{split}D=&a^{2}\nu e^{-\beta V_{b}}e^{\gamma\phi^{4}}\Bigl{\{}\cosh(\eta{\boldsymbol{\nabla}}\phi\cdot{\boldsymbol{a}})\\ &+\left[\frac{2\phi^{2}-2\phi}{{\boldsymbol{\nabla}}\phi\cdot{\boldsymbol{a}}}-\frac{1}{2}{\boldsymbol{\nabla}}\phi\cdot{\boldsymbol{a}}\right]\sinh(\eta{\boldsymbol{\nabla}}\phi\cdot{\boldsymbol{a}})\Bigr{\}},\end{split}$$

$$D=a^2\nu e^{-\beta V_b}{e}^{\gamma {\varphi }^4}\left[1+2\eta \varphi (\varphi -1)\right],$$

Full text

I Introduction

In the study of the diffusion of small molecules into polymeric materials, several different types of behavior have been observed. They are classified according to the equation: $M_{t}=kt^{n}$, where $M_{t}$ is the amount of solvent absorbed per unit area of polymer at an elapsed time $t$ after the two have been placed in contact, and $k$ is a constant. When the exponent $n=1/2$, it is the familiar Fickian diffusion; when $n=1$, it is known as Case II diffusion, and when $1/2<n<1$ it is termed anomalous diffusion. masaro Jacques et al. jac found that there is an acceleration of diffusion when the two fronts originating from the two opposite sides of a sample slab meet during Case II diffusion, and called it super Case II diffusion.

Predicting the detailed nature of the diffusion of the penetrant has turned out to be a challenging problem, and most of the proposed theoretical models apply only to a specific subinterval of the diffusion realm. For example, Thomas and Windle proposed a widely discussed model for Case II diffusion. thom Qian and Taylor made an effort to explain the crossover from Case II to Fickian behavior using two separate equations. qian Hansen has also tried to incorporate different kinds of diffusion behavior by modifying the boundary conditions of the diffusion equation and using diffusion coefficients that are exponentially dependent on the concentration. hansen Wilmers and Bargmann also used an exponentially concentration-dependent diffusion coefficient while incorporating a dual-phase-lag concept to Fick’s first law and derived a promising description of Case II behavior. wilmers

In this paper, we present a single model that exhibits all the observed types of diffusion. We begin by noting that solvent-induced swelling can affect the diffusivity in two distinct ways. One consequence of the swelling is a decrease in the potential energy of solvent molecules in the interstices of the polymer matrix. A second consequence is that the height of the potential-energy barriers that must be surmounted is reduced. We explore the respective effects of these phenomena in detail, using first a theoretical analysis and then a numerical simulation. We find the latter effect to be the principal determinant of the diffusive behavior. With an appropriate parameter set, all the different types of diffusion mentioned above can be reproduced.

II Theoretical Analysis

II.1 The background

In the simplest possible model of diffusion, the diffusing particles execute a random walk within the interstices of a fixed network of obstructions. On a macroscopic scale, in which a local concentration $\varphi({\bf r})$ of diffusant can be defined, one finds that $\varphi$ obeys Fick’s law, crank

[TABLE]

and the process is known as Fickian diffusion. The development with time of the function $\varphi({\bf r},t)$ is entirely determined by the initial conditions and the constant diffusivity $D$.

The process by which an organic solvent diffuses into a real polymer is much more complex than the simple Fickian picture. The most important difference is that the presence of the solvent modifies the diffusion process. This can occur through a variety of mechanisms, the simplest of which is that entry of solvent may cause the polymer to swell, thus widening the channels through which the solvent molecules must pass. These processes have the effect of making the diffusivity $D$ a function of solvent concentration. The non-linear diffusion equation that results is

[TABLE]

when other complications are ignored.

Although Eq. (2) is an oversimplification of a real polymer diffusion problem, it can provide good insight into the type of non-linear diffusion that is possible. An example would be the exponential dependence of $D$ on $\varphi$ first proposed by Prager, prager which has been the subject of a recent study by Marais et al. marais In another special case, $D=D_{0}/(\varphi_{0}-\varphi)$ with $D_{0}$ and $\varphi_{0}$ constants. We then have a diffusivity that is inversely proportional to the polymer density, and Eq. (2) has a solution for diffusion in one dimension of the form

[TABLE]

with $v$ a velocity and $q=\varphi_{0}v/D_{0}$. This represents a solitary wave in which the solvent advances into the polymer as a front moving with constant speed $v$, as illustrated in Fig. 1. It thus represents an example of what has been termed Case II diffusion, thomas1 ; thom ; alfrey ; rossi ; thomas3 ; lasky1 ; lasky2 ; hui1 ; hui2 ; mills in which a glassy polymer is transformed to a more porous, rubbery state by the intrusion of solvent. The model of idealized Case II diffusion, in which the speed of the advancing front is constant, can only be an approximation to the real situation, in which there will always be some slowing of the rate at which the front advances. The model does, however, serve an important function in indicating that a study of the form of the concentration-dependent diffusivity $D(\varphi)$ can be a first step in developing a realistic macroscopic theory with good predictive power. We note, however, that the speed $v$ of the diffusion front is not determined by Eq. (2) and Eq. (3), which contain no parameter having the dimensions of a velocity.

II.2 The model

We consider a model of a random walker, for which the direction of motion of a solvent molecule from one site to an adjacent site is completely independent of the direction of the previous motion by which the molecule arrived at that site. This is used because, in molecular dynamics simulations, one sees a solvent molecule oscillating many times within a cavity before making a transition to a new site. The momentum with which it arrives at a site is thus lost, and the direction of the next jump to another site is mostly uncorrelated with the direction of the previous one. There may remain a small negative correlation due to the swelling of a previously occupied site.

This model is illustrated in Fig. 2. Within the polymer network there are sites at positions ${\boldsymbol{l}}$ at which solvent molecules may reside. A particular configuration will have some sites occupied (shown as filled circles) and some empty. We write the probability of occupancy of site ${\boldsymbol{l}}$ as $\phi_{l}$. This will be related to the local concentration $\varphi$ of solvent through the equation $\varphi=\rho_{0}\phi_{l}$, with $\rho_{0}$ the concentration of sites. We then examine the current $J$ of solvent flowing between sites ${\boldsymbol{l}}$ and ${\boldsymbol{l}}+{\boldsymbol{a}}$, with ${\boldsymbol{a}}$ the displacement between a cavity and one of its nearest neighbors.

For a solvent molecule to make a transition from ${\boldsymbol{l}}$ to ${\boldsymbol{l}}+{\boldsymbol{a}}$ we must have site ${\boldsymbol{l}}$ occupied and site ${\boldsymbol{l}}+{\boldsymbol{a}}$ empty. There is thus a factor of $\phi_{l}(1-\phi_{l+a})$ in the expression for the current. There will also be an attempt frequency $\nu$ and an Arrhenius factor, $e^{-\beta V}$, for the probability of success of the transition, with $\beta=1/k_{\rm B}T$, where $k_{\rm B}$ is Boltzmann’s constant and $T$ is the absolute temperature. The quantity $V$ is an energy barrier height, which will differ for a transition in the forward direction from that for a reverse transition. The energy landscape is depicted in Fig. 3(a). The potential energy at the local sites varies with solvent concentration, and hence with position, as $V_{1}$, while the potential at the summit of the barriers varies as $V_{2}$.

We can thus write the barrier $V$ for forward and reverse transitions as $V_{2}({\boldsymbol{l}}+\frac{{\boldsymbol{a}}}{2})-V_{1}({\boldsymbol{l}})$ and $V_{2}({\boldsymbol{l}}+\frac{{\boldsymbol{a}}}{2})-V_{1}({\boldsymbol{l}}+{\boldsymbol{a}})$, respectively. We define

[TABLE]

and

[TABLE]

When we sum the contributions from forward and reverse transitions we find the current ${\boldsymbol{J}}$ of solvent flowing between sites ${\boldsymbol{l}}$ and ${\boldsymbol{l}}+{\boldsymbol{a}}$ as

[TABLE]

where $V_{0}$ and $W$ are evaluated at position ${\boldsymbol{l}}+\frac{{\boldsymbol{a}}}{2}$. The small difference $W$ in barrier height can be written as

[TABLE]

where

[TABLE]

since it can reasonably be assumed to depend linearly on the concentration gradient in all but the most extreme circumstances. In our present calculation, we assume that $K$ is a non-negative constant. For concreteness, we also take a simple model for the form of $V_{0}(\phi)$. We expect the presence of solvent to have the most marked effect on lowering the potential barrier $V_{0}$ when $\phi$ is close to unity; a small concentration of solvent will only have a marginal effect on swelling. It is only when $\phi$ is largest that the polymer matrix will be loosened sufficiently to permit a uniform expansion. This indicates that $V_{0}(\phi)$ is a monotonically decreasing concave function on $[0,1]$, such as $V_{0}(\phi)=A-B\phi^{r}$ ($A$, $B$ and $r$ are positive numbers, and $r\geq 2$). Taking an example among the variety of possible forms of $V_{0}$, the relation that we adopt to model this behavior is

[TABLE]

where $V_{b}$ is the potential barrier in the absence of solvent ($\phi=0$), while $V_{a}$ is the lowered barrier at saturation ($\phi=1$). The form of $V_{1}(\phi)$ and $V_{2}(\phi)$ are then as shown in Fig. 3(b).

At this point we introduce two dimensionless parameters, $\gamma$ and $\eta$, to characterize the dependence of the potential-energy landscape on the solvent occupancy $\phi$. The first parameter is defined as

[TABLE]

As we see from Fig. 3(a) and Eqs. (4) and (9), this measures the strength of the dependence of the barrier height on $\phi$. The second parameter is

[TABLE]

This is a measure of the dependence of the potential energy at a local minimum on $\phi$, as is also seen from Fig. 3(a).

We can now rewrite Eq. (6), which expresses the current \mbox{\boldmathJ}_{{\boldsymbol{l}},{\boldsymbol{l}}+{\boldsymbol{a}}} in terms of the function $V_{0}(\phi)$ and $W(\phi)$, into a form in which the current depends on the dimensionless parameters $\gamma$ and $\eta$. We interpret the value of $\phi$ appearing in $V_{0}(\phi)$ and $W(\phi)$ as $(\phi_{l}+\phi_{l+a})/2\simeq\phi_{l+\frac{a}{2}}$ and find

[TABLE]

Conservation of matter tells us that:

[TABLE]

of which the discretized form is (considering only diffusion in one dimension)

[TABLE]

II.3 Results

The combination of Eq. (14) with Eq. (12) was solved numerically. The initial and boundary conditions are:

[TABLE]

[TABLE]

[TABLE]

To minimize the calculation time while keeping the solution convergent, the prefactor $\frac{\Delta t}{a}$ was set to $2.5\times 10^{-6}$.

Results for various values of $\gamma$ and $\eta$ were obtained as time-varying concentration profiles and as plots of mass intake as a function of time. In one limit, when $\gamma=0$ and $\eta=0$, Eq. (14) becomes Eq. (1), which is Fick’s second law, and the mass intake varies as $t^{1/2}$. In another limit, when $\gamma\gg 1$, the factor $e^{-\beta V_{0}}$ has a shape similar to $D=D_{0}/(\varphi_{0}-\varphi)$, and we expect Eq. (14) will give us a solitary wave solution similar to Eq. (3), for which the mass intake varies as $t$. This suggests that the mass-time curves calculated for different values of $\gamma$ and $\eta$ might be fitted to a function of the form $f(t)=kt^{n}$, with $k$ and $n$ the fitting parameters, to determine the nature of the diffusion mechanism as discussed below.

The results for $\gamma=0$ are shown in Fig. 4 as concentration profiles at several times and for three different values of $\eta$ in Figs. 4(a), 4(b), and 4(c), while the time variation of the mass intake is shown in Fig. 4(d) for the same $\eta$ values. The fitting parameter $n$ varies with the elapsed time $t$; it is shown in Fig. 5. For $\eta=0,0.5,1.0,1.5,2.0$, and $2.5$ the $n$-$t$ curves are well fitted by a function of the form $n=n_{\infty}+bt^{-c}$, as can be seen in Fig. 5(a). This allows an accurate estimate of $n_{\infty}$, the long-time limit of $n$, as 0.49992, 0.49998, 0.50002, 0.50006, 0.50006 and 0.49989 respectively for these values of $\eta$. The $n$-$t$ relations shown in Fig. 5(b) for $\eta=3.0,3.5,4.0,4.5$ could not be fitted well, but we can estimate values of $n_{\infty}$ as 0.49937, 0.49933, 0.50044 and 0.50279 by extending the diffusion time. The $n_{\infty}$ for $\eta=5.0$, which is outside the range shown in Fig. 5, is calculated to be 0.49918. Generally speaking, the mass-time curves follow asymptotically a $t^{1/2}$ relation, characteristic of Fickian diffusion. The influence of $\eta$ on $n$ is significant only at an early stage in the diffusion process, and is negligible at long diffusion times.

When $\gamma=5.0$, the results are as shown in Fig. 6. We again fit the mass-time curves with the function $f(t)=kt^{n}$, giving the results shown in Fig. 7. We see that, although $n$ is decreasing with the diffusion time, it remains greater than 0.5, particularly for larger values of $\eta$. According to the classification, this is the regime of anomalous diffusion.

Finally, the results for $\gamma=9.6$ are as shown in Fig. 8. We see that a sharp diffusion front develops in this case, which is a characteristic of Case II diffusion. Moreover, as $\eta$ becomes larger, a small increase in $\phi$ is seen to occur ahead of the sharp front. This phenomenon has been described as a “Fickian precursor”. qian The mass intake is proportional to the time for all $\eta$ at small times. For $\eta=0$, some curvature is visible at long times, but this is not observable in our results for $\eta=2.5$ or $\eta=5.0$, partly as a result of the very slow diffusion in the latter case. The increasing curvature of the mass-time curve as $t$ increases can be interpreted as a decrease in the parameter $n$ for the case where $\eta=0$. This is illustrated in Fig. 9, where $n$ is seen to drop from a value close to unity all the way down to 0.6 after 100 million time steps. Because we are assuming a sample with infinite thickness and all the reported experiments used very thin samples, this leads us to conclude that Case II diffusion is only a “transient” phenomenon that happens at an early stage of diffusion under the above conditions. This has also been pointed out by Windle in a discussion of the limits of Case II sorption. windle According to Fig. 9, Case II diffusion happens in the first 3 million time steps for $\eta=0$. The duration can be extended by applying a larger $\eta$. For example, when $\eta=2.5$, the index $n$ remains as high as 0.992 after 15 million time steps; while for $\eta=5$, the index $n$ remains as high as 0.996 after 300 million time steps. The duration can also be extended by increasing the exponent in the $\phi^{4}$ term in Eq. (9) to a higher value.

Super Case II diffusion happens when the two opposite sides of a sample are in contact with solvent. To examine this phenomenon we replace the boundary condition Eq. (17) by the requirement:

[TABLE]

where $L$ is the thickness of the sample. In our model, super Case II diffusion originates from the small energy difference $W$ in barrier height expressed by Eq. (7). When the two fronts originating from opposite sides of a simple slab approach each other, the swelling due to solvent molecules in one front propagates elastically into the vicinity of the other front. The parameter $\eta$ must necessarily be positive in order for this phenomenon to occur.

Because the presence of a solvent molecule at site ${\boldsymbol{l}}$ will cause swelling that will be significant at sites more distant from ${\boldsymbol{l}}$ than the nearest-neighbor sites, it is necessary for the $W$ in our model to reflect this non-local influence. The barrier-height difference $W$ was represented in Eq. (7) as proportional to ${\boldsymbol{a}}\cdot\nabla\phi$, which was approximated in the previous calculation as $\phi_{l+a}-\phi_{l}$. For super Case II diffusion, this simple difference was replaced by a weighted average over a thickness of 1, 5, 9, and 13 layers, respectively. The weighting followed a binomial distribution, for which

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The corresponding $\phi_{l+\frac{a}{2}}$ is also calculated using the same binomial distribution:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The resulting concentration profiles are shown in Fig. 10. We see that with the increase of chosen thickness, the “precursor” becomes smaller, and the diffusion speed becomes larger.

The corresponding curves of mass intake versus time are shown in Fig. 11 for layers of thickness 5, 9, and 13 (the single-layer case is omitted, as no super Case II behavior is observable). In this figure, the trends of linear sorption are shown as dashed lines. For the three larger thicknesses considered, the acceleration in the rate of mass intake remains constantly observable in this model. This suggests that the particular phenomenon of super Case II diffusion requires the inclusion of longer-range elasticity effects for an accurate description.

Our overall conclusion in this section is that it is possible to reproduce all the observed types of diffusion by extending the well known diffusion equation. The time variation of the mass intake, however, deviates from a simple power law when the type of diffusion is anomalous or Case II.

III A stochastic model

III.1 Simulation details

In the preceding section, diffusion was modeled as an instantaneous process in which the local transport of solvent was determined solely by the concentration gradient at that instant of time. It is, however, also possible that effects that are nonlocal in time could occur. The barrier to hopping between two sites could reflect a previous history of occupation if local solvent-induced swelling persists for a brief time even in the absence of solvent molecules.

To study whether such persistence effects would change the predictions made in the preceding sections, simulations were performed on a model in which the transition probabilities included a memory factor. The polymer network is represented by a three-dimensional grid as shown in Fig. 12. The cubes in this grid denote cavities inside the polymer bulk, the walls between the cubes being polymer matrix. Transport is in the $x$-direction, and the number of cavities in each layer in the $y$-$z$ plane is set to 10$\times$10. Periodic boundary conditions are used in the $y$ and $z$ directions. During the simulation, the $x$ = 0 surface of this grid is considered to be in contact with the solvent, and accordingly all cavities on this surface are filled at all times. All other sites may be either empty or singly occupied.

The diffusion occurs through a series of events in which a molecule in a cavity can transition only to one of its six nearest neighbors. To simulate the effect of previous history on the local transport, we characterize each wall as being either intact or ruptured, according to whether there has been a previous transit through that wall.

The probability of a transition occurring in which a solvent molecule moves between two sites will depend, as in the previous section, on the height of the barrier; this height, however, will now have a value that is dependent on past history rather than on the local concentration of solvent. We retain the language used in the previous section by denoting the barrier height of an intact wall as $V_{b}$, and that of a penetrated wall as $V_{a}$. The ratio of probabilities for transitions will then be $e^{-\beta V_{b}}/e^{-\beta V_{a}}=e^{-\gamma}$, where the dimensionless quantity $\gamma$ carries the same connotation as before, i.e., a measure of the barrier height reduction due to the presence of solvent. Now, however, the solvent may have been present in the past, but be currently absent.

We must also take into account the tendency for the potential energy at a site to be lowered by the presence of solvent. This effect was characterized in the preceding section by the parameter $\eta$, which was proportional to $-dV_{1}/d\phi$. To introduce this effect into the simulation, we define $\phi_{l}\equiv\phi(x=l)$ as the average occupancy of the 100 sites within a slice of the grid in the $y$-$z$ plane. The probabilities of hopping in the positive-$x$ and negative-$x$ directions are then multiplied by a factor $e^{\eta\delta\phi_{+}}$ and $e^{-\eta\delta\phi_{-}}$, respectively, where $\delta\phi_{+}=\phi_{l+a}-\phi_{l}$ and $\delta\phi_{-}=\phi_{l}-\phi_{l-a}$.

The combined prescription for the transition probabilities is then as follows: In the forward direction the probability to make a transition is proportional to $e^{\eta\delta\phi}$ if the wall is unpenetrated, and $e^{\gamma+\eta\delta\phi}$ if it has been previously traversed; in the backward direction the probabilities are proportional to $e^{-\eta\delta\phi}$ for an unpenetrated wall and $e^{\gamma-\eta\delta\phi}$ for a penetrated one. Transitions at constant $x$ have no dependence on $\eta$, but the probability of occurrence contains the factor of $e^{\gamma}$ for travel through a ruptured wall. For each set of $\gamma$ and $\eta$, we perform simulations on 10 parallel systems, and then calculate ensemble averages to reduce the statistical error.

III.2 Results for the stochastic model

As in the preceding section, we find that the type of diffusion is mainly determined by $\gamma$ and mildly influenced by $\eta$. When $\gamma=0$, the diffusion is mainly Fickian; when $\gamma\gg 1$, for example 9.6 in the following calculation, it is mainly Case II diffusion; and when $\gamma$ takes an intermediate value, such as 5.0, it is mainly anomalous diffusion. The major role of the parameter $\eta$ is its influence in producing a sharp diffusion front and in the formation of a precursor.

III.2.1 Fickian diffusion

When $\gamma=0$, the transition probability is not history dependent, and the probability of passing through a previously penetrated wall is the same as for an unpenetrated wall.

The calculated concentration profiles at different times and different values of $\eta$ (0, 2.5, 5.0) are shown in Figs. 13(a)-(c). In these plots we see that, with increasing $\eta$, a sharp diffusion front begins to develop, as seen earlier in Fig. 4. The plots when $\eta=2.5$ in both Fig. 4 and Fig. 13 resemble Fig. 8 of Thomas and Windle’s paper, thom where, in addition to the sharp front, there is also a pronounced precursor ahead of it, and a steep concentration gradient following it.

Figure 13(d) shows the mass intake of solvent versus time at different values of $\eta$. These mass-time curves are then fitted with the function $f(t)=kt^{n}$, with $k$ and $n$ the fitting parameters. We find that, unlike the time local theory in Section II, the variation of $n$ with respect to time $t$ is greatly reduced after taking the non-local effect in time into account. The results for the exponent $n$ are shown in the second row of Table 1. The mass-time curves generally follow a $t^{1/2}$ relation, which is a fingerprint of Fickian diffusion. When $\eta$ is larger than 5.0, the uncertainty in $n$ increases rapidly as a consequence of the slow diffusion speed and relatively short simulation time.

III.2.2 Anomalous diffusion

When $\gamma=5.0$, the computed concentration profiles at different times and for different values of $\eta$ are shown in Figs. 14(a)-(c). In contrast with the Fickian diffusion discussed above, we notice the appearance of a sharp front at $\eta=0$. It is also noticed that, with the increase of $\eta$, the concentration gradient behind the front becomes less steep, while the precursor ahead of the front penetrates more deeply into the polymer matrix, although with reduced concentration range.

The mass-time curves are shown in Fig. 14(d); their fits with the function $f(t)=kt^{n}$ are shown in the third row of Table 1. We see that the exponents remain well above 0.5 with varying $\eta$. According to the present definition, this is the so-called “anomalous” diffusion.

III.2.3 Case II diffusion

When we choose $\gamma=9.6$, the calculated concentration profiles at different times and values of $\eta$ are as shown in Figs. 15(a)-(c). We see that the sharp front now becomes pronounced, and the concentration gradient behind this front becomes small, which are the characteristic features of Case II diffusion. As $\eta$ becomes larger, the concentration gradient behind the front becomes less steep.

The positions of concentration front at equally spaced time intervals in Figs. 15(a)-(c) indicate a nearly constant diffusion speed. This can be seen more clearly in Fig. 15(d), where the mass-time curves show that the amount of mass intake is proportional to time $t$, and the diffusion speed of the front is constant. The fitting parameter $n$ to the function $f(t)=kt^{n}$ is listed in the fourth row of Table 1.

III.2.4 Super Case II diffusion

As discussed in the theoretical analysis, the rate of mass uptake can rise rapidly when the precursors of the two fronts of diffusing solvent meet at the middle of a sample.

A grid of 20 layers was used to do the simulation. During the simulation, both the $x=0$ and $x=L=19a$ surfaces of the grid were filled by solvent molecules at all times. The value of $\eta$ used here is 5.0, and $\gamma=9.6$. The approximated concentration gradient at position $x=l$ was calculated in the same way as in the previous theoretical section, according to Eqs. (20), (21), and (22).

The corresponding curves of mass intake versus time are shown in Fig. 16 for layers of thickness of 5, 9, and 13. In the plot, where the trends of linear sorption are shown as dashed lines, we see an acceleration in the rate of mass intake at late times as the two fronts meet.

The conclusion from the simulations reported in this section is that our results from a model that is non-local in time are qualitatively similar to those of the time-local theory in Section II. One exception is that the slowing-down effect apparent in the theoretical results in Figs. 6(d) and 8(d) is much less pronounced in the results of the simulations.

IV The relation of our model to the Thomas-Windle model

The widely discussed Thomas-Windle (TW) model can explain reasonably well the Case II diffusion. thom To relate our model to the TW model, we develop our theoretical results a few steps further. Using the relations

[TABLE]

Eq. (12) can be rewritten as:

[TABLE]

We convert to a current density ${\boldsymbol{j}}$ by multiplying the density of sites, $\rho_{0}$, and transform ${\boldsymbol{j}}$ to take the form of Fick’s first law, with

[TABLE]

where we have dropped the subscript $l+\frac{a}{2}$ for $\phi$. We see that $D$ is a function of both $\phi$ and ${\boldsymbol{\nabla}}\phi$. When $\eta$ is small ($\eta<0.5$), expansion of Eq. (29) to terms linear in ${\boldsymbol{\nabla}}\phi$ gives

[TABLE]

where all terms in ${\boldsymbol{\nabla}}\phi$ cancel out, and $D$ returns to being a function of $\phi$ only.

According to the TW model, thom the rate of change of volume fraction of penetrant at each very thin material element is

[TABLE]

where $k_{c}$ is a constant of proportionality. Here $P$ is osmotic pressure, which can be expressed as

[TABLE]

where $k_{B}$ is the Boltzmann constant, $\Omega$ is the partial molecular volume of penetrant, and $a_{T}$ is the penetrant activity. The quantity $\eta_{T}$ is the viscosity of the polymer, and is represented as

[TABLE]

where M is a constant, and $\eta_{T0}$ is the viscosity of the unswollen polymer. (Note: the $\eta_{T}$ here has a different physical meaning from the $\eta$ in our theoretical model.) The diffusion process is governed by the equation

[TABLE]

By integrating the coupled Eqs. (31) and (34) numerically with the extra assumption that $a_{T}$ is a state function of $\phi$, hui2 i.e., $a_{T}=a_{T}\left(\phi(x,t)\right)$, Thomas and Windle were able to reproduce the main features of Case II diffusion. From Eq. (32), we have

[TABLE]

and thus

[TABLE]

We also have

[TABLE]

Substituting Eq. (36) into (37), and then substituting Eqs. (35) and (37) into Eq. (34) yields

[TABLE]

Because $a_{T}$ is a state function of $\phi$, according to Eq. (32), $P$ is also a state function of $\phi$: $P=P(\phi)$. If we define a new entity, $D$, as

[TABLE]

where we have used $\beta=1/k_{B}T$, then $D$ is also a state function of $\phi$: $D=D(\phi)$. From the comparison with Eq. (30), we have

[TABLE]

[TABLE]

Thus the parameter $\eta$ is related to the osmotic pressure $P$ in the TW model. The quantities $\eta$ and $P$ play a similar role in each model. In our model, we simply assume $\eta$ to be constant, and mainly analyze the role played by $D^{*}(\phi)$, while the TW approach calculates $P$ by combining Eq. (31) and Eq. (33) to find

[TABLE]

and fixes the form of $D^{*}(\phi)$ to an exponential function of the concentration $\phi$. We note that if $P$ is a state function of $\phi$, then $\frac{\partial\phi}{\partial t}$ is also state function of $\phi$ only, which puts some constraint on the functional form of $\phi$.

V Relation to Experimental Data

Now let us apply our model to some realistic systems. We choose the benzene-polystyrene system and the toluene-polystyrene system as examples. To calculate $\gamma$, we need to know the potential barrier in the absence of solvent ($V_{b}$) and at saturation ($V_{a})$. These effective potential barrier heights are often referred to as the “activation energy” $E_{\rm D}$ for diffusion. braune Experimental results for the activation energy at these two extreme concentrations are sparse, and so we can only make an estimate based on limited sources. zielinski ; vrentas The detailed analysis leading to this estimation is in the Appendix. It is found that for the benzene-polystyrene system at 30 ∘C, $V_{b}\simeq 4.86\times 10^{-20}$ J, and $V_{a}\simeq 2.78\times 10^{-20}$ J. Thus, $\gamma=\beta(V_{b}-V_{a})=(4.86-2.78)\times 10^{-20}\rm{J}/(1.38\times 10^{-23}\rm{J/K}\times 303\rm{K})\simeq 5.0$. According to our model, this indicates that the diffusion of benzene into polystyrene is anomalous diffusion, which is in accord with the results reported by Long and Kokes. long For the toluene-polystyrene system at 30 ∘C, we have $V_{b}\simeq 9.24\times 10^{-20}$ J, and $V_{a}\simeq 2.57\times 10^{-20}$ J. Then we have $\gamma=\beta(V_{b}-V_{a})=(9.24-2.57)\times 10^{-20}\rm{J}/(1.38\times 10^{-23}\rm{J/K}\times 303\rm{K})\simeq 15.9$, which indicates that the diffusion of toluene into polystyrene should be Case II type according to our model. Experiments by Gall and Kramer for the diffusion of deuterated toluene in polystyrene led to the same conclusion. gall

We have seen in the preceding sections that the Taylor expansion of Eq. (29) at small $\eta$ leads to a concentration-dependent diffusivity, which can also lead to all types of diffusion. Now we give a theoretical derivation of how to verify the existence of a concentration-dependent diffusivity using the diffusion profile observed in experiments. From the non-linear diffusion equation

[TABLE]

we have

[TABLE]

So

[TABLE]

Inverting $\phi=\phi(x,t)$ to $x=x(\phi,t)$, we have

[TABLE]

Ideally, from the experimentally observed concentration profiles, we could calculate $g[x(\phi,t),t]$ according to Eq. (45). If $g[x(\phi,t),t]$ is independent of $t$, then the diffusivity $D(\phi)$ is concentration dependent only, and provides a complete description of diffusion. Unfortunately, the precision of most measurements is not sufficient to achieve this goal. The relevant data can be found in the work of Hermes et. al. hermes and Ogieglo et. al. ogi1 ; ogi2

VI Summary and Conclusions

In this paper we have set up a random walker model to theoretically analyze and a three-dimensional-grid model to simulate the diffusion behavior of small solvent molecules in a glassy polymer matrix. In the theoretical analysis, we generalized the diffusion equation, and then solved it numerically. Results showed that the diffusion type was mainly determined by a dimensionless parameter $\gamma\equiv\beta(V_{b}-V_{a})$, where $V_{b}$ is the potential barrier between cavities in the polymer matrix in the absence of solvent ($\phi=0$), and $V_{a}$ is the lowered barrier at saturation ($\phi=1$). In one limit, when $\gamma=0$, the diffusion is Fickian, while in another limit, when $\gamma\gg 1$, the diffusion becomes Case II type. When $\gamma$ takes intermediate values, it is anomalous diffusion. Super Case II diffusion, which happens when the two opposite sides of a sample are in contact with solvent, originates from the swelling due to solvent molecules in one front propagating elastically into the vicinity of the other front. Using the same assumptions as that in the theoretical model, with an extra memory factor reflecting a previous occupation history of cavities, we performed simulations, and obtained similar results.

We also derived the formal expression of diffusivity by comparing with Fick’s first law. Generally speaking, the diffusivity depends on $\phi$ and ${\boldsymbol{\nabla}}\phi$; but when another dimensionless parameter $\eta$ (which is a measure of the variation of the base potential energy $V_{1}$ in the polymer matrix with respect to solvent concentration $\phi$) in our model is small, the dependence of $D$ on ${\boldsymbol{\nabla}}\phi$ vanishes, and the expression for $D$ becomes a function of $\phi$ only. At this limit the Thomas-Windle model can be related and incorporated into our model.

According to our model, the immense range of diffusion coefficients comes from polymer matrix properties (the average distance between cavities $a$, the height of potential barriers between cavities $V_{a},V_{b}$), environment (temperature $T$), and solvent properties (attempt frequency $\nu$). However, only the relative magnitudes of the diffusion coefficients determine the diffusion type, while their absolute magnitudes only influence the time scale at which diffusion takes place. The present analysis provides an overview that incorporates the many specialized models interpreting specific types of diffusion qian ; thom ; tsige1 ; tsige2 ; ferreira into one unified model.

Acknowledgements.

This work was supported by the National Science Foundation under Grant DMR-1410290 (M.T.) and by the Petroleum Research Fund of the American Chemical Society under Grant PRF# 51995-ND7 (J.M.&P.L.T.), and was made possible by use of facilities at the Case ITS High Performance Computing Cluster.

Appendix A Estimation of the activation energy

The effective potential barrier height defined in Eq. (4) is often referred to as the “activation energy” $E_{\rm D}$ for diffusion. Our work shows that the activation energy for non-Fickian diffusion is concentration-dependent. Duda et al. zielinski calculated the concentration dependence of the activation energy $E$ required to overcome attractive forces in the Vrentus-Duta free-volume diffusion model, according to which the self diffusion coefficient is given by: vrentas

[TABLE]

where $D_{0}$ is a constant pre-exponential factor; $E$ is the energy per mole that a solvent molecule needs to overcome the attractive forces that constrain it to its neighbors; $R$ is the ideal gas constant; $\omega_{i}$ is the mass fraction for component $i$ ($i$ equals to 1 or 2, 1 for solvent and 2 for polymer); $\hat{V}_{i}^{*}$ is the specific critical hole free volume required for jump for component $i$; $T_{gi}$ is the glass transition temperature for component $i$; $\xi$ is the ratio of molar volumes for the solvent and polymer jumping units; $\gamma^{*}$ is an overlap factor (between $1/2$ and 1) introduced because the same free volume is available to more than one molecule, and thus has a different physical meaning from the $\gamma$ in our present model; $K_{11}$ and $K_{21}$ are free-volume parameters for the solvent; $K_{12}$ and $K_{22}$ are free-volume parameters for the polymer. The activation energy for diffusion $E_{\rm D}$, on the other hand, can be defined according to the Arrhenius equation:

[TABLE]

which finally gives us

[TABLE]

Substituting Eq. (47) into Eq. (49), we have the relation between $E_{\rm D}$ and $E$:

[TABLE]

Now let us calculate the $E_{\rm D}$ for self diffusion of benzene in polystyrene at 30 ∘C. The values of parameters for this system are: zielinski $\hat{V}_{1}^{*}=0.901$ cm3/g, $\hat{V}_{2}^{*}=0.850$ cm3/g, $K_{11}/\gamma^{*}=1.48\times 10^{-2}$ cm3/(g$\cdot$K), $K_{12}/\gamma^{*}=5.82\times 10^{-4}$ cm3/(g$\cdot$K), $K_{21}-T_{g1}=-178.52$ K, $K_{22}-T_{g2}=-327.0$ K, $\xi=0.51$, $D_{0}=1.80\times 10^{-3}$ cm2/s. We also have $T=303.15$ K, $R=1.987$ cal/(mol$\cdot$K), and $\omega_{2}=1-\omega_{1}$. Using Eq. (50) and the $E$-$\omega_{1}$ relation in Fig. 4 of Duda et al.’s work, zielinski we are able to obtain the $E_{\rm D}$-$\omega_{1}$ relation as shown in Fig. 17 (a).

One point to be noted is that the $\phi=1$ state in our model, when all the cavities are occupied, does not correspond to the $\omega_{1}=1$ state here, which stands for the pure benzene state, but should correspond to a smaller $\omega_{1}$, say, $\omega_{1}=0.5$. The activation energy at this point is about 4000 cal/mole. This gives an estimation of the value of $V_{a}$ in our model: $V_{a}\simeq 4000\ \rm{cal/mole}=2.78\times 10^{-20}$ J. The largest possible activation energy at low concentration can also be estimated from Fig. 17 (a), and is about 7000 cal/mole, which gives an estimation of the value of $V_{b}$: $V_{b}\simeq 7000\ \rm{cal/mole}=4.86\times 10^{-20}$ J.

Using the same method, we also calculated the $E_{\rm D}$-$\omega_{1}$ relation for the toluene-polystyrene system at 30 ∘C as shown in Fig. 17 (b). The activation energy at $\phi=1$ can be estimated from the value of $E_{\rm D}$ at $\omega_{1}\simeq 0.5$, which is about 3700 cal/mole. The estimation of the value of $V_{a}$ turns out to be: $V_{a}\simeq 3700\ \rm{cal/mole}=2.57\times 10^{-20}$ J. The largest possible activation energy at low concentration is estimated to be 13300 cal/mole. Then the corresponding value of $V_{b}$ is estimated to be: $V_{b}\simeq 13300\ \rm{cal/mole}=9.24\times 10^{-20}$ J.