v0.14.0
Wood Cell (Response to an Environmental Change)

This is an example of how to run a simultation on a cell structure, similar to what is found within wood.

Wood is a naturally occurring material, which interacts strongly with water from the air in the surrounding environment. These changes in moisture content due to the external relative humidity can cause distortions that are observable on the macro-scale. However the interactions with the water happens within the cell wall polymers, where the incoming water (solvent) is chemically bound to the polymers within through hydrogen bonding.

The transportation of this bound-water through the wood cell is a very complex process. It is not only determined by simple Fickian diffusion due to the flux of the chemical potential of the air/water in the environment, but the rate can also be controlled by the rate at which the polymer substrate deforms viscoelastically. The question of whether or not the moisture transport through these polymers is viscoelastically limited has been widely debated in the literature [74]. Using the polymeric gelmodel, in the form of a wood cell, we can determine numerically whether or not viscoelastically limited solvent uptake occurs within the amorphous polymers within wood.

The following graph demonstrates the principles of viscoelastically limited solvent migration using the principle of the Deborah number. When the Deborah number (DEB_D) is close to 1 the diffusion and relaxation times are almost equal and it is around this region that the viscoelastic process controls the diffusion of bound-water through the sample. In this case the material temperatures are thought of as concentration dependent as we are dealing with isothermal conditions. TE is the temperature, below which the polymer acts as an elastic solid, and is characterised by a very stiff viscous damper. Tv is temperature above which the polymer acts as a viscous liquid and Tg is the glass transition temperature where the material behaviour switches from glassy to rubbery. The Deborah number is given as follows:

\[ DEB_D = \frac{T}{t} \]

where t is the time for diffusion to occur and T is the viscoelastic relaxation time, defined as follows:

\[ T = \frac{\hat{G}^\beta}{G^\beta} \]

Viscoelastically Limited Diffusion as Demonstrated by Vrentas et al 1975

The model used is based upon the theory of concurrent poro-elastic and viscoelastic swelling within polymeric gels by Yuhang Hu and Zhigang Suo. Viscoelasticity and poroelasticity in elastomeric gels. Acta Mechanica Solida Sinica, Vol. 25, No. 5 [32]. Further details on the mechanical model can be found in Mechanical model.

Properties

The wood cell is composed of a mixture of hemicelluloses and lignin, with the properties obtained from a paper which determines the viscoelastic properties of wood polymers using multi-scale modelling [25]. There properties are based upon a moisture content of around 10% within the wood, and as such we assign the environment to around 40% relative humidity. The shear modulus of the hemicellulose and lignin mixture is taken as 3.7 GPa, which is split evenly across the two springs, and the dashpot shear modulus is taken as 96 GPa.s. The diffusivity of bound water is taken as 1 e-05 mm2/s and the size of 1 mole of water is taken 18000 mm3. The permeability is set using the equation for diffusivity.

# Material parameters for Gel material model
[block_1]
# Spring Alpha
gAlpha = 1.85e03
vAlpha = 0.23
# Spring Beta
gBeta = 1.85e03
vBeta = 0.23
# Daspot
gBetaHat = 96e03
vBetaHat = 0.23
# Volume per solvent molecule
oMega = 18000
# Initial chemical potential
mU0 = 9143945
# Solvent transport material parameters
vIscosity = 1e-06
pErmeability = 8.77e-17
# User can set approximation order for each blocks independently
# Set approximation order of block 1 as quadratic function.
oRder = 3

Above file can be found $GEL_MODULE_DIRECTORT/gel_config_cell.in

Boundary Conditions

The wood cell is exposed to an environment of increased relative humidity, where the incoming water vapour is allowed to diffusive into the wood cell in the form of bound-water. The influx of bound-water comes through the lumens within the cell.

Wood Cell Problem
Wood Cell BC
Wood Cell Dimensions

The cell is constrained in the x and y directions, with the primary mode for adsorption and desorption being expansion into, or expansion out from, the lumen within the cell. This replicated the cell being within a large network of cells tightly packed together.

The change in relative humidity is applied slowly, with a simple 1D bar problem being solved to determined the evolution of the chemical potential through a sample observable on the macro-scale, with the diffusion coefficient being set to that of water vapour through wood, 2 e-03 mm2/s.

1D Diffusion Problem
1D Diffusion Test
1D Diffusion Test Results

The test is carried out for several lengths of diffusion, however within this example we focus on a length of 100 micrometers. The time series for chemical potential can then be applied as a dirichlet BC within the lumens using a history file. This file can be found in $GEL_MODULE_DIRECTORT/my_chemical_history_cell.in. The history file is very long so it will not appear in this documentation, however examples an example on history files can be found in Boundary Conditions.

The chemical load history contains two columns, the left hand column is the time and the right hand column is the factor which the boundary chemical potential is multiplied by. In this case the reference boundary conditions is \mu = 3989505 J mol-1.

The boundary conditions can be calculated as shown in Dirichelt boundary condition for chemical potential

Cell Mesh

The mesh is genrated using the following script:

reset
set duplicate block elements on
brick x 0.0008 y 0.0016 z 0.00015
brick x 0.0008 y 0.0016 z 0.00015
rotate volume 2 angle 15 about z include_merged
move Volume 2 y 0.0015 include_merged
move Volume 2 x -0.0004 include_merged
brick x 0.0008 y 0.0016 z 0.00015
rotate volume 3 angle 345 about z include_merged
move Volume 3 y 0.0015 include_merged
move Volume 3 x 0.0004 include_merged
subtract volume 2 from volume 1
subtract volume 3 from volume 1
brick x 0.0008 y 0.0016 z 0.00015
rotate volume 4 angle 345 about z include_merged
move Volume 4 y -0.0015 include_merged
move Volume 4 x -0.0004 include_merged
brick x 0.0008 y 0.0016 z 0.00015
rotate volume 5 angle 15 about z include_merged
move Volume 5 y -0.0015 include_merged
move Volume 5 x 0.0004 include_merged
subtract volume 4 from volume 1
subtract volume 5 from volume 1
brick x 0.001 y 0.0016 z 0.00015
move Volume 5 y -0.0015 include_merged
volume 1 copy
volume 1 copy
move Volume 1 y 0.0011342 include_merged
move Volume 7 y -0.0011342 include_merged
subtract volume 1 from volume 6
subtract volume 7 from volume 6
brick x 0.001 y 0.0016 z 0.00015
rotate volume 9 angle 15 about z include_merged
move Volume 9 y 0.00112 include_merged
brick x 0.001 y 0.0016 z 0.00015
rotate volume 10 angle 345 about z include_merged
move Volume 10 y -0.00112 include_merged
subtract volume 9 from volume 8
subtract volume 10 from volume 8
brick x 0.001 y 0.0016 z 0.00015
move Volume 11 x -0.0005 include_merged
volume 8 copy
rotate volume 12 angle 180 about z include_merged
move Volume 8 x 0.0001 include_merged
move Volume 12 x -0.0001 include_merged
volume 11 copy
move Volume 13 x 0.001 include_merged
subtract volume 11 from volume 8
subtract volume 13 from volume 12
move Volume 12 x -0.00005 include_merged
move Volume 8 x 0.00005 include_merged
subtract volume 12 from volume 6
subtract volume 8 from volume 6
brick x 0.001 y 0.0016 z 0.00015
volume 14 copy
move Volume 14 x -0.00095 include_merged
move Volume 15 x 0.00095 include_merged
subtract volume 14 from volume 6
subtract volume 15 from volume 6
brick x 0.001 y 0.0016 z 0.00015
volume 16 copy
move Volume 16 y 0.001527405 include_merged
move Volume 17 y -0.001527405 include_merged
subtract volume 16 from volume 6
subtract volume 17 from volume 6
volume all scale 30
# Set block with volume with gel material
block 1 volume all
block 1 name 'GEL_1'
# Make a mesh
volume all scheme Tetmesh
volume all size auto factor 4
mesh volume all
# Set Dirichelt boundary condition for solvent concentration
block 3 surface 81 84 83 82 138 137 139 75 74 76 73 131 130 132 182 153 184 174 146 177 209 206 201 214
block 3 name 'CHEMICAL_LOAD_1'
block 3 attribute count 1
block 3 attribute index 1 3.98950543760577e06
# Set solvent flux, i.e.
block 5 surface 179 180 172 171 173 186 178 181 185 183 210 204 202 212 208 207 199 200 213 211
block 5 name 'FLUX_CHEMICAL_LOAD_2'
block 5 attribute count 1
block 5 attribute index 1 0
# Kinematic boundary condition for mechanical field
create displacement on surface 208 207 199 200 dof 2 fix 0
# Kinematic boundary condition for mechanical field
create displacement on surface 210 204 202 212 dof 1 fix 0
# Kinematic boundary condition for mechanical field
create displacement on surface 211 dof 3 fix 0
# Set block 4 and set 10 node tetrahedrons in that block
block 4 tet all
block 4 element type TETRA10

The mesh file used in this example is located in $MODULE_DIRECTORY/meshes/wood_cell_example.cub

The journal file used in this example is located in $MODULE_DIRECTORY/meshes/wood_cell_example.jou

Executing Problem

mpirun -np 4 ./gel_analysis -my_file wood_cell_example.cub\
-my_gel_config gel_config_cell.in -my_chemical_load_history my_chemical_history_cell.in \
-ksp_type fgmres -ksp_final_residual -ksp_monitor -ksp_converged_reason \
-pc_type lu -pc_factor_mat_solver_package mumps -snes_monitor -ts_type beuler \
-ts_dt 10 -ts_final_time 1200 -ts_max_steps 10000000 -my_output_prt 1 \
-snes_linesearch_type basic -ts_monitor -snes_atol 1e-14 -snes_rtol 1e-14

Notes:

  • Option -my_file gives name of the mesh file
  • Option -my_gel_config gives the name of the material configuration file
  • Option -my_chemical_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
  • Option -snes_atol gives the absolute convergence tolerance
  • Option -snes_rtol gives the relative convergence tolerance
  • Option -snes_linesearch_type gives the linesearch type (currently can use basic or l2)

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

Results

The results of the test confirm the theory of viscoelastically limited solvent migration within wood polymers. Several different values of dashpot shear modulus are used to compare the time for equilibrium. In the case where the dashpot modulus is less than the value assigned in the config file, the uptake of moisture happens much fast, behaving in a Fickian manner. The rest of the results all display different degrees of viscoelastic solvent migration. This is evidence of the models ability to capture complex phenomena within natural cell wall polymers in wood and other cellulosic materials.

Wood Cell Animation
Evidence of Viscoelastically Limited Solvent Migration
Comparison of Dashpot Shear Moduli

Future Work

It is hoped that this gel model can be applied to more complex multi-scale problems within wood and other natural ligno-cellulosic materials. It can also be applied to sorption theory, including the analysis of sorption hysteresis etc. Future documentation will be provided for this.

To download the gel model visit: https://bitbucket.org/likask/mofem_um_gels

For reference on how to apply to multi-layered materials, as is found within wood cells, see Layered Gel Problem (Response to an Environmental Change).

surface
Definition: surface.py:1
scale
double scale
Definition: plastic.cpp:119
convert.type
type
Definition: convert.py:64
Beta
constexpr double Beta
Definition: radiation.cpp:40