Multiphysics modeling of responsive deformation of dual magnetic-pH-sensitive hydrogel

Journal Pre-proof

Multiphysics modeling of responsive deformation of dual magnetic-pH-sensitive hydrogel Qimin Liu , Muyu Liu , Hua Li , K.Y. Lam PII: DOI: Reference:

S0020-7683(19)30463-9 https://doi.org/10.1016/j.ijsolstr.2019.11.002 SAS 10529

To appear in:

International Journal of Solids and Structures

Received date: Revised date: Accepted date:

6 February 2019 30 September 2019 4 November 2019

Please cite this article as: Qimin Liu , Muyu Liu , Hua Li , K.Y. Lam , Multiphysics modeling of responsive deformation of dual magnetic-pH-sensitive hydrogel, International Journal of Solids and Structures (2019), doi: https://doi.org/10.1016/j.ijsolstr.2019.11.002

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Ltd.

Highlights  A magneto-chemo-electro-mechanical model is developed for dual magnetic-pHsensitive hydrogel in ionic solution.  Four physicochemical responsive mechanisms are characterized, such as hydrogel magnetization, diffusions of solvent and ions, ionic polarization, and large hydrogel deformation.  Multiple interactions are considered, including the interactions between (i) the fixed charges and the mobile ions, (ii) the polymeric networks and solvent, and (iii) the mobile ions.  Both the hydrogel and surrounding solution are covered in the computational domain.  The equilibrium and transient performances of the hydrogel are investigated via the validated model under various solution pH, initial fixed charge densities, hydrogelmagnet distances, maximum magnetic field, and hydrogel geometry.

1

Multiphysics modeling of responsive deformation of dual magneticpH-sensitive hydrogel Qimin Liua, b, Muyu Liua, Hua Lib, *, K. Y. Lamb a

Hubei Key Laboratory of Roadway Bridge & Structure Engineering, Wuhan University of Technology, Wuhan 430070, P. R. China

b

School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Republic of Singapore *

Corresponding author

Email: [email protected] (Hua Li) Tel: +65 67904953; Fax: +65 67924062

2

ABSTRACT: A magneto-chemo-electro-mechanical model is developed for simulation of the swelling behavior of the dual magnetic-pH-sensitive hydrogel that is placed in an ionic solution. In this work, four physicochemical responsive mechanisms are characterized, such as the magnetization of the hydrogel, the diffusions of solvent and ions, the ionic polarization, and the nonlinear large deformation of the hydrogel. Moreover, multiple interactions are considered, including the interactions between (i) the fixed charges and the mobile ions, (ii) the polymeric networks and solvent, and (iii) the mobile ions. Furthermore, both the hydrogel and surrounding solution are covered in the computational domain, in which the Maxwell stress is included over the hydrogelsolution interface as an additional mechanical boundary. After the multiphysics model is validated via both theoretical and experimental findings in the open literature, the magnetic, electrochemical, and mechanical performances of the magnetic-pH-sensitive hydrogel are investigated in detail, and the result shows that the abrupt change in magnetic intensity occurs and the edge effect is more pronounced when approaching the hydrogel-solution interface. Furthermore, the smaller maximum magnetic field, the higher pH level, and the longer hydrogel-magnet distance contribute to the larger swelling deformation of the hydrogel. These findings may be employed to systematically design and optimize the dual magnetic-pH-sensitive hydrogel and its relevant devices. Keywords: Multiphysics model; dual magnetic-pH-sensitive hydrogel; responsive deformation

3

1. Introduction The advent of the environmental-sensitive material opens new perspectives in material science, due to its ability to convert the external physical or chemical cues to the observable change of state (Chikh Alard et al., 2017). One such candidate is the stimulusresponsive hydrogel, which is composed of three-dimensional soft polymeric networks and liquid that can react to the small perturbation in surrounding with significant changes in its morphology, the surface characteristics, and the sol-to-gel transition (Jeong and Gutowska, 2002). Currently, most of the studies on hydrogel system focused on a single stimulus. However, it is anticipated that a hydrogel for a specific application may experience not only one stimulus but several. For example, the P(NIPAAm-co-MAA) hydrogel can simultaneously react to both the temperature and pH cues, making it appropriate for drug delivery and release in human body (Brazel and Peppas, 1996). When exposed to a single stimulus (e.g. pH) at a time, the pH-sensitive hydrogel is often faced with slow response rate and poor controllability (Beebe et al., 2000). Usually, the response rate can be improved by two strategies, using the macroporous hydrogel or multiple hydrogels with a smaller size, whereas the hydrogel still suffers from low mechanical stiffness (Birgersson et al., 2008; Lee et al., 2015). To overcome the limitations, great attempt was made to develop dual magnetic- and pH-sensitive hydrogel by immobilization of the magnetic particles into the polymeric chain with fixed charge groups (Galicia et al., 2009). However, it is obvious from the literature that the responsive mechanism of the hydrogel still remains unclear, such that it is worthwhile to 4

develop a multiphysics model to investigate the magneto-chemo-electro-mechanical behavior of the hydrogel subject to various physicochemical cues at different environmental conditions.

To date, the dual magnetic-pH-sensitive hydrogels have attracted considerable attention due to their peculiar characteristics, such as biocompatibility, fast response, and remote actuation (Ghadban et al., 2015; Li et al., 2013). As such, they were employed in diverse areas including heavy ions adsorption (Tiwari and Sharma, 2013), drug delivery (Mahdavinia et al., 2014), and catalysis (Sahiner and Ozay, 2011). They were also good candidates as pH chemical sensors, since the magnetic property of the hydrogel varied due to pH-dependent hydrogel deformation (Song et al., 2014; van Berkum et al., 2015). Moreover, the hydrogel-based scaffold was synthesized and characterized to provide appropriate cell microenvironment for efficient tissue promotion by delivering various biophysical and biochemical cues, when subjected to magnetic- and pH- coupled stimuli (Li et al., 2017). For preparation of the dual responsive hydrogel, most researchers opted the hydrogen bonding (Galicia et al., 2009), coordination bonds (Zhou et al., 2012), or covalent attachment (Messing et al., 2011) to incorporate the magnetic particles into pHsensitive hydrogel with fixed charge groups, including poly (methacrylic acid) (Muzzaupo et al., 2015), poly (acrylic acid) (van Berkum et al., 2015), and poly (2hydroxyethyl methacrylate) (Hao et al., 2016).

Despite numerous experimental works on dual magnetic-pH-sensitive hydrogel, a literature survey unveils that no studies were carried out on mathematical modeling and numerical simulation of the hydrogel. According to the literature review, it is found that 5

the majority of the relevant models concentrated on the magneto-active elastomer which works in air surrounding, with the assumption of constant volume during deformation (Weeber et al., 2012). Currently, the magneto-mechanical models for the magneto-active elastomer were generally classified by the length scales. At the microscale level, the interplays between the magnetic particles, and between the particle and individual polymer chain are involved in dipole-spring models (Weeber et al., 2015). At mesoscopic scale, the polymeric networks were regarded as elastic continuum, while the embedded magnetic particles were taken as dipoles (Zubarev, 2013). At macroscopic level, various macroscopic magneto-mechanical models were formulated, in which the Maxwell equations were coupled with linear or nonlinear elasticity (Brigadnov and Dorfmann, 2003; Raikher and Stolbov, 2005). Regarding the magnetic-sensitive hydrogel, a thermodynamic model (Filipcsei and Zrínyi, 2009) was formulated to describe the equilibrium deformation of the magnetic-sensitive hydrogel submerged in water and exposed to a uniform magnetic field. However, the predicted swelling deformation is independent of magnetic field direction, which is qualitatively inconsistent with experimental measurements (Backes et al., 2015; Genoveva and Miklos, 2010). To overcome the limitation, a multi-effect-coupling magnetic stimulus model with finite deformation was developed for magnetic-sensitive hydrogel in water, where the deformation and the instability coupled with hysteresis were well captured (Liu et al., 2017). However, the multiphysics model formulated in the present work emphasizes on the deformable hydrogel immersed in bath solution with mobile ionic species, in which the hydrogel magnetization is characterized by a magnetic field, the diffusions of solvent and ions are characterized by a chemical field, the ionic polarization is described by an 6

electric field, and the remarkable volume variation of hydrogel is described by a mechanical field.

In this work, a magneto-chemo-electro-mechanical model is developed for the dual magnetic-pH-sensitive hydrogel that is placed in an ionic solution with mobile ions. Four physicochemical responsive mechanisms are considered, such as magnetization of the hydrogel, diffusions of solvent and ions, ionic polarization, and the nonlinear deformation of the hydrogel. The governing equations are based on the conservation laws of mass and momentum, where the magnetic migration on ion transport is included. The constitutive relations are formulated by the second law of thermodynamics for the influences of the electromagnetic effect, chemical potential, and finite deformation. To verify the present model, three comparisons are carried out, where satisfactory agreements are achieved. Furthermore, the magnetic, electrochemical, and mechanical responses are simulated and predicted via the validated model for further understanding the kinetic and equilibrium behaviors of the hydrogel.

