21.5. V3.04.156: FORMA07 - Circular crack in an infinite medium under tension¶
This tutorial explores the 3D quasi-static problem of a circular crack in an infinite body under uniaxial tension.
We will explore two models in this tutorial:
Model A: A 3D FEM model with an explicit crack
Model B: A 3D X-FEM model with an implicit crack represented by X-FEM enrichment
21.5.1. Reference problem¶
Fig. 21.5.1 Schematic of domain with circular crack¶
Though the body is infinite, we will assume that it can be represented reasonably accurately by a cube of side \(10 a\), where \(a\) is the radius of the crack.
21.5.1.1. Geometry¶
We will consider a crack of radius \(a = 2\) m and a cube of side \(H = 20\) m.
21.5.1.2. Material properties¶
The material is linear elastic (isotropic and homogeneous):
Young’s modulus: \(E\) = 200 GPa
Poisson’s ratio: \(\nu\) = 0.3
21.5.1.3. Boundary conditions and loads¶
The body is in tension (\(P\) = 1 MPa).
21.5.2. Reference solution¶
21.5.2.1. Method used for the reference solution¶
For a circular crack of radius \(a\) in an infinite medium, subjected to a uniform traction \(P\) normal to the crack faces, the local energy release rate \(G(s)\) is independent of the curvilinear abscissa along the crack front \(s\) and is expressed as 1:
The stress intensity factor \(K_I(s)\) is given by Irwin’s formula:
21.5.2.2. Reference values¶
For the data used in this tutorial, \(K_I\) = 1.5957 MPa \(\sqrt{m}\) and \(G\) = 11.59 J \(m^{-2}\).
21.5.2.3. References¶
- 1
TADA, H., Paris, P., & Irwin, G. (2000). The analysis of cracks handbook. New York: ASME Press, 2, 1.
21.5.3. Model A: A 3D FEM model with an explicit crack¶
21.5.3.1. Creating and running the tutorial¶
21.5.3.1.1. Geometry and mesh with GIBI
¶

Fig. 21.5.2 Mesh generated with GIBI¶
The Code_Aster
installation contains an existing mesh in the share/aster/tests
directory in the file
forma07a.mmed
. This quadratic mesh of the body is provided in MED format. Taking into
account certain symmetries of the problem, only a quarter of the body is represented. The mesh
was generated with the GIBI software, and toroids are defined around the crack tip:
radius of the smallest torus: 0.12 m
radius of the largest torus: 0.53 m
The mesh contains the following face groups: FACE_SUP (TOP_FACE), FACE_INF (BOTTOM_FACE), FACE_AV (FRONT_FACE), FACE_LAT (LEFT_FACE), LEV_SUP (TOP_CRACK_FACE), LEV_INF (BOTTOM_CRACK_FACE). It also contains the edge group: LFF (CRACK_FRONT) and the node groups NFF1 (CRACK_FRONT_START), NFF2 (CRACK_FRONT_END), and D (a point that is fixed to prevent rigid body motions).
Barsoum elements are generated at the crack tip by moving the middle nodes of the edges of the elements touching the crack tip to a quarter of these edges.
21.5.3.1.2. Geometry and mesh with gmsh
¶

