Single Edge Crack

This example deals with the stability of a plate with an edge crack that is subjected to a tensile load. To analyze the stability of existing cracks, you can apply the principles of fracture mechanics.

A common parameter in fracture mechanics, the so-called stress intensity factor KI, provides a means to predict if a specific crack causes the plate to fracture. When this calculated value becomes equal to the critical fracture toughness of the material, KIc (a material property), then usually catastrophic fracture occurs.

Determining the stress intensity factor directly from the local state at the crack tip is problematic, since the stresses are singular there. Because of this, more indirect energy based methods are attractive. In this example, KI is computed using the J-integral and from the energy release rate.

In addition, the crack growth rate and number of cycles needed to propagate the crack a certain distance are computed.

Model Definition

A plate with a width w = 1.5 m and height h = 4.5 m has a single horizontal edge-crack at the middle of the left vertical edge. The length of the crack is varied from a = 0.5 m to a = 0.7 m, see Figure 1. An external load is pulling the plate such that the top and bottom edges experience tensile stress, σ, of 20 MPa.

Because of the symmetry, only half of the plate is modeled. There are three paths for computing the J-integral:

A circle around the crack tip with the radius being half the crack length.

A circle around the crack tip with the radius being 0.7 times the crack length.

The external boundaries, excluding the crack surface.

Figure 1: Plate geometry and loading (before applying symmetry conditions).

You apply a tensile load to the upper horizontal edge, while the lower horizontal edge is constrained in the y direction using a symmetry condition. One point is constrained in the horizontal direction in order to suppress rigid body motion.

Material Model

The material properties are representative for steel.

Table 1: Material Data.
Quantity	Name	Expression
Young's modulus	E	206·109 Pa
Poisson's ratio	ν	0.3
Coefficient in Paris’ law	A	1.4·10-11 (KI unit system: MN/m3/2)
Exponent in Paris’ law	m	3.1

The J-Integral

In this model, you determine the stress intensity factor KI using the so-called J-integral.

The J-integral is a two-dimensional path independent line integral along a counterclockwise contour, Γ, surrounding the crack tip. For a crack extending along the positive x-axis, the J-integral is defined as

where W is the strain energy density

and T is the traction vector defined as

σij denotes the stress components, εij the strain components, and ni the normal vector components.

The J-integral has the following relation to the Mode I stress intensity factor for a plane stress case and a linear elastic material:

(1)

where E is Young’s modulus.

Energy Release Rate

For a linear elastic material it is actually possible to compute the value of the J-integral without using the path integrals. The reason is that its value equals the value of the energy release rate, G,

(2)

Here U is the potential energy, a is the crack length, and t is the thickness. By computing the potential energy for two slightly different crack lengths, G can be estimated as

(3)

The potential energy of an elastic body is

The first term is the strain energy in the volume, and the second term is the potential of the prescribed tractions on the boundary. Because of the linearity,

Thus, it is possible to compute the potential energy using either of these terms independently.

The total strain energy density exists as a built-in variable, making the first expression attractive for determining G.

Crack Propagation

When subjected to a periodic load, the crack growth rate (in meters per load cycle) can be expressed by Paris’ law:

(4)

Here A and m are material parameters and ΔKI is the range of the stress intensity factor. It is assumed that the load varies between 0 and 20 MPa, so that ΔKI equals the computed KI. Note that the value of the coefficient A depends on the unit for the stress intensity factor in a rather complex way because of the power m, which in general is non-integer.

Results

Figure 2 shows the stress singularity at the crack tip.

Figure 2: von Mises stress and the deformed shape of the plate when the crack length is 0.6 m. The displacement is exaggerated to illustrate the deformation under the applied load.

Based on Ref. 1 an analytical solution for the stress intensity factor is

where σ = 20 MPa, and f is a correction factor given in Ref. 1.

(5)

The expression in Equation 5 assumes that the plate is long, so that the height is significantly larger than the width. The stress field in Figure 2 indicates that this assumption is fulfilled.

The calculated stress intensity factors for the three different contours, tabulated in the evaluation group , show excellent agreement with the reference value. The largest deviation for any of the studied crack lengths and integration paths is less than 0.3%.

