v0.14.0
Gel Model Creep (response due to external load)

This is an example on how to run a simulation of a gel undergoing creep and sorption due to an external load.

To examine the creep response of a gel sample, a bar section (1/8th bar) is modelled in an environment of constant relative humidity (RH) on the boundary, and with a tensile load applied in the longitudinal direction of the bar. A dirichlet BC ( \(\mu\)) is imposed on surface 4 and a Neumann symmetry condition is imposed on surface 6.

Mechanical model

Native constitutive model is implemented after Yuhang Hu and Zhigang Suo. Viscoelasticity and poroelasticity in elastomeric gels. Acta Mechanica Solida Sinica, Vol. 25, No. 5 [32].

Gel model

The model represents a spring, \(\alpha\), in parallel with a spring/dashpot, \(\beta\)/ \(\hat{\beta}\), in series. Initially the load is taken up by both springs, and then over time the load in spring \beta is relaxed as the dashpot, \(\hat{\beta}\), is compressed. This happens until the spring \beta is completely relaxed and the load is taken up entirely by spring \alpha. At this point the system is in equilibrium with the environment.

Zener Evolution

For more details about model implementation, equations and physical constants look here GelModule::Gel::ConstitutiveEquation

Dirichelt boundary condition for chemical potential

For given relative humidity chemical load is calculated form following equation,

\[ \mu = {\mu_0 + RT ln(\text{RH})) \over \text{L}} \]

where \(\mu_0\) is initial chemical potential, \(R\) is a gas constant, \(T\) is temperature, \(RH\) is relative humidity of environment (surrounding air) and \(L\) is Avogadro consistent [32].

The water vapour pressure in described in terms of relative humidity, i.e. the ratio of partial water vapour pressure to saturation pressure of the air at a given temperature. Above equation estimates the Gibb's free energy of each solvent molecule (in this case water vapour) entering the system.

Mesh

Mesh is a 1/8th section of a rectangular cube, with 3 Dirichlet boundary conditions on the external faces and 3 Neumann conditions on the inner faces to impose symmetry. All dimensions are given in meters (m) and the chemical load is Gibb's free energy per solvent molecule entering the system with the units in Joules (J).

Creep Mesh

Journal File (creep_force.jou)

Mesh is generated in Cubit using following script,

reset
set duplicate block elements on
brick x 0.0025 y 0.005 z 0.0025
# Set block with volume with gel material
block 1 volume all
block 1 name 'GEL'
# Set Dirichelt boundary condition for solvent concentration
block 2 surface 2 4 5
block 2 name 'CHEMICAL_LOAD_1'
block 2 attribute count 1
block 2 attribute index 1 1.6102573617124e-17
# Set solvent flux, i.e.
block 4 surface 1 3 6
block 4 name 'FLUX_CHEMICAL_LOAD_2'
block 4 attribute count 1
block 4 attribute index 1 0
# Static boundary condition, i.e. pressure
create pressure on surface 5 magnitude -500
# Kinematic boundary condition for mechanical field
create displacement on surface 6 dof 1 fix 0
# Kinematic boundary condition for mechanical field
create displacement on surface 3 dof 2 fix 0
# Kinematic boundary condition for mechanical field
create displacement on surface 1 dof 3 fix 0
# Make a mesh
volume all scheme Tetmesh
volume all size auto factor 7
mesh volume all
# Set block 3 and set 10 node tetrahedrons in that block
block 3 tet all
block 3 element type TETRA10

The mesh file used in this example is located in $MODULE_DIRECTORY/meshes/creep_force.jou.

Material Properties

A simple section of amorphous cellulosic material is described. This material exhibits a viscoelastic response and is highly interactive with atmospheric water vapour content and makes a good case for analysis. The volume per solvent molecule (water) in it's liquid form, after chemically binding to the gel, is given by oMega in m^3. In this case the sample is initially in equilibium before the force is applied (i.e. \(\mu = \mu_0\)).

# Material parameters for Gel material model
# Spring Alpha
gAlpha = 6e3
vAlpha = 0.232
# Spring Beta
gBeta = 6e3
vBeta = 0.232
# Daspot
gBetaHat = 9e6
vBetaHat = 0.232
# Volume per solvent molecule
oMega = 2.988970058e-20
# Initial chemical potential
mU0 = 0
# Solvent transport material parameters
vIscosity = 9e6
pErmeability = 3.4e-8
# User can set approximation order for each blocks independently
# Set approximation order of block 1 as quadratic function.
[block_1]
oRder = 2

Above file can be found in $MODULE_DIRECTORY/gel_config_creep_sorption_example.in

Load History

To model the creep behaviour, a tensile force is applied in the y direction to a sample in equilibrium with it's surroundings. The application of a tensile force will encourage solvent uptake within the sample until the sample reaches equilibrium with it's environment again (i.e. increased solvent content). The load is then released, encouraging ejection of solvent from the sample until the sample once again reaches equilibrium with it's external environment.

To implement the tensile force in this manner, a load history file is used. The file includes a left hand column (time in seconds (s)) and a right hand column of load factor. In this example, a load factor of 1.0 is applied for 50000s until it is released (load factor zero) and another 50000s is allowed for the sample to recover to it's initial state.

0 1.00
50000 1.00
50001 0.00
100000 0.00

Executing Problem

mpirun -np 8 ./gel_analysis -my_file creep_force.cub -my_gel_config \
gel_config_creep_sorption_example.in -my_load_history my_load_history.in -ksp_type fgmres \
-ksp_final_residual -ksp_monitor -ksp_converged_reason -pc_type lu \
-pc_factor_mat_solver_package mumps -my_order 3 -snes_monitor -ts_type beuler \
-ts_dt 50 -ts_final_time 100000 -my_output_prt 10 -snes_linesearch_type basic \
-ts_monitor

Notes:

  • Option -my_file gives name of the mesh file
  • Option -my_gel_config gives the name of the material configuration file
  • Option -my_load_history gives the name of the load history file
  • Option -ts_dt gives the size of each time step in seconds (s)
  • Option -ts_final_time gives the final time in seconds (s)
  • Option -my_output_prt gives the number of vtk files output (i.e. -my_output_prt 1, means each time step is processed)
  • Option -np gives the number of processors to run on

Results

VTK files for ParaView

Running dynamic analysis out_values_1.h5m, out_values_2.h5m, ... are created for each time step. Post-processing mesh in output fails is stored in native MoAB data format using standard h5m.

VTK files can be generated using script located in users_modules/nonlinear_elasticity/do_vtk.sh. For example

../nonlinear_elasticity/do_vtk.sh out_*h5m

Paraview

Creep Test

To view the creep response of the material, first calculate the displacements (SPATIAL_POSISTION - coords). Plot the outside edge, furthest from the mid-point, against time. The solvent migration can also be plotted by selecting the corner that intersects all 3 inner faces (i.e. surface 1,3 and 6) and plotting against time. The results are plotted below.

Left: Chemical Load Evolution, Right: Time Series Data

The left hand graph shows the chemical load, \(\mu\), through the sample from the top to bottom at various times given in minutes in the legend. The graph on the right gives the displacement at the top corner of the sample (left axis) vs the time in minutes. It shows the reversible creep behaviour that has been induced by the loading conditions, where immediately after the load has been released (t=33min), the sample does not go back to it's reference configuration immediately, instead taking a further 35 minutes to do so.

The right hand axis shows the chemical load at the bottom corner of the sample. Initially during the tensile loading phase, the chemical load decreases at this area of the sample, as space has been created, thus reducing the concentration. The sample responds to this by taking in more solvent in order to reach equilibrium with the environment (i.e. \( \mu \) inside the sample is equal to \( \mu \) in the environment).

After the load has been released, the extra solvent that has migrated into the sample from the environment is compressed as the sample tries to regain it's reference state, forcing the extra solvent out of the sample until equilibrium with the outside environment has been reached again.