91桃色

Modeling Development and Numerical Simulation of Transient Nonlinear Behaviors of Electric sensitive Hydrogel Membrane under an External Electric Field

Zhen Yuan1,2* and Hua Li3
1Center for Strategic Communication, Arizona State University, POB 871205, Tempe, AZ 85287-1205, USA
2Department of Biomedical Engineering, University of Florida, Gainesville, FL 32611, USA
3School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore
Corresponding Author :	Zhen Yuan Department of Biomedical Engineering University of Florida Gainesville, FL 32611, USA E-mail: 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 Electric-sensitive Hydrogel Membrane under an External Electric Field. J Biochips Tiss Chips 3:103. doi:10.4172/2153-0777.1000103
Copyright: © 2013 Yuan Z, et al. This is an open-access 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

Visit for more related articles at Journal of Bioengineering and Bioelectronics

Abstract

A multi-physics model is developed to predict the transient nonlinear behavior of electric-sensitive hydrogel membrane, based on a multi-phasic mixture theory. In the developed model involving chemo-electro-mechanics, the transient convection-diffusion 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 convection-diffusion problem and hydrogel deformation, iteratively. Unlike the mesh-based
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 electric-sensitive 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 Bio-MEMS [1-3]. Their reversible volume changes can be induced by external bio-stimuli 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 electric-stimuli responsive hydrogel membranes with fixed-charge groups can bend reversibly, when they are immersed into a bathing-solution (Figures 1 and 2). A number of experimental studies and numerical simulations were performed to investigate the electric-sensitive 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 mechano-electrochemical 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 electric-sensitive hydrogels, with large-deformation. Nemat-Nasser and Li [14] developed an electromechanical model for ionic polymer metal composite. Accordingly, Neubrand [15] and Grimshaw et al. [16] developed electro-chemical model and electro-mechanical model for ion-exchange membrane and gel, respectively. Recently, Wallmerperger et al. [17] and Zhou et al. [18] made further studies of electric-sensitive hydrogel membranes. However, the multi-physic characteristics of the hydrogels are yet not fully understood, and now under extensive investigations.

The multi-phasic mixture theory was early developed by Lai et al. [19-21] or the swelling and deformation behaviors of articular cartilage. Based on this theory, a new multi-physics model is developed, where the governing equations are composed of the continuity equation, describing the transient deformation of solid phase, the transient convection-diffusion 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 Hermite-Cloud 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 chemo-electro-mechanical coupled behaviors of the smart hydrogel membranes.

Model Development of Electro-stimulus-sensitive Hydrogel Membrane

Multi-phasic mixture theory

The classical triphasic and multi-phasic mixture models [19-22], developed for simulation of charged hydrated biological soft tissues is extended to study the transient nonlinear behaviors of charged hydrogel membrane. Based on the multi-phasic mixture theory, this investigation conducts a novel mathematical model to investigate the multi-physics behaviors of the porous membrane applied by an external electric field. The main governing equations of the multi-phasic mixture model are given by Lai et al. [19-21] 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 [20-23]

(4)

The momentum equations, without the effects of body and the inertial forces

Mixture phase (5)

Water phase (6)

The kth ion phase (7)

The electro-neutrality condition

(8)

The constitutive equations

(9)

(10)

(11)

The parameters definitions of the multi-phasic theory are indicated in table 1.

Model development based on modified multi-phasic mixture theory

As described above, the classical triphasic and multi-phasic mixture models are unable to directly calculate the distribution of electrical potentials due to the use of electro-neutrality 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 multi-phasic 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 [19-22], 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 us is the displacement of the solid phase and Bw is set to zero, since its value is very small [23,24].

In order to incorporate the effects of migration and diffusion, the Nernet-Planck type of convection-diffusion 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 convection-diffusion 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 hydrogels-solution interfaces,

(20)

(21)

where (k=1,2,…N) is the external solution concentration at the ends of bathing-solution, is the kth ionic concentration, within the hydrogels near the interfaces and , the kth 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 hydrogel-solution 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 two-level 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 fixed-charge concentration and water volume fraction can also be computed iteratively. An implicit numerical scheme is conducted to solve the transient nonlinear convection-diffusion equations (15) and (16). A recently developed meshless approach-Hermite-Cloud 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 Hermite-Cloud Method (HCM)

