PDF
Large Swelling in Polymer Hydrogels
Introduction
Polymer hydrogels consist of a crosslinked network of long-chained polymers and a large number of imbibed solvent molecules. A wide range of hydrogels can be produced by controlling the chemistry of the polymer network, its crosslinking, and the type of solvent. This class of materials is therefore used in a wide range of applications, such as medical devices, targeted drug delivery, design of tissue-like materials, and in stimuli-sensitive actuators (Ref. 1).
Polymer gels can undergo large deformations as a result of mechanical loads as well as changes in the chemical potential of the surrounding solvent. At time scales much shorter than that typical of diffusion of the solvent molecules, the mechanical response is isochoric, allowing only changes in gel shape. At longer time scales, however, the solvent molecules have sufficient time to rearrange and the gel can thus undergo large changes in both shape and volume.
This example model demonstrates how to implement a coupled nonlinear field theory of large deformation and solvent diffusion for polymer hydrogels in COMSOL Multiphysics® based on the work in Ref. 1. The implementation is validated by a comparison with two benchmark problems: inhomogeneous swelling of a gel bonded to a rigid substrate (Ref. 2) and uniaxial consolidation under creep loading (Ref. 1).
Theoretical Background
The gel is considered as a homogenized, continuous, single-phase material. The presence of the solvent is described by the molar concentration c(X, t), which denotes the number of moles of solvent per unit reference volume of the polymer gel at position X and time t. As is common, the reference configuration is taken to be the dry polymer network (that is, c(X, 0) = 0). The gel is regarded as an open system since it can exchange solvent molecules with the environment. The concentration c must therefore obey the local form of the volume balance
(1)
Herein, denotes the molar flux of solvent (SI unit: mol/(m2·s)), Q is the volumetric flux, and Vm is the molar volume (SI unit: m3/mol) of the solvent. When inertia is neglected, the local momentum balance takes the form
(2)
with P being the first Piola–Kirchhoff stress in the gel, ρS the mass density of the polymer network, Mm the molar mass of the solvent, and B the mass-specific body force. In the following, body forces are neglected.
The free energy density of polymer gels is often split additively into two contributions (Ref. 1): one related to stretching of the polymer chains and one related to the mixing between the polymers and the solvent molecules, that is,
(3)
Herein, the free energy of stretching is given as
(4)
and the free energy of mixing takes the form
(5)
In these expressions, G is the initial shear modulus of the polymer network, I1 is the first invariant of the right Cauchy–Green deformation tensor C = FTF, J = det(F) is the volume ratio, Kmix is the bulk stiffness of mixing, and χ is a dimensionless constant related to the enthalpy of mixing. Note that the free energy density is expressed per unit reference volume and is written as a function of the deformation gradient F and the solvent concentration c.
For most hydrogels, it is natural to assume that volume changes are exclusively due to changes in solvent content, whereas the polymers and the solvent themselves are incompressible. This constraint establishes a direct link between the solvent concentration and the volume ratio, that is,
(6)
Because of this constraint, it is convenient to introduce another form of the free energy density by a Legendre transform, see Ref. 2. Let
(7)
with denoting the chemical potential of the solvent in units of J/mol. By invoking thermodynamic restrictions, the expression for the first Piola–Kirchhoff stress, the solvent concentration, and the solvent flux read
(8)
where the second-order tensor M (SI unit: m4/(N s)) is called the mobility tensor, which is required to be symmetric and positive definite. For convenience, we have also introduced the chemical potential in units of pressure, . Note that the constitutive law for the solvent flux is analogous to Darcy’s law in porous media; in fact, we can identify the material permeability tensor κ = μfM (SI unit: m2) given the solvent viscosity μf. In Ref. 1, an expression for the spatial solvent mobility tensor is derived based on the solvent diffusivity D, which is assumed to be isotropic in the current configuration. After converting the units, applying a pull-back operation to the material frame, and using the molecular incompressibility constraint, the material mobility tensor analogous to Equation (30) in Ref. 1 reads
(9)
Here, NA and kB denote the Avogadro number and the Boltzmann constant, respectively. The free energy density (Equation 7) and the local form of the volume balance (Equation 1) can now be expressed entirely in terms of the new dependent variables. By using the molecular incompressibility constraint, we obtain the free energy density
(10)
Furthermore, inserting the constitutive equation for the solvent flux (Equation 83) into the local volume balance (Equation 1) and making use of the molecular incompressibility constraint yields
(11)
Finally, the first Piola–Kirchhoff stress follows by inserting Equation 10 into Equation 81
(12)
Thus, all equations have been formulated in terms of F and μ as dependent variables instead of F and c. This is the natural choice of dependent variables in terms of boundary conditions because the chemical potential is continuous across boundaries, whereas the solvent concentration is not. Instead, once the boundary value problem is solved, c can be computed from the volume ratio J using Equation 6.
Reference Configuration for the Finite Element Implementation
So far, the dry polymer network has been used as a reference configuration. In practice, however, hydrogels always contain a large amount of solvent. Moreover, the dry state is problematic in a numerical implementation because the free energy density of mixing is singular in the limit J → 1. For these reasons, a new reference state is chosen where the gel is in a state of free swelling; that is, under no external forces and in equilibrium with a solvent of chemical potential μref, see Ref. 2. Since the gel is homogeneous and isotropic, the free swelling is characterized by the deformation gradient F0 = λ0I. The deformation gradient F can thus be multiplicatively decomposed as F = F 'F0, with F ' denoting the deformation gradient between the free swelling reference configuration and the current configuration. The stretch λ0 can be computed by solving the free swelling boundary value problem, which reduces to a single algebraic equation,
(13)
With λ0 determined, the free energy density and the mobility tensor in the new reference configuration follow by Piola transforms
(14)
with J0 = det(F0) = λ03. The expressions for Ψ' and M' that are needed for the COMSOL Multiphysics implementation are therefore given by
(15)
for the free energy of stretching,
(16)
for the free energy of mixing, and
(17)
for the mobility tensor with respect to the swollen reference configuration. Note that the variables J' = det(F'), C' = F'TF', and I'1 = tr(C') are functions of the deformation gradient F' that maps material line elements from the swollen reference configuration to the current configuration, which is the deformation that we solve for. In the absence of intrinsic inelastic processes, these variables are available in the Solid Mechanics interface as solid.Jel, solid.Cel, and solid.I1Cel, respectively.
Benchmark Problems
To validate the model implementation, two axisymmetric benchmark problems are studied. First, the long-term, stationary solution to the swelling of a cylindrical gel that is immersed in a pure solvent and bonded to a rigid substrate is computed and compared with results reported in Ref. 2. The aspect ratio of the gel is varied in a parametric sweep, and the results are compared to semi-analytical calculations of homogeneous free and uniaxial swelling in the limiting cases of pillar-like and disk-like gel geometries, respectively. For wide and short (that is, disk-like) gels, the lateral edge of the gel deforms severely and comes in contact with the rigid surface. Although only the long-term, stationary solution is of interest in this benchmark, a time-dependent study is used to demonstrate how contact-dependent boundary conditions for the solvent flux can be implemented.
The second benchmark problem concerns the time-dependent uniaxial consolidation of a cylindrical gel under creep compression, see Ref. 1. A cylindrical gel is constrained radially, and a load is applied at the top while the bottom boundary is fixed. Free fluid flow is allowed through a porous filter at the top. The coupled problem is solved in a 2D axisymmetric geometry, and the results are compared with those computed in Ref. 1 where the consolidation problem is reduced to a one-dimensional diffusion-like differential equation.
Model Parameters
The model parameters used for the two benchmark problems are given in the following table.
Table 1: Model parameters for the two benchmark problems.
Parameter
Inhomogeneous swelling (Ref. 2)
Uniaxial consolidation (Ref. 1)
 G
40 kPa
40 kPa
 Kmix
40 MPa
40 MPa
 χ
0.1
0.2
 Vm 
6.022·10-5 m3/mol
6.022·10-5 m3/mol
 D
103 m2/s
8·10-10 m2/s
 μref
-2.6376·105 Pa
0 Pa
 λ0
2.0
3.215
Results and Discussion
The swelling ratio J relative to the dry polymer network is shown in Figure 1 for cylindrical gels bonded to a rigid substrate and swollen to reach equilibrium with a pure solvent. The results are displayed for gels with three different aspect ratios (D/H), and they agree well with the deformed shapes reported in Ref. 2, see Figure 7 therein. Note that, most obvious in the case D/H = 10, the free lateral edge deforms strongly by folding over and making contact with the rigid substrate. Furthermore, note that the swelling ratio (and thus the solvent concentration) is inhomogeneous within the gel and along its outer boundaries.
Figure 1: Swelling ratio in a cylindrical polymer hydrogel bonded to a rigid substrate.
Figure 2 shows the change in gel height as a function of the aspect ratio of the cylinder. For thin, pillar-like gels, the gel height approaches that predicted by homogeneous free swelling (green curve), as most of the domain is unaffected by the constraint. In contrast, for wide, disk-like gels (high D/H), the gel is constrained in the in-plane directions and the state of swelling is close to uniaxial (red curve).
Figure 2: Equilibrium swelling ratio in the thickness direction as a function of gel aspect ratio.
Next, consider the time-dependent solution under uniaxial consolidation. In Figure 3, surface plots of the normal stress and the chemical potential are shown in the spatial frame at dimensionless time . Note that the normal stress is constant through the thickness, as required by the momentum balance. At the bottom of the gel, most of the load is borne by fluid pressurization, whereas at the top of the gel, the load is borne primarily by the polymer network.
The stretch in the thickness direction is shown as a function of the normalized reference coordinate for a few selected time points in Figure 4. For easier comparison with Ref. 1 (see Figure 4 therein), the direction of the Z-axis has been reversed, now going from top to bottom. Both the short- and long-term limits as well as the temporal evolution are in excellent agreement with the results reported in Ref. 1.
Figure 3: Normal stress (left) and chemical potential and solvent flux (right) in a gel under uniaxial consolidation at dimensionless time Dt/H0 = 200. Note that the normal stress is constant through the thickness. The chemical potential is maximal at the bottom surface and minimal at the top surface where the load is applied through a porous filter and the free flow boundary condition applies.
Figure 4: Stretch in the thickness direction during uniaxial consolidation as a function of the position along the gel thickness, shown for various dimensionless times. The data are in excellent agreement with the benchmark in Ref. 1 (circles).
Notes About the COMSOL Implementation
This multiphysics model is implemented by combining the Solid Mechanics and the Darcy’s Law interfaces using the Poroelasticity multiphysics coupling for solving the coupled momentum and volume balance equations. Since both the polymer network and the solvent are assumed materially incompressible, the Biphasic poroelasticity model is used. In addition, the Global ODEs and DAEs interface is used to solve the nonlinear algebraic equations for free and uniaxial homogeneous swelling. The geometry and the displacement field are discretized using quadratic Lagrange shape functions, whereas linear shape functions are used for the chemical potential; these choices ensure a consistent interpolation order when computing the stresses.
The free energy density is implemented as a user-defined Hyperelastic material. The load contribution from the chemical potential is added automatically in the Poroelasticity coupling node. Note that the stress variables solid.sij correspond to the so-called effective stress; that is, they do not include the contribution from the pore pressure. The total Cauchy stress corresponding to Equation 12 is available from the multiphysics coupling node in the variables poro1.sij, poro1 being the Name of the Poroelasticity node.
In problems where the contact patch is changing during a simulation, the boundary condition for the solvent flux depends on the contact status. In case of a gel coming into contact with an impermeable surface, a No Flow boundary condition applies for the portion of the gel boundary in contact with the surface, whereas free flow applies elsewhere. Such complex boundary conditions can be implemented with an if-statement embedded in a Pressure node.
References
1. W. Hong, X. Zhao, J. Zhou, and Z. Suo, “A theory of coupled diffusion and large deformation in polymeric gels,” J. Mech. Phys. Solids, vol. 56, pp. 1779–1793, 2008.
2. W. Hong, Z. Liu, and Z. Suo, “Inhomogeneous swelling of a gel in equilibrium with a solvent and mechanical load,” Int. J. Solids Struct., vol. 46, pp. 2382–3289, 2009.
Application Library path: Porous_Media_Flow_Module/Poromechanics/hydrogel_swelling
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D Axisymmetric.
2
In the Select Physics tree, select Structural Mechanics > Poroelasticity > Poroelasticity, Solid.
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies > Time Dependent.
6
Click  Done.
Global Definitions
Geometrical Parameters
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, type Geometrical Parameters in the Label text field.
3
Locate the Parameters section. In the table, enter the following settings:
Name
Expression
Value
Description
height
1[mm]
0.001 m
Gel height
ar
0.2
0.2
Gel aspect ratio
radius
ar*height/2
1E-4 m
Gel radius
elsize
if(ar<=2,0.05*height,0.10*height)
5E-5 m
Element size
Material Parameters
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Material Parameters in the Label text field.
3
Locate the Parameters section. In the table, enter the following settings:
Name
Expression
Value
Description
G
4e4[Pa]
40000 Pa
Initial shear modulus
Kmix
4e7[Pa]
4E7 Pa
Bulk modulus of mixing
chi
0.1
0.1
Enthalpy of mixing
temp
293.15[K]
293.15 K
Reference temperature
Vm
1e-28[m^3]*N_A_const
6.0221E-5 m³/mol
Solvent molar volume
D
1e3[m^2/s]
1000 m²/s
Solvent diffusivity
M0
Vm*D/(N_A_const*k_B_const*temp)
2.4707E-5 m³·s/kg
Solvent mobility
rhof
1[g/cm^3]
1000 kg/m³
Solvent density
muf
1e-3[mPa*s]
1E-6 Pa·s
Solvent viscosity
kappa0
M0*muf
2.4707E-11 m²
Permeability
Model Parameters
1
In the Home toolbar, click  Parameters and choose Add > Parameters.
2
In the Settings window for Parameters, type Model Parameters in the Label text field.
3
Locate the Parameters section. In the table, enter the following settings:
Name
Expression
Value
Description
stretch0
2
2
Free swelling stretch
J0
stretch0^3
8
Free swelling volume ratio
mu_ref
G/J0*(stretch0^2-1)+Kmix*(log(1-1/J0)+1/J0+chi/(J0^2))
-2.6376E5 Pa
Reference chemical potential
mu_sw
0[Pa]
0 Pa
Chemical potential of pure solvent
poro0
1-1/J0
0.875
Initial porosity
t_end
1000[s]
1000 s
Duration
stress
0[Pa]
0 Pa
Applied stress
Interpolation 1 (int1)
1
In the Home toolbar, click  Functions and choose Global > Interpolation.
2
In the Settings window for Interpolation, locate the Definition section.
3
In the table, enter the following settings:
t
f(t)
0
mu_ref
t_end
mu_sw
4
Locate the Units section. In the Function table, enter the following settings:
Function
Unit
int1
Pa
5
In the Argument table, enter the following settings:
Argument
Unit
t
s
Geometry 1
1
In the Model Builder window, expand the Component 1 (comp1) > Geometry 1 node, then click Geometry 1.
2
In the Settings window for Geometry, locate the Units section.
3
From the Length unit list, choose mm.
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type radius.
4
In the Height text field, type height.
Rectangle 2 (r2)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type 1.4*radius.
4
In the Height text field, type 0.1*height.
5
Locate the Position section. In the z text field, type -0.1*height.
Partition Edges 1 (pare1)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Partition Edges.
2
On the object r2, select Boundary 3 only.
3
In the Settings window for Partition Edges, locate the Positions section.
4
In the table, enter the following settings:
Relative arc length parameters
1/1.4
5
In the Geometry toolbar, click  Build All.
Form Union (fin)
1
In the Model Builder window, under Component 1 (comp1) > Geometry 1 click Form Union (fin).
2
In the Settings window for Form Union/Assembly, locate the Form Union/Assembly section.
3
From the Action list, choose Form an assembly.
4
Clear the Create pairs checkbox.
5
In the Geometry toolbar, click  Build All.
Component 1 (comp1)
Set the geometry shape function to Quadratic Lagrange polynomials. We will use the same shape functions for the geometry as for the displacement field.
1
In the Model Builder window, click Component 1 (comp1).
2
In the Settings window for Component, locate the Curved Mesh Elements section.
3
From the Geometry shape function list, choose Quadratic Lagrange.
Definitions
Create a Contact Pair between the free edge of the rigid substrate (source) and the lateral edge of the gel (destination), which can come into contact when the gel swells.
Contact Pair 1 (p1)
1
In the Definitions toolbar, click  Pairs and choose Contact Pair.
2
Select Boundary 4 only.
3
In the Settings window for Pair, locate the Destination Boundaries section.
4
Click to select the  Activate Selection toggle button.
5
Select Boundary 9 only.
Variables 1
Define Variables for the free energy density, solvent mobility, and the permeability tensor. These are provided in a separate text file.
1
In the Model Builder window, right-click Definitions and choose Variables.
2
In the Settings window for Variables, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domain 2 only.
5
Locate the Variables section. Click  Load from File.
6
Browse to the model’s Application Libraries folder and double-click the file hydrogel_swelling_variables.txt.
Solid Mechanics (solid)
1
In the Model Builder window, under Component 1 (comp1) click Solid Mechanics (solid).
2
In the Settings window for Solid Mechanics, click to expand the Discretization section.
3
From the Displacement field list, choose Quadratic Lagrange.
Rigid Material 1
1
In the Physics toolbar, click  Domains and choose Rigid Material.
2
Select Domain 1 only.
Fixed Constraint 1
In the Physics toolbar, click  Attributes and choose Fixed Constraint.
Hyperelastic Material 1
1
In the Physics toolbar, click  Domains and choose Hyperelastic Material.
2
Select Domain 2 only.
3
In the Settings window for Hyperelastic Material, locate the Hyperelastic Material section.
4
From the Material model list, choose User defined.
5
In the Ws text field, type Psi_s+Psi_mix.
6
Click to expand the Advanced section. In the Eeq text field, type 3*G. When simulating contact, the equivalent Young’s modulus is used to estimate the contact pressure. For soft materials like gels, the default value is too high, leading to an ill-conditioned problem.
Contact 1
1
In the Model Builder window, click Contact 1.
2
In the Settings window for Contact, locate the Contact Pressure Penalty Factor section.
3
From the Penalty factor control list, choose Manual tuning.
4
In the fp text field, type 10.
Fixed Constraint 1
1
In the Physics toolbar, click  Boundaries and choose Fixed Constraint.
2
Select Boundary 7 only.
Darcy’s Law (dl)
Now, set up the Darcy’s Law interface in the hydrogel domain to solve the volume balance equation. Set the discretization to Linear to enable a standard P2–P1 interpolation of displacement and pressure.
1
In the Model Builder window, under Component 1 (comp1) click Darcy’s Law (dl).
2
Select Domain 2 only.
3
In the Settings window for Darcy’s Law, click to expand the Discretization section.
4
From the Pressure list, choose Linear.
5
Click to expand the Dependent Variables section. In the Pressure (Pa) text field, type mu.
Fluid 1
1
In the Model Builder window, under Component 1 (comp1) > Darcy’s Law (dl) > Porous Medium 1 click Fluid 1.
2
In the Settings window for Fluid, locate the Fluid Properties section.
3
From the Fluid type list, choose Incompressible.
4
From the ρref list, choose User defined. In the associated text field, type rhof.
5
From the μ list, choose User defined. In the associated text field, type muf.
Porous Matrix 1
1
In the Model Builder window, click Porous Matrix 1.
2
In the Settings window for Porous Matrix, locate the Matrix Properties section.
3
From the εp list, choose User defined. In the associated text field, type poro0.
4
From the κ list, choose User defined. From the list, choose Symmetric.
5
Specify the κ matrix as
KRR
KRZ
KRZ
KZZ
Initial Values 1
1
In the Model Builder window, under Component 1 (comp1) > Darcy’s Law (dl) click Initial Values 1.
2
In the Settings window for Initial Values, locate the Initial Values section.
3
In the mu text field, type mu_ref.
Continuity 1
Disable the default Continuity condition for the contact pair. In this case, the proper contact condition needs to be set up manually.
1
In the Model Builder window, click Continuity 1.
2
In the Settings window for Continuity, locate the Advanced section.
3
Select the Disconnect pair checkbox.
Contact (Swelling)
1
In the Physics toolbar, click  Boundaries and choose Pressure.
2
In the Settings window for Pressure, type Contact (Swelling) in the Label text field.
3
Select Boundary 9 only.
4
Locate the Pressure section. In the p0 text field, type if(incontact_p1,mu,int1(t)).
Free Flow (Swelling)
1
In the Physics toolbar, click  Boundaries and choose Pressure.
2
In the Settings window for Pressure, type Free Flow (Swelling) in the Label text field.
3
Locate the Pressure section. In the p0 text field, type int1(t).
4
Select Boundary 8 only.
Multiphysics
Poroelasticity 1 (poro1)
1
In the Model Builder window, under Component 1 (comp1) > Multiphysics click Poroelasticity 1 (poro1).
2
In the Settings window for Poroelasticity, locate the Poroelastic Coupling Properties section.
3
From the Poroelasticity model list, choose Biphasic.
4
From the pref list, choose Reference pressure level (dl).
Mesh 1
Mapped 1
1
In the Mesh toolbar, click  Mapped.
2
In the Settings window for Mapped, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domain 2 only.
Distribution 1
1
Right-click Mapped 1 and choose Distribution.
2
Select Boundaries 7 and 8 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
In the Number of elements text field, type floor(radius/elsize)+1.
Distribution 2
1
In the Model Builder window, right-click Mapped 1 and choose Distribution.
2
Select Boundaries 6 and 9 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
In the Number of elements text field, type floor(height/elsize)+1.
Mapped 2
1
In the Mesh toolbar, click  Mapped.
2
In the Settings window for Mapped, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domain 1 only.
Distribution 1
1
Right-click Mapped 2 and choose Distribution.
2
Select Boundary 2 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
In the Number of elements text field, type floor(1.4*radius/elsize)+1.
Distribution 2
1
In the Model Builder window, right-click Mapped 2 and choose Distribution.
2
Select Boundaries 1 and 5 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
In the Number of elements text field, type 1.
5
Click  Build All.
Inhomogeneous Swelling
Set up a parametric time-dependent study to compute the inhomogeneous swelling of the constrained hydrogel under a change in the external chemical potential for different gel aspect ratios. Note that the diffusion coefficient has been set to a large value so that solvent equilibrium is achieved within the study duration.
1
In the Model Builder window, click Study 1.
2
In the Settings window for Study, type Inhomogeneous Swelling in the Label text field.
3
Locate the Study Settings section. Clear the Generate default plots checkbox.
Parametric Sweep
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
Click  Add.
4
In the table, enter the following settings:
Parameter name
Parameter value list
Parameter unit
ar (Gel aspect ratio)
0.1 0.6 1.0 2.0 6.0 10.0 20.0 100.0
 
Step 1: Time Dependent
1
In the Model Builder window, click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type range(0,0.1,1)*t_end.
Solver Configurations
Modify the default solver suggestion to improve convergence.
Solution 1 (sol1)
1
In the Model Builder window, right-click Solver Configurations and choose Show Default Solver.
2
Expand the Solution 1 (sol1) node.
3
In the Model Builder window, expand the Inhomogeneous Swelling > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 node, then click Pressure (comp1.mu).
4
In the Settings window for Field, locate the Scaling section.
5
From the Method list, choose Manual.
6
In the Scale text field, type 1e5.
7
In the Model Builder window, click Displacement Field (comp1.u).
8
In the Settings window for Field, locate the Scaling section.
9
In the Scale text field, type height.
10
In the Model Builder window, click Time-Dependent Solver 1.
11
In the Settings window for Time-Dependent Solver, click to expand the Time Stepping section.
12
Clear the Interpolate solution at end time checkbox.
13
In the Model Builder window, expand the Inhomogeneous Swelling > Solver Configurations > Solution 1 (sol1) > Time-Dependent Solver 1 node, then click Direct.
14
In the Settings window for Direct, locate the General section.
15
From the Solver list, choose PARDISO.
16
In the Model Builder window, click Fully Coupled 1.
17
In the Settings window for Fully Coupled, click to expand the Method and Termination section.
18
From the Nonlinear method list, choose Constant (Newton).
19
From the Jacobian update list, choose Once per time step.
20
In the Maximum number of iterations text field, type 10.
21
In the Study toolbar, click  Compute.
Result Templates
Create a Surface Plot to visualize the swelling ratio of the gel in the deformed configuration for three different aspect ratios. Start from the Predefined Plot of the stress.
1
In the Results toolbar, click  Result Templates to open the Result Templates window.
2
Go to the Result Templates window.
3
In the tree, select Inhomogeneous Swelling/Parametric Solutions 1 (sol2) > Solid Mechanics > Stress (solid).
4
Click the Add Result Template button in the window toolbar.
5
In the Results toolbar, click  Result Templates to close the Result Templates window.
Results
Swelling Ratio
1
In the Settings window for 2D Plot Group, type Swelling Ratio in the Label text field.
2
Click to expand the Title section. From the Title type list, choose Manual.
3
In the Title text area, type Surface: Swelling volume ratio (1).
4
Locate the Plot Settings section. From the View list, choose New view.
5
Clear the Plot dataset edges checkbox.
Surface 1
1
In the Model Builder window, expand the Swelling Ratio node, then click Surface 1.
2
In the Settings window for Surface, locate the Expression section.
3
In the Expression text field, type comp1.J.
4
Locate the Data section. From the Dataset list, choose Inhomogeneous Swelling/Parametric Solutions 1 (sol2).
5
From the Parameter value (ar) list, choose 0.6.
6
Locate the Coloring and Style section. From the Color table list, choose Rainbow.
7
From the Color table type list, choose Discrete.
Mesh 1
1
In the Model Builder window, right-click Swelling Ratio and choose Mesh.
2
In the Settings window for Mesh, locate the Data section.
3
From the Dataset list, choose Inhomogeneous Swelling/Parametric Solutions 1 (sol2).
4
From the Parameter value (ar) list, choose 0.6.
5
Locate the Coloring and Style section. From the Element color list, choose None.
6
Click to expand the Inherit Style section. From the Plot list, choose Surface 1.
Deformation 1
Right-click Mesh 1 and choose Deformation.
Mesh 1, Surface 1
Right-click and choose Duplicate.
Surface 2
1
In the Settings window for Surface, locate the Data section.
2
From the Parameter value (ar) list, choose 2.
3
Click to expand the Inherit Style section. From the Plot list, choose Surface 1.
Transformation 1
1
Right-click Surface 2 and choose Transformation.
2
In the Settings window for Transformation, locate the Transformation section.
3
In the R text field, type 1.25.
Mesh 2
1
In the Model Builder window, under Results > Swelling Ratio click Mesh 2.
2
In the Settings window for Mesh, locate the Data section.
3
From the Parameter value (ar) list, choose 2.
Transformation 1
1
Right-click Mesh 2 and choose Transformation.
2
In the Settings window for Transformation, locate the Transformation section.
3
In the R text field, type 1.25.
Mesh 2, Surface 2
1
In the Model Builder window, under Results > Swelling Ratio, Ctrl-click to select Surface 2 and Mesh 2.
2
Right-click and choose Duplicate.
Surface 3
1
In the Settings window for Surface, locate the Data section.
2
From the Parameter value (ar) list, choose 10.
Transformation 1
1
In the Model Builder window, expand the Surface 3 node, then click Transformation 1.
2
In the Settings window for Transformation, locate the Transformation section.
3
In the R text field, type 3.5.
Mesh 3
1
In the Model Builder window, under Results > Swelling Ratio click Mesh 3.
2
In the Settings window for Mesh, locate the Data section.
3
From the Parameter value (ar) list, choose 10.
Transformation 1
1
In the Model Builder window, expand the Mesh 3 node, then click Transformation 1.
2
In the Settings window for Transformation, locate the Transformation section.
3
In the R text field, type 3.5.
Swelling Ratio
In the Model Builder window, under Results click Swelling Ratio.
Table Annotation 1
1
In the Swelling Ratio toolbar, click  More Plots and choose Table Annotation.
2
In the Settings window for Table Annotation, locate the Data section.
3
From the Source list, choose Local table.
4
In the table, enter the following settings:
x-coordinate
y-coordinate
Annotation
0.2
-0.35
D/H = 0.6
2.0
-0.35
D/H = 2
6.0
-0.35
D/H = 10
5
Locate the Coloring and Style section. Clear the Show point checkbox.
6
From the Anchor point list, choose Center.
7
In the Swelling Ratio toolbar, click  Plot.
8
Click the  Zoom Extents button in the Graphics toolbar.
Gel Height
Create a plot of the normalized gel height, similar to Fig. 6 in Ref. 2.
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Gel Height in the Label text field.
3
Locate the Data section. From the Dataset list, choose Inhomogeneous Swelling/Parametric Solutions 1 (sol2).
4
From the Time selection list, choose Last.
5
Click to expand the Title section. From the Title type list, choose Manual.
6
In the Title text area, type Equilibrium gel height.
7
Locate the Plot Settings section.
8
Select the x-axis label checkbox. In the associated text field, type Gel aspect ratio (1).
9
Select the y-axis label checkbox. In the associated text field, type Normalized gel height h/H<sub>0</sub>.
10
Locate the Axis section. Select the Manual axis limits checkbox.
11
In the x minimum text field, type 0.1.
12
In the x maximum text field, type 100.
13
In the y minimum text field, type 3.2.
14
In the y maximum text field, type 5.
15
Select the x-axis log scale checkbox.
16
Locate the Grid section. Select the Manual spacing checkbox.
17
In the y spacing text field, type 0.2.
18
Locate the Legend section. Clear the Show legends checkbox.
Point Graph 1
1
Right-click Gel Height and choose Point Graph.
2
Select Point 7 only.
3
In the Settings window for Point Graph, locate the y-Axis Data section.
4
In the Expression text field, type (height+w)/(height/stretch0).
5
Locate the x-Axis Data section. From the Axis source data list, choose ar.
6
Click to expand the Coloring and Style section. From the Width list, choose 2.
7
Find the Line markers subsection. From the Marker list, choose Circle.
Add Physics
For comparison, compute the swelling stretches when the gel is allowed to swell homogeneously. The free and uniaxial swelling stretches can be computed by solving two nonlinear global DAEs.
1
In the Home toolbar, click  Add Physics to open the Add Physics window.
2
Go to the Add Physics window.
3
In the tree, select Mathematics > ODE and DAE Interfaces > Global ODEs and DAEs (ge).
4
Click the Add to Component 1 button in the window toolbar.
5
In the Home toolbar, click  Add Physics to close the Add Physics window.
Global ODEs and DAEs (ge)
Global Equations 1 (ODE1)
1
In the Settings window for Global Equations, locate the Global Equations section.
2
In the table, enter the following settings:
Name
f(u,ut,utt,t) (1)
Initial value (u_0) (1)
Initial value (ut_0) (1/s)
Description
lam0_free
G/(lam0_free^3)*(lam0_free^2-1)+Kmix*(log(1-lam0_free^(-3))+lam0_free^(-3)+chi*lam0_free^(-6))-mu_sw
2
0
Free swelling ratio
lam0_ua
G/(stretch0^2*lam0_ua)*(lam0_ua^2-1)+Kmix*(log(1-1/(stretch0^2*lam0_ua))+1/(stretch0^2*lam0_ua)+chi/(stretch0^4*lam0_ua^2))-mu_sw
2
0
Uniaxial swelling ratio
3
Locate the Units section. Click  Define Source Term Unit.
4
In the Source term quantity table, enter the following settings:
Source term quantity
Unit
Custom unit
Pa
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies > Stationary.
4
Click the Add Study button in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Homogeneous Swelling
1
In the Settings window for Study, type Homogeneous Swelling in the Label text field.
2
Locate the Study Settings section. Clear the Generate default plots checkbox.
Step 1: Stationary
Disable all physics interfaces except for the Global ODEs and DAEs interface. Also, since only global DOFs are solved for, no mesh is required.
1
In the Model Builder window, under Homogeneous Swelling click Step 1: Stationary.
2
In the Settings window for Stationary, locate the Physics and Variables Selection section.
3
Select the Modify model configuration for study step checkbox.
4
In the tree, select Component 1 (comp1) > Definitions > Contact Pair 1 (p1).
5
Click  Disable.
6
In the tree, select Component 1 (comp1) > Solid Mechanics (solid), Controls spatial frame and Component 1 (comp1) > Darcy’s Law (dl).
7
Click  Disable in Model.
8
In the tree, select Component 1 (comp1) > Multiphysics > Poroelasticity 1 (poro1).
9
Click  Disable in Model.
10
Click to expand the Mesh Selection section. In the table, enter the following settings:
Component
Mesh
Component 1
No mesh
11
In the Study toolbar, click  Compute.
Results
Create an Evaluation Group to read out the homogeneous swelling stretches for use in the Gel Height plot.
Homogeneous Swelling Stretches
1
In the Results toolbar, click  Evaluation Group.
2
In the Settings window for Evaluation Group, type Homogeneous Swelling Stretches in the Label text field.
3
Locate the Data section. From the Dataset list, choose Homogeneous Swelling/Solution 11 (sol11).
Global Evaluation 1
1
Right-click Homogeneous Swelling Stretches and choose Global Evaluation.
2
In the Homogeneous Swelling Stretches toolbar, click  Evaluate.
Global 1
1
In the Model Builder window, right-click Gel Height and choose Global.
2
In the Settings window for Global, locate the y-Axis Data section.
3
In the table, enter the following settings:
Expression
Unit
Description
3.3900
1
Free homogeneous swelling
4.7815
1
Uniaxial homogeneous swelling
4
Locate the x-Axis Data section. From the Axis source data list, choose ar.
5
Click to expand the Coloring and Style section. From the Width list, choose 2.
Gel Height
In the Model Builder window, click Gel Height.
Table Annotation 1
1
In the Gel Height toolbar, click  More Plots and choose Table Annotation.
2
In the Settings window for Table Annotation, locate the Data section.
3
From the Source list, choose Local table.
4
In the table, enter the following settings:
x-coordinate
y-coordinate
Annotation
3.2
3.45
Free swelling
3.2
4.85
Uniaxial swelling
5
Locate the Coloring and Style section. Clear the Show point checkbox.
6
From the Anchor point list, choose Center.
7
In the Gel Height toolbar, click  Plot.
Global Definitions
Now, set up the uniaxial consolidation load case. Start by adding parameter cases where some of the material parameters will be modified, see Table 1.
Material Parameters
1
In the Home toolbar, click  Parameter Case.
2
In the Settings window for Case, type Uniaxial Consolidation Case in the Label text field.
3
Locate the Parameters section. In the table, enter the following settings:
Name
Expression
Description
chi
0.2
Enthalpy of mixing
D
8e-10[m^2/s]
Solvent diffusivity
Model Parameters
1
In the Home toolbar, click  Parameter Case.
2
In the Settings window for Case, type Uniaxial Consolidation Case in the Label text field.
3
Locate the Parameters section. In the table, enter the following settings:
Name
Expression
Description
stretch0
3.215
Free swelling stretch
mu_ref
0[Pa]
Reference chemical potential
stress
0.05*Kmix/(stretch0^2)
Applied stress
Step 1 (step1)
Add a smooth Step function, which will be used to apply the creep load.
1
In the Home toolbar, click  Functions and choose Global > Step.
2
In the Settings window for Step, locate the Parameters section.
3
In the Location text field, type 0[s].
4
Click to expand the Smoothing section. From the Location definition list, choose Beginning of step.
Solid Mechanics (solid)
Boundary Load 1
1
In the Physics toolbar, click  Boundaries and choose Boundary Load.
2
Select Boundary 8 only.
3
In the Settings window for Boundary Load, locate the Force section.
4
Specify the fA vector as
-stress*step1(t)
z
Roller 1
1
In the Physics toolbar, click  Boundaries and choose Roller.
2
Select Boundary 9 only.
Darcy’s Law (dl)
Free Flow (Consolidation)
1
In the Physics toolbar, click  Boundaries and choose Pressure.
2
In the Settings window for Pressure, type Free Flow (Consolidation) in the Label text field.
3
Select Boundary 8 only.
4
Locate the Pressure section. In the p0 text field, type mu_ref.
Mesh 2
In the Mesh toolbar, click Add Mesh and choose Add Mesh.
Mapped 1
1
In the Mesh toolbar, click  Mapped.
2
In the Settings window for Mapped, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domain 2 only.
Distribution 1
1
Right-click Mapped 1 and choose Distribution.
2
Select Boundaries 7 and 8 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
In the Number of elements text field, type 1.
Distribution 2
1
In the Model Builder window, right-click Mapped 1 and choose Distribution.
2
Select Boundaries 6 and 9 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
In the Number of elements text field, type 100.
Mapped 2
1
In the Mesh toolbar, click  Mapped.
2
In the Settings window for Mapped, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domain 1 only.
Distribution 1
1
Right-click Mapped 2 and choose Distribution.
2
Select Boundaries 1 and 5 only.
3
In the Settings window for Distribution, locate the Distribution section.
4
In the Number of elements text field, type 1.
5
Click  Build All.
Add Study
Add a new time-dependent study to compute the consolidation load case. Modify the model configuration for the study to exclude the boundary conditions specific to the swelling load case.
1
In the Study toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies > Time Dependent.
4
Click the Add Study button in the window toolbar.
5
In the Study toolbar, click  Add Study to close the Add Study window.
Uniaxial Consolidation
1
In the Settings window for Study, locate the Study Settings section.
2
Clear the Generate default plots checkbox.
3
In the Label text field, type Uniaxial Consolidation.
Parametric Sweep
Use a parameter switch to select the correct set of parameters for the uniaxial consolidation case.
1
In the Study toolbar, click  Parametric Sweep.
2
In the Settings window for Parametric Sweep, locate the Study Settings section.
3
From the Sweep type list, choose Parameter switch.
4
Click  Add twice.
5
In the table, enter the following settings:
Switch
Cases
Case numbers
Material Parameters
User defined
1
Model Parameters
User defined
1
Step 1: Time Dependent
1
In the Model Builder window, click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
In the Output times text field, type {0 1 2 5 10 20 50 100 200 1000}*(height/stretch0)^2/D.
4
Locate the Physics and Variables Selection section. Select the Modify model configuration for study step checkbox.
5
In the tree, select Component 1 (comp1) > Definitions > Contact Pair 1 (p1).
6
Click  Disable.
7
In the tree, select Component 1 (comp1) > Darcy’s Law (dl) > Contact (Swelling) and Component 1 (comp1) > Darcy’s Law (dl) > Free Flow (Swelling).
8
Click  Disable.
9
In the tree, select Component 1 (comp1) > Global ODEs and DAEs (ge).
10
Click  Disable in Model.
Solution 12 (sol12)
1
In the Model Builder window, right-click Solver Configurations and choose Show Default Solver.
2
Expand the Solution 12 (sol12) node.
3
In the Model Builder window, expand the Uniaxial Consolidation > Solver Configurations > Solution 12 (sol12) > Dependent Variables 1 node, then click Pressure (comp1.mu).
4
In the Settings window for Field, locate the Scaling section.
5
From the Method list, choose Manual.
6
In the Scale text field, type 1.0e5.
7
In the Model Builder window, click Displacement Field (comp1.u).
8
In the Settings window for Field, locate the Scaling section.
9
In the Scale text field, type height.
10
In the Model Builder window, click Time-Dependent Solver 1.
11
In the Settings window for Time-Dependent Solver, locate the Time Stepping section.
12
From the Steps taken by solver list, choose Strict.
13
Select the Initial step checkbox.
14
In the Model Builder window, expand the Uniaxial Consolidation > Solver Configurations > Solution 12 (sol12) > Time-Dependent Solver 1 node, then click Direct.
15
In the Settings window for Direct, locate the General section.
16
From the Solver list, choose PARDISO.
17
In the Model Builder window, click Fully Coupled 1.
18
In the Settings window for Fully Coupled, locate the Method and Termination section.
19
From the Jacobian update list, choose Once per time step.
20
In the Maximum number of iterations text field, type 10.
21
In the Study toolbar, click  Compute.
Results
Add a Selection to the Dataset for the Uniaxial Consolidation study to exclude the rigid substrate when creating surface plots.
Uniaxial Consolidation/Solution 12 (sol12)
In the Model Builder window, expand the Results > Datasets node, then click Uniaxial Consolidation/Solution 12 (sol12).
Selection
1
In the Results toolbar, click  Attributes and choose Selection.
2
In the Settings window for Selection, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domain 2 only.
Result Templates
1
In the Results toolbar, click  Result Templates to open the Result Templates window.
2
Go to the Result Templates window.
3
In the tree, select Uniaxial Consolidation/Solution 12 (sol12) > Solid Mechanics > Stress (solid).
4
Click the Add Result Template button in the window toolbar.
5
In the Results toolbar, click  Result Templates to close the Result Templates window.
Results
Stress and Chemical Potential
1
In the Model Builder window, expand the Results > Stress (solid) node, then click Stress (solid).
2
In the Settings window for 2D Plot Group, type Stress and Chemical Potential in the Label text field.
3
Locate the Data section. From the Time (s) list, choose 24187.
4
Locate the Plot Settings section. From the View list, choose New view.
5
Click to expand the Plot Array section. Select the Enable checkbox.
6
From the Displacement list, choose Absolute.
7
In the Cell displacement text field, type 0.4.
Mixture Stress
1
In the Model Builder window, under Results > Stress and Chemical Potential click Surface 1.
2
In the Settings window for Surface, type Mixture Stress in the Label text field.
3
Locate the Expression section. In the Expression text field, type -poro1.sGpzz.
4
From the Unit list, choose kPa.
5
Select the Description checkbox. In the associated text field, type Normal stress.
Chemical Potential
1
Right-click Mixture Stress and choose Duplicate.
2
In the Settings window for Surface, type Chemical Potential in the Label text field.
3
Locate the Expression section. In the Expression text field, type mu.
4
In the Description text field, type Chemical potential.
5
Locate the Inherit Style section. From the Plot list, choose Mixture Stress.
Arrow Surface 1
1
In the Model Builder window, right-click Stress and Chemical Potential and choose Arrow Surface.
2
In the Settings window for Arrow Surface, locate the Expression section.
3
In the R-component text field, type dl.u.
4
In the Z-component text field, type dl.w.
5
Locate the Arrow Positioning section. Find the R grid points subsection. In the Points text field, type 4.
6
Find the Z grid points subsection. In the Points text field, type 20.
7
Click to expand the Inherit Style section. From the Plot list, choose Mixture Stress.
8
Clear the Arrow scale factor checkbox.
9
Clear the Color checkbox.
10
Clear the Color and data range checkbox.
11
Click to expand the Plot Array section. Select the Manual indexing checkbox.
12
In the Index text field, type 1.
Deformation 1
Right-click Arrow Surface 1 and choose Deformation.
Table Annotation 1
1
In the Stress and Chemical Potential toolbar, click  More Plots and choose Table Annotation.
2
In the Settings window for Table Annotation, locate the Data section.
3
From the Source list, choose Local table.
4
In the table, enter the following settings:
x-coordinate
y-coordinate
Annotation
0.05
-0.015
Normal stress $-\sigma_{zz}$
0.45
-0.015
Chemical potential $\mu$
5
Select the LaTeX markup checkbox.
6
Locate the Coloring and Style section. Clear the Show point checkbox.
7
From the Anchor point list, choose Center.
Stress and Chemical Potential
1
In the Model Builder window, click Stress and Chemical Potential.
2
In the Settings window for 2D Plot Group, click to expand the Title section.
3
From the Title type list, choose Manual.
4
In the Title text area, type
Surface: Normal stress (kPa) Surface: Chemical potential (kPa)
Arrow Surface: Total Darcy velocity field (spatial frame)
on separate lines.
5
In the Stress and Chemical Potential toolbar, click  Plot.
6
Click the  Zoom Extents button in the Graphics toolbar.
Thickness Stretch
Follow the instructions below to reproduce Fig. 4 in Ref. 1.
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Thickness Stretch in the Label text field.
3
Locate the Data section. From the Dataset list, choose Uniaxial Consolidation/Solution 12 (sol12).
4
Locate the Plot Settings section.
5
Select the x-axis label checkbox. In the associated text field, type Z/H.
6
Select the y-axis label checkbox. In the associated text field, type \lambda<sub>3</sub>.
7
Locate the Axis section. Select the Manual axis limits checkbox.
8
In the x minimum text field, type 0.
9
In the y minimum text field, type 0.5.
10
In the y maximum text field, type 3.5.
11
Locate the Grid section. Select the Manual spacing checkbox.
12
In the x spacing text field, type 0.2.
13
In the y spacing text field, type 0.5.
14
Locate the Legend section. Clear the Show legends checkbox.
Line Graph 1
1
Right-click Thickness Stretch and choose Line Graph.
2
Select Boundary 6 only.
3
In the Settings window for Line Graph, locate the y-Axis Data section.
4
In the Expression text field, type solid.stchp3*stretch0.
5
Select the Description checkbox. In the associated text field, type Third principal stretch.
6
Locate the x-Axis Data section. From the Parameter list, choose Expression.
7
In the Expression text field, type 1-Z/height.
8
Click to expand the Coloring and Style section. From the Color cycle list, choose Long.
Thickness Stretch
In the Model Builder window, click Thickness Stretch.
Table Annotation 1
1
In the Thickness Stretch toolbar, click  More Plots and choose Table Annotation.
2
In the Settings window for Table Annotation, locate the Data section.
3
From the Source list, choose Local table.
4
In the table, enter the following settings:
x-coordinate
y-coordinate
Annotation
0
3.2
$0^+$
0.06
2.9
1
0.09
2.8
2
0.14
2.7
5
0.2
2.6
10
0.28
2.5
20
0.43
2.2
50
0.59
2.0
100
0.76
1.5
$Dt/H_0^2 = 200$
0.8
0.8
1000
5
Select the LaTeX markup checkbox.
6
Locate the Coloring and Style section. Clear the Show point checkbox.
7
From the Anchor point list, choose Lower left.
8
In the Thickness Stretch toolbar, click  Plot.
Reference Data, Hong and others (2008), t=2
Compare the results with data from Ref. 1 for a few time points.
1
In the Results toolbar, click  Table.
2
In the Settings window for Table, type Reference Data, Hong and others (2008), t=2 in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file hydrogel_swelling_uniaxial_t2.txt.
Reference Data, Hong and others (2008), t=20
1
In the Results toolbar, click  Table.
2
In the Settings window for Table, type Reference Data, Hong and others (2008), t=20 in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file hydrogel_swelling_uniaxial_t20.txt.
Reference Data, Hong and others (2008), t=100
1
In the Results toolbar, click  Table.
2
In the Settings window for Table, type Reference Data, Hong and others (2008), t=100 in the Label text field.
3
Locate the Data section. Click  Import.
4
Browse to the model’s Application Libraries folder and double-click the file hydrogel_swelling_uniaxial_t100.txt.
Reference Data, Hong and others (2008), t=2
1
Right-click Thickness Stretch and choose Table Graph.
2
In the Settings window for Table Graph, locate the Coloring and Style section.
3
Find the Line style subsection. From the Line list, choose None.
4
Find the Line markers subsection. From the Marker list, choose Circle.
5
From the Color list, choose From theme.
6
In the Label text field, type Reference Data, Hong and others (2008), t=2.
Reference Data, Hong and others (2008), t=20
1
Right-click Reference Data, Hong and others (2008), t=2 and choose Duplicate.
2
In the Settings window for Table Graph, locate the Data section.
3
From the Table list, choose Reference Data, Hong and others (2008), t=20.
4
In the Label text field, type Reference Data, Hong and others (2008), t=20.
Reference Data, Hong and others (2008), t=100
1
Right-click Reference Data, Hong and others (2008), t=20 and choose Duplicate.
2
In the Settings window for Table Graph, type Reference Data, Hong and others (2008), t=100 in the Label text field.
3
Locate the Data section. From the Table list, choose Reference Data, Hong and others (2008), t=100.
4
In the Thickness Stretch toolbar, click  Plot.
Inhomogeneous Swelling
Finally, disable the boundary conditions added for the consolidation study in the settings for Step 1: Time Dependent under the Inhomogeneous Swelling study node if you want to recompute the swelling load case.
Step 1: Time Dependent
1
In the Model Builder window, under Inhomogeneous Swelling click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Physics and Variables Selection section.
3
Select the Modify model configuration for study step checkbox.
4
In the tree, select Component 1 (comp1) > Global ODEs and DAEs (ge).
5
Click  Disable in Model.
6
In the tree, select Component 1 (comp1) > Solid Mechanics (solid), Controls spatial frame > Boundary Load 1 and Component 1 (comp1) > Solid Mechanics (solid), Controls spatial frame > Roller 1.
7
Click  Disable.
8
In the tree, select Component 1 (comp1) > Darcy’s Law (dl) > Free Flow (Consolidation).
9
Click  Disable.