Brittle Damage in Uniaxial Tension

Modeling crack formation in quasi-brittle materials such as concrete is associated with the phenomenon of strain localization due to material softening. In a finite element model, this can cause the solution to be mesh dependent, which is an undesirable property. This tutorial model shows how to avoid this problem by using two different regularization techniques available in COMSOL Multiphysics.

The example describes the axial stretching of a bar, where a damage model is used to account for tensile cracking.

Model Definition

The geometry of the model consists of a 10 cm long bar with height and thickness equal to 2 mm. The analysis is done under plane stress conditions. Due to the symmetry, only half the height is modeled.

The bar is considered to be made of a brittle material with the following properties.

•

Young's modulus is 30 GPa.

•

Poisson's ratio is 0.2.

•

The tensile strength is 2 MPa.

•

The fracture energy is 60 J/m2. This is the energy dissipated during the creation of a single crack. The cracking process is modeled using an isotropic damage model with a single damage variable that only considers the tensile failure of the material.

One of the ends of the bar is subjected to an incrementally increasing displacement, while a roller boundary condition is applied at the other end. A roller boundary condition is also applied to the bottom edge of the bar where symmetry is assumed.

To avoid unwanted mesh dependency of the solution during cracking, the damage model needs to be regularized to ensure that a consistent amount of energy is dissipated during mesh refinement or for different discretization orders. Two techniques are available in COMSOL Multiphysics, and these are exemplified 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 damage model.

A mapped mesh with 101x1 elements is used, either with linear or quadratic shape functions for the displacement field. To force a strain localization, a weakness is introduced in the middle of the bar by reducing the tensile strength by 2.5%.

Results and Discussion

Four studies are set up to model different combinations of discretization orders and regularization methods:

Quadratic displacement field with the crack band method

Quadratic displacement field with the implicit gradient method

Linear displacement field with the crack band method

Linear displacement field with the implicit gradient method

Contour plots of the damage variable are shown in Figure 1 for the simulations using the crack band method and in Figure 2 for those using the implicit gradient method. The difference between the two regularization methods is clearly visible by comparing the two figures. For the crack band method, damage is only nonzero in a single element. On the contrary, for the implicit gradient method the damage is distributed over several elements. The effect of this is also seen in Figure 3, where the curves of force versus average strain from each study are compared. For studies using the crack band method (Study 1 and 3), the force-deformation curve has the same shape as the exponential strain softening prescribed by the damage model. On the other hand, for studies 2 and 4 where the implicit gradient method is used, the shape of the force-deformation curves differs, which is a consequence of the evolution of the damage zone during the incremental stretching of the bar. There are also differences between the quadratic and linear solutions; these discrepancies are due to how well the respective discretization can resolve strains in the damage zone.

Figure 1: Damage distribution using crack-band regularization.

Figure 2: Damage distribution using implicit regularization.

Figure 3: Force versus prescribed end displacement.

Figure 4 shows the stress versus strain and the evolution of the damage variable in the centroid of the middle element. Even though the model is globally supplied with the same fracture energy, the different regularization methods result in significantly different stress versus strain curves. The same conclusion can be drawn from studying the damage evolution curves in Figure 5, where the damage grows significantly faster for the models using the crack band model.

Figure 4: Stress versus strain at the center of the damaged region.

Figure 5: Evolution of the damage variable at the center of the damaged region.

To study the strain localization phenomenon in more detail, the equivalent strain and the damage variable are plotted along the bar in Figure 6 for studies 1 and 2 (quadratic displacement field), and in Figure 7 for studies 3 and 4 (linear displacement field). When the crack band method is used, all damage and deformation is concentrated into a single element. Note, however, that especially for the quadratic displacement field, some unwanted damage also appeared in the adjacent elements. To minimize this unwanted effect, it is in general recommended to use linear interpolation for the displacements when using the crack band method. When comparing this localization behavior to that of the implicit gradient method, the latter exhibits a distributed damage zone. However, it can also be noticed that the strain field is much more narrow. This is in fact also visible in Figure 2 where the middle elements are much more stretched than those toward the edge of the damaged zone. When using the implicit gradient method, the distribution of damage and the localization of strains are clearly much less mesh dependent than when using the crack band method. This can be further investigated by doing a mesh refinement study.

Finally, an important consequence of the different methods are the properties of the displacement field. In Figure 8, the horizontal displacement component is plotted along the middle part of the bar for all four studies. It can be noticed that for both studies 1 and 3 (crack band method), most of the horizontal displacements are localized in the middle element and its derivative is not well defined (see also Figure 6 and Figure 7). However, for studies 2 and 4, where the implicit gradient method is used, both the displacements and its derivatives remain continuous and well defined across the crack.

Figure 6: Distribution along the bar of equivalent strain and damage when using quadratic shape functions.

Figure 7: Distribution along the bar of equivalent strain and damage when using linear shape functions.

Figure 8: Displacement along the bar

Notes About the COMSOL implementation

Because the material softens, it would not be possible to perform this analysis with a prescribed load. Using a prescribed displacement, however, solves the problem.

In this example, the material is made slightly weaker at the center of the bar. It can be extremely difficult to make a model like this to converge without this trick. In this particular example, the damage had appeared in the entire bar based on numerical roundoff, and the solution would jump back and forth between iterations. Homogeneous stress states are notoriously difficult to handle with softening material models. Fortunately, this is seldom the case in real life structures. One approach, with some physical interpretation is to modify the material data with a function with a random spatial distribution.

Geomechanics_Module/Damage/damage_test_bar

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
L	0.1[m]	0.1 m	Bar length
H	L/50	0.002 m	Bar thickness
max_average_strain	5e-4	5E-4	Maximum average strain
average_strain	0	0	Average strain (loading parameter)
lscale	0.001[m]	0.001 m	Length scale

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type L.

In the text field, type H/2.

Click

Add an integration parameter to calculate the total reaction force in postprocessing.

Definitions

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Select Boundary 4 only.

Locate the section. From the list, choose .

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

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.

Linear Elastic Material 1

In the window, click .

Damage: Implicit Gradient

In the toolbar, 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.

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