COMSOL 6.4 - Cracking of a Notched Beam

Cracking is an important failure mechanism in many materials, especially in brittle and quasi-brittle materials such as concrete. To model the failure of such materials requires not only an appropriate material model, but also modeling techniques to properly describe the size effect in order to avoid mesh-dependent results.

In this example, failure is examined by simulating a three-point bending fracture test of a notched concrete beam. The model is based on the experiments presented in Ref. 1. In that paper, several tests are presented on specimens of different sizes and configurations (different notch sizes or no notch) to study the size effect. In this example, one of those specimens is modeled, although the model is parameterized such that all specimens from the reference can be reproduced. The results are compared to measured data.

Model Definition

The geometry of a general specimen from Ref. 1 is given in Figure 1. The parameter Dn is defined by

(1)

where D0 = 400 mm and n = (1, 2, 3, 4).

Furthermore, the notch is described by its depth αm = (0.5, 0.2, 0). In this example, the specimen given by n = 2 and αm = 0.2 is modeled. Notice that the width of the notch is set to one mesh element. The depth of the beam is constant and equal to 50 mm. Plane stress conditions are used. The beam is simply supported and loaded with a ‘distributed’ point load above the notch. As in the experimental setup described in Ref. 1, the load is controlled through the crack mouth opening displacement (CMOD), which is monotonically increased for all tested specimens. Although the deflection is also monotonically increased for n = 2 and αm = 0.2, this might not be the case for larger specimens.

Figure 1: Geometry of the specimens, as specified in Ref. 1.

According to Ref. 1, measurements on the concrete used gave a Young’s modulus of 37 GPa and a Poisson’s ratio of 0.21. The measured compressive strength is 42.3 MPa and the tensile strength is 3.9 MPa. The material is described by a linear elastic material model with damage to account for tensile cracking. The damage model is set up to describe tensile material failure by using the Rankine equivalent strain definition. A damage evolution law with exponential strain softening is used. With this option, the fracture energy enters as an additional input parameter, which is here set to 85 J/m2. When modeling tensile cracking it is also necessary to add regularization to the damage model in order to ensure that a consistent amount of energy is dissipated during mesh refinement or for different discretization orders. Two techniques for this are available in COMSOL Multiphysics and both are evaluated in this model:

•

The crack band method;

•

The implicit gradient method.

The crack band method considers the current discretization and modifies the damage model locally at each material point based on the element size. A more refined approach is to use the implicit gradient method which enforces a predefined width of the damage zone through a localization limiter. This is done by adding a nonlocal strain variable and an internal length scale to the material model through an additional PDE that is solved simultaneously with the displacements. While the recommendation to use a linear displacement field when using the crack band method is used, the study with the implicit gradient method is set up to use the default discretization settings (that is, quadratic serendipity). For simplicity, the same mesh size is used for both cases. This can be considered a bit unfair when comparing results; when using first-order discretization a significantly finer mesh would be needed to obtain a similar accuracy in both methods.

Results and Discussion

Two studies are defined to compare the two regularization techniques:

•

The crack band method with a linear displacement field;

•

The implicit gradient method.

The results from both studies are compared in Figure 2 showing the force versus CMOD curves as well as the measured data from Ref. 1. As an additional comparison, the force versus deflection curves are shown in Figure 3, despite no measured data being available. The crack band method in Study 1 is in very good agreement with the measured data at the peak load, while for the post-peak behavior, the simulated force is underestimated for a given CMOD. For the implicit gradient method in Study 2, the peak load is overestimated by approximately 1 kN. Also in this case, the force is underestimated for the post-peak response. However, both solutions can describe the overall failure behavior of the beam.

Figure 2: Load versus crack mouth opening displacement.

Figure 3: Load versus vertical deflection.

The results of Study 1 can be analyzed in more detail in Figure 4 and Figure 5, which show the damage and strain fields, respectively, at the final step of the solution. As expected for the crack band method, both the damage and the strains are localized into a single row of elements. It can also be noticed that the crack almost propagates through the entire height of the beam.

Figure 4: Damage distribution when using the crack band model.

Figure 5: Equivalent strain distribution when using the crack band model.

The same plots are shown for the implicit gradient method from Study 2 in Figure 6 and Figure 7. For this modeling technique, a distributed zone of damage over several elements is obtained. Also, the strains are much less localized as compared to the crack band method. In Figure 6, it can be noticed that the damage field stretches below the tip of the notch. This is an unphysical behavior that results from the nonlocal modeling approach. This could explain the overestimation of the peak load seen in Figure 2 and is a general issue of nonlocal regularization methods for problems where a crack initiates at nontrivial physical boundaries, such as a notch. Finally, the stress versus strain response of the models is compared in Figure 8 at four points located at different heights above the notch.

Figure 6: Damage distribution when using the implicit gradient model.

Figure 7: Equivalent strain distribution when using the implicit gradient model.

Figure 8: Horizontal stress versus horizontal strain at four points above the notch.

Notes about the COMSOL implementation

Apart from being prone to tensile cracking, concrete is also highly nonlinear in compression. Accounting for this can improve the force versus CMOD curves. This can for example be done by using the Mazars damage model for concrete available with the Geomechanics Module.

The loading of the beam is controlled through the CMOD. A suitable modeling technique to accomplish this is to use an algebraic equation that controls the applied force so that the model reaches the desired CMOD. This is implemented using a , and the parametric solver steps up the desired CMOD.

Although a single specimen is analyzed, the model is parameterized so that all set-ups presented in Ref. 1 can be described. For example, the size effect can be examined by running a parametric sweep over n. However, if a specimen without a notch is to be modeled by setting i = 0, some minor modifications to the geometry and load control are necessary:

•

Remove the part of the geometry that creates the notch.

•

Redefine the position of the CMOD coupling operators so that the measuring points are further from the midsection.

For some configurations, the nonlinear solver also has to be modified to obtain a converging solution. In particular, this is the case for the specimens without a notch. For example, better convergence is often obtained if in the node the is set to with on every iteration.

Reference

1. D. Grégoire, L.B. Rojas-Solano, and G. Pijaudier-Cabot, “Failure and size effect for notched and unnotched concrete beams,” Int. J. Numer. Anal. Meth. Geomech., vol. 37, no. 10, pp. 1434–1452, 2013.

Geomechanics_Module/Damage/notched_beam_damage

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Right-click and choose .

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

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 3.5*Dn.

In the text field, type Dn.

Locate the section. From the list, choose .

In the text field, type Dn/2.

Click to expand the section. Add a number of subdivisions to the beam in order to facilitate application of loads and constraints, as well as mesh control.

Select the checkbox.

Clear the checkbox.

In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	0.5*Dn-loadSurf/2
Layer 2	loadSurf
Layer 3	0.75*Dn-loadSurf/2
Layer 4	0.5*Dn-loadSurf/2
Layer 5	loadSurf
Layer 6	0.5*Dn-loadSurf/2
Layer 7	0.75*Dn-loadSurf/2
Layer 8	loadSurf

Parameter Check 1 (pch1)

In the toolbar, click

and choose .

In the window for , locate the section.

In the text field, type am<=0.

In the text field, type Notch size must be larger than zero.

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type eSize.

In the text field, type am.

Locate the section. From the list, choose .

In the text field, type am/2.

Difference 1 (dif1)

In the toolbar, click

and choose .

Click in the window and then press Ctrl+A to select both objects.

Select the object only.

In the window, click .

In the window for , locate the section.

Click to select the

toggle button for .

Select the object only.

Click

Polygon 1 (pol1)

In the toolbar, click

Add a help line to simplify mapped meshing.

In the window for , locate the section.

From the list, choose .

In the text field, type -Dn*1.75 Dn*1.75.

In the text field, type am am.

Click

Most internal edges and vertices are only needed to facilitate a mapped mesh. Use to only make them visible during meshing.

Mesh Control Edges 1 (mce1)

In the toolbar, click

and choose .

On the object , select Boundaries 4, 6, 8, 9, 11, 13, 14, 16, 18, 19, 21, 23, 24, 27, 30, 31, 33, 34, 36, 38, 39, 41, 43, 44, 46, 48, and 49 only.

In the window for , locate the section.

Clear the checkbox.

Mesh Control Vertices 1 (mcv1)

In the toolbar, click

and choose .

On the object , select Points 2, 5, 6, 8–14, 20, 21, 23–25, 27, 28, 30, 31, and 33 only.

In the toolbar, click

Definitions

Create a variable measuring the crack mouth opening displacement (CMOD).

Crack Mouth Opening Displacement, Right

In the toolbar, click

and choose .

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

In the text field, type Crack Mouth Opening Displacement, Right.

Locate the section. From the list, choose .

Select Point 8 only.

Crack Mouth Opening Displacement, Left

In the toolbar, click

and choose .

In the window for , type Crack Mouth Opening Displacement, Left in the text field.

In the text field, type CMOD_left.

Locate the section. From the list, choose .

Select Point 6 only.

Variables 1

In the window, right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
CMOD	CMOD_right(u)-CMOD_left(u)	m	Crack mouth opening displacement

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

Linear Elastic Material 1

In the window, under > click .

Damage: Crack Band

In the toolbar, click

and choose .

In the window for , type Damage: Crack Band in the text field.

Damage: Implicit Gradient

Right-click and choose .

In the window for , type Damage: Implicit Gradient in the text field.

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

In the lint text field, type lscale.

In the hdmg text field, type 3*lscale.

Rigid Connector 1

In the toolbar, click

and choose .

Select Boundary 4 only.

In the window for , locate the section.

Select the checkbox.

Rigid Connector 2

In the toolbar, click

and choose .

Select Boundary 12 only.

In the window for , locate the section.

Select the checkbox.

Boundary Load 1

In the toolbar, click

and choose .

Select Boundary 6 only.

In the window for , locate the section.

From the list, choose .

Specify the Ftot vector as


0	x
-load	y

The applied load is controlled through the CMOD using a .

Click the

button in the toolbar.

In the dialog, in the tree, select the checkbox for the node > .

In the tree, select the checkbox for the node > .

Click .

Load Control

In the toolbar, click

and choose .

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

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


Name	f(u,ut,utt,t) (1)	Initial value (u_0) (1)	Initial value (ut_0) (1/s)	Description
load	CMOD-para*maxCMOD	0	0

Locate the section. Click

In the dialog, type force in the text field.

In the tree, select > .

Click .

In the window for , locate the section.

Click

In the dialog, type displ in the text field.

In the tree, select > .

Click .

Discretization, Linear

In the toolbar, click

and choose .

For the crack band method, linear shape order for the displacements is preferred.

In the window for , type Discretization, Linear in the text field.

Locate the section. From the list, choose .

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	37[GPa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.21	1	Young’s modulus and Poisson’s ratio
Density	rho	2400	kg/m³	Basic
Peak strength	sigmap	3.9[MPa]	N/m²	Scalar damage
Fracture energy per area	Gf	85	J/m²	Scalar damage