Zhen Yuan^{1,2*} and Hua Li^{3}  
^{1}Center for Strategic Communication, Arizona State University, POB 871205, Tempe, AZ 852871205, USA  
^{2}Department of Biomedical Engineering, University of Florida, Gainesville, FL 32611, USA  
^{3}School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore  
Corresponding Author :  Zhen Yuan Department of Biomedical Engineering University of Florida Gainesville, FL 32611, USA Email: yzhen.star@gmail.com 

Received November 20, 2012; Accepted December 07, 2012; Published December 14, 2012  
Citation: Yuan Z, Li H (2013) Modeling Development and Numerical Simulation of Transient Nonlinear Behaviors of Electricsensitive Hydrogel Membrane under an External Electric Field. J Biochips Tiss Chips 3:103. doi:10.4172/21530777.1000103  
Copyright: © 2013 Yuan Z, et al. This is an openaccess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.  
Related article at Pubmed Scholar Google 
Visit for more related articles at Journal of Bioengineering and Bioelectronics
A multiphysics model is developed to predict the transient nonlinear behavior of electricsensitive hydrogel membrane, based on a multiphasic mixture theory. In the developed model involving chemoelectromechanics, the transient convectiondiffusion equations for ionic concentrations incorporate the migration and diffusion terms; the Poisson equation is employed to compute the distribution of electric potential directly, and the transient hydrogel deformation is implemented easily by the continuity and momentum equations. To solve the present mathematical model consisting of transient nonlinear partial differential governing equations, a true meshfree, implicit numerical scheme is conducted for solution of convectiondiffusion problem and hydrogel deformation, iteratively. Unlike the meshbased
methods, the employed meshless Hermite Cloud Method (HCM) uses a fixed reproducing kernel approximation for construction of the interpolation functions, and employs the point collocation technique for discretization of partial differential boundary values and initial value problems. The transient responds of electricsensitive hydrogels, including the membrane deformation, ionic concentrations and electric potentials, interior and exterior the membranes are numerically simulated, and the parameters having important influence on the transient hydrogel deformation are also investigated.
Keywords  
Hydrogels; Membrane; Multiphasic mixture theory; Mathematical model; Diffusion; Convection; BioMEMS; Meshless method; Biomechanics  
Introduction  
Over the past decades, the actuators/sensors based on stimuliresponsive polymer hydrogels have attracted much attention for widerange biological applications, such as artificial muscles and BioMEMS [13]. Their reversible volume changes can be induced by external biostimuli including pH, light, temperature and electric field. Usually, hydrogel membranes are composed of the solid, interstitial water and ion phases. If an external electric field is applied, the electricstimuli responsive hydrogel membranes with fixedcharge groups can bend reversibly, when they are immersed into a bathingsolution (Figures 1 and 2). A number of experimental studies and numerical simulations were performed to investigate the electricsensitive behaviors of hydrogels which covered the preparation of such materials, studies of deformation mechanism and design of actuators. Tanaka et al. [4] reported a gel collapsed in an acetone/water mixture applied by electric field. Osada et al. [5] reported a polymer gel with electrically driven motility, and Kim et al. [6] studied the electric sensitive behavior of IPN hydrogel. Homma et al. [7,8] and Fei et al. [9] discussed the factors having important effects on the swelling deformation of electricsensitive hydrogels. Sun and Mak [10] studied the mechanoelectrochemical behavior of chitosan composite fibers by experiment. Doi et al. [11], Shiga et al. [12] and Shahinpoor [13] investigated the dynamics of ionic polymer gels, subject to an external electric field. Brock et al. [3] studied the dynamic model for electricsensitive hydrogels, with largedeformation. NematNasser and Li [14] developed an electromechanical model for ionic polymer metal composite. Accordingly, Neubrand [15] and Grimshaw et al. [16] developed electrochemical model and electromechanical model for ionexchange membrane and gel, respectively. Recently, Wallmerperger et al. [17] and Zhou et al. [18] made further studies of electricsensitive hydrogel membranes. However, the multiphysic characteristics of the hydrogels are yet not fully understood, and now under extensive investigations.  
The multiphasic mixture theory was early developed by Lai et al. [1921] or the swelling and deformation behaviors of articular cartilage. Based on this theory, a new multiphysics model is developed, where the governing equations are composed of the continuity equation, describing the transient deformation of solid phase, the transient convectiondiffusion equations computing the diffusive ionic concentrations, the Poisson equation calculating the electric field, and the momentum equation describing the mechanical field. To solve the governing equations with remeshing requirement in hierarchical iteration procedure, a true meshfree, implicit numerical scheme is conducted. The meshless HermiteCloud Method (HCM) based on a fixed kernel approximation is used for spatial discretization, in which the interpolation functions are constructed according to a set of points scattered in problem domain. The point collocation technique is employed in the problem domain for discretization of the governing equations, boundary and initial conditions. By the meshless implicit numerical scheme, the transient nonlinear partial differential governing equations are solved to simulate the distributions of diffusive ionic concentrations and electric potentials, interior and exterior the hydrogel membranes, as well as the swelling and bending deformations of the hydrogels. Moreover, several parameter having important influences on the transient hydrogels deformations are also investigated to enhance the understanding of the chemoelectromechanical coupled behaviors of the smart hydrogel membranes.  
Model Development of Electrostimulussensitive Hydrogel Membrane  
Multiphasic mixture theory  
The classical triphasic and multiphasic mixture models [1922], developed for simulation of charged hydrated biological soft tissues is extended to study the transient nonlinear behaviors of charged hydrogel membrane. Based on the multiphasic mixture theory, this investigation conducts a novel mathematical model to investigate the multiphysics behaviors of the porous membrane applied by an external electric field. The main governing equations of the multiphasic mixture model are given by Lai et al. [1921] and Hon et al. [22].  
The saturation condition,  
(1)  
The continuity equations  
(2)  
The conservation of fixed charge groups attached on the polymer network  
(3)  
The continuity equation of mixture phase by neglecting , when compared with and [2023]  
(4)  
The momentum equations, without the effects of body and the inertial forces  
Mixture phase (5)  
Water phase (6)  
The k^{th} ion phase (7) 

