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

../../_images/fig1_v3.04.156.svg

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:

\[G(s) = 4 P^2 a\,\frac{4(1-\nu^2)}{\pi E}\]

The stress intensity factor \(K_I(s)\) is given by Irwin’s formula:

\[G(s) = \frac{(1-\nu^2)}{E}\,K_I^2 \quad \implies \quad K_I(s) = \frac{2P\sqrt{a}}{\sqrt{\pi}}\]

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

../../_images/fig2_v3.04.156.png

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

../../_images/fig4_v3.04.156.png

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 (or forma07a.mmed if you want to run the simulation on the mesh provided with Code_Aster) mesh 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 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.

../../_images/fig5_v3.04.156.png

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’)

../../_images/forma07a_8.svg

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.

../../_images/forma07a_31.svg

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.

../../_images/fig3_v3.04.156.png

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.

../../_images/fig7_v3.04.156.svg

Fig. 21.5.4 Refined mesh.

21.5.4.1.2. Creation of the command file without postprocessing of the crack

  1. 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)

  2. 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)

  3. 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)

../../_images/fig8_v3.04.156.png ../../_images/fig9_v3.04.156.png

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.

  1. 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).

  2. 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.

../../_images/forma07b_1_8_G.svg ../../_images/forma07b_1_8_K.svg

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}\).