This is example how to run a simple problem for gel model with visocus-relaxation. Is is dedicated for beginners starting work with MoFEM.
In this model mechanical and solvent concentration field is modeled. In addition viscous relaxation is taken into account by considering model with two springs and dashpot, with additional internal state variable. Detail about implementation and equations can be found in structure GelModule::Gel
Analysis of simple gel brick
Input files
It is assumed that user parent working directory is users_modules/gel in build directory.
Journal file
Mesh file with boundary conditions is created in Cubit journal file with prismatic bar (brick_of_gel.jou),
reset
set duplicate block elements on
brick x 10 y 1 z 1
webcut volume 1 with plane yplane imprint merge
# Sete block with volume with gel material
block 1 volume 1
block 1 name 'GEL_1'
block 2 volume 2
block 2 name 'GEL_2'
# Set Dirichelt boundary condition for solvent concentration
block 6 surface 11 13
block 6 name 'CHEMICAL_LOAD_1'
block 6 attribute count 1
block 6 attribute index 1 1
# Set solvent flux, i.e. Neumann boundary condition
block 5 surface 8 15 5 14 10 3
block 5 name 'FLUX_CHEMICAL_LOAD'
block 5 attribute count 1
block 5 attribute index 1 0
# Kinematic boundary condition for mechanical field
create displacement on surface 9 16 dof 1 fix 0
create displacement on vertex 7 dof 1 dof 2 dof 3 fix 0
create displacement on vertex 4 dof 2 fix 0
create displacement on vertex 8 dof 3 fix 0
# Static boundary condition, i.e. pressure
create pressure on surface 6 magnitude 1
# Make a mesh
volume all scheme Tetmesh
volume all size auto factor 8
mesh volume all
# Set block 3 and set 10 node tetrahedrons in that block
block 3 tet all
block 3 element type TETRA10
Defining material properties
Material properties file (gel_config.in).
# Material parameters for Gel material model
[block_1]
# Spring Alpha
gAlpha = 1
vAlpha = 0
# Spring Beta
gBeta = 1
vBeta = 0
# Daspot
gBetaHat = 1
vBetaHat = 0
# Volume per solvent molecule
oMega = 10
# Initial chemical potential
mU0 = 0.5
# Solvent transport material parameters
vIscosity = 1
pErmeability = 1
# User can set approximation order for each blocks independently
# Set approximation order of block 1 as quadratic function.
oRder = 2
[block_2]
# Spring Alpha
gAlpha = 1
vAlpha = 0
# Spring Beta
gBeta = 1
vBeta = 0
# Daspot
gBetaHat = 1
vBetaHat = 0
# Volume per solvent molecule
oMega = 5
# Initial chemical potential
mU0 = 0.5
# Solvent transport material parameters
vIscosity = 1
pErmeability = 1
# User can set approximation order for each blocks independently
# Set approximation order of block 1 as quadratic function.
oRder = 2
Native constitutive model is implemented after Yuhang Hu and Zhigang Suo [32].
Gel model
In basic implementations all springs and dashpot are implemented using Hook'e equation. However generalization to large deformations and displacements is straightforward, code uses automatic differentiation (ADOL-C), and user need only implement stress-strain physical equation and all tangent matrices are calculated automatically.
For more details about model implementation, equations and physical constants look here GelModule::Gel::ConstitutiveEquation
Executing code
Example how to execute example problem:
mpirun -np 4 ./gel_analysis \
-my_file brick_of_gel.cub \
-my_gel_config gel_config.in \
-ksp_type fgmres -ksp_final_residual \
-ksp_atol 1e-10 -ksp_rtol 1e-10 \
-ksp_monitor -ksp_converged_reason \
-pc_type lu -pc_factor_mat_solver_package mumps \
-snes_atol 1e-8 -snes_rtol 1e-8 \
-snes_type newtonls -snes_linesearch_type basic \
-snes_linesearch_monitor -snes_monitor \
-ts_type beuler -ts_dt 0.1 -ts_final_time 10 \
-ts_monitor
Notes:
- Option -my_file gives name of mesh file
- Option -my_gel_config give name of material config file
- Time integration scheme is controlled by -ts_type, look to PETSc manual for more details.
- Resolution of post-processing mesh is set by -my_max_post_proc_ref_level 0. If more than 0, i.e. 1,2,... a denser post-processing mesh is generated. This results in biger size of post processing files and compromise overall efficiency.
- Option -snes_atol control absolute tolerance of nonlinear solver, should be bigger than ksp_atol 1e-10
- Option -snes_rtol control relative tolerance of nonlinear solver
- Using -my_output_prt 1, each time step is post-processed.
Result:
Example of animation obtained from this example for approximation order 4
Gel brick
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
User (hackable) constitutive equation
In file UserGelModel.hpp user can implement own constitutive equation, for details see UserGelConstitutiveEquation. Any change in UserGelModel.hpp have to follow rebuilding the code
Similiar (not exact) method is used to implement general nonlinear elastic materials, see for example NeoHookean.hpp or Hooke.hpp