The electroneutrality condition  
(8)  
The constitutive equations  
(9)  
(10)  
(11)  
The parameters definitions of the multiphasic theory are indicated in table 1.  
Model development based on modified multiphasic mixture theory  
As described above, the classical triphasic and multiphasic mixture models are unable to directly calculate the distribution of electrical potentials due to the use of electroneutrality condition. The migration and diffusion terms aren’t considered in the ionic continuity equations, and the solving domain only covers the hydrogel membranes. To overcome the drawbacks mentioned, a novel mathematical model based on modified multiphasic mixture theories is developed, and the governing equations are derived as follows.  
By summarizing Equations (6) to (7) for all N ionic species, and neglecting and , when compared with [1922], we have 

(12)  
Substituting the constitutive equations (10)(11) into equation (12), and assuming the osmotic coefficients of all ionic species equal to Φ, the equation (12) is rewritten,  
(13)  
It is noted, the continuity equation (4) can be rewritten as  
(14)  
Substituting Equation (13) into Equation (14), the continuity equation of the mixture phase is obtained,  
(15)  
Where u^{s} is the displacement of the solid phase and B_{w} is set to zero, since its value is very small [23,24].  
In order to incorporate the effects of migration and diffusion, the NernetPlanck type of convectiondiffusion equations for each diffusive ionic flux, replacing the continuity equation (2) of the ionic phase, is employed,  
(16)  
Moreover, as the velocity of propagation in the electric field is much higher than the one occurring in the convectiondiffusion equation, the Poisson equation is given by,  
(17)  
where is the fluid velocity relative to the polymer network, and can be computed directly by the equation (13), diffusive coefficient, the source term resulting from the chemical conversation of the molecules, ε the dielectric constants and ε_{0} the permittivity of the free space. As such, the developed model is composed of the momentum equation (5), the continuity equations (15) and (16) and the Poisson equation (17) with unknown variables, p, , ψ and (k=1, 2, … N). The boundary conditions consist of Dirichlet boundary conditions at the ends of bathing solution,  
(18)  
(19)  
And the interfacial conditions at hydrogelssolution interfaces,  
(20)  
(21)  
where (k=1,2,…N) is the external solution concentration at the ends of bathingsolution, is the k^{th} ionic concentration, within the hydrogels near the interfaces and , the k^{th} ionic concentrations of bathing solution near the interfaces, Ve is the applied electric voltage, and represents the osmotic pressure at reference configuration. The interfacial conditions above are based on the assumptions that, for each time step on the hydrogelsolution interface, the chemical potentials of hydrogel membrane are equal to those of bathing solution [24], and the total mixture stress is equal to zero. A twolevel hierarchical iteration technique is carried out for each time step, to solve the coupled nonlinear partial differential governing equations. The inner iterations are used for the computation of the diffusive ionic concentrations and electric potential simultaneously. Substitute the implemented numerical results into the subsequent outer iteration loop for calculation of the displacements and pressure. The fixedcharge concentration and water volume fraction can also be computed iteratively. An implicit numerical scheme is conducted to solve the transient nonlinear convectiondiffusion equations (15) and (16). A recently developed meshless approachHermiteCloud Method (HCM) [23], is employed for the spatial discretization of the governing equations.  
A Numerical Scheme for Developed Mathematical Model based on Meshless HCM  
Meshless HermiteCloud Method (HCM)  
In this section, meshless HermiteCloud Method (HCM) [23], is introduced for the spatial discretization of the above nonlinear partial differential initial value problems. As an extension of the classical Reproducing Kernel Particle Method (RKPM) [2527], this meshless approach employs Hermite theorem for the construction of the interpolation functions, and uses the point collocation technique for discretization of Partial Differential equations with Boundary Value (PDBV) problems.  
For any unknown real function , representing here the k^{th} ionic concentration c^{k} , electric potential ψ, fluid pressure p and solidphase displacement u, the approximation can be constructed by the meshless HCM,  
(22)  
Where N_{T} and are total numbers of discrete points covering the computational domain, and are the shape functions defined as, 

