Standard nonlinear finite element solution procedures apply for creating a crack model with proper solution-control settings, loadings and boundary conditions.

SMART uses higher-order tetrahedral element SOLID187, and the finite element model must be meshed with that element.

To create a finite element model with an initial crack, you can use Ansys Workbench, Ansys Mechanical, Mechanical APDL, or any third-party meshing tools that work with Mechanical APDL.

Fracture mechanics deals with cracks (defects), and a singularity always exists around the crack tip/front. The crack-tip/-front mesh is therefore of utmost importance in a crack analysis, as stress-analysis and fracture-parameters calculation accuracy depend on the crack mesh. Size and shape differences in the elements ahead of and behind the crack tip/front affect the accuracy of the fracture-parameters calculation, and therefore the crack-growth simulation.

For more information, see Understanding How Fracture-Mechanics Problems Are Solved and Fracture-Parameter Calculation Process.

The SMART crack-growth method uses either J-integral or stress-intensity factors (SIFS) as the fracture parameter (driving force) and the criteria for crack-growth calculation.

For each crack, only one fracture parameter can be specified. The parameter must be consistent with the specified crack-growth criterion (CGROW,FCOPTION).

The CINT command initiates the fracture-parameter calculation and specifies options for the calculation.

For each crack, only one fracture criterion can be specified. The criterion parameter must be consistent with the defined fracture-parameter calculation (CINT).

For fatigue crack-growth analysis, you can include a crack-closure effect (with certain crack-growth criteria).

If using J-integral as a fracture parameter, the crack is assumed to grow along the initial direction. It is therefore suited for Mode I crack-growth only.

SMART supports the following options for the calculations of equivalent stress-intensity factor and kink angle. These options are valid only when stress-intensity factors (SIFS) are specified for fracture-parameter calculations by issuing CINT,TYPE,SIFS. Note that any negative value of (mode-I stress-intensity factor) is non-physical and is treated as zero for the calculations.

References

The following references are cited:

When a crack begins to grow, a new mesh is generated. To correctly predict crack-growth behavior, any external loading associated with the remeshing volume must be mapped to the new meshes from the old meshes. Boundary conditions not in the remeshing volume are maintained as initially defined.

You can use nodes or element components to define initial boundary conditions; however, Mechanical APDL does not maintain them after remeshing. It is therefore best to avoid using those components for any purpose (such as applying new boundary conditions at later load steps) after crack-growth has begun.

 

3.3.1.5.9. Enhancing Crack Surfaces with Cohesive Zone Elements

When a crack grows, the crack surfaces can be enhanced using INTER204 3D 16-node cohesive interface elements to simulate a cohesive effect (for bonding tractions and to prevent crack-surface penetration). The cohesive elements are generated automatically during the first remeshing around the associated crack.

The process of generating the cohesive elements requires close topological matching between the upper and lower crack surfaces. Both surfaces should have the same number of vertices and edges on both sides, and each vertex/edge should have a counterpart on the opposite surface. Each pair of the topological entities should be as close as possible.

Generate INTER204 elements on a new crack surface as follows:

CGROW,CSFL,CZM,TB_CZM_ID,ESYS_ID

where the cohesive zone material model parameters are defined via TB_CZM_ID with TB,CZM,,,,BILI, TB,CZM,,,,CEXP, TB,CZM,,,,CLIN, or TB,CZM,,,,REXP.

If the parameters are not specified or are set to zero, Mechanical APDL uses linear material behavior with a penalty slope to prevent surface penetration on the cohesive interface as the default cohesive model type (TB,CZM,,,,CLIN) . The program sets reasonable CZM parameters based on the solid elements around the crack fronts/surfaces.

The penalty stiffness for CLIN is determined by:

where:

= characteristic element length around the crack (the average element size)

= effective elastic modulus, calculated by:

where subscripts and denote upper and lower crack-surface materials, respectively.

The penalty stiffness for CLIN is evaluated at the reference temperature only. For temperature-dependent materials, define the CZM parameters (TB,CZM,,,BILI / CEXP / CLIN / REXP) to obtain an accurate result.

When TB,CZM,,,,FRIC / REXP is applied and no valid local coordinate system (LCS) ID is specified (ESYS_ID in CGROW,CSFL,CZM,MAT_ID,ESYS_ID), the program creates an internal LCS for CZM on crack surfaces. Information about this internal LCS such as (ID, locations, and orientations) is provided via echo messages. The LCS ID is assigned to the ESYS attribute of the CZM elements inserted on crack surfaces; the program calculates and displays tractions and separations in CZM elements in the element coordinate system. The LCS is crucial for maintaining consistency in CZM during remeshings.

When neither FRIC or REXP is specified, a valid ESYS_ID in the CGROW,CSFL command applies, although no internal LCS is created when ESYS_ID is invalid or unavailable.