In this tutorial we will also pursue an alternative approach for mesh generation with gmsh
.
Note that we only use one torus in this model instead of the five
that were used in the GIBI
model. The torus minor radius is 0.5 m.
The mesh is saved in the file forma07a_gmsh.msh
in gmsh
2.2 format.
The Python script for creating and meshing the model is in
forma07a_gmsh.py
).
The annotated script is listed below.
import gmsh
import sys
import math
# Use OpenCASCADE kernel
gmsh.initialize()
model = gmsh.model
occ = model.occ
mesh = model.mesh
model.add("forma07a")
# Parameters
# Box size
L = 20
# initial crack opening
w = 0.005
# Initial crack length
a = 2
# Element size (at points, overwritten later when using
# transfinite interpolation)
lc = 20/10.0
# Create the top and bottom box
box_bot = occ.addBox(0, 0, -L/2, L/2, L/2, L/2)
box_top = occ.addBox(0, 0, 0, L/2, L/2, L/2)
# Join them while retaining the interace between them
box = occ.fragment([(3, box_bot)], [(3, box_top)])
# Add a torus around the crack tip location
torus = occ.addTorus(0, 0, 0, a, 0.5, tag = 1000, angle = math.pi/2)
# Partition the box with the torus
box_torus = occ.fragment(box[0], [(3, torus)])
# Make sure all new geometric entities are synchronized
occ.synchronize()
# Add the circular crack
# Crack tip
p0 = occ.addPoint(0, 0, 0, lc/10)
p1 = occ.addPoint(a, 0, 0, lc/10)
p2 = occ.addPoint(0, a, 0, lc/10)
l1 = occ.addCircleArc(p1, p0, p2)
# Crack face top/bottom points
p3 = occ.addPoint(0, 0, w/2, lc/5)
p4 = occ.addPoint(0, 0, -w/2, lc/5)
# Crack top face lines
l2 = occ.addLine(p3, p1)
l3 = occ.addLine(p2, p3)
# Crack bottom face lines
l4 = occ.addLine(p4, p1)
l5 = occ.addLine(p2, p4)
# Crack opening line
l6 = occ.addLine(p3, p4)
# Add curve loops
loop1 = occ.addCurveLoop([l2, l1, l3])
loop2 = occ.addCurveLoop([l4, l1, l5])
loop3 = occ.addCurveLoop([l6, l4, -l2])
loop4 = occ.addCurveLoop([l6, l5, -l3])
# Add faces
face1 = occ.addSurfaceFilling(loop1)
face2 = occ.addSurfaceFilling(loop2)
face3 = occ.addSurfaceFilling(loop3)
face4 = occ.addSurfaceFilling(loop4)
# Add surface loop
sloop1 = occ.addSurfaceLoop([face1, face4, face2, face3])
# Add volume
crack = occ.addVolume([sloop1])
# Subtract the crack from the box-torus partitioned geometry
box_torus_crack = occ.cut(box_torus[0], [(3, crack)], -1, True, True)
occ.removeAllDuplicates()
occ.synchronize()
# Add physical group for the mesh
box_with_crack = model.addPhysicalGroup(3, [1, 2, 3, 4])
model.setPhysicalName(3, box_with_crack, "box_with_crack")
# Add physical groups for the top and bottom faces
bot_face = model.addPhysicalGroup(2, [11])
top_face = model.addPhysicalGroup(2, [19])
model.setPhysicalName(2, bot_face, "bot_face")
model.setPhysicalName(2, top_face, "top_face")
# Add physical groups for defining the crack
crack_front = model.addPhysicalGroup(1, [32])
model.setPhysicalName(1, crack_front, "crack_front")
crack_top_surf = model.addPhysicalGroup(2, [22, 25])
crack_bot_surf = model.addPhysicalGroup(2, [8, 15])
model.setPhysicalName(2, crack_top_surf, "crack_top_surf")
model.setPhysicalName(2, crack_bot_surf, "crack_bot_surf")
# Add physical groups for the crack start/end vertices
crack_tip_xz = model.addPhysicalGroup(0, [26])
crack_tip_yz = model.addPhysicalGroup(0, [25])
model.setPhysicalName(0, crack_tip_xz, "crack_tip_xz")
model.setPhysicalName(0, crack_tip_yz, "crack_tip_yz")
# Add physical groups for symmetry faces and fixed corner
symm_face_xz = model.addPhysicalGroup(2, [7, 18, 16, 26])
symm_face_yz = model.addPhysicalGroup(2, [6, 17, 13, 24])
model.setPhysicalName(2, symm_face_xz, "symm_face_xz")
model.setPhysicalName(2, symm_face_yz, "symm_face_yz")
fixed_vert = model.addPhysicalGroup(0, [8])
model.setPhysicalName(0, fixed_vert, "fixed_vert")
# Update
occ.synchronize()
# Set element sizes near the crack (these may be overwritten by
# the transfinite interpolation definitions)
crack_pts = [(0, 26), (0, 25), (0, 9), (0, 4), (0, 20), (0, 30), (0, 23), (0, 32)]
mesh.setSize(crack_pts, lc/10)
# Set element sizes for transfinite interpolation
num_nodes = 10
radial_curves = [(1, 31), (1, 33), (1, 48), (1, 29), (1, 30), (1, 47)]
for curve in radial_curves:
mesh.setTransfiniteCurve(curve[1], num_nodes)
mesh.setSmoothing(curve[0], curve[1], 100)
num_nodes = 20
torus_curves = [(1, 22), (1, 42), (1, 18), (1, 37)]
for curve in torus_curves:
mesh.setTransfiniteCurve(curve[1], num_nodes)
mesh.setSmoothing(curve[0], curve[1], 100)
num_nodes = 10
crack_curves = [(1, 23), (1, 41)]
for curve in crack_curves:
mesh.setTransfiniteCurve(curve[1], num_nodes, coef=0.83)
mesh.setSmoothing(curve[0], curve[1], 100)
crack_curves = [(1, 17), (1, 38)]
for curve in crack_curves:
mesh.setTransfiniteCurve(curve[1], num_nodes, coef=1.2)
mesh.setSmoothing(curve[0], curve[1], 100)
num_nodes = 15
box_corner_edges = [(1, 16), (1, 34)]
for curve in box_corner_edges:
mesh.setTransfiniteCurve(curve[1], num_nodes, coef=1.2)
mesh.setSmoothing(curve[0], curve[1], 100)
box_corner_edges = [(1, 16)]
for curve in box_corner_edges:
mesh.setTransfiniteCurve(curve[1], num_nodes, coef=0.83)
mesh.setSmoothing(curve[0], curve[1], 100)
box_symm_edges = [(1, 8), (1, 4)]
for curve in box_symm_edges:
mesh.setTransfiniteCurve(curve[1], num_nodes, coef=1.5)
mesh.setSmoothing(curve[0], curve[1], 100)
# Set meshing options
gmsh.option.setNumber('Mesh.RecombinationAlgorithm', 1)
gmsh.option.setNumber('Mesh.Recombine3DLevel', 2)
gmsh.option.setNumber('Mesh.ElementOrder', 2)
gmsh.option.setNumber('Mesh.OptimizeNetgen', 1)
gmsh.option.setNumber('Mesh.MeshSizeMax', 5)
gmsh.option.setNumber('Mesh.MshFileVersion', 2.2)
# Generate mesh
mesh.generate(3)
mesh.removeDuplicateNodes()
# Save mesh
gmsh.write("forma07a_gmsh.msh")
# Finish
gmsh.finalize()
After the gmsh
mesh has been generated, you will need to convert it to MED format in
AsterStudy by running the command file listed below.
# Convert gmsh mesh to MED format
DEBUT(LANG='EN')
# Read gmsh mesh
mesh = LIRE_MAILLAGE(FORMAT = 'GMSH',
UNITE = 20)
# Read MED mesh
IMPR_RESU(RESU = _F(MAILLAGE = mesh),
UNITE = 80)
# Finish
FIN()
The resulting MED file will not have transferred the geometry groups with the names that
you had assigned to the physical groups in gmsh
. You will
have to change the names of the node, edge, and face groups from GMxx to readable ones
with the Salome-Meca
Mesh module. The vertex used to prevent rigid body motions
is fixed_vert, the crack front end points are crack_tip_xz and crack_tip_yz,
the crack front edge is crack_front, the top and bottom faces of the crack are
crack_top_surf and crack_bot_surf. The tractions will be applied on
top_face and bot_face and symmetry displacement boundary conditions will be
applied on symm_face_xz and symm_face_yz.
After the group names have been changed, save the mesh in MED format (you can
call it forma07a_alt.mmed
to distringuish it from the file provided in your Code_Aster
installation.
21.5.3.1.3. Creation of the command file without post-processing the crack¶
Start Salome_Meca
, then choose the AsterStudy module. Add the following steps to your
case:
Read the mesh (LIRE_MAILLAGE). Select the
forma07a_alt.mmed
(orforma07a.mmed
if you want to run the simulation on the mesh provided withCode_Aster
) mesh and the MED formatModify the mesh (MODI_MAILLAGE). Choose the mesh read previously and select reuse. Select the action ORIE_PEAU_3D to reorient the face normals outward for element groups symm_face_xz, symm_face_yz, top_face, bot_face, crack_top_surf, and crack_bot_surf (named FACE_AV, FACE_LAT, FACE_SUP, FACE_INF, LEV_INF and LEV_SUP in the mesh provided with
Code_Aster
).Assign finite elements (AFFE_MODELE). Choose the mechanical phenomenon and 3D modeling of continuous media (3D)
Define the material (DEFI_MATERIAU) and assign the material (AFFE_MATERIAU)
Define boundary conditions (AFFE_CHAR_MECA):
Symmetry on the symm_face_yz (‘FACE_LAT’) symmetry plane (prescribe DOF)
Symmetry on the symm_face_xz (‘FACE_AV’) symmetry plane (prescribe DOF)
Prevent rigid modes (prescribe DOF on the GROUP_NO fixed_vert (‘D’))
Apply traction on top_face (‘FACE_SUP’) and bot_face (‘FACE_INF’) with PRES_REP
Solve the elastic problem: Static mechanical analysis (MECA_STATIQUE)
For viewing with Results/Paravis:
Calculate the stress field extrapolated to the nodes (CALC_CHAMP, option ‘CONTRAINTE’ with the field ‘SIGM_NOEU’)
Calculation of the field of equivalent stresses (CALC_CHAMP, option ‘CRITERES’ with the field ‘SIEQ_NOEU’)
We will enrich the concept resulting from MECA_STATIQUE by using the same name when computing the stresses.
Print the results in MED format: (IMPR_RESU).
View the displacement and stress fields in Results/ParaViS.

21.5.3.1.4. Postprocessing for the fracture¶
In order to separate the calculation and the post-processing, you can add a new stage (New stage) to your case study.
Define the crack tip with DEFI_FOND_FISS starting with the edge group containing the crack front crack_front (‘LFF’) and the crack faces crack_top_surf and crack_bot_surf (‘LEV_INF’ and ‘LEV_SUP’).
Calculate of the energy release rate with CALC_G (OPTION = ‘CALC_K_G’). Complete the information on the THETA field:
the crack front FOND_FISS
the radii of the contours of the \(\theta\) field (R_INF, R_SUP), to be defined according to the mesh size.
Print the values of \(G\) (IMPR_TABLE).
Calculate \(K\) and \(G\) with POST_K1_K2_K3:
provide the crack tip concept created earlier
enter the parameter ABSC_CURV_MAXI
Print the results in an array (IMPR_TABLE)
Plot the values of \(G\) and \(K_I\) resulting from CALC_G and POST_K1_K2_K3 along curvilinear abscissa of the crack front (column ‘ABSC_CURV’)
21.5.3.1.5. Study of the influence of 3D smoothing¶
In CALC_G, add the keyword LISSAGE (SMOOTHING) with the keywords LISSAGE_G = ‘LAGRANGE’ and LISSAGE_THETA = ‘LAGRANGE’.
Compare the values of \(K\) and \(G\) with those obtained previously.
Note
When LISSAGE is not provided, the default value is LISSAGE_G = ‘LEGENDRE’ and LISSAGE_THETA = ‘LEGENDRE’.
Add the keyword NB_POINT_FOND = 5 (NUM_PTS_FOR_SMOOTHING) to “Lagrange” smoothing and observe the values of \(G\) or \(K_I\) along the crack front. Also compare the compute times for CALC_G with and without NB_POINT_FOND.
21.5.3.1.6. The Python command file for the mesh created with gmsh
¶
You can find the annotated command file used to run the simulation on the gmsh
mesh at
forma07a_alt.comm
). The file is listed below.
# Start
DEBUT(LANG='EN')
# Read mesh
mesh = LIRE_MAILLAGE(FORMAT = 'MED')
# Assign 3D mechanics finite elements
model = AFFE_MODELE(MAILLAGE = mesh,
AFFE = _F(TOUT = 'OUI',
PHENOMENE = 'MECANIQUE',
MODELISATION = '3D'))
# Reorient faces where tractions are applied to point outwards
mesh = MODI_MAILLAGE(reuse = mesh,
MAILLAGE = mesh,
ORIE_PEAU_3D = _F(GROUP_MA = ('top_face','bot_face',)))
# Convert to Barsoum elements at crcak tip
mesh = MODI_MAILLAGE(reuse = mesh,
MAILLAGE = mesh,
MODI_MAILLE = _F(OPTION = 'NOEUD_QUART',
GROUP_MA_FOND = 'crack_front'))
# Define material
steel = DEFI_MATERIAU(ELAS = _F(E = 2.E11,
NU = 0.3))
# Assign material
mater = AFFE_MATERIAU(MAILLAGE = mesh,
AFFE = _F(TOUT = 'OUI',
MATER = steel))
# Assign displacement bcs
disp_bc = AFFE_CHAR_MECA(MODELE = model,
DDL_IMPO = _F(GROUP_NO = 'fixed_vert',
DZ = 0.0),
FACE_IMPO = (_F(GROUP_MA = 'symm_face_xz',
DNOR = 0.0),
_F(GROUP_MA = 'symm_face_yz',
DNOR = 0.0)))
# Assign traction bcs
trac_bc = AFFE_CHAR_MECA(MODELE = model,
FORCE_FACE = (_F(GROUP_MA = 'top_face',
FZ = 1.E6),
_F(GROUP_MA = 'bot_face',
FZ = -1.E6)))
# Solve
result = MECA_STATIQUE(MODELE = model,
CHAM_MATER = mater,
EXCIT = (_F(CHARGE = trac_bc),
_F(CHARGE = disp_bc)))
# Add stresses to result
result = CALC_CHAMP(reuse = result,
RESULTAT = result,
CONTRAINTE = ('SIGM_ELNO', 'SIGM_NOEU'))
# Write the result
IMPR_RESU(FORMAT = 'MED',
UNITE = 80,
RESU = _F(RESULTAT = result))
# Define the crack front
crack = DEFI_FOND_FISS(MAILLAGE = mesh,
FOND_FISS = _F(GROUP_MA ='crack_front',
GROUP_NO_ORIG = 'crack_tip_xz',
GROUP_NO_EXTR = 'crack_tip_yz'),
SYME = 'NON',
LEVRE_SUP = _F(GROUP_MA = 'crack_top_surf'),
LEVRE_INF = _F(GROUP_MA = 'crack_bot_surf'),
)
# Calculate G
# Using Legendre interpolation (default)
r_max=0.5
r_min=0.2
G_LEG = CALC_G(RESULTAT = result,
OPTION = 'CALC_G',
THETA = _F(FOND_FISS = crack,
R_SUP = r_max,
R_INF = r_min),
LISSAGE = _F(LISSAGE_THETA = 'LEGENDRE',
LISSAGE_G = 'LEGENDRE',
DEGRE = 5))
# Write table
IMPR_TABLE(TABLE = G_LEG)
# Calculate G
# Using Lagrange interpolation (default)
G_LAG = CALC_G(RESULTAT = result,
OPTION = 'CALC_G',
THETA = _F(FOND_FISS = crack,
R_SUP = r_max,
R_INF = r_min),
LISSAGE = _F(LISSAGE_THETA = 'LAGRANGE',
LISSAGE_G = 'LAGRANGE'))
# Write table
IMPR_TABLE(TABLE = G_LAG)
# Calculate G
# Using node-to-node Lagrange interpolation (default)
G_LAGN = CALC_G(RESULTAT = result,
OPTION = 'CALC_G',
THETA = _F(FOND_FISS = crack,
R_SUP = r_max,
R_INF = r_min),
LISSAGE = _F(LISSAGE_THETA = 'LAGRANGE',
LISSAGE_G = 'LAGRANGE_NO_NO'))
# Write table
IMPR_TABLE(TABLE = G_LAGN)
# Calculate K with POST_K1_K2_K3
#-------------------------------
K = POST_K1_K2_K3(RESULTAT = result,
FOND_FISS = crack,
ABSC_CURV_MAXI = 1.5,
TYPE_MAILLAGE = 'LIBRE')
# Write table
IMPR_TABLE(TABLE = K)
# Save x-y data
C_G_LEG = RECU_FONCTION(TABLE = G_LEG,
PARA_X = 'ABSC_CURV',
PARA_Y = 'G')
C_G_LAG = RECU_FONCTION(TABLE = G_LAG,
PARA_X = 'ABSC_CURV',
PARA_Y = 'G')
C_G_LAGN=RECU_FONCTION(TABLE = G_LAGN,
PARA_X = 'ABSC_CURV',
PARA_Y = 'G')
IMPR_FONCTION(FORMAT='TABLEAU',
UNITE=31,
COURBE=(_F(FONCTION = C_G_LEG),
_F(FONCTION = C_G_LAG),
_F(FONCTION = C_G_LAGN)),
SEPARATEUR = ',')
# Verify reference values (for inifnite plate)
G_ref = 11.58648
K1_ref = 1.5957e6
TEST_TABLE(REFERENCE = 'ANALYTIQUE',
PRECISION = 8.0E-3,
VALE_CALC = 1607011.55098,
VALE_REFE = 1.595700E6,
NOM_PARA = 'K1',
TYPE_TEST = 'MAX',
TABLE=K)
# Finish
FIN();
21.5.3.2. Verification tests¶
You can verify that the value of \(K_I\) computed by POST_K1_K2_K3 is within 0.8% of 1.5957 MPa \(m^{1/2}\).
21.5.4. Model B: A 3D X-FEM model with an implicit crack represented by X-FEM enrichment¶
The problem considered is the same as that studied in Model A. However, unlike the previous model where the crack was meshed explicitly, we consider in this part a mesh that does not contain a pre-crack. The crack is thus not explicitly meshed and we need to use the X-FEM method.
21.5.4.1. Creating and running the tutorial¶
21.5.4.1.1. Refinement of the mesh¶
The objective of this first stage is to use HOMARD
to obtain a refined mesh starting from a
relatively coarse mesh.
You can create the coarse mesh with Salomé
(a quarter of the structure is modeled) or use that
provided with the Code-Aster
installation (forma07b.mmed
).
Warning
In order to refine the mesh, it is necessary to use a Python loop. The command file for this stage is
not editable graphically with AsterStudy
. We will therefore use the text mode of
AsterStudy
for the first stage of computation.
Launch Salome-Meca
, then choose AsterStudy.
Go to the preferences (File → Preferences → AsterStudy) and check External editor.
Create a new stage for your calculation (New Stage) that you can rename (Refinement) if you want. Left click on this new step and click on Edit. A window will appear on the left and you can edit your command file.
First, it will be necessary to calculate the number of iterations necessary to obtain the refined mesh. To do this, you must first import the math library and then set up parameters:
import math # Size of the largest element in the crack region init_size = 5. # Crack depth crack_depth = 2. # Target element size after refinement target_size = crack_depth/20. # Compute number of refinement iterations n = (log(init_size) - log(target_size)) / log(2.) num_refine = int(n) num_iter = num_refine + 1 # The actual element size after refine (must be close to target_size) final_size = init_size / 2**num_refine # Radius of refinement zone around crack front r_refine = 5*final_size
Initialize the set of vectors that will be used in the refinement loop:
mesh = [None]*(num_iter+1) referr = [None]*num_iter crack = [None]*num_iter
Read the initial mesh in MED format MED (LIRE_MAILLAGE):
mesh[0] = LIRE_MAILLAGE(UNITE = 20, FORMAT = 'MED')
Choose the
forma07b.mmed
mesh by adding it in the Data Files tab under unit 20.Create the refinement loop:
# Refinement loop for i_refine in range(num_iter) : # Define the crack crack[i_refine] = DEFI_FISS_XFEM(MAILLAGE = mesh[i_refine], TYPE_DISCONTINUITE = 'FISSURE', DEFI_FISS = _F(FORM_FISS = 'ELLIPSE', DEMI_GRAND_AXE = crack_depth, DEMI_PETIT_AXE = crack_depth, CENTRE = (0,0,0,), VECT_X = (1.,0.,0.), VECT_Y = (0.,1.,0.))) # Compute refinement criterion referr[i_refine] = RAFF_XFEM(FISSURE = crack[i_refine], TYPE = 'ZONE', RAYON = r_refine) # Update mesh name mesh[i_refine+1] = CO('mesh_%d' % (i_refine+1)) # Refine the mesh MACR_ADAP_MAIL(ADAPTATION = 'RAFFINEMENT', CHAM_GD = referr[i_refine], CRIT_RAFF_ABS = 0.5, DIAM_MIN = target_size, MAILLAGE_N = mesh[i_refine], MAILLAGE_NP1 = mesh[i_refine+1])
To visualize the new mesh, print it in MED format:
IMPR_RESU(FORMAT = 'MED', UNITE = 80, RESU = _F(MAILLAGE = mesh[num_iter]))
Name the output mesh
forma07b_refined.mmed
by adding it in the Data Files tab under unit 80.