2. Theory In this section, the physicochemical responsive mechanisms of the dual magnetic-pHsensitive hydrogel are characterized by the theoretical description of the hydrogel magnetization, solvent and ions diffusions, ionic polarization, as well as the nonlinear large deformation of the hydrogel.

As shown in Figure 1, the present dual magnetic-pH-sensitive hydrogel is a combination of acidic hydrogel embedded with magnetic particles at micro- or nano7

level, which can simultaneously react to the external magnetic field and the alteration of environmental solution pH (Mahdavinia et al., 2015). The exterior solution is comprised of four species, solvent molecules (i.e. water), hydrogen ion, co-ions (i.e. chloride ion) that have the same charge as the fixed charges, and the counterions (i.e. sodium ion) that possess the opposite charge to the fixed charges. If the hydrogel is submerged in bath solution, part of the functional groups dissociates in the polymeric chain, which causes the formations of the fixed charges attached to the chain and the mobile ions. Due to the fixed charges, the ionic concentration within interior hydrogel differs from that at the exterior solution. Therefore, an osmotic pressure is generated and it contributes to the swelling or deswelling of the hydrogel. Concurrently, if an external magnetic field is imposed, the embedded magnetic particles are magnetized and then attracted to the place with the higher magnetic field. Since the particles are tightly anchored onto the flexible chain, the movement of the particles may drive the morphology of the whole hydrogel, which in turn redistributes the magnetic field, the fixed charge density, the concentration of mobile ions, as well as the electric potential. Hence, it is a full-loop magneto-chemoelectro-mechanical coupled system, which requires an efficient mathematical model to characterize its fundamental mechanism.

2.1.

Governing equations

Without free current, the distribution of magnetic field in magnetostatics is described by (Dorfmann and Ogden, 2005)

 H  0,

B  0

8

(1)

where   () and   () denote the material curl and divergence operators respectively, B and H are the vectors of magnetic induction and intensity in reference configuration

respectively, and they are related via B  0 JC1 (H  M)

(2)

where  0 denotes the vacuum magnetic permeability, M is the magnetization, F is the deformation gradient with the determinant J  det(F) , and C  FT F is the right CauchyGreen tensor. Across the interface between the hydrogel and the environmental solution, the magnetic intensity H and induction B are needed to meet the jump conditions below (Dorfmann and Ogden, 2005), N  [[H]]  0 ,

N  [[B]]  0

(3)

where N is the unit outward normal, and the double square bracket represents a quantity jump across the material surface from the inside to the outside.

For description of the chemical field associated with the solvent diffusion, the conservation law of mass is employed as (Coussy et al., 1998)

C s    J s  Rs

(4)

where Cs denotes the solvent concentration, J s and Rs are the solvent flux and the rate of solvent generation by chemical reaction.

Similarly, the chemical field relating to the ionic diffusion is expressed as 9

C k    J k  Rk

(5)

where Ck , J k , Rk are ionic concentration, flux, and generation for the kth mobile species respectively. By the Gauss’s theorem (Lai and Li, 2011), the electric displacement D in the hydrogel and surrounding solution satisfies the equation as follows,

  D  Fr ( z f C f   z k C k )

(6)

k

where z f and C f are the valence and concentration of the fixed charges, z k and C k are the valence and concentrations of the mobile ions, and Fr is the Faraday constant.

To convert the magnetic, chemical, and electric coupled effects to the large nonlinear deformation of the hydrogel, the conservation law of momentum is utilized via (Liu et al., 2017)   P f 0 V h b

(7)

where  0 , Vh , P , and f b denote the mass density of the hydrogel, the deformation velocity, the nominal stress tensor, and the external force density respectively.

2.2.

Constitutive equations

To form the constitutive relations for the present hydrogel system with finite deformation at the isothermal condition, the inequality of free energy imbalance is used by (Gurtin et al., 2010) 10

d 1 2 {  (   0 Vh )dv}  W  C  E  M dt V0 2

(8)

where  is the Helmholtz free energy density, and W is the mechanical power given by

W   PN  Vh ds   fb  Vh dv S0

(9)

V0

For the chemical power C , two contributions are considered in the present study, namely, (i) the diffusions of the solvent and ions, and (ii) the mass generation Rs and Rk by chemical reaction. As a result, the power C is expressed as (Nardinocchi et al., 2011) C     s J s  Nds    s Rs dv   (   k J k  Nds   k Rk dv) S0

k

V0

S0

(10)

V0

where  s and  k represent the chemical potentials of solvent and mobile ions respectively. The rate of work E done by the electric field is given by (Drozdov and deClaville Christiansen, 2015)

  F   z C )dv E   (E  D r k k k

V0

where E and  are the electric field strength and the electric potential respectively.

11

(11)

In general, the magnetic power M is mainly contributed by two parts, one is related to magnetization of the hydrogel (Han et al., 2011), and the other is associated with the magnetic migration on the mobile ions. Consequently, we have

  F L z C )dv M   (H  B r k k

(12)

k

V0

where L is a function that satisfies the relationship v  B  L (Zhang et al., 2018).

By sum of Equations (9) to (12), the balances of mass and momentum (4), (5), and (7), and the divergence theorem, the local form of the inequality (8) is rewritten as

  P : F  s Cs   k Ck  J s s   (J k k ) k

k

 H  B  Fr L zk Ck  E  D  Fr  zk Ck k

(13)

k

If the free energy density  is specified as a function of the five independent variables, namely the magnetic induction B , the electric displacement D , the concentrations of solvent C s and ions Ck , as well as the deformation gradient F , we have

(

     P) : F  (  H)  B  (  E)  D  (  s  )Cs F B D Cs

 [  k  Fr (  L) zk ]Ck  J s s   J k k  0 Ck k k

(14)

Owing to the arbitrariness of the variables, the coefficient in any bracket vanishes, such that the constitutive equations are given below 12

P

     , H , E , s  , k   Fr (  L) zk F B D Cs Ck

(15)

To guarantee the inequality (14) held, the solvent flux J s and the ionic flux J k are related to the chemical potential gradient  s and  k respectively through (Hong et al., 2008)

Js  

Cs Ds C1  s , k BT

Jk  

Ck Dk C1  k k BT

(16)

where k BT is the temperature in the unit of energy, Ds and Dk are the diffusion coefficients of solvent and ions respectively.

2.3.

Free energy function

For simulation and analysis of the responsive behavior of the dual magnetic-pH-sensitive hydrogel, a specific free energy density  is required, in which five components are considered in the present study, including the stretching of the polymer networks  ela , the mixing of the polymers and solvent  sol , the mixing of the solvent and ions  ion , the magnetization  mag , and the polarization  ele . Based on the Flory-Rehner theory, the free energy density  ela and  sol with the effect of particle volume fraction  p are given respectively by (Ding et al., 2016; Flory, 1953; Flory and Jr., 1943; Liu et al., 2017)

ela 

1 G{tr(FT F)  3  2 ln[det(F)]} 2

13

(17)

 sol 

v C p 1 p k BT [vs Cs ln( s s )   H vs C s ] vs det(F) det(F)

(18)

where G denotes the shear modulus of hydrogel, vs is the volume per solvent species, and  H is the Flory-Huggins parameter that characterizes the interaction between solvent and polymer networks.

The free energy density due to the mixing of the solvent with mobile ions is expressed as (Drozdov and deClaville Christiansen, 2015)

ion  k BT  Ck (ln Ck  ln Cs  1)

(19)

k

The free energy densities of magnetization  mag and polarization  ele for an ideal hydrogel in reference configuration are given in the form of (Pelteret et al., 2016)

 mag 

B  CB D  CD ,  ele  2 0  r J 2 0 r J

(20)

where  0 is the dielectric permittivity in vacuum,  r and  r are the relative dielectric and magnetic permeabilities respectively.

By sum of Equations (17) to (20), and with consideration of the molecular incompressibility condition 1  vs Cs  det(F) (Hong et al., 2008; Liu et al., 2016), the free energy density  is given by

14



vs C s   p 1p 1 k T G{tr(FT F)  3  2 ln[det(F)]}  B [vs Cs ln( )   H vs C s ] 2 vs 1  vs C s 1  vs C s D  CD B  CB    k BT  Ck (ln Ck  ln Cs  1)  [1  vs Cs  det(F)] 2 0 r J 2 0  r J k

(21)

where  is the Lagrange multiplier that is also interpreted as osmotic pressure (Hong et al., 2008).

By the free energy equation (21) and the constitutive relation (15), the magnetic intensity H , the electric field strength E , the chemical potentials  k and  s , and the stress P are expressed respectively as

H  CB /( 0  r J )

(22)

E  CD /( 0 r J )

(23)

k  kBT ln

 s  k BT [ln(

vCs   p 1  vCs

)

P  G(F  F T )  JF T 

Ck  Fr (  L) zk Cs

1 p

1p 2  H ( ) ]  k BT  Ck / Cs  vs 1  vCs 1  vCs k

(24)

(25)

(FD)  D D  CD T (FB )  B B  CB T (26)  F   F  0 r J 2 0 r J 0  r J 20  r J

On the basis of the Faraday's law, the magnetic intensity H is associated with magnetic scalar potential  via H   , and the electric field strength E relates to the electric potential

 via E   . By Equations (1) and (22), the magnetic potentials

 in the hydrogel and the exterior solution are written as 15

  ( J0 r C1 )  0

(27)

and in a similar way, by Equations (6) and (23), the electric potentials

  ( J 0 r C 1 )   Fr ( z f C f   z k Ck )

 are rewritten as (28)

k

where the fixed charge density of acidic hydrogel C f is achieved through the Langmuir monolayer adsorption theory (Lai and Li, 2011) with consideration of the interaction between the fixed charges and the mobile ions

Cf 

s0C 0f JK d JK d  CH

(29)

where K d denotes the dissociation constant,  s0 is the volume fraction of the polymer at reference state, C 0f is the concentration of fixed charges at dry state, as well as CH  is the concentration of the hydrogen ion.

Based on Equations (5), (16)2 and (24), the ion transport without chemical reaction is given by    zCF Ck     Dk C1  Ck  k k r（  v  B)   kBT   