When INTER204 elements are enhanced to a crack surface with pressure loads, Mechanical APDL accounts for the contributions from both the cohesive tractions and pressure loads during fracture-parameter calculations.

When applying predefined pressures to a crack surface, and no pressures are defined for the new open crack surface (CGROW,CSFL,PRES), a fine mesh around the original crack front is recommended. Issue CGROW,RMCO,ESIZE,Par1,COMP,Par2, where Par1 is the expected mesh size and Par2 is the node-component name of the crack front. (See CGROW,RMCO,ESIZE).

If a crack has existing cohesive zone elements on its surfaces, enable crack-surface loading (CGROW,CSFL,CZM) for the associated crack. All existing cohesive zone elements are maintained during remeshing.

If you enable crack-surface loading, Mechanical APDL detects the element type, material model, and element coordinate system for cohesive zone elements attached to the crack surfaces, then applies additional cohesive zone elements on new opened crack surfaces. (If your specified values are inconsistent with those expected by the program, a message appears and your values are not used; however, the program still applies cohesive zone surface loading to the crack surface using the values that it detects.)

For more information, see:

• Exponential Cohesive Zone Material for Preventing Surface Penetration in the Material Reference

• Linear Cohesive Zone Material for Preventing Surface Penetration in the Material Reference

• CGROW,CSFL,CZM

Fatigue crack-growth under nonproportional loading is a complex phenomenon and can be difficult to analyze. The varying mode mix occurring at different loading states, and the varying stress ratio , are among the challenges. SMART-based crack-growth analysis supports nonproportional loading problems via cyclic loading (CLOAD).

CLOAD,DEFINE,BEGIN
*DIM,_cycload1,TABLE,IMAX 
*DIM,_cycload2,TABLE,IMAX
CLOAD,DEFINE,END

CLOAD,CYCNUM,100   ! Total number of cycles

 

3.3.1.6.2. Applying Nonproportional Crack-Growth

SMART-based fatigue crack-growth offers options for defining the crack-propagation direction under nonproportional cycle loading, and for defining effective stress-intensity factors for fatigue-life calculation:

• 3.3.1.6.2.1. Defining the Stress-Intensity-Factor Range
• 3.3.1.6.2.2. Defining the Crack-Growth Direction

 

3.3.1.6.2.1. Defining the Stress-Intensity-Factor Range

The crack-growth rate is typically dependent on the stress-intensity-factor range calculated from the maximum and minimum stress-intensity factor during a loading cycle.

With nonproportional loading, different fracture-modes typically have the corresponding minimum and maximum stress-intensity-factors occurring at different loading states:

Figure 3.9: Nonproportional Loading Example

The stress-intensity-factor range drives fatigue crack-growth and must therefore be calculated from the history of a loading cycle. The following methods are available to introduce the stress ratio’s effect:

Effective stress-intensity-factor range method

The simplest method for calculating the stress-intensity-factor range is to use the difference between the maximum and minimum equivalent stress-intensity-factors at each crack-front node, determined by:

where:

= equivalent stress-intensity-factor range at each crack-front node

= stress ratio at each crack-front node

and = minimum and maximum equivalent stress-intensity factors at each crack-front node, respectively

The effective stress-intensity-factor range at each crack-front node is defined as

If there is no crack closure effect, the closure parameter , and thus .

At the end of each loading cycle analysis, the effective stress-intensity-factor range is used in the fatigue law as:

where are the material parameters of fatigue crack-growth law.

This command introduces the effective stress-intensity-factor range method in nonproportional fatigue problems:

CGROW,NPLOAD,METH,DKEFF

This method overwrites any previously defined stress ratio (CGROW,FCG,SRAT). It works with any available fatigue law.

Stress-intensity-factor range method based on mode-separation

In some mixed-mode problems, mode-separation may be required to analyze complex nonproportional fatigue crack-growth. You can specify a stress-intensity-factor range based on mode-separation. The stress-intensity-factor range is formulated from the ranges of mode-separated stress-intensity factors at each crack front node as:

where:

and = stress-intensity-factor ranges for mode I and mode II at each crack-front node, respectively

and = maximum stress-intensity-factors for mode I and mode II at each crack-front node, respectively

and = minimum stress-intensity-factors for mode I and mode II at each crack-front node, respectively

The equivalent stress-intensity-factor range is formulated as a function of the mode-separated stress-intensity-factor ranges , , and as:

where:

= function for the equivalent stress-intensity-factor formulation (CGROW,SOPT,KEQV).

= nonproportional crack-growth direction (CGROW,NPLOAD,DIRECTION)

At the end of each loading-cycle analysis, the equivalent stress-intensity-factor range is used in the fatigue law as:

where:

