Polymer Gel Based Actuator:
Dynamic model of gel for real time control
by
Woojin Lee
B.S.M.E., Massachusetts Institute of Technology, 1988
M.S.M.E., Massachusetts Institute of Technology, 1990
SUBMITTED TO THE DEPARTMENT OF MECHANICAL ENGINEERING
IN PARTIAL FULFILLMENT OF THE REQUIREMENTS
FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY
AT THE
MASSACHUSETTS INSTITUTE OF TECHNOLOGY
MAY 1996
Massachusetts Institute of Technology 1996
Signature of Author: _________________________________________________________
Department of Mechanical Engineering
May 1996
Certified by: ________________________________________________________________
Dr. David L. Brock
Thesis Supervisor
Certified by: ________________________________________________________________
Dr. J. Kenneth Salisbury, Jr.
Thesis Supervisor
Accepted by: _______________________________________________________________
Professor A.A. Sonin
Chairman, Department Graduate Committee
Polymer Gel Based Actuator:
Dynamic model of gel for real time control
by
Woojin Lee
Submitted to the Department of Mechanical Engineering
on May 3, 1996 in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy
Abstract
This thesis explored the potential of polymer gel based actuator as an alternative actuation method by focusing on static/dynamic modeling of actuator system suitable for real-time control implementation. Such model must be simple enough to be implemented in real-time control, yet still able to capture the essential dynamics such as the coupling effect among chemical, electrical and mechanical energy domains of polymer gel. The developed model was tested and verified through computer simulation and experimentation employing the gel based actuator system prototype. The verified model was then incorporated into real-time control of the actuator prototype.
Thesis Supervisor: Dr. David L. Brock
Title: Research Scientist, MIT AI Laboratory
Thesis Supervisor: Dr. J. Kenneth Salisbury, Jr.
Title: Principal Research Scientist, MIT AI Laboratory
Acknowledgments
I would like to thank the members of the Artificial Muscle Laboratory for their support and encouragement. I am grateful to David L. Brock and Ken Salisbury for their continuous support, encouragement and insights.
My O-thu-gi, Meesook devoted tremendous, endless love, and has been patiently waiting for this moment. Meesook and my two kids, Calvin and Jonathan, made this process more than enjoyable, and they have never failed to put smile back on my face whenever things were going tough. Thanks for having faith in me.
Thanks also to my grandmother, parents and my brothers for their continuous love, patience, encouragement and support. My grandmother, mom and dad always believed in me, and my two brothers were always there for me. Words cannot describe how much I'm grateful.
Special thanks goes to my parents, sisters and brother in law. Their day-to-day support made this long journey possible, and without them, we would not be where we are today. Thank you.
Lastly, I would like to thank God for giving me everything I have.
TABLE OF CONTENTS
Abstract 2
Acknowledgments 3
Table of Contents 4
List of Symbols 6
Chapter
1 Introduction 8
1.1 Motivation 8
1.2 Polymer Gel Actuation 9
1.3 Objective 10
2 Background 12
3 Dynamic Model of Polymer Gel for Control System 16
3.1 Introduction 16
3.2 Dynamic Model of Non-Electrolyte Gel 20
3.2.1 Mixing/diffusion process: Flory's model 22
3.2.2 Rubber elasticity: Flory's model (statistical thermodynamics) 27
3.2.3 Rubber elasticity: simple mechanical model 30
3.2.4 Overall process: mixing/diffusion and deformation 35
3.2.5 Polymer gel swelling kinetics 35
3.3 Dynamic Model of Electrolyte Gel 39
3.3.1 Equilibrium conditions 39
3.3.2 Chemical capacitance model 43
3.3.3 Electric potential capacitance model 44
3.3.4 Kinetic model 49
3.3.5 Ion diffusion limited reaction 53
3.3.5 Overall ion transport process: binary solution case 60
3.3.6 Ionic contribution and electrolyte gel swelling 65
4 Model Verification 69
4.1 PVA-PAA in Acetate Buffer Solution 69
4.1.1 Acetate buffer chemistry 70
4.1.2 PVA-PAA gel preparation 75
4.2 Equilibrium swelling model verification 75
4.2.1 Equilibrium swelling model 75
4.2.2 Experimental procedure 77
4.2.3 Experimental result 78
4.2.4 Model parameter estimation 82
4.2.5 Model result 84
4.2.6 Discussion 92
4.3 Swelling kinetics model verification 96
4.3.1 Limiting process identification procedure 97
4.3.2 Experimental result 100
4.3.3 PVA-PAA gel dynamic model 101
4.3.4 Model simulation and discussion 104
5 Implementation 111
5.1 Demonstration actuator device 111
5.2 Results 113
6 Conclusion 117
6.1 Review 117
6.2 Contributions 123
6.3 Future 123
Appendix
A Actuator system design 126
B PAN fiber mechanical testing 140
References 148
LIST OF SYMBOLS
To ambient temperature
Po ambient pressure
Ucv internal energy of control volume
Scv entropy of control volume
Vcv volume of control volume
Fx uniaxial tensile force
L elongation
chemical potential of polymer
chemical potential of solvent
np number of moles of polymer inside the control volume
ns number of moles of solvent inside the control volume
Gcv1 Gibbs free energy of control volume at state 1
Gcv2 Gibbs free energy of control volume at state 2
Gso standard Gibbs free energy of solvent
Gpo standard Gibbs free energy of polymer
change in Gibbs free energy of control volume due to mixing
Scv1 entropy of control volume at state 1
Scv2 entropy of control volume at state 2
Sso standard entropy of solvent
Spo standard entropy of polymer
change in entropy of control volume due to mixing
R universal gas constant
is the volume fraction of solvent,
is the volume fraction of polymer, and
s is the number of chain segments of a polymer molecule such that the volume of a polymer molecule is equal to s time the volume of a solvent molecule, i.e., molar volume of polymer = s times molar volume of solvent
Flory-Huggins polymer-solvent interaction parameter
change in Gibbs free energy of control volume due to elastic deformation
change in entropy of control volume due to elastic deformation
k Boltzmann's constant
effective number of chains in the gel
elongation factor (length / dry length) in x-axis
elongation factor in y-axis
elongation factor in z-axis
Lo,s swollen, unstretched length of gel
Lo dry length of gel
ei
normal axial strains in i-axis (i= x, y, z)si
normal axial stress in i-axisE elastic modulus of elastomer
n
Poisson’s ratio (or "Poisson-like" ratio in the context of simple mechanical model of rubber elasticityLido dry, stress-free gel of linear dimension in i-axis
Li isotropically swollen and stretched dimension in i-axis
Aid cross sectional area of deformed elastomer unit aligned in i-axis
a
i elongation factors in i-axisvso molar volume of pure solvent
Pm permeability factor of the solvent in polymer gel network
U macro-scale "superficial velocity", defined as the flow rate Q divided by the total (solid plus fluid) cross-sectional area
Km permeability in the context of Darcy's law
m
vis solvent’s fluid dynamic viscosityD
z gel’s characteristic thickness solvent travels during swellinga radius of the fiber in Happel and Brenner's model
b fluid envelope radius in Happel and Brenner's model
m
i,e electrochemical potential of ion ixi,g mole fraction of ion i
zi valence of ion i
f
i,g electrical potential of the ion i inside the gelF Faraday constant, 9.648 x 104 coul/mol
CE,i electric potential coupling capacitance
dielectric constant of the electrolyte
d typical maximum separation of charges in a volume element, ~ 10-10 m
qi charge of the ion [coul]
CE,i bulk electroneutrality potential capacitance of ion species i
r
u total net charge![]()
Di diffusion coefficient [m2/sec]
ci concentration [mol/m3]
ui mobility [m2/volt*sec]
kf forward reaction constant
kr reverse reaction constant
Keq reaction's equilibrium constant
t
process time constantKa acid ionization constant of acetic acid
Chapter 1
Introduction
1.1 Motivation
The electric motor is the primary actuator used in numerous applications ranging from household appliances to robotics. However, in the fields of robotics, its high weight, limited sizes, complex transmissions and restrictive shapes have constrained the design and development of robotic systems. Although improvements in volume to weight ratio, clever transmissions, and alternatives such as pneumatics and hydraulics are being developed, actuators remain the fundamental constraint on robot performance. Nitinol, a shape memory alloy, has also been used in robotics, however the thermal to mechanical energy conversion has proven difficult in practical applications.
Over the last 10 years there has been a great deal of effort directed toward developing more dexterous robots (Salisbury 1984, Jacobsen 1986, Hess 1989) and master controllers (Bejczy and Salisbury, 1984, Agronin, 1987, Jacobsen 1989 and Iwata, 1990) capable of permitting human control of these complex devices. In the development of these systems a great deal of effort was spent in developing the mechanism of getting power to the individual finger joints and in making the system as light as possible so that the limited power available was not wasted on supporting the structure. The robot hand most similar to the human geometry, the UTAH/MIT hand, has a "remotizer" which consists of numerous pulleys (nearly 400) and tendons filling approximately a 12" by 4" by 3" area. This remotizer allows the pneumatic cylinders to be tightly packed in a remote location as it was not possible to locate them at the individual joints. It unfortunately makes integration of this device into a robot with a full range of motion at the wrist difficult. The Salisbury design has fewer degrees of freedom and uses steel cables routed inside sheaths around the wrist area. This device is more practical, but has higher friction and lower bandwidth. Others also tried placing the actuators in the base of the palm with fewer degrees of freedom, making the device easy to integrate with other robots while sacrificing bandwidth and backdrivability.
All of these groups were working around a fundamental limitation in today's actuator technology. Actuators small and strong enough to be placed proximal to the joint, as muscles are in the human system, were not available. What would be ideal for robotic actuation and rehabilitation devices would be an analog of a human muscle --- a contractile compliant "artificial muscle" material driven by chemical, thermal or electrical signals. In addition to direct applicability to robotics and rehabilitation, such material may be useful in producing smaller valves, bio-compatible end effectors for laparoscopic and arthroscopic surgery, and effective tactile displays. Recent research in polymer gels offers the hope for such an "artificial muscle" and a potential revolution in actuator design.
1.2 Polymer Gel Actuation
The polymer gel can swell or shrink considerably in response to various external stimuli such as temperature, pH, electric fields, or solvent. It has been suggested that these gels could be used as the basis for a compliant actuator or "artificial muscle". Recent research on polymer gel has made actuation based on this technology more practical. First, much research on gel based actuation has focused on the material itself, hence the physics of gel contraction has become well understood. Second, a number of innovations have made gel based actuation much more practical. For instance, thin films and fine fibers have been fabricated, allowing faster contraction rates and higher strength (Chiarelli, 1989). Finally, technological innovations such as robotics and implantable artificial biological organs have created a demand for such devices.
1.3 Objective
The design of a practical polymer gel actuator requires consideration not only of the material, but also the supporting mechanics, stimulation method, energy storage system, energy efficiency, power delivery technique, packaging method, static/dynamic model, and control system. Among these issues, availability of a simple thermodynamically correct static/dynamic model of gel/support system greatly enhances the investigations of the rest of the issues such as the stimulation method, energy efficiency, power delivery technique, etc. In addition, such a simplified model is integral to a real-time control system implementation. Therefore, I have chosen the development of a simplified model of polymer gel based actuator system suitable for real-time control as the objective.
The objective of this thesis was to develop a lumped parameter model of a gel based actuator for real-time control. The research consisted of three parts: (1) model development, (2) experimental model verification, and (3) application of the dynamic model using a real-time control system. A prototype actuator was constructed based on pH stimulated poly(vinyl Alcohol)-poly(acrylic acid) (PVA-PAA) gel since its swelling behavior is continuous with minimum hysteresis. The model was experimentally verified, and once verified, it was incorporated into real-time control system of prototype actuator system.
Chapter 2
Background
Previous work by researchers in the area of polymer gel largely fall into three categories: (1) gel fabrication and characterization, (2) theoretical modeling of swelling behavior, and (3) application (primarily drug delivery systems.) Mainly due to the limitation of the available material, much of the research on gel based actuation has focused on material fabrication and characterization, including studies of the mechanical properties (Chiarelli, 1988; Suzuki, 1991). Recent development in polymer gel fabrication have produced thinner and stronger gel fibers. For example, a gel based actuator formed by bundling 25mm diameter poly(acrylonitrile)-polypyrrole (PAN-PPY) and poly(vinylalchohol) (PVA) fibers was constructed, and demonstrated contraction rates on the order of a second and a force up to 100 N/cm2 (Chiarelli, 1989). At Hokkaido University in Japan, an electrically-driven chemomechanical poly(2-acrylamido-2-methylpropanesulfonic acid) (PAMPS) polymer gel sheet was fabricated which produced swinging motion under alternating 10V electric field. This actuation's response time was on the order of 2 seconds (Okuzaki and Osada, 1994).
The theoretical modeling of gel's swelling behavior mainly consist of equilibrium swelling volume studies and electromechanochemical kinetic studies. Much of the work in the area of equilibrium swelling volume study (Hirotsu, et. al., 1987; Li and Tanaka, 1989) are largely based on the pioneering work by Flory (Flory, 1953). In the area of gel kinetics, the need to predict the dynamic behavior for applications such as the drug delivery and the gel based actuation stimulated recent researchers (Tanaka and Fillmore, 1979; Matsuo, 1988; Grimshaw, et. al., 1990; Gehrke, et al, 1992; Segalman, et. al., 1994).
With the availability of promising materials, interest in polymer gel based actuator has increased. A number of contractile gel devices have already been constructed, including a robot gripper (Caldwell, 1989), multifingered hand (Kurauchi, et. al., 1991), artificial fish (Kurauchi, et. al., 1991), artificial urethral sphincter (Chiarelli, 1989), and artificial muscle (Brock, 1991; Brock, Lee, Segalman, Witkowski, 1994; Shahinpoor and Mojarrad, 1994). While these devices successfully demonstrated the feasibility of polymer gel based actuation, they have not yet fully addressed the engineering design and control issues necessary for these devices to be practical.

Figure 2.1: A force generating element, an artificial urethral sphincter and a robot gripper are some of the mechanical devices constructed from contractile polymer gel.
The following tables summarize some of the published research on developing polymer gels.
pH
|
Kuhn 1950: pH artificial muscle (PAA) |
|
DeRossi 1986: Isotonic and isometric force measurement (PVA-PAA) |
|
Chiarelli 1989: Artificial Urethral Sphincter (PAN) (PAN-PPY) |
|
Grimshaw 1989, 90: Dynamic membranes, selective macromolecular transport |
Electric Field
|
Hamlem 1965: Contractile gel fiber (PVA-PAA) |
|
Tanaka 1982: Gel cylinders (PAM) |
|
Chiarelli 1988, 89, Osada 1985: Gel disks and strips (PVA-PAA) |
|
DeRossi 1986: Gel strips (PVA-PAA) |
|
Grimshaw, et. al. 1990: Kinetics of gel strips (PMAA) |
|
Ozuki and Osada 1994: Gel sheet (PAMPS) |
Temperature
|
Tanaka 1978, 80, 81: N-isopropylacrylamide (NIPA) |
|
Ilavsky 1982: N-isopropylacrylamide (NIPA) |
|
Hirotsu 1987: N-isopropylacrylamide (NIPA) |
|
Matsuo 1988: Kinetics of contraction/dilation of N-IPA gel spheres (NIPA) |
Solvents
Sodium Acrylate
|
Tanaka 1987: Surface patterns on gel spheres (PAM) |
|
Nagasaws 1988: Surface patterns on gel spheres (PAM) |
Acetone/Water
|
Hrouz 1981: Anionic gels: Diethylacrylamide-sodium methacrylate copolymer) |
|
Ilavsky 1982: Anionic gels: Diethylacrylamide-sodium methacrylate copolymer) |
|
Nicoli 1983: Role of ionized groups on (PAM) |
|
Katayama 1985: Cationic gels -- Acrylamide, trimethyl(N-acryloyl-3-aminopro, N,N'-methylenebisacrylamide copolymer |
|
Caldwell 1990: Servo controlled gripper |
Salt
|
Ohmine 1982: Salt concentrations |
Light
|
Mamada 1990: Ultraviolet light ionizes gel and creates osmotic pressure |
|
Tanaka 1990: Visible light directly heats the polymer network inducing phase transition (NIPA/Ch) |
Enzymes
|
Kokufata 1991: Lectin filled gel (Concanavalin A) responses to saccharide: Saccharide dextran sulfate - dilation Methyl-D-manopyranoside - contraction (NIPA/Lectin) |
Definitions::
PAA Poly(acrylic Acid)
PVA-PAA Poly(vinyl Alcohol) Poly(acrylic Acid)
PAN Poly(acrylonitrile)
PAN/PPY Poly(acrylonitrile) Poly(pyrrole)
PAM Poly(acrylamide)
PAMPS Poly(2-acrylamido-2-methylpropanesulfonic acid)
NIPA N-isopropylacrylamide
Chapter 3
Dynamic Model of Polymer Gel for Control System
3.1 Introduction
In order for an actuator to be practical, its behavior must be predictable and controllable. In both cases, a quantitative static/dynamic model of the actuator is essential, and this model should be simple enough for the real-time control system implementation. In other words, the lumped parameter dynamic model, on the contrary to the continuum model, of both the gel and the support systems are necessary. Even though many researchers have investigated the equilibrium swelling behavior (Flory, 1953; Tanaka, 1981; Vasheghani-Farahani, et. al., 1990) and continuum kinetics of polymer gel (Grimshaw, et. al., 1980; Li and Tanaka, 1991; Quinn and Grodzinsky, 1993; Segalman, et. al., 1994), a dynamic model which is simple enough to be incorporated into the control system, yet able to capture essential dynamics was not yet available. Most of the gel kinetics researches as mentioned above were primarily motivated by applications such as the drug delivery and the selective macromolecular transports where both the detailed spatial and temporal descriptions of gel/membrane's state are essential. Subsequently, kinetic models of gel swelling developed so far were mainly based on the macroscopic continuum description. On the other hand, for applications such as the gel actuation, not only the detailed spatial description offered by the continuum approach is often unnecessary, the complexity of continuum model would be simply much too high to be useful as a simple design or real-time control model. In fact, in designing practical actuation devices, it would be often desirable to treat the gel as a uniform bulk unit with average values of gel's state properties, thereby neglecting the spatial variations all together for simplicity, and only the temporal description of gel's swelling behavior would then be of concern.
In addition, most of the previous models treated the equilibrium behavior and kinetics separately (Quinn and Grodzinsky (1993) investigated the relationship between the gel swelling kinetics and the gel's longitudinal modulus in the context of a poroelastic theory, thereby relating the swelling kinetics to the equilibrium swelling to some extent. However, tight integration of swelling kinetics and the equilibrium state was still not offered.), and therefore did not sufficiently address the electro-chemo-mechanical coupling effects of gel's static/dynamic behavior. Researchers such as Genuini, et. al. (1990) have proposed a lumped parameter model where a diffusion-reaction process was simplified as a time delay process and the mechanical readjustment process was simplified as a standard linear viscoelastic solid. The parameters of the mechanical readjustment process were somewhat artificially related to the free swelling volume of the gel at the given pH. Also, the possible effect of external mechanical forces on the gel's swelling behavior was not mentioned.
Polymer gel actuation involves interactions among multiple energy domains -- electrical, chemical and mechanical -- thus an energy based unified approach to dynamic modeling was appropriate. In this thesis, we employed the bond graph (network thermodynamic) formalism (Paynter, 1961, 1981; Oster, Perelson and Katchalsky, 1973) to aid us in developing the gel model. In the bond graph approach, each subunit of the gel was modeled as a multiport capacitance with ports to multiple energy domains which experiences quasistatic changes where the equilibrium states were defined by modified Flory's equilibrium swelling model. The subunits were then linked to each other by mechanical transmissions and the ions/solvent diffusion processes. Not only did this approach correctly account for the coupling effect between the chemical and the mechanical energy domains, but the topological layout and connectivity of the bond graph model components enhance the intuition of gel actuation dynamics. Dividing the gel structure into many fine subunits should yield results equivalent to the finite difference method, and thus may be used for checking validity. However, it is worthwhile to reemphasize that, in strict sense, the lumped parameter model of polymer gel as presented in this thesis is valid only in the case where the spatial variations of gel's properties can be safely ignored. As the gel's dimension increases, the prediction based on simple few subunits model would increasingly deviate from the actual behavior as well, and in such case, finer division of the gel system into more subunits would be necessary to capture the spatial variations. As with lumped parameter modeling of any system, one should proceed with caution before applying the model as presented in this thesis blindly.
In the following sections, the bond graph based lumped parameter models of both non-electrolyte and electrolyte gel are described. According to Flory, the swelling of gel can be viewed as the result of the reduction of interstitial solvent's chemical potential where the gel's swelling process is driven by the solvent's chemical potential gradient across the gel and bath boundary. The forces influencing the solvent's chemical potential can be largely categorized into three effects; mixing of polymer and solvent, elastic deformation of the polymer network, and the mixing with the mobile ion constituents. Assuming that the gel solution is sufficiently diluted, Flory treated these three contributions to be additive. In the case of non-electrolyte gel such as poly(vinyl) alcohol (PVA), only the first two effects are relevant since non-electrolyte gel's polymer network lacks fixed ionic groups necessary to interact with mobile ion constituents. In the case with electrolyte gel, the non-electrolyte gel model can be simply extended by adding the mobile ions mixing effect. Adopting the Flory's model on gel's swelling, I have first developed the lumped parameter model of non-electrolyte gel based on the quasistatic assumption, and then it was extended to the electrolyte gel case by adding the effect of the mobile ions' mixing.
In the development of the lumped parameter model of non-electrolyte gel, the fundamental modeling assumption employed were that the gel's swelling process is relatively slow, solvent diffusion-limited, isothermal, isobaric and quasistatic, and the mixing/diffusion effect has been modeled after the Flory's work. Although the theory of rubber elasticity of swelling gel has already been developed by Flory based on complex statistical theory, due to its complexity and difficulty in use, I have developed an alternative elasticity model based on simple mechanical analog. The resulting model of non-electrolyte gel was then extended to arrive at the lumped parameter model of electrolyte gel where the states of interstitial mobile ions and polymer network's fixed charge groups were determined by electrical and chemical effects such as the Donnan partition, bulk electroneutrality and diffusion-limited reactions.
The electrolyte gel model was then applied to the case of poly(vinyl Alcohol)-poly(acrylic acid) (PVA-PAA) gel in acetate buffer solution case where the kinetics limiting process has been identified in order to further simplify the complexity of the model. The equilibrium and kinetics predictions of the simplified model was then experimentally verified, and as a demonstration, the model was incorporated into a real-time control system of prototype actuator system.
3.2 Dynamic Model of Non-Electrolyte Gel
Consider the polymer gel network and solvent system inside a control volume (c.v.) as shown in Figure 3.1. The control volume changes as the gel network absorbs the surrounding solvent fluid. The modeling assumptions were: (1) relatively slow process, (2) solvent diffusion-limited, (3) inertial and convection effects negligible, (4) quasistatic (i.e., the properties inside the c.v. are uniform at any instant of time), (5) isothermal, (6) isobaric, (7) two sub-processes (mixing/diffusion and elastic expansion of the polymer-solvent), and (8) polymer and solvent incompressible. The c.v. is assumed to be surrounded by the solvent fluid at constant temperature, To, and constant pressure, Po. As the solvent fluid enters the c.v., the c.v. expands, and the state of the system inside the control volume at an instant of time can be thermodynamically described by
![]()
where Ucv is the internal energy,
Scv is the entropy,
Vcv is the volume,

Figure 3.1: Control volume for non-electrolyte gel model. As the gel's volume changes by either contracting or swelling, the control volume is changed together with the gel's surface.
To is the ambient temperature,
Po is the ambient pressure,
F is the uniaxial tensile force,
L is the elongation,
is the chemical potential of polymer,
is the chemical potential of solvent,
np is the number of moles of polymer, and
ns is the number of moles of solvent inside the control volume.
Since the amount of polymer inside the c.v. is fixed,
. In addition, PodVcv, represented the differential expansion work done by the control volume. Assuming both the polymer and the solvent to be incompressible, PodVcv = 0 since the expansion of c.v. alone was exactly canceled by the flow work done by the solvent flowing into the c.v.. The above expression then simplifies to
.
With the assumption of a quasistatic process, the system inside the control volume was treated as a lumped parameter multiport capacitance field with the relation
.
The first term on the right hand side (RHS) represented the effects of heat transfer and entropy transfer due to convection and the entropy generation due to irreversibility. The second term on RHS represented the chemical energy brought in by the incoming surrounding solvent, and the last term represented the effect of mechanical work on the gel. In Figure 3.2, the gel was represented as a multiport capacitance with three energy ports corresponding to thermal, chemical and mechanical interactions respectively. Constitutive relations describing the polymer-solvent gel network system have been derived by Flory, and were the base of our multiport capacitance model.

Figure 3.2: Bond graph representation of non-electrolyte gel. The gel is represented as a multiport capacitance with three energy ports corresponding to thermal, chemical and mechanical interactions.
3.2.1 Mixing/diffusion process: Flory’s model
In state 1 of Figure 3.3, pure polymer network and pure solvent inside the control volume (c.v.) are completely separated whereas in state 2, polymer network and solvent inside the control volume mix and form homogeneous polymer solution, while maintaining constant volume of c.v. The process from state 1 to state 2 is therefore entropy generating run-down mixing process, and the changes in thermodynamic properties of control volume due to such mixing process are described by Flory.
Change in Gibbs free energy of control volume due to mixing,
, may be defined as
![]()

Figure 3.3: Illustration of entropy generating mixing process. In state 1, pure polymer network and pure solvent inside the control volume are completely separated whereas in state 2, they mix and form a homogeneous solution, while maintaining constant volume.
where Gcv1 is the Gibbs free energy of c.v. at state 1,
Gcv2 is the Gibbs free energy of c.v. at state 2,
Gso is the standard Gibbs free energy of solvent,
Gpo is the standard Gibbs free energy of polymer.
Similarly, change in entropy of control volume due to mixing,
, may be defined as
![]()
where Scv1 is the entropy of c.v. at state 1,
Scv2 is the entropy of c.v. at state 2,
Sso is the standard entropy of solvent,
Spo is the standard entropy of polymer.
According to Flory, the Gibbs free energy change due to mixing of polymer and solvent is considered to be of two effects; combinatorial effect and polymer-solvent contact effect. The combinatorial effect accounts for the many different conformations the amorphous polymer may adopt, and the polymer-solvent contact effect accounts for the effects of non-randomness induced by intermolecular interactions. The combinatorial effects on entropy and Gibbs free energy are
![]()
![]()
![]()
where R is the universal gas constant,
ns is the molar quantity of solvent,
np is the molar quantity of polymer,
is the volume fraction of solvent,
is the volume fraction of polymer, and
s is the number of chain segments of a polymer molecule such that the volume of a polymer molecule is equal to s time the volume of a solvent molecule, i.e., molar volume of polymer = s times molar volume of solvent.
The polymer-solvent contact effect involves enthalpy change, and its effect on Gibbs free energy is
![]()
where
is the Flory-Huggins polymer-solvent interaction parameter.
is a temperature dependent dimensionless quantity which characterizes polymer-solvent interactions, and can also be expressed in form
![]()
where a and b are temperature independent quantities.
can be then written as
![]()
Since

we then have
![]()
It is interesting to note that
represents the thermal energy generated due to the mixing process, and thus can be modeled as the heat transferred to the isothermal surrounding solvent environment.
Combining both the combinatorial and polymer-solvent contact effects, we have
![]()
![]()
![]()
For quasistatic mixing/diffusion process, the rate of change of
with respect to time can be modeled as
![]()
and furthermore, the rate of change of
may be expressed as

With the assumptions of isothermal and isobaric conditions, the quasistatic mixing/diffusion process implies

and the above expressions simplify to
![]()
Since the entropy and Gibbs free energy of gel system inside the control volume are defined as
![]()
where
are constant during the process, the rate of change of Scv, and Gcv with respect to time are
![]()
Combining the above results, we then arrived at
![]()
This process can be graphically represented by a simple bond graph structure as shown in Figure 3.4. The entropy generating mixing/diffusion process is modeled by a resistance field driven by solvent’s chemical potential difference across the gel’s control volume boundary, and the total entropy change of control volume is sum of the entropy generated and the entropy transferred by the heat/mass transfer processes.

Figure 3.4: Bond graph representation of entropy generating mixing/diffusion process. The resistance field models the mixing/diffusion process driven by solvent's chemical potential difference across the gel's control volume boundary. The total entropy change is sum of the entropy generation and entropy transfer.
3.2.2 Rubber elasticity: Flory’s model (statistical thermodynamics)
In addition to the mixing model, Flory also developed a model for the rubber elasticity effect of polymer gel, and extensive research has followed upon his work. Here we will briefly describe his model, and investigate how his model can fit into a lumped parameter frame work.
According to Flory's work, deformation of polymer gel network may be modeled as a two step processes; first, expansion due to free swelling and second, stretching at constant volume without changing the enthalpy of the gel system. Based on the statistical theory of rubber elasticity where the end-to-end distances of the polymer chains is described by the Gaussian distribution, the enthalpy change due to the deformation is assumed to be negligible, and the Gibb's free energy change,
, and the entropy change,
, are described by
![]()
![]()
where k is the Boltzmann's constant,
is the effective number of chains in the gel,
represents the elongation factor (length / dry length) in x-axis,
represents the elongation factor in y-axis,
represents the elongation factor in z-axis.
Flory assumed that the polymer network formed at volume Vp is subsequently swollen isotropically by a diluent to a volume V such that the volume fraction of polymer is
. In the subsequent deformation due to stretching, the volume is assumed to stay constant. Letting
,
, and
represent the changes in dimensions resulting from the combination of swelling and elongation,
is therefore constant during the elongation. Considering only a simple elongation in the x-direction, let
represent the length in this direction relative to the swollen, unstretched length Lo,s = (V/Vo)1/3 Lo. Since the swelling is isotropic,
![]()
Then

![]()
Substituting the above results, Flory then arrived at

For the deformation process during swelling and stretching, Flory assumed the heat of deformation
to be negligible, and therefore the Gibbs free energy change due to elastic deformation according to Flory is

With the constant temperature and pressure assumption,
, which is equal to the deformation work input, is a function of the solvent quantity in the c.v. and the elongation only, i.e.,
. In addition, according to the statistical thermodynamic theory of rubber elasticity, the change in internal energy of elastomer due to the deformation is small. Thus, for constant volume deformation process, we have
![]()
![]()

where F and Lx are tensile force and length in x-axis.
Taking the above results by Flory, the quasistatic description of the process can be represented as

where

This further simplifies to
![]()
and
![]()
which are now in lumped parameterized form.
3.2.3 Rubber elasticity: simple mechanical model
The above theory on rubber elasticity of swelling gel has been developed by Flory based on complex statistical theory which remarkably relates how the molecular structure of elastomer responds to an applied strain to the macroscopic deformation behavior. However, this theory still offers only a qualitative approximation to the actual behavior of elastomer, mainly due to the assumption that end-to-end distances of the chains can be described by the Gaussian distribution (Young and Lovell). In addition, parameters such as effective number of chains in the gel are difficult to obtain. Given the complexity of the statistical model, the purpose of this section is to approximate the gel's elastic behavior with a simpler mechanical analog.
Let’s start our analysis with general equations of elasticity for the case with normal stresses only;
![]()
![]()
![]()
where ei is normal axial strains in i-axis (i= x, y, z),
si is normal axial stress in i-axis,
E is elastic modulus of elastomer, and
n is Poisson’s ratio.
(Note the orientation independent common values of E and n for all three axes implies that the material is homogeneous and isotropic.) First, we need to define our control volume, whether it should be just the elastomer or the whole gel containing both the elastomer and solvent. In the case where the whole gel is the system, the above relations may take rather simple forms since the osmotic pressure (swelling pressure) can directly replace the axial stresses, and the strains are proportional to length changes. In addition, the applied stress in x-axis is simply the force divided by the cross-sectional area facing x-axis. However, in this approach, the material changes continuously as the gel undergoes swelling along with its system parameters. Although this approach may be appropriate for analysis of small differential swelling effect, it seems more suitable to take the "elastomer only" approach for analysis of finite large swelling effect.
Consider a dry, stress-free gel of linear dimensions Lxdo, Lydo and Lzdo, isotropically swollen and stretched to dimensions Lx, Ly and Lz as shown in Figure 3.5. Assume the elastomer is made up of elastomer units aligned in principal axes, x, y, z, and its volume, Vp, is assumed to be constant. Assuming the volume occupied by the intersecting joint region of elastomer units is relatively small, the constant elastomer volume assumption may be expressed as
![]()
where Aid is cross-sectional area of deformed elastomer unit aligned in i-axis. Since
Vp = AxdoLxdo = AydoLydo = AzdoLzdo
above expression may also be written as
![]()
In addition, define the elongation factors ax, ay and az such that ax = Lx/Lxdo, ay = Ly/Lydo, az = Lz/Lzdo. Furthermore, it has been observed with actual gel under uniaxial tensions that ax does not vary significantly from ay and az. Therefore, we approximated the length ratios and the cross-sectional area ratios in all three axes to be similar in magnitudes, i.e.,

Substituting this to the previous expression, we arrived at
and
(i=x,y,z)
where Ai is cross-sectional area of deformed gel in i-axis. (Note that the deformed elastomer units’ lengths are also equal to the linear dimensions of the polymer gel as well.) The normal

Figure 3.5: Simple mechanical model of gel's rubber elasticity effect. The gel is consisted of elastomer units aligned in principal axes and the interstitial solvent. The elastomer units are assumed to maintain constant volume, and deformation effects in three axes are coupled by Poisson effect at the intersecting joints. Uniaxial tension in x-axis only is consider, and the net effect of the rubber elasticity is to change the osmotic pressure of solvent within the control volume.
stresses and strains of the elastomer units in x, y and z directions due to swelling and uniaxial tension in x-axis are then
,
, ![]()
,
, ![]()
where P is solvent pressure inside the gel, Po is ambient bath pressure, and Fx is tensile force in x-axis. Substituting above results to the general elasticity equations, we have
.
Treating the effect of rubber elasticity as purely mechanical, the change in chemical potential of solvent due to the rubber elasticity is then
. In addition,
where vso is molar volume of pure solvent. Combining above results, we have
.
It should be noted that n here is not same as the Poisson's ratio in conventional sense. It represents the overall mechanical coupling effect due to the coupling only at the intersecting joint region (not of solid elastomer for which n = 0.5), and we will refer to it as the "Poisson-like" ratio. In addition, one should note the limitations of this simple mechanical model, namely the lack of thermo-mechanical coupling --- the entropy spring effect. With real elastomer, there is thermal-mechanical energy conversion, as Flory's model describes in terms of entropy change due to deformation. This energy coupling behavior is expected to be unimportant in our isothermal, isobaric system. However, the importance of this coupling effect needs to be investigated for other systems.
3.2.4 Overall process: mixing/diffusion and deformation
Combining the results on mixing/diffusion process and elastic deformation process, i.e., Dms=(Dms,M+Dms,el), the bond graph representation of overall process is shown in Figure 3.6. Note that for the surrounding solvent, hso = mso + Tsso.
3.2.5 Polymer gel swelling kinetics
Just as the fluid flow rate through a restricted pores as shown in Figure 3.7 can be characterized in terms of the net pressure difference between the inlet and the outlet, the solvent flow rate into and out of the gel during swelling and deswelling processes was

Figure 3.7: The gel's solvent's flow kinetic driven by chemical potential difference is analogous to the fluid flow through restrictive pores driven by hydrostatic pressure difference.
phenomenologically described in terms of solvent chemical potential difference across the system boundary;
![]()
where Pm represents the permeability of the solvent in polymer gel network. Recognizing that it may also be expressed in terms of the osmotic pressure difference as
, dependencies of Pm on critical parameters such as gel thickness, cross-sectional area, and swelling degree can be explicitly defined by comparing with Darcy’s law for flow through porous media,
![]()
where U is the macro-scale "superficial velocity", defined as the flow rate Q divided by the total (solid plus fluid) cross-sectional area,
P is the macro-scale pressure,
Km is the permeability
mvis is the solvent’s fluid dynamic viscosity.
For simplified case where the pressure is a function of z-axis only, Darcy’s law may also be expressed in terms of solvent’s molar flow rate as
![]()
where Dz is the gel’s characteristic thickness solvent travels during swelling, and subsequently, we arrived at
![]()
Furthermore, the gel can be treated as a collection of essentially cylindrical polymer fibers aligned either parallel or perpendicular to the zero Reynolds number solvent flow as shown on Figure 3.8. It is assumed the fiber matrix consists of one-third parallel fibers and two-third perpendicular fibers, and an approximate analytic expression was obtained for the permeability of a periodic porous matrix (Happel and Brenner, Low Reynolds Number Hydrodynamics),


where a is the radius of the fiber and b is the fluid envelope radius (which is half the distance between fibers). The net effect on Km is then
![]()
Since

the dependence of permeability coefficient Km on swelling degree can be easily defined. Note that the only adjustable parameter is the cylindrical fiber radius, a, and according to sedimentation experiments (Ogston et. al.; 1973), the fiber radius of linear polymer is expected to lie between 0.3 nm and 0.8 nm.

Figure 3.8: Happel and Brenner model of flow through porous media. The gel can be treated as a matrix of essentially cylindrical polymer fibers aligned either perpendicular or parallel to the zero Reynolds number solvent flow (parallel cylinders not shown in above diagram), and Happel and Brenner's model can be used to predict the permeability of interstitial solvent.
3.3 Dynamic Model of Electrolyte Gel

Figure 3.9: Control volume for electrolyte gel model. As with the non-electrolyte gel case, the gel's control volume changes with the material volume where the gel is consisted of polymer network, interstitial solvent fluid and mobile ions.
In addition to the modeling assumptions for the non-electrolyte gel model, let us also assume: (1) the characteristic length of gel element is much greater than the Debye length, (2) ionic concentrations inside the gel are uniform at any instant of time, and (3) ions behave as in ideal dilute solution;
where mi,e is the electrochemical potential, xi,g is the mole fraction, zi is the valence and fi,g is the electrical potential of the ion i inside the gel.
3.3.1 Equilibrium conditions
The Donnan membrane equilibrium phenomena is due to the balance between the diffusion tendencies toward uniform concentrations of each individual chemical species and the electrical attraction/repulsion among ions/fixed charge groups. Since the equilibrium condition requires the ions in both inside and outside the gel/membrane to have equal potential while satisfying the electroneutrality condition, electric potential develops across the gel/membrane interface. For example, diffusion tendency of the positively charged ion in Figure 3.10 tries to maintain equal concentration of the positive ion between the bath and the gel. Meanwhile the fixed negatively charged side groups of the polymer network attracts the positive ions from the bath, causing the positive ion's concentration to increase in the gel. As the result of this tug-of-war between the two forces, an electric potential is developed across the interface. It should be noted that the macroscopic electroneutrality is valid for dimensions much greater than the Debye length (Grodzinsky; lecture note).
Consider the ionic gel immersed in ionic solvent bath. Safely assuming the concentration of the mobile ions in both the gel and the solvent bath to be much smaller than the concentration of polymer or solvent, the mobile ion concentration can be treated as in dilute solution, and the electrochemical potential of the ions can be described by
![]()
where
is the electrochemical potential of the ion specie i [J/mol],
is the chemical potential of the pure ion specie i,
R is the universal gas constant, 8.314 J/(K*mol),
T is the temperature [K],
ai is the activity coefficient of the ion,
zi is the ionic valence,
F is the Faraday constant, 9.648 x 104 coul/mol,
is the electrical potential of the ion [volt].
When the electrolyte in both the gel and the bath is in true equilibrium, the electrochemical potential of the ion in the bath will be equal to the electrochemical potential of the ion inside the gel,
![]()
where
is the electrochemical potential of the ion inside the gel, and
Since
is same for both the gel and the bath phases, we have
![]()
or
![]()
where
is
, the electrical potential difference across the gel-bath interface. Since this electrical potential difference is imposed to all ions, the activity coefficient ratio across the interface of different ions, say i and j, are related by

For the 1:1 electrolyte where zi = 1 and zj = -1, the above expression reduces to
![]()
For the case where the activity coefficients in the gel and bath phases are equal or for the ideal dilute solution case, we can also write in terms of the mole fractions or the concentrations
![]()
or
![]()
where xig, xib , cig, cib are the mole fractions and concentrations [mol/m3] of the ion i in the gel and the bath, and xjg, xjb , cjg, cjb are the mole fractions and concentrations [mol/m3] of the ion j in the gel and the bath, respectively. In addition, the bulk electroneutrality for the bath and the gel are
for external bath
for inside the gel/membrane
where
is the fixed charge group density [coul/m3], and
is the ion concentration inside the gel/membrane.
The Donnan equilibrium and electroneutrality conditions, along with charge conservation relations, can be used to compute individual ion concentrations. Strictly speaking, the above result applies only to true equilibrium case. However, in general, it can also be applied to non-equilibrium case as long as the net transport is much slower than the diffusion or the drift components of the transport process (Grodzinsky). In the following section, lumped parametrized dynamic elements representing the chemical and electric potential capacitances are developed.
3.3.2 Chemical capacitance model

Figure 3.11: Ideal dilute solution. Since the concentrations of ions in the gel are much less than the effective concentrations of either the solvent or the polymer network, the ions are modeled to be in ideal dilute solution.
With the assumption of dilute solution, the chemical potential of ion in the gel may be described by
.
Treating the polymer molecule as solvent/ion size molecule with equivalent scaling factor, s, the effective mole fraction of ion, xi is then
![]()
and

Then, the rate of change of ion’s chemical potential with respect to time may be approximated by

Since the chemical potential of ion can be changed by either the ion transport or the solvent flow into the control volume, the ion's chemical potential needs to be modeled as a coupled two port capacitance as shown in Figure 3.11 where Cc,i represents the chemical potential capacitance of ion specie i.
3.3.3 Electric potential capacitance model
Bulk electroneutrality
The ions with greater diffusivity will tend to move across the medium at faster rate than the ions with lower diffusivity. However, such tendency for charge separation induces local electric field among ions, and it in turn acts to retard the faster ions and to accelerate the slower ions. For instance, if the faster ions are positively charged whereas the slower ions are negatively charged, then positive local electric potential gradient will build up to retard the positive ions and to accelerate the negative ions. Such a "local" electric field, which acts at scale of a few Debye length, is quite large and thus it is rather difficult to induce charge separation in macroscopic scale. Therefore, at macroscopic scale, the electroneutrality prevails.
The electrical potential coupling among ions due to such bulk electroneutrality effect can be modeled as a simple ideal parallel plate capacitor (Oster)
, ![]()
![]()
where
is the induced electric potential,
CE,i is the electric potential coupling capacitance,
is the dielectric constant of the electrolyte,
A is the cross-sectional area of the ion transport [m2],
d is the typical maximum separation of charges in a volume element, ~ 10-10 m,
qi is the charge of the ion [coul].
The bond graph representation is shown in Figure 3.13 where CE,i represents the bulk electroneutrality potential capacitance of ion species i. Note that the "local" electric potential acts in the scale of a few Debye length in strict sense, and for mediums with much greater scale,
here only models its approximate net effect on the bulk’s electric potential. Furthermore,
is modeled as a transient effect which diminishes at an equilibrium even though the local electric field among ions would continue to exist in microscopic scale.
Macroscopic electric potential gradient

Figure 3.14: Macroscopic electric potential gradient due to difference in fixed charge concentrations. As with the Donnan potential across gel/membrane interface, an electric potential gradient can build up even between two regions of gel if their fixed charge concentrations differ from each other.
Although the bulk electroneutrality requires zero net background charge density in each region of electrolytic medium, this does not require total net charge, ru, in a region to be nonexistent as well. As an illustration, consider a monoprotic binary ion case as explained by Grodzinsky. Defining the total concentrations of charges c+ and c- as sums of background and perturbation components,
![]()
the total net charge, ru, in a region can be written as
.
The bulk electroneutrality requires zero net background charge density
![]()
and therefore, the net charge, ru, is
![]()
which is much less in magnitude than that of background densities of either species. According to Gauss’ law which relates the electric field gradient to the net charge at a point,
![]()
existence of such net charge thus induces electric potential gradient at a point while maintaining bulk electroneutrality in the region. (Although one needs to solve the Gauss’ law in order to describe the electric potential profile exactly, it will not be analyzed in detail here due to complexity as well as lack of necessity for the exact electric potential profile inside the gel for our purpose.) For example, the Donnan equilibrium membrane potential is the electric potential gradient that arises across the membrane interface due to net charges induced by the balance between the individual ion diffusion tendencies and the bulk electroneutrality in each region.
In the case of true equilibrium, the electrochemical potentials of the ions are equal across the two regions, and as discussed previously, the resulting electric potential difference between the regions is described by
![]()
for any ion species, i or j. Although the system of our interest is not at true equilibrium during the process, we will assume that the condition is near equilibrium, and we will model the electric potential difference between two uniform media regions (Figure 3.14) as the average of the individual ion specie’s near equilibrium contribution to the electric potential buildup. Defining the reference electrical potential to be at the electrically neutral surrounding bath, the individual ion contribution to the electric potential,
and the average electric potential
at a region are modeled as
![]()

Overall effect
Combining above results, the total electric potential at a region may be expressed as

or

The rate of change of electric potential is then

Note that applying above result to the surrounding bath, we get

which becomes zero, i.e. reference potential, since it is defined to be electrically neutral.
3.3.4 Kinetic model
The transport of ions due to diffusion and migration in the electrolytic medium is described by the following relation (neglecting the convective contribution)
![]()
where
is the molar flux of the ith species ion [mol/m2*sec],
Di is the diffusion coefficient [m2/sec],
ci is the concentration [mol/m3],
ui is the mobility [m2/volt*sec],
zi is the valence of the ith species,
is the electrical potential [volt].

Figure 3.15: Bond graph representation of ion transport. The ion transport is driven by the combination of chemical and potential gradients. The entropy generation is neglected here.
Modeling both the concentration and electrical potential gradient inside the gel medium as linear and considering only one dimensional case for simplification, the average ion transport relation of ions into the gel can be expressed as
![]()
where
,
are the ion concentration and electrical potential net differences between inside the gel and outside bath, i.e.
, and
is the separating distance. The first term describes the diffusion component due to the concentration gradient, and the second term describes the migration component due to the electric potential gradient. Since the electric potential gradient is commonly experienced by all ions, it is the source of coupling effect among ions. Studying the expression for the ion transport above, one can easily represent this using a bond graph resistive structure as shown in Figure 3.15 where the ion transport is a function of the net electrochemical potential difference, i.e.:
. Note that the above bond graph structure suggests a entropy generation port at the resistance field. However, we assumed it is negligible for simplicity since its effect on the gel's swelling would be minimal.
The ion transport kinetics of the resistive element suggested by the bond graph structure is simply
![]()
Comparing this expression to the generally used ion transport equation

we get the expression for the ion's transport coefficient Ki

For our model, we used the average value of the above Ki's.

Combining above results, the ion transport process between two regions is represented in bond graph formalism in Figure 3.16(a) where CDi, CEi, CDi represent the Donnan membrane, bulk electroneutrality and chemical potential capacitances of ion at a region i. Lumping the individual capacitance components CDi and Cei together into single
capacitance field produces less cluttered views in Figure 3.16 (b) and (c). Note the solvent port at the chemical capacitances were not shown in the figure to reduce clutter.

Figure 3.16: Bond graph representations of ion transport process between two regions. CDi, CEi, CDi represent the Donnan membrane, bulk electroneutrality and chemical potential capacitances of ion at a region i. Lumping the individual capacitance components CDi and Cei in (a) together into single capacitance field produces less cluttered views in (b) and (c).
3.3.5 Ion diffusion limited reaction
As the ion moves into the gel/membrane network, it may also react with the fixed charged side groups as illustrated in Figure 3.17. Namely diffusing H+ or OH- in solvent can react with corresponding conjugate basic or acidic structures on the network, and it was shown that such diffusion-reaction can significantly slow down the overall ion transport process (Grodzinsky). As the ions are moving into the gel network while reacting with the fixed charges of network simultaneously, not only is the overall kinetic of the ion’s transport affected, also the concentrations of both the mobile ions and the fixed charge groups are changed as well, subsequently changing the overall equilibrium behavior. Therefore, proper integration of the dynamic model of ion diffusion limited reaction process is necessary.
Chemical reaction
Assuming the predominant chemical reaction is the acid/base reaction between the H+/OH- and the fixed charge groups, we will limit the discussion to acid/base reaction only. However, the results can be translated to other types of reactions as well.
Consider the case of gel/membrane with negatively charged fixed side groups with pKa below 7. The side group is negatively charged when the internal pH is above its pKa, whereas it is neutral when the internal pH is below its pKa. When the internal pH is equal to the pKa, this implies exactly half of the total side group is dissociated. In acidic environment, there are two reactions occurring simultaneously. First, the neutral side group of gel/membrane network can dissociate into anion and hydronium ion
![]()

Figure 3.17: Diffusion limited reaction of H+ ion. As the H+ ion moves into the gel network, it may react with corresponding conjugate basic structures on the network, and such diffusion-limited reaction can significantly slow down the overall ion transport process.
where AHb is the neutral side group of gel network,
A- is the network-bound anion,
H3O+b is the hydronium ion originally bound to AHb.
The bond graph representation shown below.

The chemical reaction is represented by a two-port resistance field where the thermal (entropy generation) port is not included. In the case where the thermal port is of importance, it can be modified as

Second, the hydrogen ion diffusing in from outside can associate with the network-bound anion, A- to form the neutral side group.
![]()
where H3O+d is the hydronium ion diffusing into the gel network,
AH d is the neutral side group formed.
The corresponding bond graph representation (neglecting the entropy generation) is

and since H3O+d is identical to H3O+b, the above two reactions are identical.
It should be noted that the reaction model by itself assumes that there is no additional creation or depletion of the chemical species. In other words, the specified reactions are the only processes through which the concentration of chemical species can change, and thus the total molar quantity of all species must stay constant. Otherwise, the individual species’ chemical capacitances cannot be specified by the reaction coordinate alone, and they must be modeled as more complicated coupled multiport capacitance fields. Unfortunately, such chemical reactions are rather common, and our case is not an exception; the amounts of H2O and H3O+ within the control volume can vary independently of the acid/base reaction through swelling and H3O+ ion transport, and the bond graph representation of gel system with multiport capacitance fields then would be

In this representation, the additional ports at the capacitance fields are not drawn explicitly due to complexity, but they are denoted as flow variables in parenthesis.
However, we can simplify above model since the molar quantity of H2O is orders of magnitude greater than any other chemical species, including that of H3O+, and thus the effects of other species on the chemical potential are negligible. In this case, the additional ports for H2O flow would be significant at the capacitance fields. In addition, note that since the net charge produced by the chemical reaction is always zero for this particular case, it has no effect on macroscopic electroneutrality. However, it changes H3O+ concentration, and thus the Donnan potential would be affected.
Chemical reaction kinetics
According to the simple mass-action law for elementary chemical reaction expressed as
![]()
where kf is the forward reaction constant, and
kr is the reverse reaction constant,
the reaction kinetic is commonly expressed as
![]()
where the equilibrium constant Keq in this case would then be
![]()
We assumed that the electric potential has negligible effect on the chemical reaction, and for simplicity, modeled both the fixed side groups and the hydrogen ions as dilute solutions, i.e.,

where xAH, xA- and xH+ are mole fractions of AH, A- and H+ respectively. Since concentrations of AH, A- and H+ are much smaller than the concentration of solvent, the mole fractions can be related to corresponding concentrations as
![]()
where vs is the molar volume of the solvent. In addition, the rate of concentration changes of AH is related to the rate of mole of AH change by
![]()
The reaction kinetic expression can then be expressed as

where




![]()
![]()
Quasi-equilibrium chemical reaction
In many cases, simplifications can be made with assumption that the reaction occurs at much faster rate than the diffusion, and thus can be considered as "instantaneous" and "always" at equilibrium (quasi-equilibrium). This does not imply that the concentrations of A-, H+ or AH do not vary with time, but they change at a rate much slower than the reaction rate, limited by the diffusion of the hydrogen ion into the gel/membrane unit. The quasi-equilibrium reaction model says at every instant of time
![]()
or

Since the total possible number of the binding site -A is fixed
![]()
the equilibrium constant Keq can also be written as

Differentiating this with respect to time, we get

For such quasi-equilibrium reaction, the bond graph representation of the diffusion-reaction (electric potential port not shown) are shown in Figure 3.18. Note the bond graph model shows that only three out of total four capacitances have integral causality implying only three capacitances are independent energy storage elements whereas the other one is dependent. This result is in agreement with the above analysis since the AH and A- concentrations are dependent of each other. As before, the entropy generation is assumed to be unimportant in the overall dynamics of the ion’s transport for the quasi-equilibrium chemical reaction model.
3.3.5 Overall ion transport process: binary solution case
When an ionic gel with thickness that is much greater than the Debye length is immersed in an electrolytic solution, the gel’s surface interface region, which is only a few Debye length thick, reaches the equilibrium state almost instantaneously while the inner bulk regions equilibrate at a much slower rate. Therefore, the interface region needs to be treated separately from the bulk units, and modeling of gel’s ion transport process thus involves at least 3 phases: bath, gel interface and gel’s bulk unit. (For simplification, the external solvent bath may be modeled as a ideal ions/solvent source as well as an electric potential reference ground.)
As an illustration of how the model elements developed so far can be integrated together, the simplest case where a hypothetical negatively charged gel is immersed in the 1:1 binary electrolyte solution (H+, Cl-) is considered, and the bond graph representation of the resulting model is shown in Figure 3.19. The model consists of the gel’s interface and the bulk unit connected in series to the solvent bath. The chloride ion is a completely dissociated ion of the strong acid, HCl, and is assumed not to react or bind with other ions or fixed charges. The dissociation of gel bound groups AH produces fixed charges A- and hydronium ion H3O+ which in turn influence the ion transport in the gel through changing both the H+ concentration and the Donnan electric potential. As mentioned previously, the electroneutrality is not affected by the dissociation chemistry of AH since the valences of produced charges are exactly opposites, but the Donnan potential is affected.
As before, the solvent and ions' chemical capacitances should be modeled as coupled multiport capacitance field in strict sense. However, since the amounts of ions are much

Figure 3.18: Bond graph representation of quasi-equilibrium diffusion-limited reaction of H+ ion. Assuming the reaction to occur at much faster rate than the diffusion, the reaction is considered to be "instantaneous" and "always" at equilibrium. Note only three out of total four capacitances have integral causality implying only three capacitances are independent energy storage elements. This is in agreement with the quasi-equilibrium assumption that the AH and A- concentrations are dependent of each other.
smaller than that of solvent and polymer, we can safely assume that inclusion of only the solvent flow port should suffice in this case as well.
According to the bond graph model, it is interesting to note that the electric potential differences can exist even between the interface region and the bulk unit (or even between bulk units if there are more than one) during the transport process since the internal pH and subsequently the fixed charge density of each region may well differ from each other. Unlike the Donnan potential difference between the bath and the gel’s interface, this electric potential differences would vanish at equilibrium when the uniform fixed charge density is achieved between interface and bulk unit. Note that the governing dynamics of interface unit and the bulk unit are identical.
Furthermore, the region situated just inside of the gel’s surface by a few Debye lengths reaches the equilibrium states almost instantaneously when compared to the rest bulk of the gel, and its equilibrium states may be treated as the boundary conditions for the bulk

Figure 3.19: Bond graph representation of ionic gel immersed in aqueous bath with H+ and Cl- ions. The gel is consisted of the surface interface region which equilibrates almost instantaneously and the much slower bulk unit connected in series. The H+ ion also undergoes chemical reaction with the fixed charge groups of network.
unit. In such case, the above representation further simplifies to representation in Figure 3.20.
Governing equations
Corresponding governing dynamic equations are
![]()

where

![]()
![]()



According to above dynamic equations, the time constants of hydronium and chloride ions’ transport processes are
.

Figure 3.20: Bond graph representation of ionic gel immersed in aqueous bath with H+ and Cl- ions. Since the interface region equilibrates at several orders of magnitude faster than the bulk unit, the equilibrium states of the interface unit can be treated as the boundary conditions for the bulk unit.
When the diffusing hydronium ion is reacting with the fixed charges of the gel whose pKa is within the gel's active pH ranges, the transport of hydronium ion will be the limiting process of overall ion transport kinetics since the term in the bracket of hydronium ion’s time constant would be much greater than unity. For the complete model of polyelectrolyte gel, above model needs to be incorporated into the model representing the swelling behavior of the gel as described in the following section.
3.3.6 Ionic contribution and electrolyte gel swelling
Mixing/diffusion process
Similarly to the non-electrolyte gel model, the quasistatic differential energy change of the control volume due to the mixing/diffusion within the gel itself can be expressed as
![]()
Since the ion concentration inside the gel is much less than the concentrations of either the polymer or the solvent, the contribution of ions to the change in free energy is assumed to be additive (Flory), i.e.,
![]()
![]()
![]()
Note the first term in the bracket of right hand side is the
for the non-electrolyte gel, and since
, above expression reduces to
Similarly, the entropy change due to mixing is
![]()
or
![]()
where the ionic contribution to enthalpy change is assumed to be zero and the second term on right hand side represents the entropy generated due to mixing of ion within the gel. In addition, the ionic contribution to the chemical potential change of the solvent is also assumed to be additive (due to sufficient dilution of ions)
where
(van’t Hoff equation)
The solvent flow kinetic will be identical to the non-electrolyte gel case.
Bond graph representation
The bond graph representation of polyelectrolyte gel’s mixing process derived above can be expressed as in Figure 3.21, and incorporating with the ion transport model and applying for the case of binary electrolyte solution, we arrive at the bond graph as shown in Figure 3.22.
Limiting process
Even with the simple binary ion case and a single bulk unit, we have four coupled independent dynamic equations describing the overall dynamic process: two for each ion specie’s chemical capacitances, one for the electric potential capacitance and one for the solvent’s chemical capacitance. In order to be useful for real-time control, it is evident the model needs to be simplified further, and fortunately, further simplification can be made. Grimshaw et. al. has shown that for some gels, the swelling rate is limited by the diffusion-limited reaction of H+ ions with the gel’s charge, and according to De Rossi, such is the case for the gel of our interest, poly(vinyl Alcohol)-poly(acrylic acid) (PVA-PAA), when the pH of the solution (modulated by HCl in H2O) in equilibrium with the gel is higher than 2.9. As expected, De Rossi showed the overall dynamics of the gel is of first order with time constant,
.
Therefore, our strategy is to first identify the limiting process for the gel and solution of our interest, and reduce the overall dynamics to the first order coupled dynamic equation.

Figure 3.21: Bond graph representation of polyelectrolyte gel's mixing process. Since the ion concentration inside gel is small, the contribution of ions to the change in free energy of gel is assumed to be additive to the solvent mixing/diffusion and the rubber elasticity effects.

Figure 3.22: Bond graph representation of negatively charged electrolyte gel in binary electrolyte solution. The interface unit's equilibrium conditions are treated as the boundary conditions for the bulk unit.
Chapter 4
Model Verification
4.1 PVA-PAA in Acetate Buffer Solution
In order to verify the lumped parameter model of the polymer gel, a poly(vinyl Alcohol)-poly(acrylic acid) (PVA-PAA) gel in acetate buffer bath solution was used. PVA-PAA gel was selected for several reasons. First, the PVA gel’s polymeric characteristics (from which the PVA-PAA gel's characteristics were estimated for our purpose) have been investigated extensively, and the parameters such as the polymer-solvent interaction parameter has been well characterized. Second, the PVA-PAA gel is a univalent gel with carboxylic charge groups, and it exhibits continuous swelling with respect to pH within the active range. On the other hand, gels such as the poly(acrylonitrile) (PAN) gel fibers possess excellent mechanical properties, but exhibit severe hysteresis on swelling behavior. Such hysteresis is believed to be due to the gel’s amphoteric nature. Third, the PVA-PAA gel is relatively easy to prepare and cast into desired shape. The gel preparation procedure will be described in detail later. Lastly, researchers such as De Rossi have already studied the PVA-PAA gel, and parameters such as the gel’s pKa have already been found (Chiarelli and De Rossi, 1988).
Although most of the studies on the PVA-PAA gel have been conducted in aqueous solution, where pH was modulated by addition of hydrochloric acid (HCl), such solution environment is not appropriate for controlling the swelling degree of the gel. Because the HCl is a strong acid, it completely dissociates instantaneously, and the pH level of the solution will be extremely sensitive to the amount of HCl added to the water, especially at near the neutral range. Modulation of the solution's pH to the desired level precisely will be difficult, and subsequently, controlling the swelling degree of the gel would be difficult as a result. Therefore, we have opted to use the acetate buffer solution with its pH adjusted to neutrality (pH 7) with sodium acetate, and during the experiment, the solution's pH is subsequently modulated to desired value with HCl. The acetate buffer solution resists changes in pH by the ability to combine with H+ ions, and the desired pH level can be therefore accurately achieved. In this case, note that four monovalent ions (rather than two of HCl aqueous solution) are involved in gel swelling process, and they are hydronium ion (H+), sodium ion (Na+), chloride ion (Cl-), and acetate ion (C2H3O2-).
4.1.1 Acetate buffer chemistry
A buffer is a solution characterized by the ability to withstand changes in pH when limited amounts of acid or base are added to it, and in this section, the acid/base chemistry of the acetate buffer solution is reviewed with understanding that the knowledge of the ion concentrations in the bath solution is critical in predicting the swelling of the ionic PVA-PAA gel.
Acetate buffer solution preparation
Given pure acetic acid and sodium acetate, we first need to figure out what concentrations of each chemical should be mixed together to produce the acetate buffer solution with desired pH of 7. Consider the acid ionization of the acetic acid, HC2H3O2
![]()
where the acetic acid and the sodium acetate with initial concentrations A and C are mixed in water;
![]()
Upon mixing, suppose x mole of the acetic acid is dissociated. As a result, the concentration of the acetic acid would decrease by x mole, and the concentrations of the hydronium and acetate ions would increase by x mole respectively.
![]()
Substituting the final concentrations into the equilibrium equation

we get
![]()
Assuming x << A, C for simplicity (which should be checked), above expression further reduces to
![]()
Acid ionization constant, Ka, of acetic acid is 1.7x10-5. Therefore, in order to produce the acetate buffer solution with pH = 7 ([H3O+] = 1x10-7), the initial concentration ratio of sodium acetate to acetic acid should be 170. In our case, we have chosen the initial concentrations of acetic acid and sodium acetate to be 0.001 M and 0.17 M respectively. (Note that both A and C are much greater than x, thus satisfying the assumption.) At equilibrium, the resulting ion concentrations were

pH modulation: lowering pH
As HCl is added to the acetate buffer solution to modulate the bath pH, the equilibrium ion concentrations and pH can be calculated given the added amount of HCl following an analysis similar to the one illustrated above.
Suppose VH in volume of Ho molar HCl is added to the pH 7 buffer solution of volume Vo with equilibrium ion concentrations of A and C for acetic acid and acetate ion, respectively
![]()
For the purpose of analysis, the overall chemical reaction can be broken into two separate steps: complete reaction between strong acid and buffer ion, and acid ionization equilibrium. The added HCl is a strong acid, and it will completely dissociate instantaneously. Therefore, the newly formed hydronium ions will then completely react with the acetate ions to form acetic acid, and this initial reaction between H+ and C2H3O2- is
.
The new total volume of the solution is Vo+VH, and consequently, the concentration of the H+, C2H3O2- and HC2H3O2 are changed to
![]()
Following the complete reaction between H+ and C2H3O2-, the acid ionization process occurs. As before, suppose x mole of the acetic acid is dissociated due to ionization. As a result, the concentration of the acetic acid would decrease by x mole, and the concentrations of the hydronium and acetate ions would increase by x mole respectively.
![]()
Substituting the final concentrations into the equilibrium equation

we get

Assuming
![]()
for simplicity, above expression then reduces to

Therefore, in order to achieve the desired hydronium concentration of x, the amount of HCl necessary is
![]()
pH modulation: raising pH
Driving the polymer gel based actuator to contract and expand requires changing the pH level of bath solution, and so far we have only discussed adding HCl to the buffer solution which lowers the pH level. In order to raise the pH level to cause the gel to expand, we can replace a portion of the bath solution with the neutral buffer solution, thereby effectively lowering the H+ concentration (and thus raising pH). Such method, however, wastes excessive fluid, and definitely is not a viable option for a practical actuator. Instead, a more effective method would be to neutralize the acid with a strong base such as sodium hydroxide (NaOH). However, this option has its own disadvantage. In the process of neutralization, salt (NaCl) and water are produced,
![]()
and the ion concentrations in the bath solution are changed. Therefore, the effect of the acid-base neutralization process on the bath's ion concentrations need to be considered. Fortunately, its effect can be easily accounted for since the sodium hydroxide completely reacts with the acid, and overall effect is simply reduction of H+ and increase of H2O and Na+ by same number of mole of NaOH added. With an appropriate adjustment in initial reactant concentrations, analysis described in the previous section can be applied to compute the equilibrium ion concentrations.
4.1.2 PVA-PAA gel preparation
The PVA-PAA gel strips were prepared by drying the PVA-PAA solution on a Lexan plate into a thin film and baking in an oven at 120 C for 30 minutes. Molecular weights of the PVA and PAA used were 31,000-50,000 and 90,000, respectively, and the PVA-PAA solution was prepared by mixing 60% and 40 % by weight of PVA and PAA in distilled water. Lexan has been chosen as the drying plate because dried PVA-PAA film was found to peel off easily. On the other hand, PVA-PAA film tends to stick to glass and Plexiglass. The dried film was cut into desired shapes and sizes before they were baked, since baking causes the film to become brittle. The baked PVA-PAA strips were soaked in distilled water for 48 hours, and subsequently soaked in the acetate buffer solution for 24 hours. It should be noted that the stiffness of the gel was extremely sensitive to the baking temperature and duration, and even the difference of 15 seconds in baking time yielded noticeable difference in the degree of swelling.
4.2 Equilibrium swelling model verification
4.2.1 Equilibrium swelling model
The equilibrium swelling of the polymer gel is such that the chemical potential of the solvent (which is water in this case) in the surrounding bath is identical to that of the solvent in the gel
![]()
while also maintaining the Donnan equilibrium and bulk electroneutrality conditions of the gel.
According to our previously developed model of the electrolyte polymer gel, the chemical potential of the solvent in the gel is
![]()
where


![]()
and the chemical potential of the bath's solvent may also be approximated by
![]()
The Donnan equilibrium condition and the bulk electroneutrality condition impose

![]()
and the diffusion reaction chemistry of the fixed charge also requires

The ion concentrations of bath solution are computed based on the analysis presented in the previous section, and since the concentrations of the mobile ions and the fixed charge inside the gel also depend on the swelling degree of the gel, the final equilibrium state of the gel is such that the above expressions are simultaneously satisfied.
4.2.2 Experimental procedure
Constant tension test
Baked thin PVA-PAA gel films of dry dimensions Lydo = 4.00 mm, Lzdo = 0.103 mm were first washed in distilled water for 48 hours, followed by immersion in the buffer solution of pH 7 for 24 hours to ensure that they have reached their equilibrium states. Once the equilibrium state had been reached, the terminations were mounted at each end of the gels as shown in the diagram below. The termination blocks were made of Lexan, and the PVA-PAA gel film ends were clamped into the slots of termination blocks.

A constant tension of 6.098 x 10-2 N was applied by hanging a weight of 6.207 gram. Bath solutions with various degrees of HCl concentrations were prepared, and the PVA-PAA gel film samples in constant tension were immersed in each bath solution for 24 hours, at which point the dimensions in x and y directions were measured with caliper. The bath temperature was periodically monitored to ensure that the isothermal assumption were met. The acetate buffer to one molar HCl volume ratios tested were 100:2, 100:4, 100:6, 100:8, 100:10, 100:13 and 100:16, respectively. The corresponding pH's of the bath were 5.62, 5.27, 5.03, 4.82, 4.61, 4.25 and 3.56, respectively.
Salt content test
To test the effectiveness of the electrolyte gel model in the presence of high salt content in the bath solution, NaCl has been added to bath solutions of pH units 5.27 and 4.82 such that the salt concentrations were 0.2 M and 0.5 M for each pH unit. The PVA-PAA gel film samples were immersed in each solution for 24 hours with no tension applied, and the dimensions were measured with caliper.
4.2.3 Experimental result
The averaged results from the constant tension test and the salt content test are tabulated below:
Buffer : HCl Ratio
|
100:2 |
100:4 |
100:6 |
100:8 |
100:10 |
100:13 |
100:16 |
|
|
a x |
2.05 |
2.00 |
1.97 |
1.94 |
1.88 |
1.77 |
1.66 |
|
a y |
1.95 |
1.88 |
1.83 |
1.79 |
1.72 |
1.62 |
1.53 |
Constant tension experiment (Tension = 6.098 x 10-2 N)
Buffer : HCl Ratio
|
100:4 |
100:8 |
|
|
0.2 M NaCl |
1.87 |
1.78 |
|
0.5 M NaCl |
1.81 |
1.73 |
Length ratios at various salt and HCl content.
The length ratios in x and y directions from the constant tension experiment were plotted against the bath pH in Figure 4.1. The standard deviations of the experimental data ranged from 0.007 to 0.019 for x direction and from 0.012 to 0.035 for y direction, respectively, and the instrumentational measurement errors in length ratios for x and <