(30)

where the diffusion of the mobile ions is dependent on three driving forces, namely the ion concentration gradient, electric force, and magnetic effect. In addition, the diffusion equation developed by the present model is consistent with the theoretical work (Qin and Bau, 2011) for magnetohydrodynamic flow in microfluidic device under the static 16

electric and magnetic fields. If the present hydrogel is immersed in an unstirred solution, the ionic diffusion contributed by convection can be neglected (Li et al., 2009), such that the induction term v  B is not considered in the subsequent analysis.

For the mechanical equilibrium of the hydrogel, we assume the equilibrium always maintains at all time during the transient evolution (Bouklas et al., 2015). Thereby, Equation (7) is rewritten as

P  0

(31)

if the external body force f b is neglected.

So far the mathematical model has been theoretically developed for the dual magnetic-pH-sensitive hydrogel immersed in solution with mobile ions when subjected to magneto-chemo-electro-mechanical coupled fields. The governing equations are comprised of the magnetostatic equation (27), the electrostatic equation (28), the ionic diffusion equation (30), and the mechanical equilibrium equation (31). Furthermore, the density of free energy is specified via Equation (21), and the constitutive relations given by Equations (22) to (26), and (29).

To solve the nonlinear coupled partial differential equations (PDEs) with the relevant constitutive relations, the boundary and initial conditions are required and thus given in the following subsection.

17

2.4.

Boundary and initial conditions

In this paper, a two-dimensional (2D) simulation is conducted with the multiphysics model for performance of a cylindrical magnetic hydrogel placed in an unstirred bath solution, as shown in Figure 2. As a result of the axisymmetry, the Neumann magnetic and electrochemical boundary conditions are used at the hydrogel center, namely

C k    0,  0 at R  0  0, R R R

(32)

and at the edge of solution domain S 0 , the concentrations of the ions, the electric and magnetic field are fixed, such that the Dirichlet boundary conditions are employed, namely C H   C H  , C  C , C  C ,   0 ,    ( R, Z ) at

S 0

(33)

where    ( R, Z ) is constructed from the experimentally measured magnetic field along the axial direction by Equation (27). The overbar denotes the parameter at the solution edge under the effect of bath solution only. At the hydrogel-solution interface  0 , we have the following mechanical boundary condition,

P  N  Tm at  0

(34)

where Tm is the magnetic tractions acting on the boundary  0 and Tm  Pm N , here the Maxwell stress Pm is given by 18

Pm  FT (H  B)  (H  B)F T /2

(35)

Furthermore, at the interface  0 , the chemical potentials of ions and solvent at the interior hydrogel near the interface equals to those at the exterior hydrogel, namely

sin  sex , Hin  Hex , +in  ex , in  ex

(36)

Insertion of Equation (25) into Equation (36)1, the osmotic pressure at the interface is given below

 int erface

kT  B vs

kT  B vs

 CHin  

 vCsin   p 1   p in 1  p 2  )  C   ( )   ln( s H 1  vCsin 1  vCsin 1  vCsin    Cin  Cin CHex  Cex  Cex    Csin Csex 

(37)

where Ckin and Ckex represent the ion concentration at the interior and exterior hydrogel near the interface, respectively.

To predict the kinetic response of the magnetic hydrogel, the initial conditions are required. In the present study, it is assumed that the hydrogel is initially in the equilibrium state under the effect of bath solution with a specific pH value (pH1) and subject to an external magnetic field. This equilibrium state is taken as the initial condition for transient simulation, and the equilibrium state is broken by varying the pH1 in the bath solution to another one (pH2). Therefore, the initial conditions used for the numerical calculation are given as follows Trans Equ Trans Equ Trans Equ Equ CpH2 (X, 0)  CpH1 (X, 0)  pH1 (X, 0)   pH1 , pH2 ,  pH2 , uTrans pH2 ( X, 0)  u pH1 (38)

19

Equ Equ Equ where  pH1 , CpH1 , pH1 and u Equ pH1 denote the magnetic potential, concentration, electric

potential, and hydrogel displacement in equilibrium state when solution pH = pH1.

2.5.

Model implementation

The solution to the initial/boundary value problem includes the vector field of location x(X, t ) and the scalar fields of magnetic potential  ( X, t ) , electric potential  ( X, t ) ,

solvent chemical potential s ( X, t ) , ionic chemical potential k ( X, t ) , which are coupled and evolve concurrently with time. The weak forms of the problem are obtained by using a pair of test functions,  x ,  ,  ,  s , and k , which meet the necessary integrability conditions. Multiplying Equation (27) by  , integrating over V0 , and applying the divergence theorem, we have



0

( J 0 r C1 )( )dV0 = 

0

BN ( )dS0

(39)

where  0 is the computational domain covering both the hydrogel and surrounding, and BN  B  N is the magnetic induction along the unit normal.

Similarly, the weak forms are obtained respectively for the distributive electric field, diffusions of solvent and ions, and the mechanical deformation



0

( J  0 r C1 )( )dV0 = 

0



0

DN ( )dS0   Fr ( z f C f   zk Ck )dV0 0

Css dV0   J s (s )dV0   0

0

20

(40)

k

rs dS0   Rss dV0 0

(41)



0



0

Cks dV0   J k (k )dV0  

ik dS0   Rk k dV0

(42)

0 Vh xdV0   P :  FdV0   T   xdS0   fb   xdV0

(43)

0

0

0

0

0

0

where r  J s  N and i  J k  N denote the surface fluxes of solvent and ion respectively, and T  P  N is the nominal traction on the surface.

The weak forms (39) ~ (43) constitute a magneto-chemo-electro-mechanical coupled system sufficient to describe the field evolution, and can be directly implemented into the finite element package COMSOL MULTIPHYSICS 5.3 operating on a 64bit 12 CPUs (2.6GHz) workstation with 32GB RAM. For implementation of the PDEs and the constitutive relations, the following interfaces are utilized: (i) PDE interfaces used for the ion transport, (ii) AC/DC module involved for the distributions of the magnetic and electric fields, and (iii) Solid mechanics module used for large deformation of the hydrogel. The computational domain covers both the hydrogel and its surrounding solution, which is discretized into the free quadrilateral elements. To solve the initial/boundary value problem, usually two solution approaches are employed, fully coupled or segregated. For the fully coupled approach, all the equations are coupled together, i.e. the distributive magnetic and electric fields, ion concentration, and hydrogel deformation are calculated at the same time. For the segregated approach, the solution process is separated into multiple steps, where the results from the previous step are used as the input for the subsequent step. Since the fully coupled approach converges more robust and in less iterations compared to the segregated one (Xu et al., 2013), the fully coupled approach is employed in the present simulation. A relative convergence error of 21

10-3 is employed for all the variables. The solutions for different mesh sizes are compared to ensure that the numerical solutions are convergent, i.e. independent of the element size. In order to verify the code, the numerical results are further compared with three analytical and experimental results in the open literature.

3. Results and discussion 3.1.

Validation of the multiphysics model

To verify the model, three comparisons are conducted, such as (i) validation with the analytical solution for the distributive magnetic field with varying magnetic permeabilities, (ii) validation with the experimental finding in the open literature for equilibrium deformation of hydrogel with varying solution pH levels, and (iii) validation with the experimental results in the open literature for equilibrium deformation of magnetic hydrogel.

3.1.1. Validation for distributive magnetic field As well known, the analytical solution of the distributive magnetic field merely exists for magnetic hydrogel with an ellipsoidal shape or other shapes with infinite geometry (Liu et al., 2018). Herein, the first comparison is made with the theoretical result for distribution of magnetic field across the free space and the cylindrical hydrogel with infinite size along the axial direction under a uniform magnetic field, as shown in Figure 3(a). The hydrogel with radius R0  2 mm is centered in a surrounding with length Lx and height Ly , and Lx  Ly =16R0 . At the surrounding edge S0 , the magnetic intensity

H 0 remains unperturbed or uniform, namely 22

 ( x,  Ly / 2)  0 ,  ( x, Ly / 2)   H 0 Ly ,

 x

(  Lx /2, y )

 0,

 x

( Lx /2, y )

0

(44)

and at the hydrogel-surrounding interface  0 , we have the following boundary conditions according to the boundary condition (3), namely

[[ ]]  0 , 2

 N

0

 1

 N

0

(45)

where the magnetic permeabilities 1  0 and 2  r 0 . By Equation (27), the analytical magnetic field H is calculated as follows

r  R0   A(er sin   e cos  ) H 2 2 （  H 0 +D /r )sin  er  ( D /r  H 0 ) cos  e r  R0