= stress ratio that you define (CGROW,FCG,SRAT). Default = 0.

= material parameters for the fatigue crack-growth law.

This command introduces the stress-intensity-factor method based on mode-separation in nonproportional fatigue problems:

CGROW,NPLOAD,METH,DKMOD

This method works with any available fatigue law.

Total effective stress-intensity-factor method

In special cases, the crack-growth rate can depend on a superposition of the minimum and maximum equivalent stress-intensity factors. You can use the summation of the maximum and minimum equivalent stress-intensity-factors at each crack front node as:

where:

= superposed stress-intensity factor

and = minimum and maximum equivalent stress-intensity factors, respectively

At the end of each loading-cycle analysis, the superposed stress-intensity factor is used in the fatigue law as:

where:

= stress ratio that you define (CGROW,FCG,SRAT). Default = 0.

= material parameters for the fatigue crack-growth law.

This command introduces the total effective stress-intensity-factor method in nonproportional fatigue problems:

CGROW,NPLOAD,METH,TKEFF

This method works with any available fatigue law.

Total stress-intensity-factor method based on mode-separation

Also in special cases, the crack-growth rate can depend on a superposition of the minimum and maximum stress-intensity-factors. You can use the summation of the maximum and minimum mode-separated stress-intensity factors at each crack-front node as:

where and are the superposed stress-intensity factors for mode I and mode II, respectively.

The equivalent stress-intensity-factor range is formulated as a function of the superposed stress-intensity-factor ranges as:

where:

= function for the equivalent stress-intensity-factor formulation. (CGROW,SOPT,KEQV).

= nonproportional crack-growth direction (CGROW,NPLOAD,DIRECTION).

At the end of each loading-cycle analysis, the superposed stress-intensity factor is used in the fatigue law as:

where:

= stress ratio that you define (CGROW,FCG,SRAT). Default = 0.

= material parameters for the fatigue crack-growth law.

This command introduces the total stress-intensity-factor method based on mode-separation in nonproportional fatigue problems:

CGROW,NPLOAD,METH,TKMOD

You can combine this method with your own stress-ratio input (CGROW,FCG,SRAT).^[2]It works with any available any fatigue law.

 

3.3.1.6.2.2. Defining the Crack-Growth Direction

The varying mode mixity occurring during a nonproportional loading cycle makes it difficult to predict crack-growth direction. SMART-based fatigue crack-growth offers the following methods for determining the crack-growth direction:

Maximum-load method

This method determines the angle based on the mode mix at the instant when the equivalent stress-intensity factor is maximum (). The crack-propagation direction is determined by the angle formulation, a function of the stress-intensity factors , , and at each crack-front node:

where is the function for the crack-propagation-direction angle formulation (CGROW,SOPT,KANG).

This command determines the crack-growth angle using the maximum-load method in nonproportional fatigue problems:

CGROW,NPLOAD,DIRECTION,KMAX

Minimum-load method

This method determines the angle based on the angle calculation at the instant when the equivalent stress-intensity factor is minimum (). The crack-propagation direction is determined by the angle formulation, a function of the stress-intensity factors , , and at each crack-front node:

where is the function for the crack-propagation-direction angle formulation (CGROW,SOPT,KANG).

This command determines the crack-growth angle using the minimum-load method in nonproportional fatigue problems:

CGROW,NPLOAD,DIRECTION,KMIN

Local kink-tip stress-intensity-factors method

This method uses the kink-tip stress-intensity factors:[14]

where:

and = kink-tip stress intensity-factors at each crack-front node

= local crack angles at each crack-front node

= solution time

and the different local crack angles of and the solution time, at each crack-front node.

The kink-tip stress-intensity-factor ranged is determined via:

where:

= nodal kink-tip stress-intensity-factor range at each crack-front node

and = kink-tip stress-intensity factors at the minimum and maximum load points at each crack-front node, respectively, where and are located in the loading cycle at each crack-front node (as shown in Figure 3.9: Nonproportional Loading Example)

Although maximization of is a good agreement with simple experimental test sets, a better approach is the maximization of a coupled function:

where is a fitting parameter controlling the weight of the function (default = 0.5). The angle at which the crack-propagation direction is determined, for which becomes the maximum, is:

This command determines the crack-growth angle using the local kink-tip stress-intensity-factor method in nonproportional fatigue problems:

CGROW,NPLOAD,DIRECTION,KLOC,ParW

Note that this option is only supported when the kink angle calculation is based on MTS criterion (CGROW,SOPT,KANG,MTS).

Fatigue crack growth under spectrum loading is a common occurrence, and structures operating under mission conditions are typically subjected to complex spectrum loading. It is crucial to simulate such scenarios accurately to predict fatigue life. To simplify spectrum loading, cycle-counting techniques (such as rainflow counting, range-pair counting, and racetrack counting) are often used.