(23)  
In which is the kernel function defined as,  
(24)  
is the linearly independent basis. For example, in onedimensional or twodimensional quadratic PDBV problem, it is given respectively.  
(25)  
(26)  
And is a symmetric matrix associated with the fixed kernel centered point , and is expressed as  
(27)  
In the equation (25), is called the window function, and defined here as a cubic spline function,  
(28)  
Where , or for the x or ycomponent, Δx and Δy denote the cloud sizes of the fixed kernel point in the x and ydirections, respectively.  
Assuming , the discrete approximations of may be given by,  
(29)  
It is noted that, compared with the shape functions in the equation (26) constructed on β–order basis, the present shape functions are constructed on (β1)order basis.  
The auxiliary conditions are required for the additional unknown functions . By imposing the firstorder partial differentials with respect to the variables x and y on the approximation expressed by the equation (22), and using the equation (29), we have the following auxiliary conditions 

(30)  
(31)  
The meshless method employs the Hermite interpolation theorem to construct the approximate unknown function , in form of the equation (22), and couples the firstorder differential functions , by the equation (29) and the auxiliary conditions equations (30) and (31). 

Developed numerical scheme for convectiondiffusion problems  
An implicit numerical scheme is developed to solve the transient nonlinear convectiondiffusion equations (15) and (16), and a Newton iteration process is also implemented for each time step. In consideration of the general case, the convectiondiffusion equation is given by  
(32)  
Together with the general boundary and initial conditions  
(33)  
(34)  
In which k is the diffusion coefficient, v is the convection coefficient,c_{1}, c_{2} are known constants, and f(x,t) and are specified known functions. According to θweighted numerical scheme (θ≥ 0.5) [28], equation (31) is written,  
(35)  
If we define , respectively, equation (35) can be given by,  
(36)  
Where . It should be noted if a 2dimensional matrix A is specified,  
(37)  
Equation (36) can be rewritten,  
(38)  
Where  
are the nodes in the inner computational domain and on the boundary, respectively.  
In which  
(39)  
And  
(40)  
Numerical Validations  
To validate the presently developed model, several onedimensional numerical simulations are implemented for the transient membrane deformation in the thickness h direction of striplike electricsensitive hydrogels, immersed into bathing solution under an external electric field (Figure 2), where it is assumed that the diffusive coefficients (k=1, 2, … N), and the strain here is isotropic, i.e.  
(41)  
Where is the strain in thickness direction. An approximately average curvature Ka, resulting from the pressure difference between the two ends of the thickness h, is defined geometrically at the middle point of the thickness [18],  
(42)  
Where e_{1} and e_{2} indicate the swelling strains at the two ends of hydrogels thickness. To investigate the transient nonlinear behaviors of electricstimulus responsive hydrogels, onedimensional numerical simulations are carried out for a hydrogel strip immersed into NaCl solution applied by an externally applied electric field. In the solution domain x (m) ∈[0, 0.005] ∩ [0.010, 0.015], the Dirichlet boundary conditions for ionic concentrations and electrical potential at anode and cathode are imposed at x=0 and 15×10^{3}(m), respectively. The specified parameters are =0.8,ε_{0}=8.854×10^{12} , ε=80, (3λ+2μ)=1.2×10^{5}(Pa), Φ=1, γ_{k}=1 (k=1,2,…N), T=293(K), R=8.314(J/mol K), D=10^{7 }(m^{2}/s), (Ns/m^{4}), c/mol. 

In this validation, a timeconstant external electric voltage Ve=0.1(V) is applied, c_{0}^{f}=2(mM), Z^{f}=−1, h=5×10^{−3} (m), =1(mM), and no effect of mechanical deformation is considered. The steadystate solutions without an externally applied electric field, as shown in figures 3 and 4 are taken as initial condition for numerical simulations. The transient variations of diffusive ionic concentrations and electrical potential are numerically simulated and shown in figures 5 and 6. It is seen from figures 5 and 6 that with progressing time, the concentration of Cl^{} ion and Na^{+} ion are increasing near the gel edge, close to the cathode and decreasing near the gel edge, close to the anode side. Additionally, it is observed from figure 7 that the electric potential is also increasing on the cathode side and decreasing on the anode side of the gel. The constant distribution of concentration of fixed charges is provided in figure 8. The change of electric potential within the hydrogel is smaller than that in the exterior solution due to the higher conductivity of the mobile ions in the hydrogel strip. The interface concentration difference between the interior hydrogel membrane and exterior bathing solution on the anode side is larger than that on the cathode side. These simulated phenomena were in good agreement with the published experimental results [29]. The larger interface differences of concentrations and electric potential on the anode side result in the larger pressure, compared with those on the cathode side. This pressure difference makes the hydrogel strip have a bending deformation. Furthermore, the presently simulated electricpotential distribution agrees well with Wallmersperger’s FEM simulations [17,29].  
Figure 9 displays the transient bending deformation of the hydrogel strip subject to different electric voltages, where h=5×10^{−3}(m), =1(mM), c_{0}^{f}=2(mM), Z^{f}=−1. It is seen from figure 9 that the bending curvature versus time became steeper and steeper with the increase of electric voltage, till leveled off at a steady state. The bending curvature increases with the increase of the applied voltage, which is in good agreement with the experimental phenomena [810]. This can be explained by the fact that there is an enhancement in the transport rate of counterions in the domain, as the increase of electric voltage. A further study indicates the electrical voltage is proportional to the bending curvature, as shown in figure 10.  
Figure 11 demonstrates the influence of fixedcharge concentration c_{0}^{f} on the transient average bending curvature Ka, under an externally applied electric field, where Z^{f}=−1, h=5×10^{−3}(m), =1(mM), Ve=0.1 and c_{0}^{f} =0.5, 1(mM), respectively. It is concluded that for a given applied voltage Ve, the average curvature Ka increases with increasing fixedcharge concentration. The observation is examined by the experiment [810]. Furthermore, the influence of bathingsolution concentration on the transient average curvature Ka subject to an externally applied electric field is discussed, where c_{0}^{f}=2 (Mm), Z^{f}=−1, h=5×10^{−3}(m), Ve=0.1 and 0.5, 1(mM). Figure 12 indicates the bending deformation of hydrogel strip decreases with the increase of , where the simulations also agree well with the experimental phenomena [810].  
Numerical simulation for a steadystate case is conducted to investigate the hydrogel membrane deformation and compared with experimental data [18,30]. The parameters used are =0.8, ε_{0}=8.854×10^{12} , ε=80, (3λ+2μ)=1.2×10^{5} (Pa), h=1.0×10^{−3}(m), (mM), =5.5 (mM), Φ=1, =1 (k=1, 2,…N), T=278(K), R=8.314(J/mol K) and Z^{f}=+1. The distance between two electrodes is 2.0×10^{−2}(m), and only univalent ions are considered. It is seen from figure 13 that the average curvatures Ka increase with the increasing external electric voltage, and the experimental results are in good agreement with numerical ones. 

Conclusions  
In this study, a multiphysics mathematical model is developed for simulation of the transient nonlinear behaviors of electricsensitive hydrogel membranes immersed into bathing solution, subject to an externally applied electric field. The formulation is capable to describe the bending deformation of the hydrogel, the distributions of diffusive ionic concentrations, and electric potential in both the interior membrane and exterior solution. Onedimensional numerical simulations are carried out for a hydrogel strip. The numerically simulated results are in good agreements with experimental data and published FEM solutions, which validates the presently developed models.  
References  

Table 1 
Figure 1  Figure 2  Figure 3  Figure 4  Figure 5  
Figure 6  Figure 7  Figure 8  Figure 9  Figure 10 
Figure 11  Figure 12  Figure 13 
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals