COMSOL 6.4 - Brittle Fracture of a Holed Plate

This example studies the brittle fracture of a holed plate made of a cement mortar. The setup of the model, including dimensions and material properties, is based on the experimental data reported in Ref. 1. As the plate is loaded, a mixed-mode fracture is induced with a crack propagating from the predefined notch to the unsymmetrically placed hole in the center of the plate.

Fracture is modeled using a damage model that regularizes the sharp geometry of the crack by the phase-field approximation. This means that the crack is described in the domain material, so the nonlocal phase field makes the crack path independent of the mesh elements. The example shows how to define an efficient and stable solver configuration for the phase-field damage method, which is often unstable for this class of problems.

Model Definition

The geometry of the specimen (Ref. 1) is shown in Figure 1. The overall dimensions of the plate are equal to 65 mm in width and 120 mm in height. The thickness of the plate is constant and equal to 1 mm. A notch is placed on the left boundary, 65 mm from the bottom of the plate, in order to control where the crack initiates during loading. A mixed-mode fracture is induced by offsetting the large hole and the notch from the center of the plate. Loading is applied through displacement-controlled metal pins inserted into the two smaller holes. The model assumes a plane stress condition.

Figure 1: Geometry of the holed plate model.

The plate is made of a cement mortar composed of 22% cement, 66% sand with a grain size smaller than 1 mm, and 12% water. Material properties (from Ref. 1) are given in Table 1. In order to account for tensile cracking, the mortar is described as a linear elastic material model with damage. A phase-field approximation is made after the sharp crack geometry, and cracking is incorporated into the domain material. The setup of the phase-field damage model does not include a damage threshold, meaning that all material points subjected to tension are damaged. The only material input to the damage model is the critical energy release rate Gc (sometimes called fracture energy), whereas the tensile strength is determined by the damage evolution function and the evolution of the phase field. The latter can be modified by the internal length scale lint that controls the width of the localized phase field. In this example, the parameter lint is set to 0.25 mm.

Table 1: Material properties of the cement mortar.
Material property	Value
Young’s modulus	6 GPa
Poisson’s ratio	0.22
Critical energy release rate	2280 J/m2

To properly resolve the phase field and to achieve a stable material behavior, a high mesh density is required in the vicinity of the propagating crack. Since the expected crack trajectory is known in advance, the mesh is locally refined with a maximum mesh element size equal to lint. This length can be considered as being close to the upper limit of the appropriate mesh size. If extra computational cost is acceptable, a finer mesh would resolve the phase field better and also yield a more stable material behavior.

The two metal pins in the model are represented by two rigid connectors attached to the boundaries of the two smaller holes. Loading is applied by prescribed displacements to the two rigid connectors while they are free to rotate. No other loads nor constraints are considered.

Results and Discussion

Figure 2 shows the reaction force in the upper rigid connector versus its prescribed displacement. A first peak in the curve can be identified at a load of 0.63 kN and a lateral displacement of 0.33 mm. As shown in Figure 3, this peak corresponds to the formation of a crack at the tip of the notch. Actually, the crack starts propagating slightly before the observed peak, see the first green circle in Figure 2 and the leftmost plot in Figure 3. Once in the post-peak regime of the curve, the crack propagates rapidly and the load capacity of the plate reduces significantly. As the crack approaches the large hole, it diverts from the initially straight path; a new stable configuration is obtained once it reaches the hole. At this point the load can be increased, but with a much lower stiffness than prior to the formation of the crack.

A second peak is reached at a load of 0.15 kN and 1.7 mm displacement. This second crack propagates toward the right boundary of the plate, and it would eventually split the specimen into two pieces. However, due to local compressive stress fields and residual stiffness in the damaged mesh elements in the model, the crack propagation is slowed down before reaching the exterior boundary. By increasing the prescribed displacement of the rigid connector, the crack would eventually reach the exterior also in the simulation when the bending action of the plate changes to pure tensile loading.

The distribution of the phase field at the last step of the simulation is shown in Figure 4. The phase field has a smooth variation, which indicates that the mesh resolution is sufficiently fine. For the second crack, the phase field is slightly wider close to the hole, since it is localized without the aid of a geometric notch. Localization of strains is an unstable event, which explains the small bleeding of the phase field observed in the second crack.

Figure 2: Load versus displacement curve with green circles indicating the snapshots presented in Figure 3.

Figure 3: Crack trajectory and maximum principal stress field for different parameter steps. The red arrows indicate the relative size of the reaction force in the rigid connectors.

Figure 4: Crack phase field at the last parameter step.

Notes About the COMSOL Implementation

The loading of the specimen is done by prescribing the displacement of each . In this way, the loading is displacement-controlled, which is necessary to track the peak loads in the global response. The simulation would have stopped before reaching the first peak if a load-controlled setup had been used.

While it is possible to solve this type of brittle fracture problem using a fully coupled strategy, the phase-field damage model can often exhibit poor or slow convergence. Therefore, in this example it is shown how to use a segregated strategy to improve the convergence and stability of the numerical solution by splitting the evolution of the crack phase field and the displacement field into two groups. This type of algorithmic operator split for phase-field fracture models was originally suggested in Ref. 2, and it can conceptually be summarized as follows for step n + 1:

Initialization. At step n, the crack phase field, the displacement field, and other state variables are known.

Update state variables. Update internal state variables with the values from step n.

Solve for the Crack Phase Field. Compute the crack phase-field variable in a Newton step, with the displacement field frozen at step n.

Solve for the Displacement field. Compute the displacement field variable in a Newton step with the updated crack phase field.

These steps lead to a single-pass algorithm that corresponds to the default solver sequence that is, however, accurate only for sufficiently small parameter steps. An improvement to this scheme is to add a multipass correction, by iterating in each increment over steps 3 and 4 until either a tolerance-based convergence criterion is fulfilled or for a predefined number of iterations. The latter strategy is used in this example by setting the to 3 in the subnode under the .

Notice that the suggested solver configuration does not require convergence of the outer problem, and the solution is always accepted after 3 segregated iterations. However, each subgroup locally fulfills the defined convergence criterion. Hence, the displacement field can be considered as a converged solution given the current crack phase field.

The accuracy of the scheme can be improved if the extra computational cost of increasing the number of iterations (or the required convergence of the outer problem) is acceptable.

References

1. M. Ambati, T. Gerasimov, and L. De Lorenzis, “A review on phase-field models of brittle fracture and a new fast hybrid formulation,” Comput. Mech., vol. 55, pp. 383–405, 2015.

2. C. Miehe, M. Hofhacker, and F. Welschinger, “A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits,” Comput. Methods Appl. Mech. Eng., vol. 199, pp. 2765–2778, 2010.

Geomechanics_Module/Damage/holed_plate_fracture

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.

In the table, enter the following settings:


Name	Expression	Value	Description
height	120[mm]	0.12 m	Plate height
width	65[mm]	0.065 m	Plate width
notchHeight	0.5[mm]	5E-4 m	Notch height
notchWidth	10[mm]	0.01 m	Notch width
notchLocation	65[mm]	0.065 m	Notch location
holeRadius	10[mm]	0.01 m	Hole radius
holeX	36.5[mm]	0.0365 m	Hole center, x-coordinate
holeY	51[mm]	0.051 m	Hole center, y-coordinate
elemSize	0.25[mm]	2.5E-4 m	Mesh element size
para	0	0	Load parameter

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type width.

In the text field, type height.

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type holeRadius.

Locate the section. In the text field, type holeX.

In the text field, type holeY.

Circle 2 (c2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 5[mm].

Locate the section. In the text field, type 20[mm].

In the text field, type 20[mm].

Circle 3 (c3)

In the toolbar, click

In the window for , locate the section.

In the text field, type 5[mm].

Locate the section. In the text field, type 20[mm].

In the text field, type height-20[mm].

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type notchWidth.

In the text field, type notchHeight.

Locate the section. In the text field, type notchLocation-notchHeight/2.

Difference 1 (dif1)

In the toolbar, click

and choose .

Drag and drop below .

Select the object only.

In the window for , locate the section.

Click to select the

toggle button for .

Select the objects , , , and only.

Partition the plate to facilitate a finer mesh around the expected crack path.

Line Segment 1 (ls1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type width.

In the text field, type holeY+elemSize*12.

Locate the section. From the list, choose .

In the text field, type holeX.

In the text field, type holeY+elemSize*5.

Locate the section. Find the subsection. Click .

In the dialog, type Mesh control edges in the text field.

Click .

Line Segment 2 (ls2)

Right-click and choose .

In the window for , locate the section.

In the text field, type holeY-elemSize*12.

Locate the section. In the text field, type holeY-elemSize*5.

Polygon 1 (pol1)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

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


x (m)	y (m)
holeX	holeY
holeX	notchLocation+elemSize*12
notchWidth*3/4	notchLocation+elemSize*12
notchWidth*3/4	notchLocation+notchHeight/2

Locate the section. Find the subsection. From the list, choose .

Polygon 2 (pol2)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

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


x (m)	y (m)
notchWidth*3/4	notchLocation-notchHeight/2
notchWidth*3/4	notchLocation-elemSize*5
holeX+holeRadius	holeY

Locate the section. Find the subsection. From the list, choose .

Rectangle 3 (r3)

In the toolbar, click

In the window for , locate the section.

In the text field, type notchWidth.

In the text field, type notchHeight.

Locate the section. In the text field, type notchWidth.

In the text field, type notchLocation-notchHeight/2.

Locate the section. Find the subsection. From the list, choose .

Partition Objects 1 (par1)

In the toolbar, click

and choose .

Select the object only.

In the window for , locate the section.

Click to select the

toggle button for .

Select the objects , , , , and only.

Click

Mesh Control Edges 1 (mce1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Click the

button in the toolbar.

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 1[mm].

Linear Elastic Material 1

In the window, under > click .

Damage 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the lint text field, type 0.25[mm].

Rigid Connector 1

In the toolbar, click

and choose .

Select Boundaries 9, 10, 13, and 14 only.

In the window for , locate the section.

Select the checkbox.

Click to expand the section. Select the checkbox.

Rigid Connector 2

Right-click and choose .

In the window for , locate the section.

Click

Select Boundaries 11, 12, 15, and 16 only.

Locate the section. In the u0y text field, type para*2[mm].

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	6[GPa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.22	1	Young’s modulus and Poisson’s ratio
Density	rho	2000	kg/m³	Basic
Critical energy release rate	Gc	2280[J/m^2]	J/m²	Phase-field damage