This figure illustrates loading in one pass (primary cycle) of a complex spectrum and its simplification via some cycle-counting method:

Figure 3.10: Original Loading History (a) and Alternating Loading After Applying Cycle-Counting (b)

You can use the SMART method to analyze fatigue crack growth for a given spectrum loading. SMART requires simplified spectrum loading (shown by (b) in the figure) as input. SMART does not perform any conversion of a complex load-time history data to a simplified spectrum.

CLOAD,DEFINE,BEGIN
*DIM,_cycload1,TABLE,4 ,2
CLOAD,DEFINE,END

CLOAD,MSUB,2          ! 2 subcycles

! Time values
_cycload1(1,0) = 0.   ! time must start at 0
_cycload1(2,0) = t1
_cycload1(3,0) = 0    ! time must start at 0 for every subcycle
_cycload1(4,0) = t2
! Load values
_cycload1(1,1) = 0.   
_cycload1(2,1) = p1   ! first proportional load value
_cycload1(3,1) = 0
_cycload1(4,1) = p2   ! second proportional load value
! Number of cycles
_cycload1(1,2) = n1   ! number of cycles of the first load
_cycload1(3,2) = n2   ! number of cycles of the second load

*DIM,_tab_srat,TABLE,3,1,,TIME
 
! Time values
_tab_srat (1,0) = 0.   
_tab_srat (2,0) = t1
_tab_srat (3,0) = t1+t2 

! Stress ratio values
_tab_srat (1,1) = 0.   
_tab_srat (2,1) = R1  ! stress ratio for first subcycle
_tab_srat (3,1) = R2  ! stress ratio for second subcycle 

! CGROW command for stress ratio
CGROW,FCG,SRAT,%_tab_srat%

CLOAD,DEFINE,BEGIN
*DIM,_cycload1,TABLE,10 ,2
CLOAD,DEFINE,END

CLOAD,MSUB,4

! Time values
_cycload1(1,0) = 0.   ! time must start at 0
_cycload1(2,0) = t11
_cycload1(3,0) = t12
_cycload1(4,0) = 0.   ! time must start at 0 for every subcycle
_cycload1(5,0) = t21
_cycload1(6,0) = 0    ! time must start at 0 for every subcycle
_cycload1(7,0) = t31
_cycload1(8,0) = t32
_cycload1(9,0) = 0.   ! time must start at 0 for every subcycle
_cycload1(10,0) = t41


! Load values
_cycload1(1,1) = 0 
_cycload1(2,1) = p11   ! first subcycle load value   
_cycload1(3,1) = 0      
_cycload1(4,1) = 0     ! ramped load begins 
_cycload1(5,1) = p21   ! ramped load ends 
_cycload1(6,1) = p31   ! second subcycle load value, p31=p21
_cycload1(7,1) = p32   ! second subcycle load value
_cycload1(8,1) = p33   ! second subcycle load value, p33=p31
_cycload1(9,1) = p41   ! ramped load begins, p41=p33
_cycload1(10,1) = 0    ! ramped load ends       

! Number of cycles
_cycload1(1,2) = n1   ! number of cycles of the first load
_cycload1(4,2) = -1   ! ramped load, ignored in SMART
_cycload1(6,2) = n2   ! number of cycles of the second load
_cycload1(9,2) = -1   ! ramped load, ignored in SMART

 

3.3.1.7.2. Applying Spectrum Loading in a SMART Crack-Growth Analysis

In SMART spectrum loading, individual cycles are not repeated in the subcycles in order to reduce computational cost. Instead, the program calculates damage accumulation from cyclic fatigue loading in each subcycle using the Palmgren-Miner linear damage rule:

where:

= index of a crack-front node

= nodal damage from all of the subcycles of one primary cycle

= number of subcycles

= repeat number of an individual cycle for the ith subcycle

At the end of each subcycle, is calculated at each crack-front node (via a proportional or nonproportional solution) using the specified fatigue law. represents the number of cycles of the ith subcycle required to cause damage of unit crack extension:

After is evaluated via the Palmgren-Miner rule, the number of primary cycles (where one primary cycle is shown by (b) in the figure) required to propagate the crack by units is determined by:

The crack increment at the crack-front is determined for the cycle increment via the crack-growth law:

The default method for obtaining the crack kink angle at a node at the end of each primary cycle involves calculating a weighted average of the kink angles from each individual subcycle (CGROW,SOPT,SPLO,WEIG):

Another approach for specifying the kinking angle in crack propagation involves using the kink angle from the dominant step, where the maximum damage is created (CGROW,SOPT,SPLO,DMAX):

Similarly, the average stress ratio for a primary cycle is calculated as a weighted average of stress ratios of every individual subcycle:

The average stress-intensity-factor range at a crack-front node for each primary cycle can be obtained via a Newton-Raphson iteration of specified fatigue law function (Paris’ Law, Walker Equation, Forman Equation, NASGRO Equation). The average quantities and should satisfy the fatigue law as:

where are the fatigue law parameters.

The following output quantities are available (PRCINT or *GET) at the end of each primary cycle of spectrum loading:

Label	Quantity	Meaning
DLTK		Average stress-intensity-factor range
R		Average stress ratio
DLTN		Primary cycle increment number
DLTA		Crack-growth increment

Fracture parameter values (such as , , and ) are accessible at each solution substep (via PRCINT or *GET). The crack-growth solution (DLTK, R, DLTN, DLTA and KANG) for each primary cycle are available at the substep corresponding to the last time-point of that primary cycle.

When a crack begins to grow, SMART refines the mesh around the new crack front for better accuracy to accommodate the singularity of the crack tip. Without suitable coarsening of the mesh behind the crack front, the mesh around the crack front may be excessively fine.

Figure 3.11: Excessively Fine Mesh Around a Crack Front

In such cases, three mesh-coarsening options (specified via CGROW) are available:

When multiple cracks exist, the specified mesh-coarsening option applies to the last identified crack.

When a crack begins to grow, SMART uses the mesh around the crack front to determine the crack-increment size for the robustness and accuracy of the crack-growth analysis.

If your model has a relatively small crack (and a correspondingly fine mesh) at the start of the solution, SMART may require a very small crack increment. Using that initial crack-growth increment throughout the analysis, however, can result in an excessive number of crack-growth increments and solution steps. As the crack propagates during the solution process, a very small crack-growth increment may become unnecessary and a larger increment may be more desirable.

Conversely, suppose that a crack front approaches a critical regime (such as a boundary or interface), requiring a finer mesh around the crack front for the sake of parameter-calculation accuracy. In that case, a smaller crack-growth increment may again be necessary.

For such scenarios, two command options (specified via CGROW) are available for controlling the crack-increment size:

The control method can be a function of the solution time, accumulated maximum crack extension, or number of remeshings:

If the maximum multiplier value > 1, the crack increments increase; if < 1, the crack increments decrease. If the current variable < beginning variable of the function, Mechanical APDL uses 1; if the current variable > end variable of the function, the program uses the maximum multiplier value:

If the input is tabular data (CGROW,RMCONT,CMULT, TABLE,%tablename%), the table containing the data controls the parameter directly. You can specify one table per crack-growth set, and time is the only valid primary variable.

Accurate parameter calculation in a SMART crack-growth simulation requires a sufficiently fine mesh around the crack front.

Because a very fine mesh for the entire remeshing zone leads is computationally expensive, SMART remeshing uses a fine mesh around the crack front and a coarse mesh in the far field. This approach can, however, lead to accuracy loss in special zones such as boundaries, interfaces, voids, and holes. If your model has such areas where the mesh must be conserved or controlled, options (specified via CGROW,RMCONT,ESIZE) are available to control the element size in the remeshing zone:

Figure 3.12: Element-Size Control to Avoid Coarsening on the Boundary

Fracture mechanics deals with cracks (defects), and a singularity always exists around a crack front. The finite element solution therefore always depends on the element sizes, especially around the crack front. To ensure consistent results throughout the crack-growth process, maintain a reasonable crack-front element size and crack-growth increment. The crack-growth increments are usually of the order of the average element size at the crack front.

You can specify the maximum crack-growth increment (damax) as an absolute value or a multiple of the reference-element size (ESIZE) at the crack front: damax = m * ESIZE (Default: m = 1.5)

The minimum crack-growth (damin) is specified in a similar way. (Default: damin = 0). To explore options for specifying damin and damax, see CGROW,SOPT.

You can also control ESIZE at the crack front as the crack grows via CGROW,RMCONT,CMULT and CGROW,RMCONT,ESIZE.

The crack-growth calculation occurs in the solution phase after stress calculation. The fracture parameter is calculated first, followed by the crack extension according to the crack-growth method:

To specify the maximum crack-growth increment (damax):

To specify the minimum crack-extension increment (damin):

When the fracture criterion is reached for the specified crack-growth type, the program determines the crack-growth increment and generates the new crack surfaces.

For an analysis involving multiple cracks, you can set different and limits for individual cracks.

For a fatigue crack-growth analysis for multiple cracks, the program maintains the cycle increment at a substep for all cracks (CGROW,SOPT,SCN) by default. The program selects such that values at all cracks are within their respective limits. Crack-extension at a node is determined via:

For a static crack-growth analysis involving multiple cracks, the crack node with the maximum stress-intensity factor (among all cracks) is identified. Crack extension at that node, at crack j for example, is assigned the corresponding value (that is, ). Crack extensions at other nodes related to any crack are determined as:

To control how the program handles small crack extensions at any crack, use the options (CGROW,SOPT,DAMN).

Large Time Steps and Fracture Parameters  —  A crack-growth condition is based on whether the fracture criterion is met along the crack-front nodes; therefore, a large time step may result in significant over-prediction of the fracture parameters (and therefore the load-carrying capacity of the structures). A large time step can also cause over-prediction in the solution when crack-growth becomes unstable. In both cases, try using a small minimum time DTMIN.

Also see Example: Fatigue Crack-Growth Analysis Using SMART.

A critical component in the SMART framework is remeshing as a result of crack-growth. When a crack is growing, however, the crack profile can become very complex even for a simple initial crack geometry; therefore, be aware that remeshing is not always successful. Try the following if you encounter remeshing failure:

Following is the input file for the example fatigue crack-growth analysis using SMART:

/prep7

A = 46.6       	! crack length, mm
W = 100        	! width, mm
W1 = 125       	! total width, mm
H = 60         	! half height, mm
R = 12.5       	! radius of load circle, mm
E = 27.5       	! pin height, mm
S = 3          	! half width of notch, mm
D1 = 80        	! depth of notch, mm
D2 = 75        	! depth of notch, mm
B  = 12       	! thickness, mm

PF = 450 		! pin load, N	

k,1,A
k,2,W
k,3,W,H
k,4,,H
k,5,(W-W1),H
K,6,(W-W1),S
k,7,,S
k,8,(W-D1),S
k,9,(W-D2)
k,10,,E
k,11,,E,E

k,20,W,-H
k,21,,-H
K,22,(W-W1),-H
k,23,(W-W1),-S
k,24,,-S
k,25,(W-D1),-S
k,26,(W-D2)
k,27,,-E
k,28,,-E,E    

k,40,(W-D1),H
k,41,(W-D1),-H

circle,10,r,11,4,,8
l,1,2
l,2,3
l,3,40
l,40,4
l,4,5
l,5,6
l,6,7
l,7,8
l,8,9
l,9,1
l,4,12
l,16,7

circle,27,r,28,21,,8
l,20,41
l,41,21
l,21,22
l,22,23
l,23,24
l,24,25
l,25,26
l,26,1
l,2,20
l,24,33
l,29,21
l,8,40
l,25,41

al,9,10,11,40,17,18
al,41,35,36,9,37,29
al,12,19,8,7,6,5,20,16,40
al,34,38,24,23,22,21,39,30,41
al,38,33,32,31,39,28,27,26,25
al,13,19,1,2,3,4,20,15,14


et,1,187
type,1
mp,ex,1,200000 		! elastic modulus, MPa
mp,nuxy,1,0.33 		! Poisson's ratio

! Paris' Law Constants (units of delta-K in MPa.mm0.5, da/dN in mm/cycle) 
C=2.29E-10
M=2

! Fatigue crack-growth law specification
tb,cgcr,2,,,PARIS
tbdata,1,C,M

esize,4
vext,all,,,0,0,-B
allsel,all

vmesh,all

asel,s,,,8
asel,a,,,18
asel,a,,,13
arefine,all,,,,1

asel,s,,,23
asel,a,,,24
asel,a,,,25
asel,a,,,26
asel,a,,,50
asel,a,,,51
asel,a,,,52
asel,a,,,53
nsla,s,1

*get,numnode,node,0,count,,,
f,all,fy,PF/numnode
allsel,all

asel,s,,,33
asel,a,,,34
asel,a,,,35
asel,a,,,36
asel,a,,,39
asel,a,,,40
asel,a,,,41
asel,a,,,42
nsla,s,1

*get,numnode,node,0,count, , , 
f,all,fy,-PF/numnode
allsel,all

lsel,s,,,49
nsll,s,1
d,all,uy,0
allsel,all

lsel,s,,,74
lsel,a,,,90
nsll,s,1
d,all,ux,0
allsel,all

! plane strain, Uz constrained at front and back faces
nsel,s,loc,z,0
nsel,a,loc,z,-B
d,all,uz,0
allsel,all

! crack-front-nodes component
nsel,s,loc,x,a
nsel,r,loc,y,0
nlist
cm,crack1,node
allsel,all

! upper crack-surface-nodes component
asel,s,,,12,13
nsla,s,1
cm,crackt_sur_01,node
allsel,all

! lower crack-surface-nodes component
asel,s,,,18,19
nsla,s,1
cm,crackt_sur_02,node
allsel,all

finish

/com, *********************************************************
/com,             SOLUTION
/com, *********************************************************

/solu

antype,static
kbc,1 			! step load, fatigue