In this section, meshless Hermite-Cloud 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) [25-27], 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 kth ionic concentration ck , electric potential ψ, fluid pressure p and solid-phase displacement u, the approximation can be constructed by the meshless HCM,

(22)

Where NT 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 one-dimensional or two-dimensional 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 y-component, Δx and Δy denote the cloud sizes of the fixed kernel point in the x- and y-directions, 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 first-order 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 first-order differential functions , by the equation (29) and the auxiliary conditions equations (30) and (31).

Developed numerical scheme for convection-diffusion problems

An implicit numerical scheme is developed to solve the transient nonlinear convection-diffusion equations (15) and (16), and a Newton iteration process is also implemented for each time step. In consideration of the general case, the convection-diffusion 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,c1, c2 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 2-dimensional 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 one-dimensional numerical simulations are implemented for the transient membrane deformation in the thickness h direction of strip-like electric-sensitive 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 e1 and e2 indicate the swelling strains at the two ends of hydrogels thickness. To investigate the transient nonlinear behaviors of electric-stimulus responsive hydrogels, one-dimensional 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×105(Pa), Φ=1, γk=1 (k=1,2,…N), T=293(K), R=8.314(J/mol K), D=10-7 (m2/s), (Ns/m4), c/mol.

In this validation, a time-constant external electric voltage Ve=0.1(V) is applied, c0f=2(mM), Zf=−1, h=5×10−3 (m), =1(mM), and no effect of mechanical deformation is considered. The steady-state 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 electric-potential 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), c0f=2(mM), Zf=−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 [8-10]. 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 fixed-charge concentration c0f on the transient average bending curvature Ka, under an externally applied electric field, where Zf=−1, h=5×10−3(m), =1(mM), Ve=0.1 and c0f =0.5, 1(mM), respectively. It is concluded that for a given applied voltage Ve, the average curvature Ka increases with increasing fixed-charge concentration. The observation is examined by the experiment [8-10]. Furthermore, the influence of bathing-solution concentration on the transient average curvature Ka subject to an externally applied electric field is discussed, where c0f=2 (Mm), Zf=−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 [8-10].

Numerical simulation for a steady-state 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×105 (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 Zf=+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 multi-physics mathematical model is developed for simulation of the transient nonlinear behaviors of electric-sensitive 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. One-dimensional 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

Shahinpoor M (1994) Kim B, Peppas NA (2002) Brock D, Lee W, Segalman D, Witkkowski W (1994) Tanaka T, Nishio I, Sun ST, Ueno-Nishio S (1982) Osada Y, Okuzaki H, Hori H (1992) Kim SY, Shin HS, Lee YM, Jeong CN (1999) Homma M, Seida Y, Nakano Y (2000) Homma M, Seida Y, Nakano Y (2001) Fei J, Zhang Z, Gu L (2002) Sun S, Mak AFT (2001) Doi M, Matsumoto M, Hirose Y (1992) Shiga T, Hirose Y, Okada A, Kurauchi T (1993) Shahinpoor M (1995) Nemat-Nasser S, Li JY (2000) Neubrand W (1999) Modellbildung und Simulation von Elektronen membranverfahren. PhD thesis, Universitat Stuttgart. Grimshaw PE, Nussbaum JH, Grodzinsky AJ, Yarmush ML (1990) Wallmersperger T, Kroeplin B, Holdenried J, Guelch RW (2001) Zhou X, Hon YC, Sun S, Mak AFT (2002) Lai WM, Hou JS, Mow VC (1991) Mow VC, Ateshian GA, Lai WM, Gu WY (1998) Lai WM, Mow VC, Sun DD, Ateshian GA (2000) Hon YC, Lu MW, Xue WM, Zhou X (1999) Li H, Ng TY, Cheng JQ, Lam KY (2003) Sun DN, Gu WY, Guo XE, Lai WM, Mow VC (1999) Liu WK, Jun S, Zhang YF (1995) Liu WK, Jun S, Li S, Adee J, Belytschko T (1995) Liu WK, Chen Y, Chang CT, Belytschko T (1996) Tannehill JC, Anderson DDA, Pletcher RH (1997) Bar-Cohen Y (2001) Li H, Yuan Z, Lam KY, Lee HP, Chen J, et al. (2004)

Tables and Figures at a glance

Table 1

Figures at a glance

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