When using the built-in J-integral computation, you should ascertain that none of the contours used for the integration passes outside the computational domain, or encloses another crack. In the default plot , all integration paths are visualized. In Figure 3 and Figure 4 the contours are shown for the shortest and longest cracks, respectively.

Figure 3: The contours used for integration when a = 0.5 m.

Figure 4: The contours used for integration when a = 0.7 m.

The three different ways of computing the energy release rate, and thus KI, are compared in Figure 5. As can be seen, all three methods give essentially the same values.

Figure 5: J-integral compared with energy release rates computed using numerical differentiation.

Finally, the crack growth speed can be investigated. In Figure 6, the crack growth speed is shown as function of the crack length. The dependence is quite strong: an increase in crack length from 0.5 m to 0.7 m (40%) increases the crack growth rate by a factor of 5. According to the constants used in Paris’ law, the crack growth rate is proportional to the stress intensity factor raised to the power of 3.1. As can be seen from the previous results, the stress intensity factor increases strongly with the crack length, and this combination results in the increase in crack growth rate.

In practice, Paris’ law may not be applicable when KI approaches the critical value KIc.

Figure 6: Crack propagation rate as function of the crack length.

Notes about the COMSOL Implementation

In this analysis you compute the J-integral for three different contours traversing three different regions around the crack tip. Note that the boundaries along the crack are not included in the J-integral because they do not give any contribution to the J-integral. This is due to the following facts: for an ideal crack, dy is zero along the crack faces, and all traction components are also zero (Ti = 0) as the crack faces are not loaded.

To evaluate the J-integral along three different paths, you need to add one node with three subnodes. Since the option is used in the node, the values of the J-integrals and stress intensity factors automatically take the full structure into account.

When you add a node, the KI, KII, and (in 3D) KIII stress intensity factors are also computed.

When computing the energy release rates, the derivative of the potential energy is computed using a difference approximation. In order to access different solutions in a single expression, the withsol() operator is used.

Reference

1. H Tada, P.C. Paris, and G.R. Irwin, The Stress Analysis of Cracks Handbook, Del Research Corp, Press, 1973.

Structural_Mechanics_Module/Fracture_Mechanics/single_edge_crack

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select .

Click .

Click

In the tree, select .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file single_edge_crack_parameters.txt.

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type Wp.

In the text field, type Hp.

Click

Add a point at the crack tip.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type Xa.

Form Union (fin)

In the window, click .

In the window for , click

Materials

Steel

In the window, under right-click and choose .

In the window for , type Steel in the text field.

Locate the section. In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	2.06e11[Pa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.3	1	Young’s modulus and Poisson’s ratio
Density	rho	7850	kg/m³	Basic

Solid Mechanics (solid)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. In the d text field, type Th.

Symmetry 1

In the toolbar, click

and choose .

Select Boundaries 2 and 4 only.

It does not matter whether you select boundary 2 as a symmetry boundary or not. If selected, it will be overridden when the node is added.

Boundary Load 1

In the toolbar, click

and choose .

Select Boundary 3 only.

In the window for , locate the section.

Specify the FA vector as


0	x
20[MPa]	y

Prescribed Displacement 1

In the toolbar, click

and choose .

Suppress rigid body motion.

Select Point 4 only.

In the window for , locate the section.

Select the check box.

Crack 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Boundary 2 only.

Click to expand the section. Click

Select Point 3 only.

J-Integral 1

In the toolbar, click

and choose .

J-Integral 2

Right-click and choose .

In the window for , locate the section.

In the rΓ text field, type solid.crack1.crackSize*0.7.

Crack 1

In the window, click .

J-Integral 3

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Point 3 only.

Locate the section. Click to select the

toggle button.

Select Boundaries 1, 3, and 5 only.

Definitions

Loaded edge integration

In the toolbar, click

and choose .

In the window for , type Loaded edge integration in the text field.

In the text field, type LoadEdgeInt.

Locate the section. From the list, choose .

Select Boundary 3 only.

Use a fine mesh close to the crack tip where the stress gradients are large. The node will automatically provide that, but you may want to adjust the settings manually.

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. From the list, choose .

Size 1

In the window, under click .

In the window for , locate the section.

Click the button.

Locate the section.

Select the check box. In the associated text field, type 0.01.

Click

Set up a parametric sweep over the crack length.

Study 1

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
Xa (Crack length)	range(0.5,da,0.7)	m

In the toolbar, click

Results

Mirror 2D 1

In the toolbar, click

and choose .

In the window for , locate the section.

In row , set to 1.

In row , set to 0.

Locate the section. From the list, choose .

Stress (solid)

In the window, under click .

In the window for , locate the section.

From the list, choose .

Locate the section. Clear the check box.

Surface 1

In the window, expand the node, then click .

In the window for , locate the section.

From the list, choose .

Click to expand the section. Select the check box.

In the text field, type 0.

In the text field, type 200.

Stress (solid)

Right-click and choose .

Boundary Load

In the window for , type Boundary Load in the text field.

Click in the upper-right corner of the section. From the menu, choose .

Click to expand the section. From the list, choose .

Locate the section. From the list, choose .

Locate the section.

Select the check box. In the associated text field, type 1e-8.

Click to expand the section. From the list, choose .

Clear the check box.

Color Expression 1

Right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose .

Locate the section. From the list, choose .

From the list, choose .

Clear the check box.

Deformation 1

Right-click and choose .

Line 1

In the toolbar, click

In the window for , locate the section.

In the text field, type 1.

Click to expand the section. From the list, choose .

Locate the section. From the list, choose .

From the list, choose .

Click to expand the section. From the list, choose .

Deformation 1

In the toolbar, click

Line 1

In the window, click .

Click to expand the section. From the list, choose .

Clear the check box.

Deformation

In the window, expand the node, then click .

In the window for , locate the section.

Select the check box. In the associated text field, type 200.

In the toolbar, click

Click the

button in the toolbar.

Check the placement of the integration contours. In particular, it must be ensured that the circular paths are inside the domain.

In the toolbar, click

Add Predefined Plot

Go to the window.

In the tree, select .

Click in the window toolbar.

Results

Cracks (solid)

In the toolbar, click

In the window for , click

Check if there are any differences between the results for the three contours by adding predefined evaluation group. That could indicate a too coarse mesh.

Add Predefined Plot

Go to the window.

In the tree, select .

Click in the window toolbar.

In the toolbar, click

Since this is a pure Mode I case, remove the evaluation of the Mode II stress intensity factors.

Results

Stress Intensity Factors, Mode II

In the window, expand the node.

Right-click and choose .

Stress Intensity Factors, Mode I

Add a comparison between the values from the second integration path and the reference stress intensity factor.

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Expression	Unit	Description
K1r	1	Reference stress intensity factor
(solid.crack1.jint2.KI-K1r)/K1r*100	1	Percent difference from reference value

In the toolbar, click

Table

Go to the window.

Compare the J-integral with an energy release rate based on numerical differentiation of the strain energy with respect to the crack length.

Definitions

Variables 1

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
PE	-LoadEdgeInt(solid.bndl1.F_Ayv)Th	N·m	Potential energy, computed from load potential
dadN	AP*(solid.crack1.jint1.KI/1e6)^mP		Crack growth rate (m/cycle)

After having added a new variable, you must update the solution to make it accessible for result presentation.

Study 1

In the toolbar, click

Results

J-Integral and G

In the toolbar, click

and choose .

In the window for , type J-Integral and G in the text field.

Locate the section. From the list, choose .

From the list, choose .

In the text field, type range(2,20).

Global 1

Right-click and choose .

In the window for , click in the upper-right corner of the section. From the menu, choose .

Locate the section. In the table, enter the following settings:


Expression	Unit	Description
solid.crack1.jint1.J	J/m^2	J-integral, contour 1
((withsol('sol2',2solid.Ws_tot,setval(Xa,Xa+da))-withsol('sol2',2solid.Ws_tot,setval(Xa,Xa-da))))/(2daTh)	J/m^2	Energy release rate, strain energy
-((withsol('sol2',PE,setval(Xa,Xa+da))-withsol('sol2',PE,setval(Xa,Xa-da))))/(2daTh)	J/m^2	Energy release rate, load potential