! fracture parameter, J-integral
cint,new,1
cint,type,JINT
cint,ctnc,crack1 
cint,edir,cs,0,x
cint,norm,0,2
cint,surf,crackt_sur_01,crackt_sur_02
cint,ncon,5

! crack-growth calculations
cgrow,new, 1
cgrow,cid, 1
cgrow,method,smart
cgrow,fcg,meth,LC         ! life-cycle method
cgrow,sopt,damx,0.5        ! maximum crack-growth increment, mm
cgrow,fcg,srat,0          ! stress ratio 
cgrow,fcoption,mtab,2     ! material table data (Paris law)

nsubst,6                  ! number of substeps
outres,all,all

solve
finish

/com, *********************************************************
/com,             POST-PROCESSING
/com, *********************************************************

/post1

  *get,nstep,active,0,set,nset
  crkid = 1
  
  maxnumnd = 0
  set,first
  *do,istep,1,nstep
  ! get number of crack front nodes
     *get,pval,cint,1,nnod
     *if,pval,gt,maxnumnd,then
        maxnumnd = pval
     *endif
  set,next
  *enddo
  
  *dim,numnd,array,nstep,1
  *dim,nodenum,array,maxnumnd,nstep
  *dim,DN,array,2,nstep
  *dim,DA,array,2,nstep
  *dim,DK,array,2,nstep
  *dim,SR,array,2,nstep
  *dim,ND_X,array,2,nstep
  *dim,ND_Y,array,2,nstep
  *dim,ND_Z,array,2,nstep
  *dim,ND_A,array,2,nstep
  
  ! DN_R and DK_R store the reference results
  *dim,DN_R,array,2,nstep
  *dim,DK_R,array,2,nstep
  
  ! out0 stores the results of the first node
  ! out1 stores the results of the last node
  *dim,out0_dn,array,nstep,2
  *dim,out0_da,array,nstep,2
  *dim,out0_dk,array,nstep,2
  *dim,out1_dn,array,nstep,2
  *dim,out1_da,array,nstep,2
  *dim,out1_dk,array,nstep,2
  
  set,first  
  ! get number of crack front nodes
  *get,pval,cint,1,nnod
  istep = 1
  numnd(istep,1) = pval
  !  nodes at the two ends of crack front
  *do,inode,1,2
     *if,inode,eq,1,then
        *get,ndnum,cint,crkid,node,inode
     *else
        *get,ndnum,cint,crkid,node,numnd(istep,1)
     *endif
 
     nodenum(inode,istep) = ndnum
  
     *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,dltn
     DN(inode,istep) = pval
     *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,dlta
     DA(inode,istep) = pval
     *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,dltk
     DK(inode,istep) = pval
     *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,R
     SR(inode,istep) = pval
  
     *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,crdx
     ND_X(inode,istep) = pval
     *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,crdy
     ND_Y(inode,istep) = pval
     *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,crdz
     ND_Z(inode,istep) = pval
     *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,apos
     ND_A(inode,istep) = pval
	 

	 ! Analytical solution (reference result)
     pi=3.1416
     a0 = 46.6 	! initial crack size, mm
     aw=a0/W
     z=16.7*aw**0.5-104.7*aw**1.5+369.9*aw**2.5
     zz=-573.8*aw**3.5+360.5*aw**4.5
     zzz=PF/B*sqrt(pi/W)
     deltk= zzz*(zz+z)
     DK_R(inode,istep) = deltk
     deltn = DA(inode,1)/(C*(deltk**M))
     DN_R(inode,istep) = deltn
	 
    *if,inode,eq,1,then
      out0_dn(istep,1)=DN(1,istep)
      out0_da(istep,1)=DA(1,istep)
      out0_dk(istep,1)=DK(1,istep)
      out0_dn(istep,2)=DN_R(1,istep)
      out0_dk(istep,2)=DK_R(1,istep)
    *else
      out1_dn(istep,1)=DN(inode,istep)
      out1_da(istep,1)=DA(inode,istep)
      out1_dk(istep,1)=DK(inode,istep)
      out1_dn(istep,2)=DN_R(inode,istep)
      out1_dk(istep,2)=DK_R(inode,istep)
    *endif
  *enddo
   
  set,next

  *do,istep,2,nstep
  
  !  get number of crack front nodes
	*get, pval, cint,1, nnod
     numnd(istep,1) = pval
  !  nodes at the two ends of crack front
	*do, inode, 1, 2
        *if,inode,eq,1,then
          *get,ndnum,cint,crkid,node,inode
         *else
          *get,ndnum,cint,crkid,node,numnd(istep,1)
         *endif
        
       *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,dltn
       DN(inode,istep) = DN(inode,istep-1) + pval
       *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,dlta
       DA(inode,istep) = DA(inode,istep-1) + pval
       *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,dltk
       DK(inode,istep) = pval
       *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,R
       SR(inode,istep) = pval
  
       *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,crdx
       ND_X(inode,istep) = pval
       *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,crdy
       ND_Y(inode,istep) = pval
       *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,crdz
       ND_Z(inode,istep) = pval
       *get,pval,cint,crkid,ctip,ndnum,contour,1,dtype,apos
       ND_A(inode,istep) = pval
	   
	   ! Analytical solution (reference result)
       da0 = DA(inode,istep-1)
       a0 = 46.6 + da0
       aw = a0/W
       z = 16.7*aw**0.5-104.7*aw**1.5+369.9*aw**2.5
       zz = -573.8*aw**3.5+360.5*aw**4.5
       zzz = PF/B*sqrt(pi/W)
       deltk = zzz*(zz+z)
       DK_R(inode,istep) = deltk
       nst = (DA(inode,istep)-da0)/(C*(deltk**M))
       deltn = DN_R(inode,istep-1) + nst
       DN_R(inode,istep) = deltn
	   
       *if,inode,eq,1,then
         out0_dn(istep,1)=DN(1,istep)
         out0_da(istep,1)=DA(1,istep)
         out0_dk(istep,1)=DK(1,istep)
         out0_dn(istep,2)=DN_R(1,istep)
         out0_dk(istep,2)=DK_R(1,istep)
       *else
         out1_dn(istep,1)=DN(inode,istep)
         out1_dk(istep,1)=DK(inode,istep)
         out1_da(istep,1)=DA(inode,istep)
         out1_dn(istep,2)=DN_R(inode,istep)
         out1_dk(istep,2)=DK_R(inode,istep)
       *endif
     *enddo
   
   set,next
  
  *enddo
 
  *dim,label,array,nstep,1
  *dim,DNtab,table,nstep,2
  *dim,DAtab,table,nstep,1
  *dim,DKtab,table,nstep,2

  *do,i,1,nstep
     label(i)=i
    *vfill,DNtab(i,1),data,DN(1,i)
    *vfill,DAtab(i,1),data,DA(1,i)
    *vfill,DKtab(i,1),data,DK(1,i)
    *vfill,DNtab(i,2),data,DN_R(1,i)
    *vfill,DKtab(i,2),data,DK_R(1,i)	
  *enddo   

  ! plot results (save as png) 

  /show,png
  /view,1,1,1,1
  /dscale,,0
  /graphics,power

  /axlab,x,crack extension (mm)
  /axlab,y,number of cycles
  /gcol,1,MAPDL
  /gcol,2,REF
  *vplot,datab,DNtab,2

  /axlab,x,crack extension (mm)
  /axlab,y, delta-K (MPa.mm**0.5)
  /yrange,20,40
  /gcol,1,MAPDL
  /gcol,2,REF
  *vplot,datab,DKtab,2

/out
/com, ----------------------solver results comparison------------------------
/com,
/com,                    DK OF THE FIRST NODE
/com,
/com,SUBSTEP |    DELTA-A    |   REFERENCE  |    MECHANICAL APDL
/com,
/com,  --------------------
/com,
/out
*vwrite,label(1),out0_da(1,1),out0_dk(1,2),out0_dk(1,1) 
(F3.0,'   ',F14.8,'  ',F14.5,'   ',F14.5) 
/com,
/com, ----------------------solver results comparison------------------------
/com,
/com,                    DN OF THE FIRST NODE
/com,
/com,SUBSTEP |    DELTA-A    |   REFERENCE  |    MECHANICAL APDL
/com,
/com,  --------------------
/com,
*vwrite,label(1),out0_da(1,1),out0_dn(1,2),out0_dn(1,1) 
(F3.0,'   ',F14.8,'  ',F14.5,'   ',F14.5) 
/com,
/com, ----------------------solver results comparison------------------------
/com,
/com,                    DK OF THE LAST NODE
/com,
/com,SUBSTEP |    DELTA-A    |   REFERENCE  |    MECHANICAL APDL
/com,
/com,  --------------------
/com,
*vwrite,label(1),out1_da(1,1),out1_dk(1,2),out1_dk(1,1) 
(F3.0,'   ',F14.8,'  ',F14.5,'   ',F14.5) 
/com,
/com, ----------------------solver results comparison------------------------
/com,
/com,                    DN OF THE LAST NODE
/com,
/com,SUBSTEP |    DELTA-A    |   REFERENCE  |    MECHANICAL APDL
/com,
/com,  --------------------
/com,
*vwrite,label(1),out1_da(1,1),out1_dn(1,2),out1_dn(1,1) 
(F3.0,'   ',F14.8,'  ',F14.5,'   ',F14.5) 
/com,
/out,scratch

parsav,all

finish

The following reference is cited in this example analysis: