The understanding of hydrogel swelling response in different environments is essential for its use in different practical applications. This necessitates its simulation in steady state and transient conditions. This paper mainly deals with the details of the numerical simulation performed by developing coupled formulation of chemo-electro-mechanical behavior of the hydrogel in response to changing pH of the surrounding solution. Simulations were performed to determine the response of hydrogel with varying pH of the surrounding solution in a wide range of pH (2-12), mainly representing its response for different solutions. The investigation of the responsiveness of the hydrogels is focused mainly on the chemical composition of the immersed solution and to change in the solvent. The methodology used for this finite element based simulation is presented. The swelling characteristics of the hydrogel (of arbitrary geometry) with different functional group (PKa) are compared in steady state and transient conditions. This analysis is carried out using COMSOL and the effects of fixed charge density, buffer solution concentration and solution pH on the swelling were studied in different simulations. The mobile ion and fixed charge concentrations showed to take around an hour before they stabilize in response to a given change in pH. Swelling as well as deswelling responses showed similar variation in the transient conditions. These simulation results are compared with available experimental evidence to show the accuracy of the model.
Keywords: Hydrogel, multiphysics models, pH sensitive hydrogel
Introduction
Hydrogels are polymers that respond to environmental factors such as temperature, pH, elecic potential and, light by either swelling or deswelling. Their ability to absorb water is attributed to presence of hydrophilic functional groups which are attached to the polymeric network, while the crosslinking prevent complete mixing of the hydrogel from dissolving in the solvent by producing an elastic restoring force that counters the expansion of the network. In this paper, we modeled a hydrogel as well as the dynamics that enables it to swell and Deswell in response to surrounding pH changes which either swells or shrinks the hydrogel size.We used chemo-electro-mechanical properties of hydrogel to model the equilibrium swollen state of the hydrogel.The response of the hydrogel to pH was determined by immersing it in different acidic solutions. The Nernst-Planck equation was coupled with Poisson's equation for electric potential to calculate the ionic concentration inside the hydrogel which was used in the mechanical field equation to determine the swelling due to hydration [1]. In this simulation of hydrogel swelling three Partially Differential Equations were used: the Nernst-Planck equation, Poisson's equation for electric potential and mechanical field equation. Then simulations were done using COMSOL Multiphysics and employing moving mesh for 2-dimensional. In this model, we were able to simulate the distributions of different ions, as well as the electric potential which were used to calculate the mechanical deformation.
Governing Equations
In hydrogel simulation, various equations are coupled in order to describe the chemo-electro-mechanical behavior of the hydrogel in response to buffer solution pH surrounding it. The se equations are Nernst Planck, Poisson and mechanical deformation equations.
Nernst Planck Equation
The Nernst-Planck equation defines the relation between the concentrations of the various mobile species in the buffer solution. Applying continuity equation, the change in concentration flux with respect to space is equated with the rate of change of concentration which is given by:
If we modify this equation to include flux due to diffusion of ionswhich is mainly due to two factors: diffusion due to concentration gradient and migration flux because of electric potential and hence it can be modified as:
Since we are considering steady state conditions, we can neglect the concentration gradient with respect to time because there is no change with time, hence the concentration flux is written as:
(1)
In this equation, the first term represents the diffusive flux due to the concentration gradient and the second term which is coupled with the Poisson's equation is the migration flux which is due to the electric potential gradient.
Where Di, ci,zi,F,R,T and ψ are the diffusion co-efficient of the ithion, the concentration of the ithion, valence of the ith ion, the faraday constant, the universal gas constant, temperature and the electric potential respectively.
Since we are considering steady state, theNernst-Planck equation is written as:
Where μi is the mobility, which is given by:
Poisson's Equation
The Poisson's equation is used to satisfy the electroneutrality condition. It is given by the equation:
(2)
Where Ï is the charge density in the hydrogel, is the dielectric constant of the vacuum and is the relative dielectric constant of the solvent.
And zf is the fixed charge ion valence and cf is the fixed ion charge concentration. When simulating the hydrogel swelling due to changes in pH the hydrogel concentration was got by coupling the concentration flux and Poisson's equationswith cfwhich is given by the equation:
Where,K,cH and H are the ionizable charge concentration, dissociation constant, hydrogen ion concentration and hydration, respectively. The hydration state of the hydrogel is the ratio of the volume of the fluid to the volume of the solid in the gel [3].
In this equation cf is represented as a function of the surrounding pH, where both H and cH are used in defining changing ionic conditions within the hydrogel which is due to the mobile ions diffusing into the hydrogel. Due to the change in the fixed charge concentration with changes in pH, the cf is updated in every iteration.
Mechanical Field Equation
For this model, in order to understand swelling of the hydrogel, the above equation was used:
Where Ï is the effective density of the gel, u represents the vector of the xdisplacements,f represents the viscous damping parameter between the solvent and the polymer-network, σ is the stress tensor and B is the vector of the body forces. Since there are no body forces present in this simulation, we can neglect frictional factors and hence the equation can be written as:
The chemical field can be modeled using a 1-dimensional hydrogel simulation, a 2-dimensional model is necessary to model the mechanical field equations because the hydrogel does not swell in only one direction but in two directions.
From Biot's theory there is an elastic property that restores the changes in the internal pressure which is due to osmotic pressure, which is given as:
Where [C],E and I are the material elasticity matrix, Green strain tensor and the identity matrix respectively.
After solving NP and the Poisson's equations to find the concentration of both the fixed and the mobile ions are known which is used to calculate the osmotic pressure by the expression:
Where n is the number of ions, ci is the concentration of the ith ion in the hydrogel and is the ion concentration in the hydrogel at stress-free state for the ith ion.
Boundary conditions for the pH sensitive hydrogel
In this simulation, a hydrogel radius of 300µm was used in the simulation. Only a quarter of the hydrogel was used since it is a circular hydrogel. The following are the boundary conditions used in the simulation.
Subdomain 1:
Material: Hydrogel
Physics Equation
Chemical Diffusion Nernst-Planck
Electrostatics Poisson's equation
Swelling Mechanical field with moving mesh
Subdomain 2:
Material: Buffer
Physics Equation
Chemical Diffusion Nernst Planck
Moving boundaries Mechanical field equation, moving mesh frame
Boundary Type: Insulation/symmetric
The moving mesh frame equilibrium equations are:
Boundary
Type: Interface between the hydrogel and buffer
Boundary
Type: Buffer far-field
Simulation
The chemo-electro-mechanical behavior of the hydrogels was simulated in response to the pH of the buffer solution surrounding it using COMSOL 3.5a. The following modules were used:
Nernst-Planck without electro-neutrality (Chemical Engineering Module)
Conductive Media DC(AC/DC Module) for Poisson's Equation
Plane Strain (Structural Mechanics Module) for Mechanical Field Equation
Moving Mesh (ALE)
In this simulation two frames were used, namely: the fixed frame and the moving mesh frame. The chemical diffusion and electrostatic physics are considered in the moving mesh to evaluate swelling at different conditions, while the mechanical equilibrium physics is considered in the fixed frame with large deformation to calculate the hydrogel changes in size with change in pH.
Figure 1:Flow chart of the algorithm used to solve the hydrogel response to pH variation in steady state.
The algorithm used for pH simulation is shown in the figure 1.The fixed charge density together with the Poisson's dependent variable(ψ) were used to couple the NP and Poisson's equations.The fixed charge density has two parameters,hydration (H) and hydrogen ion concentration (CH), while the hydrogen ion concentration is calculated from the NP equation for steady state.
The mechanical field equation uses the osmotic pressure to in the calculation of the displacement, which is found from the moving mesh in the x and y direction which represent the displacements. Due to the change in size at every iteration, the hydration too is also updated at every iteration from the moving mesh.
The hydrogel simulationis carried out using COMSOL 3.5a.The Na+,Cl- and H+are considered in the simulationconsisting of 10 dependent variables.In this simulation, three dependent variables from Nernst-Planck equation (CNa,CCl,CH),the electric potential (ψ) from the Poisson's equation, the displacements (u, v) from the mechanical field equation, and the x and y coordinates (X & Y) and two weak constraint variables from the moving mesh module [1]. The whole domain consisted of 1810 mesh elements with 22742 degrees of freedom.
pH Sensitive Response
The diameter of the hydrogel in this simulation is fixed at 300µm. The hydrogel is assumed to be immersed in the buffer solution and electroneutrality condition was satisfied by Poisson's equation with the hydrogel taken as an isotropic material. The pH was then varied from 1-12 with a step size of 0.025 with the error convergence criterion was fixed at 1x10-4.The modulus of elasticity is 0.29 Mpa for pH<5.5 and 0.23 for pH>7.5, with a linear variation profile assumed between these two pH values. A Poisson's ratio of 0.43 was assumed for the entire range of pH. The pH simulation for various pH values is solved using three fully coupled partially differential equations with static equilibrium as the final convergence criterion. The solver used was the stationary Direct-PARADISO linear system solver and Newton iterative methods were used since the equations were highly non-linear.
RESULTS AND DISCUSSIONS:
Effect of fixed charge density on the swelling of hydrogel:
Figure 3illustrates the effect of fixed charge density on hydrogel swelling with change in pH. Figure.3 shows the effect of fixed charge density on the hydrogel swelling at various pH values. It is observed that as the fixed charge density at dry-state strongly affects the swelling equilibrium of the hydrogels at high pH values whereas a decrease in the fixed charge density from the initial fixed-charge concentration, cmos will dramatically decrease the degree of swelling at high pH. As the pH increases, diffusion of mobile ions from the buffer solution to the hydrogel is promoted hence the increase in hydrogel size with increase in pH till pH 8, where all the fixed charge sites have been occupied by the mobile ions and hence the hydrogel is said to be saturated. As the fixed charge density is increased, the hydrogel expansion increases, which can be explained by the fact that as it increases the availability of fixed charges for the mobile ions to associate with increases and hence the increase in hydrogel swelling.
Effect of buffer solution strength on hydrogel swelling:
Figure 4illustrates the effect of varying buffer solution strength on the hydrogel swelling.
Figure 4 demonstrates the effect of varying buffer solution concentration on swelling response of the poly-HEMA hydrogels,with a fixed dry-state fixed-charge concentration and Young's modulus.The highest curve represents 50mM ionic solution,whereas subsequent lower curves represent bath solutions with higher ionic strength, which can be explained in terms of the available mobile ions in the buffer solution. When the buffer solution concentration is 50mM, The mobile ions are more than when the concentration is 300mM.Since the difference in the concentration between the ith ion inside the hydrogel and the concentration of the ith ion in the buffer solution is used to calculate osmotic pressure, when the buffer solution concentration is 50mM the osmotic pressure will be greater than when the concentration I 300mM.Thus the hydrogel expands more in 50mM than in subsequent higher buffer solution concentration.
Effect of pKa on the swelling of hydrogel:
Figure 5 illustrates the effect of pKa on hydrogel swelling with varying pH.
From figure 5, it can be seen that as the dissociation constant increases, the time taken for the mobile ions to associate with the fixed charge ions decreases. This can be explained from the equation for fixed charge concentration. Theoretically, hydrogels with higher dissociation constants take long to reach the maximum swelling than hydrogels with lower dissociation constants. This can be attributed to the fact that hydrogels with high dissociation constants have mobile ions associating with the fixed charges slower than hydrogels with lower dissociation constants. Dissociation constant is used in the calculation of the fixed charge concentration.
Effect of Young's modulus on the hydrogel swelling:
Figure 6 illustrates the effect of Young's modulus on hydrogel swelling,with the modulus being a function ofpH.
Theoretically,hydrogels with high Young's modulus have decreased swelling at higher pH solutions. This can be attributed to the fact that hydrogels with high Young's modulus has low strain and hence decreased swelling. Young's modulus is used in the calculation of the elastic force that is due to the osmotic pressure change between the hydrogel and the solution, thus the degree of swelling of the hydrogel is highly dependent on the Young's modulus as evident from fig.6 ,whereas if the Young's modulus is increased the hydrogel swelling decreases.
Figure 7 illustrates comparison between simulation results and experiment date results.
Finally,the simulation results and the experimental results were compared.Figure 7shows the comparison of the plots between the experimental values and the simulation values .It should be noted that the simulation was done at 300 µm diameter HEMA hydrogel to match with the experimental gel dimensions,with the buffer concentration fixed at 300 mM [2, 4].
Conclusions:
