Mathematical Modeling of Hydrogels Swelling Based on the Finite Element Method ()
1. Introduction
Hydrogels are cross-linked polymeric networks that have the capability to hold water within their porous structure. This property allows using them in applications that require either to absorb or to expel liquids from or towards the environment in which the hydrogel is immersed. Hydrogel applications cover areas such as biomedical engineering, pharmaceutical industry, agriculture, etc. [1]. Some applications involve capturing or delivering liquids or mixtures in a dynamically controlled way, at predetermined rates for predefined periods of time. During the water uptake or delivery process, hydrogels could experience volume changes of several orders of magnitude. According to Singh et al. [2], hydrogels can absorb amounts water nearly 10 - 20 times their molecular weight. Consequently, when a dry solid hydrogel sample is immersed in water, it suffers a continuous volume transformation until the ultimate, and completely swollen state is reached [1].
Depending on the application considered, dynamical description of the swelling process could be important. In some applications, only equilibrium hydrogel states need to be considered and the dynamic transient swelling process of the hydrogel does not need to be described [3-5].
However, some applications require prior knowledge of the transient behavior of the process of uptake or delivery of drugs or solvents. For example, the description of the dynamic swelling process of the hydrogels is very important in applications that require precise control of the absorption of solvents inside the hydrogel structure.
Qualitative descriptions of the polymer relaxation process are reported by several authors [6,7]. Figure 1 shows a diagram of the water sorption process into a hydrogel slab sample.
Initially, a macroscopically homogeneous dry hydrogel sample is immersed in water. Depending on the particular characteristics of the hydrogel sample (material components, polymer structure and porosity), it can retain a maximum value of water in its structure (water equilib-
Figure 1. Qualitative description of polymer relaxation during water sorption into a slab geometry (adapted from [6]). White areas represent glassy polymer regions. Solvent absorption corresponding to state (a) to (d) are detailed in the text.
rium concentration value). Therefore, water penetrates the dry polymer sample travelling from the hydrogel surface to its interior starting the hydrogel swelling process (Figure 1(a)). Thus, simultaneously, the surface of the hydrogel sample (swelling front) moves away from its core, while a mobile interface that separates dry and hydrated polymer is displaced inside of the dry sample (Figure 1(b)). During this process, depending on the specific geometrical characteristic of the hydrogel sample, water concentration can reach its equilibrium value in the swelling front, or the dry-hydrated polymer interface disappears because water is present in all places of the hydrogel sample (Figure 1(c)). Finally, water local concentration reaches its equilibrium value all around the hydrogel sample (Figure 1(d)).
To describe the process of hydrogel volume change driven by water diffusion, different mathematical models have been developed. One approach is the so-called power-law equation where a simple empirical equation, developed by Peppas and Korsmeyer [8], assumes a time-dependent power law function which is used to determine the mechanism of diffusion in polymeric networks
(1)
where Mt and M∞ are the cumulative amounts of water at time t and infinite time (equilibrium), respectively, k is a structural/geometric constant for a particular system and n is designated as “release exponent” representing the release mechanism [9].
Variations to the power-law function approach have been made in order to improve its capabilities to take into account some weaknesses. In particular, to take into account both the solvent diffusion and the polymer relaxation, Equation (1) was modified by Peppas and Sahlin [10]:
(2)
where k1, k2, and n are constants. Terms on the right side represent the diffusion and polymer relaxation contribution to the release profile, respectively.
Despite their simplicity, two drawbacks of the powerlaw function approach limit its use. First, it is necessary to determine the k’s and n coefficients, for each geometrical shape and size of the hydrogel sample. Second, no details of internal diffusion process such as position of the dry-hydrated polymer interface can be provided.
Another family of models have considered the combination of the kinematics of the deformation, the presservation of solvent molecules, the local equilibrium conditions and the kinetics of the molecule migration, in order to simultaneously calculate the temporal evolution of the chemical potential of the solvent and the deformation of the hydrogel [11]. This approach can be implemented using a variational formulation, which leads to a set of coupled governing equations that describe mechanical and chemical equilibrium conditions, along with adequate boundary conditions [12,13]. The finite element method is frequently used to discretize the resulting system of differential partial Equations.
Even if such implementations enable numerical simulations of hydrogel swelling under various constraints, the later approach requires knowing or estimating a large number of physical parameters, such as, for example, values of the volume per solvent molecules, the environmental temperature, two dimensionless material parameters to represent the Flory-Rehner free-energy density function, besides of the coefficient of diffusion of the solvent in the hydrogel [11].
Between the two approaches described, a third class of models is found. In this approach, mathematical models based on Fick’s laws of diffusion [10] overcome limitations of the power-law Equation model and they do not require a very large number of parameters. In these models, change in solvent concentration inside the hydrogel matrix are computed from [14]
(3)
where C is the mass fraction of solvent in the hydrogel matrix, and D is the diffusivity. If the hydrogel deformation due to osmotic stress is taken into account, Equation (3) can be modified to consider an elastic stress diffusion coefficient and the internal stress state of the polymer [9]. Nevertheless, for simple systems, authors as Rossi and Mazich [15,16] and Mazich et al. [17] have established that Fick’s law of diffusion is completely adequate to describe the swelling process.
However, mathematical modeling of the hydrogel swelling process using Equation (3) needs to consider two different elements: the mathematical expression of the diffusivity coefficient and the volume change of the hydrogel matrix during the swelling process. First, to describe the diffusion process in the hydrogel matrix, the diffusivity coefficient cannot be considered constant. Good results [18,19] have been obtained by expressing the diffusion coefficients with a Fujita-type exponential dependence [20].
In addition, due to the intrinsic nature of the diffusion and penetration of the solvent in the matrix, the boundary between the hydrogel and the environment in which is immersed, the swelling front, moves during the swelling process (Figure 1). Consequently, the solution domain of the diffusion Equation (3) changes in time and a moving boundary condition needs to be included.
As a consequence of these last considerations, there are no analytical solutions to the mathematical model (3) with mobile boundary conditions, and these equations must be solved using a numerical method.
To date, most of the numerical models developed are based on the finite difference methods. Even though numerical models based on in finite differences can be easily implemented, their use is limited to samples whose shape is very simple. Consequently, most applications are developed considering axial and radial diffusion within cylindrical tablets [18].
Figure 2 shows the typical representation of the cylindrical tablets and the symmetry plane and the symmetry axis, indicating the rectangular domain that represents the overall cylindrical tablet. The plane z = 0 cuts just in the top half of the tablet. The z-axis is a symmetry axis. Thus, a reduced computational domain, the rectangle OABC, can be used to model the solvent inflow into the
Figure 2. Rectangular domain OABC representing overall cylindrical tablet for finite difference model.
cylindrical tablet. In the domain OABC, a rectangular grid can be generated. Therefore, the diffusion equation is easily discretized using a finite differences scheme.
During the swelling process, volume variations of the hydrogel sample are represented in a dynamically deformable grid that is computed each time step. Mathematical models based on this approach have been successfully used in a vast range of applications [18,19]. Siepmann et al. [18] proposed a mathematical model for modeling the water transport into and drug release from hydroxypropyl methylcellulose (HPMC) cylindrical tablets. This model [18] considers that the tablet can swell in both axial and radial directions. The swelling is assumed to be ideal; the total volume of the tablet at any instant is given by the sum of the volumes of different components: polymer, water, and drugs. Thus, as the water penetrates the tablet, the total volume changes. In this manner, the tablet radius Rt and the half-height Ht change with time. A rectangular grid is defined in the region q = 0, 0 £ r £ Rt and 0 £ z £ Ht. Once the new volume is computed in each time interval, the new values of Rt and the Zt are calculated. Due to that, [18] assumed that water imbibing in the axial direction leads to a volume increase in the axial direction, whereas water imbibing in the radial direction leads to a volume increase in the radial direction; the increase in volume in each direction is proportional to the surface area in this direction. In this approach, the shape of the tablet remains cylindrical during all the swelling process.
However, as it has been shown by Achilleos et al. [21], swelling is not a continual process [1]. Achilleos et al. [21] have developed a technique for the real-time visualization of dynamic deformation profiles during gel swelling processes. Achilleos et al. [21] visualizations allow appreciating that the evolution of the dry tablet shape does not necessarily preserve the initial cylindrical shape during all stages of the swelling process. When a hydrogel in its initial state is in contact with solvent molecules, the latter penetrate into the polymeric network. Obviously, the region close to the intersection between the plane surface (z = Zt) and the cylindrical surface (r = Rt) is exposed to a solvent flow on both surfaces. Consequently, a more substantial increase in volume is expected in this region compared to the rest of the surface of the tablet.
On the other hand, the use of finite difference schemes in hydrogel samples different from cylindrical tablets can be restrictive because of the difficulties to generate a computational grid where diffusion equations can be easily discretized. When non-rectangular grids are considered, the interpolation process at the boundaries requires special considerations.
To overcome modeling difficulties associated to the sample shape, the finite element method can be used.
Finite element methods allow the use of adaptive grids non-restricted to rectangular elements. Therefore, non a priori restrictions to the sample shape during the swelling process need to be imposed. Additionally, initial sample shapes not limited to cylindrical tablets can be modeled. Even though the flexibility of the finite element method allows modeling arbitrary tridimensional sample shapes, for practical purposes, it is convenient, from a computational point of view, to restrict sample shapes to bodies with a symmetry axis. Fortunately, this is the preferred shape in most practical applications. Figure 3 shows typical shapes that can be easily modeled with this approach.
In this study, the swelling dynamics of initially dry hydrogel samples is numerically modeled using the finite element method. Due to the numerical approach chosen, the sample shapes to be considered can correspond to those of any solid of revolution. The mathematical model is based on Fick’s second law of diffusion. The diffusion coefficient is expressed by a Fujita-type exponential dependence. The resulting equations are solved by using a deformable grid. Grid displacements are computed at the end of each time-step considering local variations in water concentration as a result of the diffusion process. The new positions of the grid nodes are computed by solving a set of non-linear equations that is solved for each timestep. This iterative process is done until the hydrogel sample reaches the final equilibrium state. Model validation was performed by the comparison between numerical results and laboratory tests, considering the swelling of hydrogel samples of different shapes (cylindrical tablets, conical and capsule samples).
This paper is organized in the following way. First, the mathematical model is presented. Equations for modeling the water diffusion process and sample shape evolution
Figure 3. Typical sample shapes and computational domain (in white) (a) cylindrical; (b) pastille; (c) capsule; (d) composed cone-cylinder and (e) general axisymmetry.
are developed. Then, a description of the numerical model, in particular, the specific implementation of the finite element model used to solve the diffusion equation and the finite difference scheme used to solve the grid equations, is carried out. Later, the experimental section describes the materials and methods used to prepare different hydrogel samples. Once the swelling laboratory tests for all hydrogel samples were finished, comparisons with numerical predictions were performed. Finally, following the analysis of results, model assumptions are evaluated and suggestions for further numerical and experimental work are provided.
2. Mathematical Model
In this study, the dynamics of the hydrogel swelling process is modeled. The following assumptions have been considered: 1) Dissolution processes are neglected, i.e., all solid material remains attached to the original sample; 2) Swelling can be modeled as a pure diffusion process; 3) Water concentration at the surface of the tablet is at its equilibrium value.
As a consequence of assumption 2) and following [18, 19], the mathematical model for modeling the water transport into the hydrogel samples is based on Fick’s second law of diffusion. A general sample and the cylindrical coordinate system are represented in Figure 4, just before the swelling process starts.
Each point P in the hydrogel sample is represented by the coordinates r, q and z. The radial coordinate r, represents the distance between the projected point of P, P’, in the plane xy and the origin O of the coordinate system; the axial coordinate, z, is the coordinate along the symmetry axis and q is the angle between the segment OP’ and the x axis.
Figure 4. General axisymmetrical sample and the coordinate system. S is the sample surface.
As there is no concentration gradient of any component with respect to q, in cylindrical coordinates, Fick’s second law (3) is expressed as:
(4)
where C and D are the concentration and the diffusion coefficient of water, respectively, r is the radial coordinate, z is the axial coordinate and t represents time.
To describe the diffusion process in the hydrogel matrix, the diffusivity coefficient is modeled following [18]. Thus, the diffusion coefficients is expressed with a Fujita-type [20] exponential dependence
(5)
where Deq and Ceq are the diffusion coefficient and solvent concentration in the equilibrium swollen state of the system, respectively, and b is a dimensionless constant characterizing the concentration dependence of D.
The hydrogel sample is allowed to swell, independently, in both axial and radial directions. The swelling is assumed to be ideal. Hence, the total volume of the sample at any instant is given by the sum of the volumes of polymer and water. However, unlike the approach of [18], the hydrogel shape is not necessarily a cylindrical tablet. In particular, any form that corresponds to a solid of revolution can be modeled, as it is shown in Figure 4.
At t = 0, the sample is dry and thus the water concentration at any position inside the sample is equal to zero while that at the sample surface, S, the water concentration is Ceq. Thus, the initial conditions in each point P of the hydrogel sample are expressed as:
(6)
During the swelling process, at the surface of the tablet, the concentration of water is assumed to be at its equilibrium value, Ceq. This boundary condition is written as follow:
(7)
During the swelling process, as a consequence of water inflow inside the hydrogel sample, the volume of the sample will change. Therefore, solution of the diffusion equation requires considering that the physical domain is changing along the swelling process. In addition, the sample shape at each time is not known in advance and it is part of the solution.
3. Numerical Model
In this work, no assumption regarding the sample shape evolution was considered. In addition, because sample shapes are not limited to cylindrical tablets, the selected numerical method should be able to manage irregular meshes. For these reasons, the finite element was chosen to discretize the mathematical model.
3.1. Finite Element Formulation
In order to obtain the finite element formulation of the problem, a weighted residual method, using the Galerkin formulation [22], was employed. The residual is defined from Equation (4) as:
(8)
This residual is minimized by imposing:
(9)
where 
is an arbitrary weight function that vanishes wherever a Dirichlet condition is specified, and dv = 2πrdA since the problem is axisymmetric.
Introducing (8) in (9) and manipulating the integrals, the weak formulation of the problem is obtained:
(10)
Using a Galerkin approach, the concentration and the weight function are interpolated from the nodal values by mean of the same matrix of interpolating functions:
(11)
With these interpolations, Equation (10) is transformed to:
(12)
where 
and 
are element matrices and vectors:
(13)
(14)
(15)
where 
represents the matrix of the derivatives of the interpolation function.
Finally, the finite element formulation of the problem is obtained by assembling all the elements in one equation:
(16)
In this work, the finite element model was constructed using axisymmetric linear triangular elements of three nodes and one degree of freedom by node, corresponding to the concentration at that point.
3.2. Grid Generation
A grid made of triangular elements was generated in several steps for t = 0. Initially, the boundary profile of the symmetry plane of the hydrogel sample is represented in the r-z plane (Step 0). Then, a first set of nr segments was generated by building a number of straight lines between the origin of the coordinate system and the external boundary of the sample (Step 1). Thus, a second set of nz lines is generated. Each line of this second set intersects those of the first set (Step 2). At this stage, the grid includes triangular and quadrilateral elements. Finally, quadrilateral elements are divided to form triangular elements (Step 4). Figure 5 shows the different steps for the generation of the initial grid.
Figure 6 shows the grid and the associated computational domain.
A set of indices i and j, 1 £ i £ nr and 1 £ j £ nz, represents each node (i, j) in the physical domain. Each node (i, j) has its representation in the computational domain. The origin in the physical domain is a “degenerated node”. It means that all nodes with j = 1 represent the origin O.
The total number of nodes is nr ´ nz, whereas the total number of elements is
.
Once the hydrogel sample starts the swelling process, the physical domain, and consequently, the computational domain changes and a new mesh must be generated each time-step. Then, it is necessary to establish nr ´ nz equations that allow computing the new position of each node. The ensemble of equations is derived from: 1) Restricting the grid nodes to move along the nr original lines (first set of lines); 2) symmetry condition along the z-axis and 3) volume of each element must contain the initial polymer mass and the water imbibed in.
Figure 6. k-element in the physical and computational domain.
Condition 1) requires that, if qi is the angle between the i-line and the r-axis, then
(17)
The symmetry condition 2) imposes that grid lines are orthogonal to the z-axis. Thus,
(18)
The resulting algebraic equation set from expression (18) can be expressed by using a forward or backward second order expression in finite differences. For example, for i=1, using forward finite differences
(19)
Finally, condition 3) needs to consider spatial local variation of water in each elementary grid element in time. In consequence, equations that express the relation between local water concentration and the volume of each element must be established.
The strategy to generate the new grid considers that the local mass of each mesh elementary volume at time t, mk, is given by
(20)
where mpk and mwk are the polymer mass and water mass at element k. The mass mpk is constant and it corresponds to the original mass of dry polymer at t = 0 at element k.
In terms of the local water concentration, the mass in the k-element can be expressed as
(21)
where rp and rw are the polymer and water volumetric density, respectively, Vk(t) and Ck(t) are the total volume and the water concentration in the k-element.
A expression such as Equation (21) can be written at t = t + Dt
(22)
Because the polymer mass is constant in the k-element, the first term on the right side of Equations (21) and (22) are equal and
(23)
Taking advantage of the geometrical properties of the hydrogel sample, the elementary volume for a solid of revolution is calculated from:
(24)
where rk(t) and Ak(t) are the centroid position and the area of the triangle associated to the k-element at the instant t.
With (24) in (23) after some simplifications it gives
(25)
The non-linear Equations (17), (19) and (25) are a complete set of algebraic equations that permits generating a new mesh each time-step.
However, when all non-linear nr ´ nz algebraic equations are solved simultaneously, some convergence problems were found. In order to avoid these problems, a numerical solution solving layer by layer was implemented. In this approach, the origin coordinates (corresponding to j = 1) are constant. The implemented procedure starts by finding the position of the grid nodes corresponding to j = 2 for t = t + Dt. Hence, only nr nonlinear algebraic equations must be solved. Once these nodes coordinates are found, the next level, j = 3, is solved. This procedure continues until the node coordinates corresponding to the level j = nz are found. In this way, a set of nr non-linear algebraic equations are solved nz times instead of a single nz ´ nr system. This procedure converged in all simulated cases.
An additional simplification, with high impact on computational time, can be introduced when the physical domain has a symmetry plane, as in the case of a cylindrical tablet. In such cases, a symmetry condition as Equation (18) can be implemented. Then, if the symmetry plane corresponds to z = 0, the symmetry condition will be
(26)
4. Materials and Methods
To validate the numerical model, a number of laboratory tests were conducted considering different sample shapes such as cylindrical tablets, capsules, pastilles and composed cone-cylinder (Figure 3). The hydrogel samples were prepared using the following materials and methods.
4.1. Materials
The following compounds were obtained from Sigma Aldrich: Acrylamide (AAm, +99% purity), N’N’metilenebisacrylamide (NMBAAm, +99% purity) and Ammonium Persulfate (PSA, +98% purity).
4.2. Methods
Polyacrylamide hydrogel was synthesized by chemical initiation, using AAm (1 g) as monomer, and N’N’Methylenebisacrylamide (11 mg) as crosslinking agent; these were mixed in 0.7 ml of distilled water for 5 minutes; then PSA was added as initiator of the reaction. The synthesis was carried out in a thermal bath, at a temperature of 30˚C ± 1˚C, under constant stirring of 200 rpm. The resultant monolith was cut in discs of different lengths. The measurements of absorption were obtained introducing the dry hydrogel discs (previously weighed, mo) in distilled water at 37˚C. The swollen hydrogels, for predetermined times (up to equilibrium swelling), were removed and gently dried with a filter paper to remove excess water and finally weighed (mt) in a scale (Adventure-OAUS with an accuracy of ±0.0001 g). The swelling percentage (%S) of the hydrogels was determined from the following relation:
(27)
During the swelling process, special care was taken to assure that the hydrogel sample was fully immersed in the distilled water avoiding contact with the walls of the recipient.
5. Numerical Results
Applications of the numerical model require advanced knowledge of three unknown parameters: the diffusion coefficient Deq, the dimensionless constant b, and the solvent concentration in the equilibrium swollen state of the system, Ceq.
The constant b is characteristic of the diffusing molecules and, for water, a value of 2.5 was calculated by Siepmann et al. [18]. The diffusion coefficient of water Deq within the fully swollen hydrogel must be determined by fitting numerical predictions with experimental water uptake data of hydrogel samples. The other parameterCeq, depends on the particular capability of the xerogel to absorb water and is determined for measuring the total amount of water in fully swollen hydrogel samples.
Since the main objective of this work is to develop a numerical method that considers the actual shape evolution of the hydrogel sample during the solvent absorption process using few parameters, the diffusion coefficient of water Deq was computed by numerical regression in order to fit the experimental data. Consequently, influence of variations in Deq due to changes in experimental conditions (size of recipient, position of the hydrogel particular sample in the original bulk, etc.) was not considered.
To validate the capabilities of the numerical model, three different kinds of samples were considered: 1) tablets; 2) capsules and 3) truncated cones.
Six different cylindrical tablets were considered. The particular characteristic of each sample are reported in Table 1.
Figure 7 shows the transient water uptake for each cylindrical tablet.
The match between experimental data and numerical model predictions is remarkable. In all analyzed cases, the proposed model was able to reproduce the dynamic water uptake.
A typical sample shape evolution and water penetration is shown in Figure 8. Initially, water uptake is more pronounced in the intersection of the top and bottom of the sample and the cylindrical wall. Consequently, this region swells more than the other surfaces of the sample.
This trend is in good agreement with experimentally observed behaviour, as shown in Figure 9.
Qualitatively, it is clear from Figure 9 that numerical predictions of shape evolution in time correspond to experimental behaviour.
Other sample with shape corresponding to a truncated cone was modelled. Figure 10 shows the initial shape of the truncated cone. This kind of sample shape cannot be modelled using the classical finite difference models.
Figure 11 shows the water uptake process for the truncated cone. The finite element model was able to reproduce accurately the swelling process.