Brief description
Three-point-bending crack propagation test in the presence of contact conditions.
1. Prepare the <em>med</em> mesh file in <em>Salome</em>
NOTE: In order to use the matching meshes contact the solids need to touch each other along the contact surface in the input mesh. Presented below is a rather general approach which permits to mesh the contacting solids separately, refine these meshes around the contact surface, and, finally, merge the meshes together. You can find additional rules on definition of contact interfaces with matching meshes in the Readme file of the simple_contact directory (Basic User Module).
Geometry
- Create the contacting solids separately
- Use Intersection to find common geometrical entities
- Create the contact surface using Build -> Face from the intersection
- Partition solids with the intersection (one by one)
- Create initial crack face
- Create all necessary groups in the partitioned solids and on the crack face:
- Groups of volumes for solid blocks
- Groups of edges for fixed BCs
- Group of faces for load BCs
- Group of faces for springs
- Group of faces for the contact interface (same name per contact interface for both solids)
- Group of all physical vertices (same name for both solids)
- Group of all physical edges (same name for both solids)
- Group of faces for the initial crack surface
NOTE: Springs can be used to eliminate rigid body motions caused by the interplay between the contact conditions and the load (pressure) control. It is recommended to append the springs to the same surface where the pressure is applied.
Mesh
- Mesh each solid separately
- Create mesh groups from the geometry groups (for each solid separately)
- Mesh the crack surface
- Create group of faces for the crack surface (from the corresponding geometry group)
- Mesh the contact surface (sufficiently fine)
- Create a group of faces from the meshed contact surface
- Create a Submesh for each solid coming to contact:
- Geometry: contact group from the solid's geometry
- Algorithm: Import 1D-2D Elements from Another Mesh
- Hypothesis -> Source Faces: group of faces from the meshed contact surface
- Recompute meshes
- Build Compound mesh of solids with node merging on
- Build Compound mesh of compounded solids and the crack without merging nodes
- Export the latter Compound mesh to the MED file
2. Prepare the config file
- To list IDs of all blocks in the MED file:
$HOME/mofem_install/um/build_release/tools/read_med -med_file three_point_bending.med
- Check BLOCKSET for the contact interface in the config file
three_point_bending.cfg # Contact interface
[block_2] # Block ID in MED file
id=2004 # Block ID in the output *.h5m file
add=BLOCKSET # Block type
name=INT_CONTACT # Block name (starts exactly like this)
- Check BLOCKSET for springs in the same config file
# Springs on the loading frame and on the brick slice
[block_8] # Block ID in MED file
id=2005 # Block ID in the output *.h5m file
add=BLOCKSET # Block type
name=SPRING_BC # Block name (starts exactly like this)
user1=0 # Spring stiffness in normal direction [MPa]
user2=1e-2 # Spring stiffness in tangential directions [MPa]
NOTE: For the considered example the normal stiffness can be set initially to 0 (will be verified below), while the tangential one can be set to 1e-6 of the the Young's modulus of the solid to which these springs are attached. The calibration of this parameter will be discussed below.
3. Generate the <em>h5m</em> file
- Generate
three_point_bending.h5m file: $HOME/mofem_install/um/build_release/tools/read_med -med_file three_point_bending.med \
-meshsets_config three_point_bending.cfg \
-output_file three_point_bending.h5m
- The correct definition of all blocks can be verified by generating
vtk files for each one of them for visualisation in Paraview: $HOME/mofem_install/um/build_release/tools/meshset_to_vtk -my_file three_point_bending.h5m
4. Prepare the param file
Check following parameters in the param file param_file.petsc:
-my_contact_order 2 # default value 2
-my_lambda_contact_order 1 # default value 1 (not to be changed)
-my_cn_value 1.e4 # default value 1.0
-my_r_value 1.0 # default value 1.0 (not to be changed)
-my_cracked_body_block_set 1001 # Block ID of the fractured solid
NOTE: Definition of the Block ID of the fractured solid in the param file is required for the crack propagation with contact conditions.
NOTE: The initial recommended value of the contact augmentation parameter cn is the Young's modulus of the contacting solids (or the harmonic mean in case of different values). The further calibration of this parameter will be discussed below.
5. Calibrate the spring stiffness (tangential)
Springs can be used to eliminate the rigid body motion in a contact simulation under the pressure control. However, their stiffness is artificial and may affect stiffness/energy of the modelled system and, consequently, the resulting load factor for the crack propagation. Therefore, it is recommended to calibrate the stiffness in order to find the minimal value which permits to avoid the unwanted motion.
- Run elasticity simulation for the mesh with the initial crack inserted:
mpirun -np 2 $HOME/mofem_install/um/build_release/mofem_um_fracture_mechanics/crack_propagation \
-my_file three_point_bending.h5m \
-my_propagate_crack 0 \
2>&1 | tee log
- Correct identification of the contact interface can be verified by observing the following line in the output (
log file): [0] <inform> [petsc] Insert INT_CONTACT (id: 2004)
- The values of the material properties and spring stiffness can also be checked in the output (
log file): [0] <inform> [petsc] Material block 1001
[0] <inform> [petsc] Young's modulus 1.0900e+04
[0] <inform> [petsc] Poisson's ratio 2.0000e-01
[0] <inform> [petsc]
[0] <inform> [petsc]
[0] <inform> [petsc] Material block 1002
[0] <inform> [petsc] Young's modulus 1.0900e+04
[0] <inform> [petsc] Poisson's ratio 2.0000e-01
[0] <inform> [petsc]
[0] <inform> [petsc]
[0] <inform> [petsc] Spring block 2005
[0] <inform> [petsc] Normal stiffness 0
[0] <inform> [petsc] Tangent stiffness 0.01
- Check output for the number of SNES and KSP iterations:
[0] <inform> [petsc] 0 SNES Function norm 1.000000000000e+00
[0] <inform> [petsc] Residual norms for elastic_ solve.
[0] <inform> [petsc] 0 KSP Residual norm 6.593276968231e+00
[0] <inform> [petsc] 1 KSP Residual norm 2.265280970427e-12
[0] <inform> [Arc] F_lambda2 = 4838 lambda = 1
[0] <inform> [petsc] 1 SNES Function norm 5.018029619260e+01
[0] <inform> [petsc] Residual norms for elastic_ solve.
[0] <inform> [petsc] 0 KSP Residual norm 1.207850396738e+00
[0] <inform> [petsc] 1 KSP Residual norm 2.005935409209e-14
[0] <inform> [Arc] F_lambda2 = 4838 lambda = 1
[0] <inform> [petsc] 2 SNES Function norm 8.692812278752e-02
[0] <inform> [petsc] Residual norms for elastic_ solve.
[0] <inform> [petsc] 0 KSP Residual norm 1.398009933371e-02
[0] <inform> [petsc] 1 KSP Residual norm 1.078518671675e-16
[0] <inform> [Arc] F_lambda2 = 4838 lambda = 1
[0] <inform> [petsc] 3 SNES Function norm 1.799626823532e-09
[0] <inform> [petsc] Residual norms for elastic_ solve.
[0] <inform> [petsc] 0 KSP Residual norm 1.137713763876e-10
[0] <inform> [petsc] 1 KSP Residual norm 1.763018416865e-24
[0] <inform> [Arc] F_lambda2 = 4838 lambda = 1
[0] <inform> [petsc] 4 SNES Function norm 1.123275702526e-09
[0] <inform> [petsc] Nonlinear elastic_ solve converged due to CONVERGED_SNORM_RELATIVE iterations 4
[0] <inform> [petsc] solution fnorm = 1.96154083e+03
- Convert output
h5m files to vtk $HOME/mofem_install/um/build_release/tools/convert.py -np 2 out_spatial* out_contact*
- Visualize solution with Paraview (see Postprocessing for more details):
open out_spatial_0.vtk out_contact_0.vtk
If at least one of the following issues is observed:
- convergence is not achieved,
- more than 1-2 KSP iterations are required per one SNES iteration,
- an unwanted motion is present (in the visualisation of
out_spatial_0.vtk),
then the spring tangential stiffness may be increased (e.g. by a factor of 10) in the three_point_bending.cfg file. After that the three_point_bending.h5m needs to be regenerated and the convergence and rigid body motions are to be checked again.
6. Calibrate contact parameter <em>cn</em>
Contact augmentation parameter cn practically does not have effect on the solution of the contact problem (in the range of its values that permits convergence), but can have significant effect on the convergence itself. Therefore, it is recommended to calibrate this parameter to reduce the number of SNES iterations:
- Run elasticity simulation for the mesh with the crack for a different value of
cn: mpirun -np 2 $HOME/mofem_install/um/build_release/mofem_um_fracture_mechanics/crack_propagation \
-my_file three_point_bending.h5m \
-my_propagate_crack 0 \
-my_cn_value 1e3 \
2>&1 | tee log
- The value of the
cn parameter can be verified in the output (log file): [0] <inform> [petsc] ### Input parameter: -my_cn_value 1.0000e+03
- Check the number of SNES iterations
[0] <inform> [petsc] 0 SNES Function norm 1.000000000000e+00
[0] <inform> [petsc] Residual norms for elastic_ solve.
[0] <inform> [petsc] 0 KSP Residual norm 5.036499506247e+00
[0] <inform> [petsc] 1 KSP Residual norm 4.457975392822e-12
[0] <inform> [Arc] F_lambda2 = 4838 lambda = 1
[0] <inform> [petsc] 1 SNES Function norm 4.965527994591e+01
[0] <inform> [petsc] Residual norms for elastic_ solve.
[0] <inform> [petsc] 0 KSP Residual norm 3.500850747304e+00
[0] <inform> [petsc] 1 KSP Residual norm 1.334586168641e-13
[0] <inform> [Arc] F_lambda2 = 4838 lambda = 1
[0] <inform> [petsc] 2 SNES Function norm 3.040087225472e-10
[0] <inform> [petsc] Residual norms for elastic_ solve.
[0] <inform> [petsc] 0 KSP Residual norm 1.266813656664e-10
[0] <inform> [petsc] 1 KSP Residual norm 6.276204766206e-24
[0] <inform> [Arc] F_lambda2 = 4838 lambda = 1
[0] <inform> [petsc] 3 SNES Function norm 2.634697341395e-10
[0] <inform> [petsc] Nonlinear elastic_ solve converged due to CONVERGED_SNORM_RELATIVE iterations 3
[0] <inform> [petsc] solution fnorm = 1.96154083e+03
- Store the optimal value in the param file
param_file.petsc:
7. Run the crack propagation analysis
Once the following points are verified:
- Contact interface is correctly identified;
- Spring stiffness is calibrated, its value is stored in the
three_point_bending.cfg file and the three_point_bending.h5m file is updated;
- Contact parameter
cn is calibrated and the value is stored in theparam_file.petsc file;
one can kick-start the crack propagation analysis (using the usual ways of supervising the computation):
mpirun -np 2 $HOME/mofem_install/um/build_release/mofem_um_fracture_mechanics/crack_propagation \
-my_file three_point_bending.h5m \
2>&1 | tee log
8. Postprocessing
All the usual output files of the crack_propagation analysis are created and can be postprocessed in the standard way. Furthermore, the nodal interpolation of Lagrange multipliers (equivalent to contact pressure) is output to files out_contact_N.h5m, where N is the number of the step. Moreover, the values of the Lagrange multipliers and the normal gap at the gauss points of the contact interface are output to files out_contact_integ_pts_N.h5m.
- Convert output
h5m files to vtk using multiprocessing script convert.py: $HOME/mofem_install/um/build_release/tools/convert.py -np 2 out_spatial* out_contact*
- Files
out_spatial_N.vtk contain material coordinates (tag MESH_NODE_POSITIONS) and current coordinates (tag SPATIAL_POSITION), which can be used to compute the displacement field with the Calculator filter as DISPLACEMENT=SPATIAL_POSITION-MESH_NODE_POSITIONS
- Files
out_contact_N.vtk contain the nodal interpolation of the Lagrange multipliers equivalent to the contact pressure (tag LAMBDA_CONTACT)
- Files
out_contact_integ_pts_N.vtk contain values of Lagrange multipliers (tag LAMBDA_CONTACT) and the normal gap (tag GAP) at gauss points of the contact interface. Note that the Point Gaussian representation or alternatively the Glyph filter should be used for the visualisation of the gauss point data.
