COMSOL 6.4 - 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 is the so-called stress intensity factor KI, which provides a means to predict if a specific crack may lead to a fracture. When KI 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. Therefore, 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 GPa
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 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. The potential energy of an elastic body is

The first term is the strain energy in the volume, and the second and third terms are the potential of the body load and prescribed tractions, respectively.

Crack Propagation

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

(3)

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

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.

(4)

The expression in Equation 4 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 two different ways of computing the energy release rate, and thus KI, are compared in Figure 5. As can be seen, these two 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 an increasing 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.

To evaluate the energy release rate, you need to add subnode under the node. The subnode has an action button to set up the deformed geometry and stationary study with sensitivity node. Since the option is used in the node, the values of the energy release rate automatically take the full structure into account.

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

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.

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
dadN	AP*(solid.crack1.jint1.KI/1e6)^mP		Crack growth rate (m/cycle)

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.

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, right-click and choose .

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	206[GPa]	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