(46)

where  is the polar angle, A  22 H 0 / (1  2 ) and D = (1 -2 ) R02 H 0 / (1  2 ) .

To select an appropriate mesh size, a mesh independent study is conducted, where the mesh size (denoted by mesh scale in COMSOL) increases continuously. Subsequently, the corresponding simulated results for the magnetic field distribution along the coordinate x  0 are compared with the analytical solution by Equation (46), which is quantified by the global error (Mukherjee and Mukherjee, 1997). Table 1 presents the mesh independent study for the simulation under different mesh scales. As seen from the table, with the decrease of the mesh size, the number of elements increases gradually while the global error decreases between the simulated results and the analytical solution. In particular, when the mesh scale is extra fine, the global error of 0.64% is found. If the mesh scale varies from extra to extremely fine, the global error changes from 0.64 to 0.61. 23

However, the number of computational elements increases significantly by about 246%. Therefore, the mesh scale of extra fine is used for the present simulation.

Apart from the 2D geometry, a relevant 3D geometry with mesh scale of extra fine is considered in this simulation. The number of the 3D elements is 147300. The computational times for the 2D and 3D simulations are 1 and 168s respectively. After computing, the global error is 0.38% for the magnetic field distribution between the between the 2D and 3D geometries, meaning no significant change is found. As such, the 2D geometry is used for the subsequent analysis.

Figure 3(b) demonstrates the comparison of the simulation with the theoretical result, for effect of the relative magnetic permeability  r on the distribution of non-dimensional magnetic intensity H y / H 0 along the coordinate x  0 . The magnetic intensity H y / H 0 remains almost 1 when the normalized coordinate is far away from the hydrogelsurrounding interface, demonstrating that the original magnetic field is not disturbed. However, a sharp decrease occurs in the magnetic intensity H y / H 0 due to the magnetic boundary condition (3), when the normalized coordinate y / R0 increases to approach the interface. In addition, the intensity H y / H 0 keeps unchanged within the hydrogel and it decreases with relative magnetic permeability  r due to the following two reasons, (i) constant demagnetizing factor for the infinite hydrogel geometry, and (ii) the increasing demagnetizing effect with the increase of the permeability  r (Liu et al., 2017). Furthermore, it is clearly found that good agreement is achieved between the simulated

24

and theoretical results, showing the satisfactory ability of the present model for predicting the magnetic response behavior with desired accuracy.

3.1.2. Validation for solution pH-actuated deformation The second comparison is conducted with the experimental data (De et al., 2002), where the acidic hydrogel was synthesized by the homogeneous copolymers of acrylic acid and hydroxyl-ethyl-methacrylate acid with the diacrylate crosslinker. The cylindrical hydrogel was initially at dry state and constrained along the axial direction by the cover glasses, such that the hydrogel swelled merely in the radial direction. Thereby, it is reasonable to employ a 2D model to predict the equilibrium hydrogel swelling that is placed in buffer solution with ionic strength of 300 mM and varying pH levels from 2 to 12. For the present simulation, the input parameters are given below (De et al., 2002):

C  C  300mM

,

CH   10(3-pH) mM

,

C 0f  1800mM

,

Kd  0.01mM

,

F  1.603 1019 C ,  0  8.8542 1012 F/m , and  r  80 . It is noted that the experimental

shear modulus G varies in a trilinear form (De et al., 2002). At low pH  5.5 , the shear modulus almost remains constant G  101.4 kPa , whereas it falls down to another constant G  80.4 kPa if pH  7.5 . At the intermediate pH range from 5.5 to 7.5, the linear decrease of shear modulus G is observed.

Figure 4 illustrates the comparison of the present simulation with the experimental result in the open literature for the equilibrium deformation of the hydrogel, where the hydrogel radius R0  200 and 350 μm respectively. The hydrogel swelling is quite dependent on the environmental pH. As seen from Figure 4(a), the diameters of the 25

hydrogel remain around 500 and 950 μm when pH  4 and pH  7 respectively, while the hydrogel undergoes a sudden change in the pH level ranging from 4 to 7. The pHactuated hydrogel swelling is highly associated with the pendant acidic groups attached to the backbone of the polymeric network. The acidic group is slightly ionized if the solution pH is far below the pKa value of the hydrogel (i.e. pKa  4.7 ), which leads to a minimal swelling of the hydrogel. However, if the environmental pH increases and approximates the pKa value of the hydrogel, the acidic functional group may be remarkably ionized, which results in a sharp increase in the fixed charges within the hydrogel. It is well known that the increase of fixed charges in the acidic hydrogel initiates the diffusion of mobile ions with positive charges. It then enlarges the difference of the ionic concentration between the interior hydrogel and exterior surrounding solution, which causes a large osmotic pressure to stimulate the hydrogel to swell more. When the ionization degree of the hydrogel reaches its saturation point (i.e. pH  7 ), no obvious change is found in the swelling of the hydrogel, if the solution pH increases further. It is also noted from the two comparisons that the predicted hydrogel deformation coincides with the experimental results qualitatively and quantitatively, such that the validated model can capture well the pH-actuated responsive behavior.

3.1.3. Validation for magnetic-actuated deformation For further examination of the present work, the simulated deformation of the magnetic hydrogel in equilibrium state is compared with the experimental results in the published work (Zrínyi et al., 1996), where the cylindrical magnetic hydrogel was immersed in water subject to a nonuniform magnetic field, as shown in Figure 2. To describe the hydrogel responsive behavior, the boundary condition on the magnetic scalar potential 26

   ( R, Z ) is needed at the exterior solution edge S 0 . In this validation, it is reconstructed from the experimental work (Zrínyi et al., 1996) on the axial magnetic field. Before inserting the magnetic hydrogel in water, the magnetic induction divided by the current intensity B(0, Z ) / I was measured as a function of the distance Z from the top surface of the electromagnet, which is fitted by a fifth-degree polynomial,

B / I  60.4934 Z 5  4.8576 Z 4  9.6825 Z 3  2.9074 Z 2  0.3136 Z  0.0124 , where the maximum magnetic field H max at top surface of the electromagnet is directly related to current intensity I by H max  9946.8I A/m . Due to the relation H   , the magnetic scalar potential  is written as a sixth-degree polynomial  (0, Z )  p6 (Z ) . In addition, the distribution of the magnetic scalar potential  in bath solution satisfies the Laplace’s equation  2  0 (27). Therefore, the analytical solution of the magnetic potential  is given by  ( R, Z )  p6 (Z )  R 2 p6( 2) (Z ) / 4  R 4 p6( 4) (Z ) / 64  R 6 p6(6) (Z ) / 2304 (Afkhami et al., 2008), which can be exerted on the solution edge S 0 . The inputs of the parameters for the validation include: LR  80 mm , LH  200 mm , r0  7 mm , h0  122.1 mm ,

z0  28 mm , G  103 Pa , and r  1.14 . The comparison of the present numerical simulation with the experimental data (Zrínyi et al., 1996) is illustrated in Figure 5 for the magnetic field-induced hydrogel length change h . As seen from the figure, the deformation h increases with the increase of the current intensity I, due to the stronger magneto-elastic effect. Furthermore, the simulation results coincide with the experimental findings, it is then indicated the validity of the developed model, such that it is rational to employ the model for the studies below. 27

3.2.

Analysis of the responsive characteristics of the dual magnetic-pH-sensitive

hydrogel

In this section, the equilibrium and transient performances of the dual magnetic-pHsensitive hydrogel are investigated by the validated multiphysics model.

3.2.1. Equilibrium Performance To gain an insight into the performance of the dual magnetic-pH-sensitive hydrogel with the magneto-chemo-electro-mechanical coupled field, it is essential to characterize the distributive profiles of the magnetic field, the ionic concentration, the electric potential, and the nonlinear hydrogel deformation. To the best of our knowledge, so far no effort has been made to investigate the multi-response of the magnetic-pH-sensitive hydrogel. Herein, the validated magneto-chemo-electro-mechanical model is employed for the first time to study the performance of the cylindrical hydrogel that is submerged in a buffer solution with the mobile ions of Na+, H+, Cl-, and OH- at physiological condition (Lai et al., 1991), as shown in Figure 2. The input parameters include: C 0f  100 mM , CNa  CCl  150 mM , r0  10 mm , h0  20 mm , G  40 kPa ,

 H  0.6 , Kd  0.01 mM ,

and  p  0.15 . The relative magnetic permeability  r is calculated via the well-known Maxwell-Garnett formula (Garnett, 1906). The original nonuniform magnetic field applied follows that in Section 3.1.3, and unless specified otherwise, the present hydrogel is submerged in ionic solution at physiological pH 7.4.

Figure 6(a) and (b) investigate the impact of the maximum magnetic intensity H max at the top surface of electromagnet on the non-dimensional magnetic intensity H z / H max , 28

and induction Bz / ( 0 H max ) along the vertical line located at the coordinate R  0 across the cylindrical magnetic-pH-sensitive hydrogel and its surrounding solution. As seen from Figure 6(a), the magnetic intensity H z / H max decreases linearly with nondimensional coordinate z / LZ , then it experiences a sharp decrease at the hydrogelsolution interface. After that, it drops to a minimum until it reaches the interface again, where it rises abruptly. Afterward, the intensity H z / H max decreases linearly again with

z / LZ . The distribution of magnetic intensity is more obviously shown in the inset of Figure 6(a), in which the dash and solid lines represent the hydrogel geometry at dry and current state respectively. Such distributive profile of the magnetic field may be explained through the two following reasons: (i) the different magnetic permeabilities between the hydrogel and the surrounding solution, and (ii) the boundary condition by Equation (3) for the magnetic intensity at the hydrogel-solution interface, which demonstrates that normal component of the magnetic intensity H discontinuously goes through the material boundary (Bustamante et al., 2011). In addition, no visible discrepancy is found between the original magnetic field and the magnetic field in the presence of hydrogel, if the coordinate z / LZ keeps far away from the hydrogel, as shown in Figure 6(a). Nevertheless, they diverge significantly from each other if the coordinate z / LZ gets close to the hydrogel-solution interface, due to the pronounced edge effect in the vicinity of the hydrogel (Liu et al., 2018). Furthermore, the distributive magnetic intensity is also varied with different magnetic intensities H max , probably due to the varying field-induced mechanical deformation of the hydrogel. Figure 6(b) illustrates the distributive profile of the non-dimensional magnetic induction 29

Bz / (0 H max ) . Similarly to the distribution of H z / H max , the edge effect is also observed and then disappear with the increase of the normalized distance z / LZ . However, unlike the discontinuous distribution of H z / H max at hydrogel-solution interface, Bz / ( 0 H max ) continuously goes through the interface to meet the boundary condition (3).

Figure 6(c) and (d) are plotted for analysis of the maximum magnetic intensity H max on the variation of the non-dimensional magnetic intensity H z / H max , and magnetic induction Bz / ( 0 H max ) along the horizontal line located at Z  22.5mm across the magnetic-pH-sensitive hydrogel and its surrounding solution. As observed from Figure 6(c), the original magnetic field applied decreases gradually with increasing coordinate

r / LR . However, it is interesting that the magnetic intensity H z / H max increases first and then it drops down. Moreover, the intensity H z / H max decreases with increasing H max , probably due to the different field-induced deformations. The non-dimensional axial magnetic intensity H z / H max is tangential to the lateral surface, on the basis of the magnetic boundary condition (3), it thus goes continuously across the hydrogel-solution interface, as shown in Figure 6(c). In terms of the distributive pattern of the nondimensional magnetic induction Bz / (0 H max ) , it jumps across the interface to satisfy the boundary condition (3), as shown in Figure 6(d).

It is commonly known that the electrochemical property of the dual magnetic-pHsensitive hydrogel is much dependent on the chemical cue in solution, e.g. pH. In general, the physiological pH level of the human body approximates 7.4 and shows distinct 30

differences among the organs, such as the bladder with urine (i.e. pH 4.6~8) and tumor tissue (i.e. pH 4.5~6) (Hao et al., 2016). Herein, the distribution patterns of ionic concentration including the sodium (Na+), chlorine (Cl-), hydrogen (H+) ions, and of the electric potential  along the coordinate R  0 , are visualized in Figure 7, as a function of environmental pH values of 4, 6, and 8, where H max  200kA/m . The ionic concentrations of Na+ and Cl- are much higher than that of H+ ion. The concentrations of H+ and Na+ ions are relatively higher within the acidic hydrogel and lower in bath solution, whereas the concentration of Cl- ion exhibits the inverse variation. In addition, the concentration of Na+ and Cl- ions are the same in the bath solution, while the former is higher than the latter in the interior of the hydrogel. The concentration difference may be attributed to the fixed charges with the valence of -1 within the hydrogel, in which the counterions are attracted into the hydrogel and the co-ions with the same valence as the fixed charges are repelled out of the hydrogel to maintain the electroneutrality inside the hydrogel (Lai and Li, 2011). The distributive electric potential  along the axial direction is shown in Figure 7(d), where it remains zero in the bath solution, due to the electroneutrality in the net total concentration of the ions. However, a very small (in millivolts) constant electric potential is visualized within the hydrogel due to the net difference between the concentrations of various positive and negative ionic species (Li et al., 2005). As also observed from the figures, it is noted that the distributions of ionic concentrations and the electric potential are also greatly affected by the ambient solution pH. For example, if the environmental pH rises up from 4 to 6, the concentration of Na+ ion increases around 27.5% from about 160 to 204 mM. However, the impact of solution pH becomes less significant with an increase of about 2% when the pH level increases 31

further from 6 to 8. Probably the majority of the fixed charge groups are ionized if pH=6, after which a further increase in pH causes an insignificant variation of the ionic concentration.

As a typical chemical component, the initial fixed charge density C 0f also plays a vital role in the electrochemical response of the hydrogel system (Lesperance et al., 1992). As a result, Figure 8 is plotted for the impact of the initial fixed charge density C 0f at dry state on variation of distributive patterns of the ionic concentrations and the

electric potential  along the vertical line located at coordinate R  0 , where

H max  200kA/m . It is observed that the hydrogels behave similarly to those in Figure 7. It is clearly seen that the ionic concentration distribution is highly dependent on the change of the initial fixed charge density C 0f . With the increase of the initial fixed charge density C 0f , the concentrations of the counterions (i.e. Na+ and H+) within the hydrogel increases, whereas the co-ion (i.e. Cl-) concentration decreases accordingly. This is 0 because the increase of C f may enlarge the probability of the protonization of the fixed

charge groups to attract more counterions into the hydrogel and repulse more co-ions out of the hydrogel. Furthermore, a sharp increase or decrease is found in the distributive ionic concentration at the hydrogel-solution interface. This is because the fixed charge groups merely exist in the anionic hydrogel instead of external solution.

Figure 9 theoretically demonstrates the influence of external magnetic field H max at the top surface of the electromagnet on the maximum hydrogel displacement wmax along the axial direction as a function of the environmental pH level subject to different 32

hydrogel-magnet distances of 12.5 and 25 mm respectively. In the present study, the swelling deformation wmax is defined as the height difference for hydrogel at the current and dry state. At the low pH level (i.e. pH  4 ), the hydrogel behaves like a hydrophobic polymeric network, since majority of the fixed charge groups still remains unionized. Subsequently, the hydrogel undergoes an abrupt transition on the displacement due to the increase of the ionization of the fixed charge groups. After reaching the ultimate value in the vicinity of pH 6.5, a continuous increase of the solution pH cannot further stimulate the hydrogel deformation, as a result of the fully ionized fixed charge groups. It is noticeable that this phenomenon is consistent with experimental findings by Siegel et al (Siegel and Firestone, 1988) and also by Naficy et al (Naficy et al., 2012). As also observed from the figures, the highest curve corresponds to the magnetic intensity H max of zero, while subsequently the lower curve related to the higher magnetic field H max . Therefore, the hydrogel shrinks with the external magnetic field. This is because the magnetic particles are firmly attached to the polymeric networks. Once a nonuniform magnetic field is applied, two different interactions occur, field-particle and particleparticle interactions, where the former results from the magnetic field gradient, and the latter from the particle magnetization. In general, the field-particle interaction dominates when the magnetic field applied is nonuniform (Liu et al., 2017). As such, the incorporated magnetic particles are attracted to the point with higher magnetic strength, which stimulates the contraction of whole hydrogel (Faidley et al., 2010). As also seen from the figures, the swelling deformation of hydrogel is significantly affected by the hydrogel-magnet distance z0 , where a relatively longer distance z0 leads to a smaller

33

contraction. Probably the magnetic field decays away from the top surface of the electromagnet. A longer distance z0 leads to a weak magnetic force imposed on the magnetic particles, which diminishes the contraction and thus results in a larger hydrogel deformation.

Figure 10 unveils the effect of maximum magnetic field H max on change of the maximum hydrogel displacement wmax with the initial fixed charge density of the hydrogel at dry state C 0f , where z0  12.5 and 25 mm respectively. The displacement

wmax increases gradually with increasing initial fixed charge density C 0f , and it further increases with the decrease of the external magnetic field H max . Probably increasing the fixed charge density C 0f enlarges the degree of ionization, which attracts more counterions into the hydrogel and repulses more co-ions out of the hydrogel. It then enhances the concentration difference of the ions between the interior hydrogel and the exterior solution, as illustrated in Figure 8. This enlarges the osmotic pressure and accordingly actuates the hydrogel to swell more. Furthermore, it is shown that the hydrogel-magnet distance z0 affects the swelling deformation dramatically, where the longer distance z0 weakens the magnetic field-induced deformation.

The variation of the hydrogel displacement wmax is examined numerically in Figure 11 against the shear modulus G, when subjected to varying maximum magnetic field

H max , where z0  12.5 and 25 mm respectively. In the absence of external magnetic field, the responsive swelling deformation of the hydrogel wmax falls linearly with the shear 34

modulus G. This is because the larger shear modulus G reinforces the hydrogel to make the polymer matrix more compacted. Furthermore, since the osmotic pressure as the driving force is balanced by the elastic restoring force, the increase of shear modulus G hinders the hydrogel to stretch further. However, the displacement wmax increases first and then it drops slightly with the shear modulus G if an external magnetic field is imposed. Probably magnetic field-induced contraction decreases, and its decrease rate drops gradually with increasing G. As such, when the shear modulus G is small (i.e. G  150 kPa under H max  200 kA/m ), the final hydrogel deformation may increase due

to the decreasing contraction. However, when the shear modulus G is large (i.e.

G  150 kPa under H max  200 kA/m ), the final hydrogel deformation decreases due to the insignificant magnetic effect and large stiffness of the hydrogel.

Figure 12 is plotted to investigate the variation of the maximum hydrogel displacement wmax against the initial hydrogel height h0 subject to varying magnetic field H max . It is found that the hydrogel geometry shows a significant impact on the swelling deformation of the hydrogel. In particular, when H max varies from 100 to 200 kA/m, the change in hydrogel displacement wmax is much more distinct than the change when H max varies from 0 to 100 kA/m. This phenomenon was also demonstrated in experiment work (Genoveva and Miklos, 2010), in which the hydrogel deformation is proportional to the squared of the external field strength and the increase rate rises with the magnetic field. Furthermore, it is found that the swelling deformation of the hydrogel

35

increases linearly with the hydrogel height h0 , where the slope decreases with increasing

H max .

3.2.2. Kinetic performance To study the transient characteristics of the dual magnetic-pH-sensitive hydrogel in response to the environmental solution, the simulations are conducted for the distributive profiles of the concentrations of Na+, Cl-, H+ ions, the electric potential  , magnetic intensity H z / H max along the coordinate R  0 , and the maximum hydrogel displacement

wmax , respectively, where the diffusivity of the mobile ions Dk  5 107 m/s2 . In the present study, two types of hydrogel are considered with different initial fixed charge groups C 0f of 100 and 500 mM, as illustrated in Figures 13 and 14 respectively. The hydrogels are initially in the equilibrium state under the effect of bath solution with pH 4 and the external magnetic field H max of 200 kA/m. After that, the hydrogel is moved to the bath solution with pH 8. Figure 13 illustrates the variation of the concentrations of Na+, Cl-, H+ ions, the electric potential  , magnetic intensity H z / H max , and the maximum hydrogel displacement wmax with time t, respectively, where C 0f  100 mM . It is found that the step change of the solution pH from 4 to 8 breaks the initial equilibrium state and makes the concentrations of the mobile ions redistributed with time inside the hydrogel. Furthermore, the concentrations of Na+, Cl-, H+ ions, and the electric potential  have the similar distribution patterns, where sharp gradients are developed at the hydrogel-solution interface at the beginning of the swelling, then the gradient becomes smaller, and finally 36

it vanishes with a new equilibrium state. As seen in Figure 13(a), the concentration of Na+ ion CNa increases gradually inside the hydrogel with time t, since the increase of the bath solution pH leads to the decreasing CH  in the solution, promoting the dissociation of the fixed charge groups to attract more Na+ from the solution for keeping the electroneutrality inside the hydrogel (Lai and Li, 2011). In addition, it is found that CNa is not distributed uniformly during the diffusion process. In particular, the concentration CNa at the center of hydrogel is lower than that in other locations. This is because the

Na+ ions initially gather at the hydrogel-solution interface, resulting in a large gradient of CNa near the interface. Therefore, CNa at the hydrogel center is lower than that at the

interface. Subsequently, when more and more Na+ ion enters into the hydrogel, the gradient decreases and the driving force for the ion diffusion decreases accordingly based on based on Equation (30), which causes a decreasing diffusion rate of the Na+ ion. It is also seen from Figure 13(f), no significant variation is found in the maximum hydrogel displacement wmax . It is because that the initial fixed charge density C 0f of 100 mM is small, the concentration difference between the hydrogel is relatively low. Thereby, the osmotic pressure due to the concentration difference is insignificant to drive the hydrogel deformation. Figure 14 is plotted to study the time evolution of the concentrations of Na+, Cl-, H+ ions, the electric potential  , magnetic intensity H z / H max , and the maximum hydrogel displacement wmax respectively, where the initial fixed charge density C 0f =500 mM . Their distributive profiles are similar to those in Figure 13. However, the deformation of 37

the hydrogel wmax rises significantly from 4.6 to 5.6 mm compared to merely 0.06 mm change in Figure 13(f) due to the variation of the initial fixed charge density C 0f . Furthermore, it is seen that the 800 s is required for the hydrogel to reach its new equilibrium state, while 1000 s required for hydrogel with C 0f =100 mM .

4. Conclusions A magneto-chemo-electro-mechanical model has been developed in this paper to predict the performance of the dual magnetic-pH-sensitive hydrogel that is placed in solution with ionic species. The physicochemical responsive mechanisms are included, including hydrogel magnetization, diffusions of solvent and ions, ionic polarization, and the nonlinear hydrogel deformation. The governing equations are based on the conservation laws of mass and momentum, and the constitutive relations are formulated by the second law of thermodynamics. For examination of the present mathematic model, three comparisons are conducted between the simulation results and both the theoretical and experimental results in the open literature, where satisfactory agreements are achieved. Furthermore, the magnetic, electrochemical, and mechanical behaviors are predicted via the validated model for further understanding the fundamental mechanism of the hydrogel. It is found that the magnetic field in the presence of hydrogel experiences a sudden change and the edge effect is more pronounced when approaching the hydrogelsolution interface. It is also shown that the smaller maximum magnetic field, the higher pH level, and the longer hydrogel-magnet distance contribute to the larger swelling deformation of the hydrogel. Furthermore, the ion concentrations are not distributed uniformly during the diffusion process by the transient analysis, and the more initial fixed 38

charge density can stimulate the larger deformation. Therefore, these results can be used for design and optimization of the dual magnetic-pH-sensitive hydrogel and its relevant devices.

Conflicts of interest The authors have no conflicts of interest to declare.

Acknowledgments This work was supported by Nanyang Technological University through the project [M4081151.050] and NTU Research Scholarships, and by Wuhan University of Technology through the Doctoral Startup Fund [315/40122081].

39

References Afkhami, S., Renardy, Y., Renardy, M., Riffle, J.S., St Pierre, T., 2008. Field-induced motion of ferrofluid droplets through immiscible viscous media. J. Fluid Mech. 610, 363380. Backes, S., Witt, M.U., Roeben, E., Kuhrts, L., Aleed, S., Schmidt, A.M., Von Klitzing, R., 2015. Loading of pnipam based microgels with cofe2o4 nanoparticles and their magnetic response in bulk and at surfaces. J. Phys. Chem. B 119 (36), 12129-12137. Beebe, D.J., Moore, J.S., Bauer, J.M., Yu, Q., Liu, R.H., Devadoss, C., Jo, B.-H., 2000. Functional hydrogel structures for autonomous flow control inside microfluidic channels. Nature 404, 588-590. Birgersson, E., Li, H., Wu, S., 2008. Transient analysis of temperature-sensitive neutral hydrogels. J. Mech. Phys. Solids 56 (2), 444-466. Bouklas, N., Landis, C.M., Huang, R., 2015. A nonlinear, transient finite element method for coupled solvent diffusion and large deformation of hydrogels. J. Mech. Phys. Solids 79, 21-43. Brazel, C.S., Peppas, N.A., 1996. Pulsatile local delivery of thrombolytic and antithrombotic agents using poly(n-isopropylacrylamide-co-methacrylic acid) hydrogels. J. Controlled Release 39 (1), 57-64. Brigadnov, I.A., Dorfmann, A., 2003. Mathematical modeling of magneto-sensitive elastomers. Int. J. Solids Struct. 40 (18), 4659-4674. Bustamante, R., Dorfmann, A., Ogden, R.W., 2011. Numerical solution of finite geometry boundary-value problems in nonlinear magnetoelasticity. Int. J. Solids Struct. 48 (6), 874-883. Chikh Alard, I., Soubhye, J., Berger, G., Gelbcke, M., Spassov, S., Amighi, K., Goole, J., Meyer, F., 2017. Triple-stimuli responsive polymers with fine tuneable magnetic responses. Polymer Chemistry 8 (16), 2450-2456. Coussy, O., Dormieux, L., Detournay, E., 1998. From mixture theory to biot's approach for porous media. Int. J. Solids Struct. 35 (34–35), 4619-4635. De, S.K., Aluru, N.R., Johnson, B., Crone, W.C., Beebe, D.J., Moore, J., 2002. Equilibrium swelling and kinetics of ph-responsive hydrogels: Models, experiments, and simulations. J. Microelectromech. Syst. 11 (5), 544-555.

40