Fig. 21.5.3 Groups defined on the provided mesh¶
The mesh produced by these commands may contain unwanted SEG3
and over-constrained face
elements which have to be removed manually before the mesh can be used in the next stage
for mechanics calculations. Procedures for cleaning the mesh are discussed in this
blog post.
If you want to skip this step, you can also directly use the refined mesh provided in the
Code_Aster
repository: forma07b.41
.
In the provided meshes (visualize them in Salomé
), the names of the mesh groups are the following:
group of nodes: D,
groups of faces: FACE_AV, FACE_LAT, FACE_SUP, FACE_INF
volume group: CUBE
In the mesh generated with Salome-Meca
these groups are called:
group of nodes: fixed_vert,
groups of faces: symm_face_xz, symm_face_yz, top_face, bot_face
The annotated command file used to run the refinement process is
forma07b_alt.comm
). The file is listed below.
#---------------------------------------------------------------
# For XFEM: Stage 1 - Create refined mesh near crack front
# You will need to run these commands with astk because the
# AsterStudy GUI does not yet support full Python functionality
#---------------------------------------------------------------
# Start
DEBUT(LANG='EN')
import math
# Size of the largest element in the crack region
init_size = 5.
# Crack depth
crack_depth = 2.
# Target element size after refinement
target_size = crack_depth/20.
# Compute number of refinement iterations
n = (log(init_size) - log(target_size)) / log(2.)
num_refine = int(n)
num_iter = num_refine + 1
# The actual element size after refine (must be close to target_size)
final_size = init_size / 2**num_refine
# Radius of refinement zone around crack front
r_refine = 5*final_size
# Initialization
mesh = [None]*(num_iter+1)
referr = [None]*num_iter
crack = [None]*num_iter
# Read the initial unrefined mesh
mesh[0] = LIRE_MAILLAGE(UNITE=20,
FORMAT='MED')
# Refinement loop
for i_refine in range(num_iter) :
# Define the crack
crack[i_refine] = DEFI_FISS_XFEM(MAILLAGE = mesh[i_refine],
TYPE_DISCONTINUITE = 'FISSURE',
DEFI_FISS = _F(FORM_FISS = 'ELLIPSE',
DEMI_GRAND_AXE = crack_depth,
DEMI_PETIT_AXE = crack_depth,
CENTRE = (0,0,0,),
VECT_X = (1.,0.,0.),
VECT_Y = (0.,1.,0.)))
# Compute refinement criterion
referr[i_refine] = RAFF_XFEM(FISSURE = crack[i_refine],
TYPE = 'ZONE',
RAYON = r_refine)
# Update mesh name
mesh[i_refine+1] = CO('mesh_%d' % (i_refine+1))
# Refine the mesh
MACR_ADAP_MAIL(ADAPTATION = 'RAFFINEMENT',
CHAM_GD = referr[i_refine],
CRIT_RAFF_ABS = 0.5,
DIAM_MIN = target_size,
MAILLAGE_N = mesh[i_refine],
MAILLAGE_NP1 = mesh[i_refine+1])
# Write the final mesh
IMPR_RESU(FORMAT = 'MED',
UNITE = 80,
RESU = _F(MAILLAGE = mesh[num_iter]))
# Finish
FIN();
These commands produce the mesh shown on the right in the figure below.
Fig. 21.5.4 Refined mesh.¶
21.5.4.1.2. Creation of the command file without postprocessing of the crack¶
Reading of the refined mesh and definition of the non-enriched model:
Launch
Salome-Meca
, then choose the AsterStudy module. Add the following steps to your case.Read the mesh (LIRE_MAILLAGE). Select the refined mesh you created and the MED format
Modify the mesh (MODI_MAILLAGE). Choose the mesh read previously and select reuse. Select the action ORIE_PEAU_3D to reorient the normals for the face groups symm_face_xz, symm_face_yz, top_face, bot_face (FACE_AV, FACE_LAT, FACE_SUP, FACE_INF)
Assign finite element type (AFFE_MODELE). Choose the mechanical phenomenon and 3D)
Definition of the crack and X-FEM elements in the Fracture and Fatigue tab:
Define the crack DEFI_FISS_XFEM. It is preferable to use the catalog of cracks (FORM_FISS = ‘ELLIPSE’)
Modify the model to incorporate X-FEM elements (MODI_MODELE_XFEM)
Definition of the material, the boundary conditions and solution of the mechanical problem:
Define a material (DEFI_MATERIAU) and assign the material (AFFE_MATERIAU)
Define the boundary conditions on the enriched model (AFFE_CHAR_MECA):
Symmetry on the symm_face_yz (‘FACE_LAT’) symmetry plane (prescribe DOF)
Symmetry on the symm_face_xz (‘FACE_AV’) symmetry plane (prescribe DOF)
Prevent rigid modes (prescribe DOF on fixed_vert (‘D’))
Apply traction on top_face (‘FACE_SUP’) and bot_face (‘FACE_INF’) with PRES_REP
Solve the elastic problem: Static mechanical analysis (MECA_STATIQUE) with the enriched model
21.5.4.1.3. Postprocessing of displacements and stresses with X-FEM and visualization¶
In order to separate the calculation and the post-processing, you can add a new stage (New stage) to your case study.
In the Fracture and Fatigue tab, create a visualization mesh (POST_MAIL_XFEM)
In the Model Definition tab, create a model for visualization (AFFE_MODELE) on the mesh created for visualization
In the Fracture and Fatigue tab, create a results field on the visualization mesh for X-FEM (POST_CHAM_XFEM)
For viewing with Results or ParaViS:
Calculate the stress field extrapolated to the nodes (CALC_CHAMP, option ‘CONTRAINTE’ with the field ‘SIGM_NOEU’)
Calculate the field of equivalent stresses (CALC_CHAMP, option ‘CRITERES’ with the field ‘SIEQ_NOEU’). For that, you will need to enrich the concept resulting from POST_CHAM_XFEM by reusing its name.
Print the results in MED format (IMPR_RESU)
![]() |
![]() |
21.5.4.1.4. Additional stage for fracture mechanics postprocessing¶
For post-processing the crack, you can also add a new stage (New stage) to your current case study.
For all of the following steps, you will have to go to the Fracture and Fatigue tab.
Calculation of \(G\) with CALC_G:
Calculate the rate of energy release \(G\) and the stress intensity factors with CALC_G (OPTION = ‘CALC_K_G’) and Legendre smoothing of degree 5.
Use the result of the static calculation (not to be confused with the result created for visualization) with POST_CHAM_XFEM (RESULTAT).
Complete the information on the THETA field:
the crack tip (FISSURE)
the radii of the integration contours of the \(\theta\) field (R_INF, R_SUP), to be defined according to the mesh size.
Print the values of \(G\) (IMPR_TABLE in the Output tab).
Calculation of \(K\) and \(G\) with POST_K1_K2_K3:
Calculate \(K\) and \(G\) with POST_K1_K2_K3:
use the result of MECA_STATIQUE (RESULTAT)
choose the crack
enter the parameter ABSC_CURV_MAXI
Print the results in a table.
Compare the results to the analytical solution. Compare the compute times for the two commands. Try LAGRANGE smoothing. What do you observe?
In order to improve the results, test the uniform distribution of \(n\) points along the crack front (CALC_G / THETA / NB_POINT_FOND) with various values of \(n\).
You can also explore other avenues, e.g., refinement/enrichment on several layers.
21.5.4.1.5. The Python command files¶
The annotated command file used to run the mechanical simulation process is
forma07b.com1
. The file is listed below.
#---------------------------------------------------------------
# For XFEM: Stage 2 - Run mechanics simulation
#---------------------------------------------------------------
# Continue
#POURSUITE(LANG = 'EN')
DEBUT(LANG = 'EN')
# Read the refined mesh (created in Stage 1)
mesh = LIRE_MAILLAGE(FORMAT = 'MED',
UNITE = 41)
# Change the face orientations
mesh = MODI_MAILLAGE(reuse = mesh,
MAILLAGE = mesh,
ORIE_PEAU_3D = _F(GROUP_MA=('symm_face_xz', 'top_face',
'symm_face_yz', 'bot_face')))
# Define material
steel = DEFI_MATERIAU(ELAS = _F(E = 2e+11,
NU = 0.3))
# Assign material
mater = AFFE_MATERIAU(AFFE = _F(MATER = steel,
TOUT = 'OUI'),
MAILLAGE = mesh)
# Assign finite elements
model = AFFE_MODELE(AFFE = _F(MODELISATION = '3D',
PHENOMENE = 'MECANIQUE',
TOUT = 'OUI'),
MAILLAGE = mesh)
# Define crack
crack_depth = 2.
crack = DEFI_FISS_XFEM(DEFI_FISS =_F(CENTRE = (0.0, 0.0, 0.0),
DEMI_GRAND_AXE = crack_depth,
DEMI_PETIT_AXE = crack_depth,
FORM_FISS = 'ELLIPSE',
VECT_X = (1.0, 0.0, 0.0),
VECT_Y = (0.0, 1.0, 0.0)),
MAILLAGE = mesh)
# Enrich the model with the crack
model_X = MODI_MODELE_XFEM(FISSURE = crack,
MODELE_IN = model)
# Boundary conditions
bcs = AFFE_CHAR_MECA(DDL_IMPO =(_F(DZ = 0.0,
GROUP_NO = 'fixed_vert'),
_F(DY = 0.0,
GROUP_MA = 'symm_face_xz'),
_F(DX = 0.0,
GROUP_MA = 'symm_face_yz')),
MODELE = model_X,
PRES_REP = _F(GROUP_MA = ('top_face', 'bot_face'),
PRES = -1000000.0))
# Solve
result = MECA_STATIQUE(CHAM_MATER = mater,
EXCIT = _F(CHARGE = bcs),
MODELE = model_X)
# Post-processing for visualization
mesh_V = POST_MAIL_XFEM(MODELE = model_X)
model_V = AFFE_MODELE(AFFE = _F(MODELISATION = '3D',
PHENOMENE = 'MECANIQUE',
TOUT = 'OUI'),
MAILLAGE = mesh_V)
result_V = POST_CHAM_XFEM(MODELE_VISU = model_V,
RESULTAT = result)
result_V = CALC_CHAMP(reuse = result_V,
CONTRAINTE = 'SIGM_NOEU',
CRITERES = 'SIEQ_NOEU',
RESULTAT = result_V)
# Write result for viz
IMPR_RESU(FORMAT = 'MED',
RESU = _F(RESULTAT = result_V),
UNITE=81)
# Finish
FIN();
The post-processing stage is in forma07b.com2
,
listed below.
#---------------------------------------------------------------
# For XFEM: Stage 3 - Fracture mechanics calcs
#---------------------------------------------------------------
# Continue
POURSUITE(LANG = 'EN')
# Parameters of contour integration
h = 0.078
r_min = 2.0 * h
r_max = 5.0 * h
# Calculate G
G_LEG = CALC_G(LISSAGE = _F(DEGRE = 5,
LISSAGE_G = 'LEGENDRE',
LISSAGE_THETA = 'LEGENDRE'),
OPTION = 'CALC_K_G',
RESULTAT = result,
THETA = _F(FISSURE = crack,
R_INF = r_min,
R_SUP = r_max))
G_LAG = CALC_G(LISSAGE = _F(LISSAGE_G = 'LAGRANGE',
LISSAGE_THETA = 'LAGRANGE'),
OPTION = 'CALC_K_G',
RESULTAT = result,
THETA = _F(FISSURE = crack,
R_INF = r_min,
R_SUP = r_max))
G_LAG_20 = CALC_G(LISSAGE = _F(LISSAGE_G = 'LAGRANGE',
LISSAGE_THETA = 'LAGRANGE'),
OPTION = 'CALC_K_G',
RESULTAT = result,
THETA = _F(FISSURE = crack,
NB_POINT_FOND = 20,
R_INF = r_min,
R_SUP = r_max))
# Save
IMPR_TABLE(TABLE = G_LEG)
IMPR_TABLE(TABLE = G_LAG)
IMPR_TABLE(TABLE = G_LAG_20);
# Calculate K
PK = POST_K1_K2_K3(ABSC_CURV_MAXI = r_max,
FISSURE = crack,
RESULTAT = result)
PK_20 = POST_K1_K2_K3(ABSC_CURV_MAXI = r_max,
FISSURE = crack,
NB_POINT_FOND = 20,
RESULTAT = result)
# Save
IMPR_TABLE(TABLE = PK)
IMPR_TABLE(TABLE = PK_20)
# Verification tests
K1_ref = 1.5957e6
TEST_TABLE(REFERENCE = 'ANALYTIQUE',
PRECISION = 0.06,
VALE_CALC = 1682482.3636702264,
VALE_REFE = 1.595700E6,
NOM_PARA = 'K1',
TYPE_TEST = 'MAX',
TABLE = G_LEG)
TEST_TABLE(REFERENCE = 'ANALYTIQUE',
PRECISION = 0.05,
VALE_CALC = 1671075.8789433436,
VALE_REFE = 1.595700E6,
NOM_PARA = 'K1',
TYPE_TEST = 'MAX',
TABLE = G_LAG_20)
TEST_TABLE(REFERENCE = 'ANALYTIQUE',
PRECISION = 0.01,
VALE_CALC = 1607458.569834848,
VALE_REFE = 1.595700E6,
NOM_PARA = 'K1',
TYPE_TEST = 'MAX',
TABLE = PK_20)
# Finish
FIN();
21.5.4.2. Verification tests¶
You can verify that the value of \(K_I\) computed by POST_K1_K2_K3 is within 1% of 1.5957 MPa \(m^{1/2}\).