Biomacromolecules 2005, 6, 109-120
109
Modeling and Simulation of the Swelling Behavior of pH-Stimulus-Responsive Hydrogels Hua Li,*,† Teng Yong Ng,†,‡ Yong Kin Yew,† and Khin Yong Lam† Institute of High Performance Computing, National University of Singapore, 1 Science Park Road, #01-01 The Capricorn, Singapore Science Park II, Singapore 117528, and School of Mechanical and Production Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798 Received June 17, 2004; Revised Manuscript Received August 31, 2004
The modulation of the swelling ability of the hydrogel matrix by pH-stimulus enables the dynamic control of the swelling forces, thereby obtaining effective diffusivity and permeability of the solutes, or mechanical energy from the hydrogel. In this work, a chemo-electro-mechanical model describing hydrogel behavior, based on multi-field effects, is developed to simulate the swelling and shrinking of these fascinating biomaterials, and it is termed the multi-effect-coupling pH-stimulus (MECpH) model. This model accounts for the ionic fluxes within both the hydrogel and solution, the coupling between the electric field, ionic fluxes, and mechanical deformations of the hydrogel. The main contribution of this model is to incorporate the relationship between the concentrations of the ionized fixed-charge groups and the diffusive hydrogen ion, which follows a Langmuir isotherm, into the Poisson-Nernst-Planck system. To validate this MECpH model, one-dimensional steady-state simulations under varying pH solution are carried out via a meshless Hermite-Cloud methodology, and the numerical results are compared with available experimental data. It is shown that the presently developed MECpH model is accurate, efficient, and numerically stable. Introduction Hydrogels are widely employed in bio-material related fields not only because of their excellent bio-compatibility and bio-stability characteristics, which can be tailored, but also due to their being relatively inexpensive.1,2 Lately, the applications of stimuli-responsive hydrogels as active elements have been explored intensively in areas, such as controlled drug release, microscale actuators/sensors, microfluidic flow control, and filtration/separation processes.3-9 Even though there are several research works detailing the experimental and theoretical studies concerning hydrogels, the responsive behaviors characteristics of these hydrogels remain poorly understood. Siegel10 provided an extensive review of thermodynamic models for describing the swelling behavior of ionic hydrogels.11-18 Despite its simple approach, the general thermodynamic formulation only provides for restricted practical use because, in most cases, the quantities involved cannot be determined by independent measurements. Therefore, models with particular properties have to be introduced for deriving equations which will describe the action of various physical forces pertinent to the particular properties of the actual system.12 Lai’s group19-21 developed a triphasic chemo-electromechanical model to describe the behavior of soft tissues, such as charged-hydrated tissues. This theory was verified for the one-dimensional equilibrium results of swelling, while * Corresponding author. Tel: (65) 64191249. Fax: (65) 64191380. E-mail:
[email protected]. † National University of Singapore. ‡ Nanyang Technological University.
neglecting geometrical nonlinearities. In this model, an assumption of “electroneutrality” condition is made thereby constraining the application range to a few particular cases. In the present study, the hydrogel is considered as a triphasic hydrophilic mixture exhibiting swelling or shrinking behaviors when immersed in a solution of varying pH. The three phases involved are the solid matrix phase, interstitial fluid phase, and the ion phase. From an experimental viewpoint, diffusion of ion species occurs between the porous hydrogel and surrounding solution when the hydrogel is immersed in a buffer solution. Due to the presence of the fixedcharge groups, there are differences in ionic concentrations between the interior hydrogels and exterior solution, generating an osmotic pressure. This pressure is the driving source of the swelling mechanism of the hydrogel. The resulting swelling or shrinking subsequently redistributes the ionic concentrations within the hydrogel. Time will elapse before the recurrent processes reach the steady-equilibrium state. A chemo-electro-mechanical model for describing hydrogel behavior is developed, and it is termed the multi-effectcoupling pH-stimulus (MECpH) model. It is presently employed for the simulation of the behavior of pH-responsive hydrogels and is developed based on the Poisson-NernstPlanck (PNP) equations. The Nernst-Planck flux system is used to describe the transport mechanisms of ionic species in a solution.12,22-27 The Poisson equation which describes the spatial distribution of electric potential completes the PNP system. This Poisson equation relates the electric potential distribution to the various ionic fluxes. More importantly, the presently developed MECpH model also couples the
10.1021/bm0496458 CCC: $30.25 © 2005 American Chemical Society Published on Web 10/29/2004
110
Biomacromolecules, Vol. 6, No. 1, 2005
Li et al.
mechanical equilibrium equations, with the PNP equations, for simulation of the hydrogel deformation. Another important contribution of the presently developed MECpH model is the incorporation of the relationship between the fixed charges attached to the hydrogel network and the diffusive hydrogen ions, into the PNP system. The contributions of the fixed-charge groups are emphasized in this study, as the discrete or continuous volume change of the hydrogel with changes in ambient bath pH, temperature, or externally applied electric field, are associated with the ionized charge groups bounded to the hydrogel network.28 The presently developed MECpH model for pH-responsive hydrogels is able to simulate the concentration distributions of all diffusive ionic species, as well as the electric-potential distribution and mechanical deformation of the hydrogels, in response to the changes of pH of the bath solution. Advanced techniques of preparing polymer gels make it possible for the fabrication of the hydrogels with a wide range of properties such as structure, Young’s modulus, degree of cross-linking, softening point, etc. It is thus also rational to investigate separately the effects of the various physical factors by means of systematically and independently varying the various properties, and these investigations will be carried out in the present study. Theoretical Formulation MECpH Model. Based on the works of Wallmersperger et al.,26,27 one possible way of describing the swelling/ shrinking behavior of a hydrogel is to use a coupled chemoelectro-mechanical multi-field formulation. The generalized Nernst-Planck equation of transport, describing the flux of the ionic species k in the solution, is given as follows:
(
Jk ) - [Dk] grad(ck) +
)
zk F c grad(ψ) + ckgrad(lnγk) + RT k ckVk (k)1, 2, 3, ..., N) (1)
where Jk (mM/s) is the flux of the kth species, N is total number of diffusive ionic species, [Dk] (m2/s) is the diffusivity tensor of the kth species, ck (mM) is the concentration of the kth diffusive ionic species, zk is the kth-ionic valence number, ψ (V) is the electrostatic potential, γk is the kthionic species chemical activity coefficient, and Vk is the kthionic species convective velocity. F, R, and T are the Faraday’s constant (9.6487 × 104 C/mol), universal gas constant (8.314 J/mol K), and absolute temperature (K), respectively. In the eq 1, the first term on the right-hand side characterizes the diffusion flux due to the gradient of the concentration and the second term the migration flux arising from the gradient of the electrical potential. The third term is related to the gradient of the chemical activity coefficient. The fourth term refers to the convective flux produced from the gradient in fluid velocity. It should be noted that chemical activity or reaction can be perceived as the ion-ion interaction, which includes association and dissociation activities. There are a number of models used to compute the chemical activity coefficient, but the more established one is the
Debye-Huckel model, ln γi ) - Az2i xI, where A is a temperature-dependent parameter and I is the ionic strength of the solution.18 Specially, when an external electric field is imposed, the rate of the ionic diffusion is much faster than the kinetics of chemical activity. Then the contribution of chemical activity coefficient becomes much small and can be neglected. According to the law of mass conservation, the change in the amount of the diffusive ionic species contained in the volume with respect to time t is given by the difference between the fluxes entering and leaving the reference volume. As is well-known, for an unstirred solution in vibration-free experimental device, the fluid movement is eliminated and then the convective flux should be neglected.29 In addition, for an ideal solution that is infinite dilution (very low concentration), the chemical activity coefficients of ionic species are always equal to 1. With increasing the ionic concentration of the solution, the chemical activity coefficient will vary from unity. A study on the influence of chemical activity coefficient is conducted by Samson et al.30 Therefore, if an ideal solution unstirred is considered here, the mass conservation of the Nernst-Planck type can be derived for the kth diffusive species as ∂ck ∂ck + div(Jk) ) + ∂t ∂t
{
(
diV - [Dk] grad(ck) +
)}
z kF c grad(ψ) ) 0 RT k (k)1, 2,...N) (2)
where the diffusion coefficient Dk is determined by the Einstein relationship Dk ) µkRT
(3)
Here, µk is defined as the mobility of the kth species ion. In eq 2, the number of equations is equal to the number of unknown diffusive ionic species ck. In addition, there is an unknown electrical potential ψ, and therefore, we have the Poisson equation to describe the spatial distribution of the electric potential in the domain and is given as ∇2ψ ) -
F 0
N
∑k zkck + zfcf)
(
(4)
where is the relative dielectric constant of the surrounding medium and 0 the vacuum permittivity or dielectric constant (8.85418 × 10-12 C2/Nm2). MacGillivray and Hare22,31 have shown that the constant field hypothesis and electroneutrality assumption are in fact the limiting cases of the Poisson equation. In the case of low salt concentration, the present PNP system is reduced to the constant field problem, by assuming a constant electric field across the considered system, i.e., a linear variation of electrical potential across the system. On the other hand, when the salt concentration is high, the electroneutrality assumption is applicable. Samson et al.25 have made a good discussion for these two limiting cases. The velocity of propagation in the electric field is much higher than that due to the effects of convection and diffusion. Hence, the quasi-static form of the Poisson
Biomacromolecules, Vol. 6, No. 1, 2005 111
Swelling Behavior of Hydrogels
equation is sufficient for the present model. The Poisson equation derived in eq 4 is based on the second and fourth Maxwell equations. The charge distribution can be summarized to be governed by the field equations as given below Eel ) - ∇ψ, Del ) 0E, ∇Del ) Fel, N
Fel ) F(
∑k zkck + zfcf)
(5)
where Del is the dielectric displacement tensor and Eel the electric field strength tensor. Fel represents the electrical charge density and it is a function of the concentration of all ions. However, the Poisson equation in eq 4 is a general form derivation. The fixed-charge group cf is one of the hydrogel properties that will constitute the electrochemical reaction of the pertaining hydrogel. Since the hydrogel of present interest is a pH responsive type of hydrogel, cf should be represented as a function of the pH of the surrounding environment. From a comprehensive literature review, it was found that the Langmuir isotherm is most suitable for computing the instantaneous concentration of the fixedcharge groups. In the present study, the diffusing ions can react with the reactive sites which refer to those fixed-charge groups bound to the hydrogel network. Following the approach outlined by Grimshaw et al.,28 the Langmuir isotherm equation, which is usually used to describe the heterogeneous binding reactions at the interface of two phases, is incorporated into the numerical simulation to compute the concentration of the fixed-charge. For the particular case of a binding site where only one ion that can bind to it, the Langmuir isotherm is derived as zfcf )
(
)
s 1 zfcm0K H K + c H+
(6)
where cf is the concentration of the fixed charge groups on the polymer chains and zf the valence of the fixed charges. For example, zf is -1 if we use the carboxylic acid groups as the fixed charges on the polymer chains. K is then the s dissociation constant of the carboxylic acid groups. cm0 is the total concentration of the ionizable groups in the hydrogel at zero state, and cH+ the concentration of hydrogen ions H+ within the hydrogel. H is the local hydration of the hydrogel and defined as the ratio of fluid volume Vf to solid volume Vs, i.e., H ) Vf/Vs. Therefore, the complete Poisson equation will take the form ∇ ψ)2
(∑ ( ))
F 0
N
zkck +
k
s 1 zfcm0K
H K + c H+
(7)
For the simulation of the mechanical deformation of the hydrogels stimulated by the pH of surrounding solution, we first look at the general equation of motion ∂u ∂2u F 2 |i + f |i ) σij,j + Fbi ∂t ∂t
(8)
where F is the mass density, t is the time variable, and u is
Figure 1. Computational domain and the boundaries conditions required for the solution of numerical simulation.
the displacement of the solid phase. f represents the viscous damping coefficient of solid-phase network, σij is the total stress tensor of the hydrogel, and bi is the body forces that is not of interest in present work. In this study, we are dealing with a steady-state equilibrium deformation, and based on the Biot’s elasticity and consolidation theory for porous material,32 the equation is given as ∇σ ) 0
(9)
Following the main assumptions of Gregor’s model,12 the stress contributions in the hydrogels are idealized to arise from two main sources, namely the mechanical elastic restoring stress tensor exerted by the solid-phase network Felastic-restoring (Pa), and the osmotic pressure Posmotic (Pa). For the present study, the simulation results are compared with the experimental works of Beebe et al.9 In the experiment, a cylindrical hydrogel in a microchannel is considered. The hydrogel in the microchannel is constrained from axial displacement by the cover glasses placed at the two ends. For that reasons, the hydrogel is assumed to undergo displacement in the radial direction as shown in Figure 1. For a one-dimensional steady-state case, the mechanical equilibrium equation, which describes the geometric nonlinear deformation of the hydrogels, is written as
{
}
∂u 1 ∂u 2 ∂ ∂σ ) + (2µ + λ) - Posmotic ) 0 ∂x ∂x ∂x 2 ∂x
(
( ))
(10)
where λ and µ are the Lame´’s coefficient33 of the solid matrix. The osmotic pressure is calculated according to the concentration difference between the stress-free and swelling states14 Posmotic ) RT
∑k (ck - c0k)
(11)
where c0k is the concentration of the kth ion species in the stress-free state in bathing solution and ck is the concentration of kth ion species within the hydrogel. For the purposes of computation, the equations are rewritten in nondimensional form by introducing a set of dimensionless parameters for the present problem
Biomacromolecules, Vol. 6, No. 1, 2005
112
ζ)
x , Lref
uj )
u , Lref
cjk )
ck , cref
Li et al.
cf , cref
cjf )
ψ h)
ψ Fψ (12) ) ψref ηRT
where ζ, uj, cjk, cjf, and ψ h are the dimensionless variables of coordinates, displacement, diffusive ionic concentration, fixed charge density, and electric potential, respectively. Lref, cref, and ψref are the characteristic length, concentration and electric potential, respectively. At steady state, the time derivative c˘ k in eq 2 is zero. This significantly reduces the computational effort when solving the system of coupled nonlinear partial differential equations. Thus, the final form of the nondimensional system of partial differential equations to be solved can be written as ∂2cjk ∂ζ2
+ ηzk h ∂ 2ψ ∂ζ2
∂cjk ∂ψ h ∂2ψ h + ηzkcjk 2 ) 0 ∂ζ ∂ζ ∂ζ
(∑ ( ))
2 F2 Lref cref
)-
0RT
η
( ( ))
∂ ∂uj
1 ∂uj
(2µ + λ) + ∂ζ ∂ζ 2 ∂ζ
2
(k ) 1,2,3,...N) (13)
N
zkcjk -
k)1
∂
1
s jcm0 K
H K + cjH+
(14)
N
function, an ideal window functions should fulfill two conditions: (1) it is orthogonal in nature and (2) its integral over the Ω domain should be unity. However, a suitable window function that simultaneously satisfies both conditions does not normally exist. Hence, only an approximation can be obtained in most cases. Based on the key features of the classical reproducing kernel particle method (RKPM), by employing a correction function C(x,y,p,q) and a kernel function K(xk - p,yk - q), the approximate unknown function can be written as ˜f (x,y) )
(
)
∑
(16)
It is noted that, in the present model, there are a few limitations in which the chemical reaction and the convection transport of ionic species have been neglected. Meshless Methodology Development: Hermite-Cloud Method. To solve the present MECpH model with difficulties such as multi-field and multiscale problems, domain remeshing, moving boundaries, and coupled nonlinear partial differential equations, a novel meshless numerical technique, termed the Hermite-Cloud method34 has been developed for the simulation of these pH-stimulus-responsive hydrogels. The Hermite-Cloud method uses the Hermite interpolation theorem for the construction technique for discretization of the partial differential equations. This true meshless technique is based on the classical reproducing kernel particle method35-39 except that a fixed reproducing kernel approximation is employed instead. Through the constructed Hermite-type interpolation functions, we are able to generate the expressions of approximate solutions of both the unknown functions and the first-order derivatives, in a direct manner. A set of auxiliary conditions has to be developed so as to construct a complete set of PDEs with mixed Dirichlet and Neumann boundary conditions. The main characteristic of the reproducing kernel function is its ability to reproduce itself by integration transform over the domain of interest. For example, if we take a twodimensional case, the reproducing kernel expression can be written as the product of the window function, φwindow(xp,y-q), and the real function of f(x,y), and integrated over the defined domain, Ω. To exactly reproduce the real
(17)
where the kernels are centered at points (xk,yk), and they are therefore known as fixed kernel functions. Usually, the correction function can be expressed as a sum of linearly independent basis functions. The highest order polynomial terms in this correction function will depend on the highest order derivative terms contained in the governing PDE of interest. For instance, by choosing β as the order of polynomial to be included in the basis function, we can express the correction function as
- crefRT ( (cjk - cjk0)) ) 0 ∂ζ k)1 (15)
s 1 cjm0K cjf ) H K + cjH+
∫Ω C(x,y,p,q)K(xk - p,yk - q)f(p,q) dp dq
C(x,y,p,q) ) B(p,q)C*(x,y)
(18)
where B(p,q) is a βth-order row vector containing the basis functions and C*(x,y) is a βth-order column vector containing the coefficients. A linearly independent basis function for a one-dimensional quadratic PDE problem is given as B(p) ) {b1(p),b2(p),...,bβ(p)} ) {1,p,p2} (β ) 3)
(19)
The correction function coefficient vector can be written as C*T(x,y) ) {c1,c2, ‚‚‚,cβ} (i ) 1,2,..., β)
(20)
where ci are the unknown coefficients which can be determined by satisfying the following consistency conditions of the reproducing kernel technique bi(x,y) )
∫Ω B(p,q)C*(x,y)K(xk - p,yk - q)bi(p,q) dp dq (i ) 1,2,..., β ) (21)
The kernel function K(xk - p,yk - q) in eq 21 may be constructed by different forms of weighted functions depending on the governing equations of the problem being considered. For the present work, a cubic spline function is employed K(xk - p,yk - q) ) W*((xk - p)/∆x) × W*((yk - q)/∆y)/(∆x∆y) (22) By introducing z ) (xk - p)/∆x for the x component and z ) (yk - p)/∆y for the y component, the cubic-spline window function W*(z) is defined as
{
Biomacromolecules, Vol. 6, No. 1, 2005 113
Swelling Behavior of Hydrogels
0 (2 + z)3 6 2 z - z2 1 + 3 2 W*(z) ) 2 z - z2 1 3 2 (2 - z)3 6 0
( (
z