Ding, Z., Toh, W., Hu, J., Liu, Z., Ng, T.Y., 2016. A simplified coupled thermomechanical model for the transient analysis of temperature-sensitive hydrogels. Mech. Mater. 97, 212-227. Dorfmann, A., Ogden, R.W., 2005. Some problems in nonlinear magnetoelasticity. Zeitschrift für angewandte Mathematik und Physik 56 (4), 718-745. Drozdov, A.D., deClaville Christiansen, J., 2015. Swelling of ph-sensitive hydrogels. Phys. Rev. E 91 (2), 022305. Faidley, L.E., Han, Y., Tucker, K., Timmons, S., Hong, W., 2010. Axial strain of ferrogels under cyclic magnetic fields. Smart Mater. Struct. 19 (7), 075001. Filipcsei, G., Zrínyi, M., 2009. Swelling of ferrogels in uniform magnetic field – a theoretical approach. Period. Polytech., Chem. Eng. 53 (2), 93-96. Flory, P.J., 1953. Principles of polymer chemistry. Cornell University, New York. Flory, P.J., Jr., J.R., 1943. Statistical mechanics of cross‐linked polymer networks ii. Swelling. J. Chem. Phys 11 (11), 521-526. Galicia, J.A., Cousin, F., Dubois, E., Sandre, O., Cabuil, V., Perzynski, R., 2009. Static and dynamic structural probing of swollen polyacrylamide ferrogels. Soft Matter 5, 26142624. Garnett, J.M., 1906. Colours in metal glasses, in metallic films, and in metallic solutions. Ii. Phil. Trans. R. Soc. A 205, 237-288. Genoveva, F., Miklos, Z., 2010. Magnetodeformation effects and the swelling of ferrogels in a uniform magnetic field. J. Phys.: Condens. Matter 22 (27), 276001. Ghadban, A., Ahmed, A.S., Ping, Y., Ramos, R., Arfin, N., Cantaert, B., Ramanujan, R.V., Miserez, A., 2015. Bioinspired ph and magnetic responsive catechol-functionalized chitosan hydrogels with tunable elastic properties. Chem. Commun. 52 (4), 697-700. Gurtin, M.E., Fried, E., Anand, L., 2010. The mechanics and thermodynamics of continua. Cambridge University Press, Cambridge, U.K. Han, Y., Hong, W., Faidley, L., 2011. Coupled magnetic field and viscoelasticity of ferrogel. Int. J. Appl. Mech. 03 (02), 259-278. Hao, L., Gwangjun, G., Seong Yong, K., Jong-Oh, P., Sukho, P., 2016. Magnetic actuated ph-responsive hydrogel-based soft micro-robot for targeted drug delivery. Smart Mater. Struct. 25 (2), 027001.

41

Hong, W., Zhao, X., Zhou, J., Suo, Z., 2008. A theory of coupled diffusion and large deformation in polymeric gels. J. Mech. Phys. Solids 56 (5), 1779-1793. Jeong, B., Gutowska, A., 2002. Lessons from nature: Stimuli-responsive polymers and their biomedical applications. Trends Biotechnol. 20 (7), 305-311. Lai, F., Li, H., 2011. Transient modeling of the reversible response of the hydrogel to the change in the ionic strength of solutions. Mech. Mater. 43 (6), 287-298. Lai, W.M., Hou, J.S., Mow, V.C., 1991. A triphasic theory for the swelling and deformation behaviors of articular cartilage. J. Biomech. Eng. 113 (3), 245-258. Lee, E., Kim, D., Kim, H., Yoon, J., 2015. Photothermally driven fast responding photoactuators fabricated with comb-type hydrogels and magnetite nanoparticles. Sci. Rep. 5, 15124. Lesperance, L.M., Gray, M.L., Burstein, D., 1992. Determination of fixed charge density in cartilage using nuclear magnetic resonance. J. Orth. Res. 10 (1), 1-13. Li, H., Luo, R., Birgersson, E., Lam, K.Y., 2009. A chemo-electro-mechanical model for simulation of responsive deformation of glucose-sensitive hydrogels with the effect of enzyme catalysis. J. Mech. Phys. Solids 57 (2), 369-382. Li, H., Ng, T.Y., Yew, Y.K., Lam, K.Y., 2005. Modeling and simulation of the swelling behavior of ph-stimulus-responsive hydrogels. Biomacromolecules 6 (1), 109-120. Li, K., Zhang, C., Qiu, L., Gao, L., Zhang, X., 2017. Advances in application of mechanical stimuli in bioreactors for cartilage tissue engineering. Tissue. Eng. Part B Rev. 23 (4), 399-411. Li, Y., Huang, G., Zhang, X., Li, B., Chen, Y., Lu, T., Lu, T.J., Xu, F., 2013. Magnetic hydrogels and their potential biomedical applications. Adv. Funct. Mater. 23 (6), 660-672. Liu, Q., Li, H., Lam, K.Y., 2017. Development of a multiphysics model to characterize the responsive behavior of magnetic-sensitive hydrogels with finite deformation. J. Phys. Chem. B 121 (22), 5633-5646. Liu, Q., Li, H., Lam, K.Y., 2018. Transition of magnetic field due to geometry of magneto-active elastomer microactuator with nonlinear deformation. J. Microelectromech. Syst. 27 (2), 127-136. Liu, Y., Zhang, H., Zhang, J., Zheng, Y., 2016. Transient swelling of polymeric hydrogels: A new finite element solution framework. Int. J. Solids Struct. 80, 246-260.

42

Mahdavinia, G.R., Etemadi, H., Soleymani, F., 2015. Magnetic/ph-responsive beads based on caboxymethyl chitosan and κ-carrageenan and controlled drug release. Carbohydr. Polym. 128, 112-121. Mahdavinia, G.R., Rahmani, Z., Karami, S., Pourjavadi, A., 2014. Magnetic/ph-sensitive κ-carrageenan/sodium alginate hydrogel nanocomposite beads: Preparation, swelling behavior, and drug delivery. J. Biomater. Sci. Polym. Ed. 25 (17), 1891-1906. Messing, R., Frickel, N., Belkoura, L., Strey, R., Rahn, H., Odenbach, S., Schmidt, A.M., 2011. Cobalt ferrite nanoparticles as multifunctional cross-linkers in paam ferrohydrogels. Macromolecules 44 (8), 2990-2999. Mukherjee, X.Y., Mukherjee, S., 1997. On boundary conditions in the element-free galerkin method. Comput. Mech. 19 (4), 264-270. Muzzaupo, R., Tavano, L., Rossi, C.O., Picci, N., Ranieri, G.A., 2015. Novel ph sensitive ferrogels as new approach in cancer treatment: Effect of the magnetic field on swelling and drug delivery. Colloids Surf. B. Biointerfaces 134, 273-278. Naficy, S., Razal, J.M., Spinks, G.M., Wallace, G.G., Whitten, P.G., 2012. Electrically conductive, tough hydrogels with ph sensitivity. Chem. Mater. 24 (17), 3425-3433. Nardinocchi, P., Pezzulla, M., Placidi, L., 2011. Thermodynamically based multiphysic modeling of ionic polymer metal composites. J. Intell. Mater. Syst. Struct. 22 (16), 18871897. Pelteret, J.-P., Davydov, D., McBride, A., Vu, D.K., Steinmann, P., 2016. Computational electro- and magneto-elasticity for quasi-incompressible media immersed in free space. Int. J. Numer. Meth. Eng 108, 1307–1342. Qin, M., Bau, H.H., 2011. When mhd-based microfluidics is equivalent to pressuredriven flow. Microfluid. Nanofluid. 10 (2), 287-300. Raikher, Y.L., Stolbov, O.V., 2005. Deformation of an ellipsoidal ferrogel sample in a uniform magnetic field. J. Appl. Mech. Tech. Phys. 46 (3), 434-443. Sahiner, N., Ozay, O., 2011. Highly charged p(4-vinylpyridine-co-vinylimidazole) particles for versatile applications: Biomedical, catalysis and environmental. React. Funct. Polym 71 (6), 607-615. Siegel, R.A., Firestone, B.A., 1988. Ph-dependent equilibrium swelling properties of hydrophobic polyelectrolyte copolymer gels. Macromolecules 21 (11), 3254-3259. Song, S.H., Park, J.H., Chitnis, G., Siegel, R.A., Ziaie, B., 2014. A wireless chemical sensor featuring iron oxide nanoparticle-embedded hydrogels. Sensors Actuators B: Chem. 193, 925-930. 43

