Modeling Development and Numerical Simulation of Transient ...

Yuan and Li, Biochip Tissue Chip 2013, 3:1 http://dx.doi.org/10.4172/2153-0777.1000103

Biochips & Tissue Chips Research Article

Open Access

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 Center for Strategic Communication, Arizona State University, POB 871205, Tempe, AZ 85287-1205, USA Department of Biomedical Engineering, University of Florida, Gainesville, FL 32611, USA 3 School of Mechanical and Aerospace Engineering, Nanyang Technological University, Singapore 1 2

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. J Biochip Tissue Chip ISSN: 2153-0777 JBTC, an open access journal

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

*Corresponding author: Zhen Yuan, Department of Biomedical Engineering, University of Florida, Gainesville, FL 32611, USA, E-mail: [email protected] 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/21530777.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.

Volume 3 • Issue 1 • 1000103

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

Page 2 of 13

Mobile Ion

Polymer Chain

Fluid Filled Region

Crosslink Undissociated

Fixed Charge

Ionizable Group Figure 1: Micro view of the hydrogel structure.

y(m) 0.0075 solution

0

0.015 x(m)

hydrogel

−0.0075 Figure 2: Bending deformation of a hydrogel strip in bathing solution under external electric field.

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, N

φ s + φ w + ∑ φ k = 1

(1)

∂ρ α + ∇ ⋅ ( ρ α v α ) = 0 (α=s,w,k) and ρ α = ρ Tα φ α (α=w,s), ρ k = ρTk φ k = c k M k φ w ∂t

(2)

k =1

The continuity equations

The conservation of fixed charge groups attached on the polymer network

∂c f + ∇ ⋅ (c f v s ) = 0 , ∂t

J Biochip Tissue chip ISSN: 2153-0777 JBTC, an open access journal

and

c f = c0f / (1 + tr (E) / ϕ0w )

(3)

Volume 3 • Issue 1 • 1000103


Page 3 of 13 The continuity equation of mixture phase by neglecting φ

k

, when compared with

φs

and

φ w [20-23]

∇ ⋅ (φ s v s + φ w v w ) =0 and ϕ w = 1 − ϕ0s / (1 + tr (E))

(4)

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

∇⋅ = 0 − ρ w ∇µ w + f ws ( v s − v w ) +

Water phase

N

∑f

wk ( v

k

− v w ) = 0

k =1

− ρ k ∇µ k + f ks ( v s − v k ) + f kw ( v w − v k ) +

The kth ion phase

N

∑

f kj ( v j − v k ) = 0

(k=1, 2, 3, … N)

(5) (6) (7)

j =1( j ≠ k )

The electro-neutrality condition N

z f c f + ∑ z k c k = 0

(8)

k =1

The constitutive equations

ó = − pI − Tc I + λstr (E)I + 2 µs E µ w = µ0w + [ p − RT

N

∑Φ c

k k

+ BwtrE]/ ρTw

(9) (10)

k =1

µ k = µ0k + (RT / M k )ln(γ k c k ) + z k Fcψ / M k

(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

f ws ( v s − v w ) = ρ w ∇µ w +

N

∑ ρ ∇µ k

k

f ks and f kj , when compared with f kw [19-22], we have

(12)

k =1 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, N N s w w k (13) ( v − v ) = (ϕ / f ws )[∇p − RT (Φ − 1)∇( c ) + Bw ∇tr (E) + Fc z k c k ∇ψ ] k =1 k =1 It is noted, the continuity equation (4) can be rewritten as

∑

∑

∇ ⋅ v s + ∇ ⋅ φ w (v w − v s ) = 0

(14)

Substituting Equation (13) into Equation (14), the continuity equation of the mixture phase is obtained,

N

∑

∇ ⋅ u,st = ∇ ⋅ {((ϕ w )2 / f ws )[∇p − RT (Φ − 1)∇( Where u is the displacement of the solid phase and s


Bw

k =1

c k ) + Fc

N

∑ z c ∇ψ ]} k k

(15)

k =1

is set to zero, since its value is very small [23,24]. Volume 3 • Issue 1 • 1000103


Page 4 of 13 Quantities and Symbols

φα

ó

Definitions Volume fractions of solid network, water and ionic species

(α=s,w,k) N

Total number of ionic species

E

Elastic strain vector of solid phase

t

Time

φ 0s

Solid-phase volume fraction at a reference configuration

ck

Concentration of the kth ionic species

Mk

Molar weight of the kth ionic species

ρTα

True mass density of phase α

vα

Velocity vector of phase α.

c0f

Fixed charge density at reference configuration.

φ 0w

Water volume fraction at reference configuration Mixture stress

µ α (α=s,w,k)

Chemical potential of phase α

fαβ = f βα

Frictional coefficients between phases α and β

z k (k=1,2,3,…N)

Fluid pressure

p I

Identity tensor Chemical expansion stress

Tc Φk

Osmotic coefficient of kth ionic species

λs, μs

µ0α

Valence of ion α

Lame coefficients of solid matrix

γk

Activity coefficient of the kth ionic species

(α=w, k)

Chemical potential of phase α at reference configuration

R T

Universal gas constant Absolute temperature

Fc

Faraday constant

ψ

Electrical potential

Bw

Coupling coefficient. Table 1: The parameter definitions of the multi-phasic mixture theory.

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,

c,kt = ∇ ⋅ (Dk ∇c k ) + (Fc / RT )z k ∇ ⋅ (c k Dk ∇ψ ) − ∇ ⋅ (c k V) + rk (c k )

(k=1,2,…N)

(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,

2

∇ ψ =−

F

εε 0

N

∑ (z c

k k

+ z f c f )

(17)

k =1


Volume 3 • Issue 1 • 1000103


Page 5 of 13 w

s

where V = (v − v ) is the fluid velocity relative to the polymer network, and can be computed directly by the equation (13), Dα diffusive coefficient, rα 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 s k equation (17) with unknown variables, p, u , ψ and c (k=1, 2, … N). The boundary conditions consist of Dirichlet boundary conditions at the ends of bathing solution,

ψ (x , t ) Anode = +0.5Ve , and ψ (x , t ) Cathode = −0.5Ve

(18)

c k ( x, t )

(19)

= c k ( x, t )

Anode

Cathode

= ck∗ ,

And the interfacial conditions at hydrogels-solution interfaces, N

∑(c

k k in −int erface ( x , t ) − cout −int erface ( x , t )) − p0 ( x , t )

(20)

pint erface (x , t )I = λstr (Eint erface ( x , t ))I + 2 µs Eint erface ( x , t )

(21)

pinterface (x , t ) = RT

k =1

k in −interface

∗ k

where c (k=1,2,…N) is the external solution concentration at the ends of bathing-solution, c is the k ionic concentration, within k the hydrogels near the interfaces and cout −interface , the kth ionic concentrations of bathing solution near the interfaces, Ve is the applied electric voltage, and p0 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. th

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.

f ( x, y ) , representing here the kth ionic concentration c k , electric potential ψ, fluid pressure p and solid~ phase displacement u, the approximation f ( x, y ) can be constructed by the meshless HCM, For any unknown real function

f (x , y ) =

NS

NT

NT

∑ N ( x , y ) f + ∑ ( x − ∑ N ( x , y )x ) M n =1

+

n

n

NS

NT

m =1 n

n

m ( x , y )Gxm

n =1

∑ ( y − ∑ N (x , y ) y )M m =1

Where NT and defined as,

n

n

n =1

N S (≤ N T )

are total numbers of discrete points covering the computational domain, and

K ( xk − pn , y k − qn )

K ( x k − p, y k − q ) =

(22)

m ( x , y )Gym

N n ( x, y ) = B ( p n , q n ) A −1 ( x k , y k ) B T ( x, y ) K ( x k − p n , y k − q n )∆S n In which

+

(n=1, 2,… N T )

N n ( x, y ) are the shape functions

is the kernel function defined as,

x − p * yk − q 1 W *( k )W ( ) ∆x ∆y (∆x∆y )


(23)

(24)

Volume 3 • Issue 1 • 1000103


Page 6 of 13

B ( p, q )

respectively.

is the linearly independent basis. For example, in one-dimensional or two-dimensional quadratic PDBV problem, it is given

B( p) = {b1( p), b2 ( p), , bβ ( p)} = {1, p, p2 }

(β=3)

(25)

B( p, q) = {b1( p, q), b2 ( p, q),, bβ ( p, q)} = {1, p, q, p2 , pq, q 2 }

(β=6)

(26)

And

A( xk , y k )

Aij (x k , yk ) =

is a symmetric matrix associated with the fixed kernel centered point ( xk , y k ) , and is expressed as

NT

∑ b ( p , q )K ( x i

n

n

k

− pn , yk − qn )b j ( pn , qn )∆Sn

( i,

j = 1,2,, β

)

(27)

2< z 1 ≤ z ≤ 2 0 ≤ z ≤1

(28)

n =1

In the equation (25),

W * ( z)

is called the window function, and defined here as a cubic spline function,

 0  W ( z ) = (2 + z ) 3 / 6  (2 / 3) − z 2 (1 − 0.5 z )  *

Where z = ( x k − p ) / ∆x , or z = ( y k − q ) / ∆y 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.

( xk , y k )

Assuming Gx ( x , y ) = ∂f ( x , y ) / ∂x and Gy ( x , y ) = ∂f ( x , y ) / ∂y , the discrete approximations of Gx ( x , y ) and may be given by, Ns Ns

 (x , y ) = Gx

∑M

m ( x , y )Gxm

m =1

It is noted that, compared with the shape functions are constructed on (β-1)-order basis.

M m ( x, y )

 (x , y ) = Gy

∑M

m ( x , y )Gym

Gy(x , y ) (29)

m =1

N n ( x, y ) in the equation (26) constructed on β–order basis, the present shape functions Gx(x , y ) and Gy(x , y ) . By imposing the first-order partial ~ f ( x, y ) expressed by the equation (22), and using the equation (29), we

The auxiliary conditions are required for the additional unknown functions differentials with respect to the variables x and y on the approximation have the following auxiliary conditions NT

∑N

n, x ( x , y ) f n

−

n =1 NT

∑N

NS

NT

∑(∑(N

n, x ( x , y )xn )) M m ( x , y )Gxm

−

m =1 n =1

n, y ( x , y ) f n

−

NS

NT

∑(∑(N

NS

NT

∑(∑(N

n, x ( x , y ) yn )) M m ( x , y )Gym

=0

(30)

n, y ( x , y )xn )) M m ( x , y )Gxm

=0

(31)

m =1 n =1

n, y ( x , y ) yn )) Mm ( x , y )Gym

−

NS

NT

∑(∑(N

n =1 m =1 n =1 m =1 n =1 The meshless method employs the Hermite interpolation theorem to construct the approximate unknown function

equation (22), and couples the first-order differential functions equations (30) and (31).

~ Gx ( x , y )

~

~ f ( x, y ) , in form of the

and Gy ( x, y ) , by the equation (29) and the auxiliary conditions

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

x ∈ Ω , t > 0

(32)

x ∈ ∂Ω , t > 0

(33)

∂u ( x, t ) / ∂t = k∇ 2 u ( x, t ) + v ⋅ ∇u ( x , t ) Together with the general boundary and initial conditions

c1u ( x, t ) + c2∂u ( x , t ) / ∂n = f ( x , t ) ,

u ( x , t ) = u 0 ( x ) , J Biochip Tissue chip ISSN: 2153-0777 JBTC, an open access journal

t = 0

(34)

Volume 3 • Issue 1 • 1000103


Page 7 of 13 In which k is the diffusion coefficient, v is the convection coefficient, c1 , c2 are known constants, and f(x,t) and functions. According to θ -weighted numerical scheme ( θ ≥ 0.5) [28], equation (31) is written,

u(x , t + ∆t ) − u(x , t ) = ∆tθ {k∇2u If we define

u n = u ( x, t n )

and t

n

t +∆t

+ v ⋅∇u

t +∆t } + ∆t (1 − θ ){k∇

2

u 0 ( x ) are specified known

u t + v ⋅∇u t }

(35)

= t n−1 + ∆t , respectively, equation (35) can be given by,

u n+1 + α∇ 2 u n+1 + β ⋅ ∇u n+1 = u n + η∇ 2 u n + ξ .∇u n Where α

= −kθ∆t , β = −θ∆tv ,η = k∆t (1 − θ )

and

(36)

ξ = ∆t (1 − θ )v . It should be noted if a 2-dimensional matrix A is specified,

NT NT   2 σ [∇2 N (x , y )]  N n ( xi , yi )xn ) M j ( xi , yi ))]N ×N + [ M j ( xi , yi )]N ×N ρ x σ [∇2 (( yi − N n (xi , yi ) yn )M j (xi , yi ))]N ×N + j i i N Ω ×NT σ [∇ (( xi − Ω Ω Ω S S S   n =1 n =1   NT NT    +[N j (xi , yi )]N ×N +[((xi − N n ( xi , yi )xn ) M j ( xi , yi ))]N ×N [(( yi − N n (xi , yi ) yn )M j (xi , yi ))]N ×N + [ M j (xi , yi )]N ×N ρ y  Ω Ω T S Ω Ω S S   n =1 n =1   c2[ M j ( xi , yi )]N ×N ix c2[c j ( xi , yi )]N ×N i y  c1[[N j (xi , yi )]]N N ×NT  N S N S   NT NT   [−( N n, x ( xi , yi )xn ) M j ( xi , yi )]N ×N [−( N n, x ( xi , yi ) yn ) M j ( xi , yi )]N ×N  [N j , x ( xi , yi )]N S ×NT  S S S S   n =1 n =1   NT NT   [−( N n, y (xi , yi )xn ) M j ( xi , yi )]N ×N [−( N n, y (xi , yi ) yn )M j (xi , yi )]N ×N  [N j , y (xi , yi )]N S ×NT  S S S S   n =1 n =1

∑

∑

∑

∑

∑

∑

∑

∑

(37)

Equation (36) can be rewritten,

Aαβ {λ }n +1 = Aηξ {λ }n + {F}n +1 Where

−

−

σ = (α ,η ) , ρ = ( β , ξ ) , {F }n+1 = {F }n+1 − {F }n

boundary, respectively.

and

(38)

N Ω , N N are the nodes in the inner computational domain and on the

In which −

{F }n = [0 . . . 0

f Ωn +1

f Ωn + 2 . . .

f Ωn + N

0 . . . 0]

T

(39)

And

{λ } = {u1 u2 . . . un Gx1 Gx2 . . . Gxm Gy1 Gy2 . . . Gym ]T

(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 Dk = D (k=1, 2, … N), and the strain here is isotropic, i.e.

tr (E) = 3e11 = 3

∂u ∂x

(41)

Where e11 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) Ka = 2(e1 − e2 ) / (h(2 + e1 + e2 )) 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 −3 ionic concentrations and electrical potential at anode and cathode are imposed at x=0 and 15×10 (m), respectively. The specified parameters are

φ0w =0.8, ε 0 =8.854×10 −12 ( C 2 / Nm2 ), ε=80, (3λ+2μ)=1.2×10 5 (Pa), Φ=1, γ k =1 (k=1,2,…N), T=293(K), R=8.314(J/mol K),

(m2/s),

f ws = 0.7 × 10−15 (Ns/m4), F = 9.6487 × 104 c/mol.

In this validation, a time-constant external electric voltage Ve=0.1(V) is applied,


c0f

=2(mM),

D = 10−7

z f =−1, h=5×10 −3 (m), ck∗ =1(mM), and Volume 3 • Issue 1 • 1000103


Page 8 of 13 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 f 0

f

−3

(m),

ck∗

=1(mM), c =2(mM), z =−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 0.000

Ionic concentration (mM)

2.6

0.002

0.004

0.006

0.008

0.010

0.012

0.014

2.6

2.4

2.4

2.2

2.2

2.0

2.0

1.8

cNa cCl

1.8

+

1.6

1.6

_

1.4

1.4

1.2

1.2

1.0

1.0

0.8

0.8

0.6

0.6

0.4 0.000

0.002

0.004

0.006

0.008

0.010

0.012

0.014

0.4

x (m) Figure 3: Distribution of ionic concentrations at steady-state without an external electric field.

Electric potential (V)

0.000

0.002

0.004

0.006

0.008

0.010

0.012

0.014

0.000

0.000

-0.005

-0.005

-0.010

-0.010

-0.015

-0.015

-0.020

-0.020

-0.025 0.000

0.002

0.004

0.006

0.008

0.010

0.012

0.014

-0.025

x (m) Figure 4: Distribution of electric potential at steady-state without an external electric field.


Volume 3 • Issue 1 • 1000103


Page 9 of 13 0.000

+

CNa (mM)

3.0

0.002

0.004

0.006

0.008

0.010

0.012

0.014

3.0

2.8

2.8

2.6

2.6

2.4

2.4

2.2

2.2

2.0

2.0

t =10 (s) t =20 (s) t =50 (s) t =100 (s) t =500 (s) t =800 (s)

1.8 1.6 1.4 1.2

1.8 1.6 1.4 1.2

1.0

1.0

0.8

0.8

0.6 0.000

0.002

0.004

0.006

0.008

0.010

0.012

0.014

0.6

x (m)

Figure 5: Distribution of Na+ concentrations under an external electric field. 0.000

1.5

0.002

0.004

0.006

0.008

0.010

0.012

0.014

1.4

1.4

1.3

1.3

1.2 1.0

_

1.2

t =10 (s) t =20 (s) t =50 (s) t =100 (s) t =500 (s) t =800 (s)

1.1

CCl (mM)

1.5

0.9 0.8

1.1 1.0 0.9 0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2 0.000

0.002

0.004

0.006

0.008

0.010

0.012

0.014

0.2

x (m)

Figure 6: Distribution of Cl- concentrations under an external electric field. 0.000

0.05

0.002

0.004

0.006

0.008

0.010

0.012

0.014

0.04

0.04

0.03

Electric potential (V)

0.05

0.03

t = 10 (s) t = 20 (s) t = 50 (s) t = 100 (s) t = 500 (s) t = 800 (s)

0.02 0.01 0.00

0.02 0.01 0.00

-0.01

-0.01

-0.02

-0.02

-0.03

-0.03

-0.04

-0.04

-0.05 0.000

0.002

0.004

0.006

0.008

0.010

0.012

0.014

-0.05

x (m) Figure 7: Distribution of electric potentials under an external electric field.


Volume 3 • Issue 1 • 1000103


Page 10 of 13 0.000

Concentration of fixed charge (mM)

2.25

0.002

0.004

0.006

0.008

0.010

0.012

0.014

2.25

2.00

2.00

1.75

1.75

1.50

1.50

1.25

1.25

1.00

1.00

0.75

0.75

0.50

0.50

0.25

0.25

0.00 0.00

0.02

0.04

0.06

0.08

0.10

0.12

0.00

0.14

x (m)

Figure 8: Distribution of concentration of fixed charges.

Curvature Ka (1/m)

0

8.0 7.5 7.0 6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0

5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 8.0 7.5 7.0

Ve= 0.02 (V) Ve= 0.10 (V) Ve= 0.20 (V)

6.5 6.0 5.5 5.0 4.5 4.0 3.5 3.0 2.5 2.0 1.5 1.0 0.5 0.0

0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 t (s)

Figure 9: Bending kinetics of hydrogels as a function of the applied electric field.

Curvature Ka (1/m)

0.00

0.05

0.10

0.15

0.20

0.25

0.30

14

14

12

12

10

10

8

8

6

6

4

4

2

2

0 0.00

0.05

0.10

0.15

0.20

0.25

0.30

0

Electric potential (V) Figure 10: Effect of externally applied electric field on the average curvature Ka.


Volume 3 • Issue 1 • 1000103


Page 11 of 13 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. f

Figure 11 demonstrates the influence of fixed-charge concentration c0 on the transient average bending curvature Ka, under an externally f −3 ∗ f applied electric field, where z =−1, h=5×10 (m), ck =1(mM), Ve=0.1 and c0 =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 ck on the transient average curvature Ka subject to an externally applied electric f f −3 ∗ field is discussed, where c0 = 2 (mM), z =−1, h=5×10 (m), Ve=0.1 and ck =0.5, 1(mM). Figure 12 indicates the bending deformation of ∗ hydrogel strip decreases with the increase of ck , 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 f w 5 −3 data [18,30]. The parameters used are φ 0 =0.8, ε 0 =8.854×10 −12 ( C 2 / Nm2 ), ε=80, (3λ+2μ)=1.2×10 (Pa), h=1.0×10 (m), c0 = 20 ∗ f −2 (mM), ck =5.5 (mM), Φ=1, γ k =1 (k=1, 2,…N), T=278(K), R=8.314(J/mol K) and z =+1. The distance between two electrodes is 2.0×10 (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 1.2

0

2

4

6

8

10

12

14

16

18

20 1.2

1.0

1.0 f

Curvature Ka (1/m)

co = 2.0(mM) f

c0 =0.5 (mM)

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.0

0

2

4

6

8

10

12

14

16

18

0.0

20

t (s) Figure 11: Effect of fixed charge concentration on the average curvature Ka.

1.0

0

2

4

6

8

10

12

14

16

18

0.9

0.9 *

c =1.0 (mM) * c =0.5 (mM)

0.8 0.7

Curvature Ka (1/m)

20 1.0

0.8 0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0.0

0

2

4

6

8

10 t (s)

12

14

16

18

0.0

20

Figure 12: Effect of bathing-solution concentration on the average curvature Ka.


Volume 3 • Issue 1 • 1000103


Page 12 of 13 120

0

1

2

3

4

5

6

7

8

9

110

110

100

100

Experimental results Numerical results

90 Curvature Ka (1/m)

120

80

90 80

70

70

60

60

50

50

40

40

30

30

20

20

10

10

0

0

1

2

3

4

5

6

7

8

9

0

Externally applied electric field Ve (V) Figure 13: Comparison of simulated results with experimental data for hydrogel membrane deformation.

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 1. Shahinpoor M (1994) Continuum electromechanics of ionic polymeric gels as artificial muscles for robotic applications. Smart Mater Struct 3: 367. 2. Kim B, Peppas NA (2002) Synthesis and characterization of pH-sensitive glycopolymers for oral drug delivery systems. J Biomater Sci Polym Ed 13: 1271-1281. 3. Brock D, Lee W, Segalman D, Witkkowski W (1994) A dynamic model of a linear actuator based on polymer hydrogel. J Intell Mater Syst Struct 5: 764-771. 4. Tanaka T, Nishio I, Sun ST, Ueno-Nishio S (1982) Collapse of gels in an electric field. Science 218: 467-469. 5. Osada Y, Okuzaki H, Hori H (1992) A polymer gel with electrically driven motility. Nature 355: 242-244. 6. Kim SY, Shin HS, Lee YM, Jeong CN (1999) Properties of electroresponsive poly(vinyl alcohol)/poly(acrylic acid) IPN hydrogels under an electric stimulus. J Appl Polym Sci 73: 1675-1683. 7. Homma M, Seida Y, Nakano Y (2000) Evaluation of optimum condition for designing high-performance electro-driven polymer hydrogel systems. J Appl Polym Sci 75: 111-118. 8. Homma M, Seida Y, Nakano Y (2001) Effect of ions on the dynamic behavior of an electrodriven ionic polymer hydrogel membrane. J Appl Polym Sci 82: 76-80. 9. Fei J, Zhang Z, Gu L (2002) Bending behavior of electroresponsive poly(vinyl alcohol)/poly(acrylic acid) semi-interpenetrating network hydrogel fibers under an electric stimulus. Polym Int 51: 502-509. 10. Sun S, Mak AFT (2001) The dynamical response of a hydrogel fiber to electrochemical stimulation. J Polym Sci B Polym Phys 39: 236-246. 11. Doi M, Matsumoto M, Hirose Y (1992) Deformation of ionic polymer gels by electric fields. Macromolecules 25: 5504-5511. 12. Shiga T, Hirose Y, Okada A, Kurauchi T (1993) Bending of ionic polymer gel caused by swelling under sinusoidally varying electric fields. J Appl Polym Sci 47: 113-119. 13. Shahinpoor M (1995) Micro-electro-mechanics of ionic polymeric gels as electrically controllable artificially muscles. J Intell Mater Syst Struct 6: 307-314. 14. Nemat-Nasser S, Li JY (2000) Electromechanical response of ionic polymer-metal composites. J Appl Phys 87: 3321. 15. Neubrand W (1999) Modellbildung und Simulation von Elektronen membranverfahren. PhD thesis, Universitat Stuttgart. 16. Grimshaw PE, Nussbaum JH, Grodzinsky AJ, Yarmush ML (1990) Kinetics of electrically and chemically induced swelling in polyelectrolyte gels. J Chem Phys 93: 4462. 17. Wallmersperger T, Kroeplin B, Holdenried J, Guelch RW (2001) Coupled multifield formulation for ionic polymer gels in electric fields. SPIE 8th Annu Symp Smart Struct Mats 4329: 264. 18. Zhou X, Hon YC, Sun S, Mak AFT (2002) Numerical simulation of the steady-state deformation of a smart hydrogel under an external electric field. Smart Mater Struct 11: 459. 19. Lai WM, Hou JS, Mow VC (1991) A triphasic theory for the swelling and deformation behaviors of articular cartilage. J Biomech Eng 113: 245-258. 20. Mow VC, Ateshian GA, Lai WM, Gu WY (1998) Effects of fixed charges on the stress-relaxation behavior of hydrated soft tissues in a confined compression problem. Int J Solids Struct 35: 4945-4962. 21. Lai WM, Mow VC, Sun DD, Ateshian GA (2000) On the electric potentials inside a charged soft hydrated biological tissue: streaming potential versus diffusion potential. J Biomech Eng 122: 336-346. 22. Hon YC, Lu MW, Xue WM, Zhou X (1999) A new formulation and computation of the triphasic model for mechano-electrochemical mixtures. Computational Mechanics 24: 155-165.


Volume 3 • Issue 1 • 1000103


Page 13 of 13 23. Li H, Ng TY, Cheng JQ, Lam KY (2003) Hermite-Cloud: a novel true meshless method. Computational Mechanics 33: 30-41. 24. Sun DN, Gu WY, Guo XE, Lai WM, Mow VC (1999) A mixed finite element formulation of triphasic mechano-electrochemical theory for charged, hydrated biological soft tissues. Int J Numer Meth Engng 45: 1375-1402. 25. Liu WK, Jun S, Zhang YF (1995) Reproducing kernel particle methods. Int J Numer Meth Fluids 20: 1081-1106. 26. Liu WK, Jun S, Li S, Adee J, Belytschko T (1995) Reproducing kernel particle methods for structural dynamics. Int J Numer Methods Eng 38: 1665-1679. 27. Liu WK, Chen Y, Chang CT, Belytschko T (1996) Advances in multiple scale kernel particle methods. Computational Mechanics 18: 73-111. 28. Tannehill JC, Anderson DDA, Pletcher RH (1997) Computational fluid mechanics and heat transfer. (2nd edn), Taylor & Francis, Washington DC, USA. 29. Bar-Cohen Y (2001) Electroactive polymer (EAP) actuators as artificial muscles: reality, potential, and challenges. (1st edn), SPIE Press. 30. Li H, Yuan Z, Lam KY, Lee HP, Chen J, et al. (2004) Model development and numerical simulation of electric-stimulus-responsive hydrogels subject to an externally applied electric field. Biosens Bioelectron 19: 1097-1107.

Submit your next manuscript and get advantages of OMICS Group submissions Unique features: • • •

User friendly/feasible website-translation of your paper to 50 world’s leading languages Audio Version of published paper Digital articles to share and explore

Special features:



• • • • • • • •

250 Open Access Journals 20,000 editorial team 21 days rapid review process Quality and quick editorial, review and publication processing Indexing at PubMed (partial), Scopus, DOAJ, EBSCO, Index Copernicus and Google Scholar etc Sharing Option: Social Networking Enabled Authors, Reviewers and Editors rewarded with online Scientific Credits Better discount for your subsequent articles

Submit your manuscript at: www.omicsonline.org/submission

Volume 3 • Issue 1 • 1000103