Tiwari, A., Sharma, N., 2013. Efficiency of superparamagnetic nano iron oxide loaded poly(acrylamide-co-maleic acid) hydrogel in uptaking cu2+ ions from water. J. Dispersion Sci. Technol. 34 (10), 1437-1446. van Berkum, S., Biewenga, P.D., Verkleij, S.P., van Zon, J.B.A., Boere, K.W.M., Pal, A., Philipse, A.P., Erné, B.H., 2015. Swelling enhanced remanent magnetization of hydrogels cross-linked with magnetic nanoparticles. Langmuir 31 (1), 442-450. Weeber, R., Kantorovich, S., Holm, C., 2012. Deformation mechanisms in 2d magnetic gels studied by computer simulations. Soft Matter 8 (38), 9923-9932. Weeber, R., Kantorovich, S., Holm, C., 2015. Ferrogels cross-linked by magnetic nanoparticles—deformation mechanisms in two and three dimensions studied by means of computer simulations. J. Magn. Magn. Mater. 383 (0), 262-266. Xu, X., Li, Z., Nehorai, A., 2013. Finite element simulations of hydrodynamic trapping in microfluidic particle-trap array systems. Biomicrofluidics 7 (5), 054108. Zhang, X., Wang, Q.J., Shen, H., 2018. A multi-field coupled mechanical-electricmagnetic-chemical-thermal (memct) theory for material systems. Comput. Methods Appl. Mech. Eng. 341, 133-162. Zhou, Y., Sharma, N., Deshmukh, P., Lakhman, R.K., Jain, M., Kasi, R.M., 2012. Hierarchically structured free-standing hydrogels with liquid crystalline domains and magnetic nanoparticles as dual physical cross-linkers. J. Am. Chem. Soc. 134 (3), 16301641. Zrínyi, M., Barsi, L., Büki, A., 1996. Deformation of ferrogels induced by nonuniform magnetic fields. J. Chem. Phys 104 (21), 8750-8756. Zubarev, A., 2013. Magnetodeformation of ferrogels and ferroelastomers. Effect of microstructure of the particles' spatial disposition. Phys. A 392 (20), 4824-4836.

44

List of Tables Table 1. Mesh independent study for the simulation ....................................................... 48

45

List of Figures Figure 1. Schematic of a dual magnetic-pH-sensitive hydrogel placed in ionic solution subject to externally applied magnetic field. .................................................... 49 Figure 2. Initial and boundary conditions for the multiphysics model. ........................... 50 Figure 3. (a) Schematic of a magnetic hydrogel immersed in surrounding medium under a uniform magnetic field. (b) Comparison of the magnetic field distribution between the simulation results and theoretical counterparts when subjected to different relative magnetic permeabilities. ....................................................... 51 Figure 4. Comparison of the present simulation and published experiment for deformation of the acidic hydrogel with the solution pH ranging from 2 to 12, when subjected to the hydrogel radius R0 of 200 (a) and 350μm (b). ............ 52 Figure 5. Comparison of the present theoretical simulation (solid line) with the experiment (square brick) for the magnetic-actuated deformation of the magnetic hydrogel. ........................................................................................... 53 Figure 6. Influence of the maximum magnetic field H max on variation of the distributive profiles of non-dimensional (a) magnetic intensity H z / H max , and (b) magnetic induction Bz / ( 0 H max ) along the coordinate R  0 , as well as nondimensional (c) magnetic intensity H z / H max , and (d) magnetic induction

Bz / (0 H max ) along the coordinate Z  22.5mm . ........................................... 54 Figure 7．Influence of ambient solution pH on variation of the concentrations of (a) Na+, (b) Cl-, (c) H+ ions, and (d) the electric potential  along the coordinate R  0 . .......................................................................................................................... 55

46

Figure 8 ． Influence of the initial fixed charge density C 0f on variation of the concentrations of (a) Na+, (b) Cl-, (c) H+ ions, and (d) the electric potential  along the coordinate R  0 ............................................................................... 56 Figure 9． Influence of the maximum magnetic field H max on variation of the maximum hydrogel displacement wmax against the varying pH level ranging from 2 to 10, when the hydrogel-magnet distance z0 is 12.5 (a) and 25 mm (b) respectively. .......................................................................................................................... 57 Figure 10. Influence of the maximum magnetic field H max on variation of the maximum hydrogel displacement wmax against the varying initial fixed charge density C 0f , when the hydrogel-magnet distance

z0 is 12.5 (a) and 25 mm (b)

respectively. ...................................................................................................... 57 Figure 11. Influence of the maximum magnetic field H max on variation of the maximum hydrogel displacement wmax with different shear modulus G, when the hydrogel-magnet distance z0 is 12.5 (a) and 25 mm (b) respectively. ............. 58 Figure 12. Influence of the maximum magnetic field H max on variation of the maximum hydrogel displacement wmax with different hydrogel height h0 , when the hydrogel-magnet distance z0 is 12.5 (a) and 25 mm (b) respectively. ............. 58 Figure 13． Time evolution of the concentrations of (a) Na+, (b) Cl-, (c) H+ ions, (d) the electric potential  , (e) magnetic intensity H z / H max along the coordinate

R  0 , and (f) the maximum hydrogel displacement wmax , when the solution pH changes from 4 to 8 and the initial fixed charge density C 0f =100 mM ...... 59

47

Figure 14． Time evolution of the concentrations of (a) Na+, (b) Cl-, (c) H+ ions, (d) the electric potential

 , (e) magnetic intensity H z / H max along the coordinate

R  0 , and (f) the maximum hydrogel displacement , when the solution pH changes from 4 to 8 and the initial fixed charge density C 0f =500 mM . .......... 60

48

Table 1. Mesh independent study for the simulation Mesh scale

Max. element size (mm)

Min. element size (mm)

Max. element growth rate

Curvature factor

No. of elements

Global error (%)

Coarse

3.2

0.064

1.4

0.4

260

12.1

Fine

1.7

0.010

1.3

0.3

699

8.7

Finer

1.2

0.004

1.25

0.5

1257

8.1

Extra fine

0.64

0.002

1.2

0.25

3180

0.64

Extremely fine

0.32

0.00064

1.1

0.2

11030

0.61

49

Figure 1. Schematic of a dual magnetic-pH-sensitive hydrogel placed in ionic solution subject to externally applied magnetic field.

50

Figure 2. Initial and boundary conditions for the multiphysics model.

51

(b)

(a)

Figure 3. (a) Schematic of a magnetic hydrogel immersed in surrounding medium under a uniform magnetic field. (b) Comparison of the magnetic field distribution between the simulation results and theoretical counterparts when subjected to different relative magnetic permeabilities.

52

(a)

(b)

Figure 4. Comparison of the present simulation and published experiment for deformation of the acidic hydrogel with the solution pH ranging from 2 to 12, when subjected to the hydrogel radius R0 of 200 (a) and 350μm (b).

53

Figure 5. Comparison of the present theoretical simulation (solid line) with the experiment (square brick) for the magnetic-actuated deformation of the magnetic hydrogel.

54

(a)

(b)

(c)

(d)

Figure 6. Influence of the maximum magnetic field H max on variation of the distributive profiles of non-dimensional (a) magnetic intensity H z / H max , and (b) magnetic induction

Bz / (0 H max ) along the coordinate R  0 , as well as non-dimensional (c) magnetic intensity H z / H max , and (d) magnetic induction Bz / ( 0 H max ) along the coordinate

Z  22.5mm .

55

(a)

(b)

(c)

(d)

Figure 7．Influence of ambient solution pH on variation of the concentrations of (a) Na+, (b) Cl-, (c) H+ ions, and (d) the electric potential  along the coordinate R  0 .

56

(a)

(b)

(c)

(d)

Figure 8 ． Influence of the initial fixed charge density C 0f on variation of the concentrations of (a) Na+, (b) Cl-, (c) H+ ions, and (d) the electric potential  along the coordinate R  0 .

57

(a)

(b)

Figure 9． Influence of the maximum magnetic field H max on variation of the maximum hydrogel displacement wmax against the varying pH level ranging from 2 to 10, when the hydrogel-magnet distance z0 is 12.5 (a) and 25 mm (b) respectively.

(a)

(b)

Figure 10. Influence of the maximum magnetic field H max on variation of the maximum hydrogel displacement wmax against the varying initial fixed charge density C 0f , when the hydrogel-magnet distance z0 is 12.5 (a) and 25 mm (b) respectively.

58

(a)

(b)

Figure 11. Influence of the maximum magnetic field H max on variation of the maximum hydrogel displacement wmax with different shear modulus G, when the hydrogel-magnet distance z0 is 12.5 (a) and 25 mm (b) respectively.

(a)

(b)

Figure 12. Influence of the maximum magnetic field H max on variation of the maximum hydrogel displacement wmax with different hydrogel height h0 , when the hydrogelmagnet distance z0 is 12.5 (a) and 25 mm (b) respectively.

59

(a)

(b)

(c)

(d)

(e)

(f)

Figure 13． Time evolution of the concentrations of (a) Na+, (b) Cl-, (c) H+ ions, (d) the electric potential  , (e) magnetic intensity H z / H max along the coordinate R  0 , and (f)

60

the maximum hydrogel displacement wmax , when the solution pH changes from 4 to 8 and the initial fixed charge density C 0f =100 mM .

61

(a)

(b)

(c)

(d)

(e)

(f)

Figure 14． Time evolution of the concentrations of (a) Na+, (b) Cl-, (c) H+ ions, (d) the electric potential  , (e) magnetic intensity H z / H max along the coordinate R  0 , and (f) the maximum hydrogel displacement , when the solution pH changes from 4 to 8 and the initial fixed charge density C 0f =500 mM . 62

Declaration of Interest Statement The authors have no conflicts of interest to declare.

63