PDF

Phase-Field Modeling of Dynamic Crack Branching
Introduction
This example considers a classical benchmark for phase-field modeling of dynamic fracture, see for example Ref. 1. An instantaneous tensile load is applied to a planar specimen with a pre-existing crack. Initially, the crack propagates perpendicular to the loading direction, after which the crack branches off symmetrically until catastrophic failure occurs. The problem is solved using the AT1 phase-field model.
Model Definition
The AT1 phase-field model is described in terms of the crack surface density function
(1)
where w(ϕ) = ϕ, cw = 8/3 is a normalization constant, and lint is an internal length scale for the phase-field regularization. The fracture energy density is then written as
(2)
where Gc is the critical energy release rate.
The governing equations for the phase-field damage model include the standard momentum balance
(3)
together with the phase-field evolution equation
(4)
Herein, g(ϕ) is the damage degradation function and is the crack driving force, which is based on a tension–compression split of the strain energy density
(5)
For the AT1 model, Equation 4 does not allow a zero solution for the phase field at zero crack driving force H, which implies that the phase field is unbounded.
A strain energy threshold can be introduced in order to ensure that ϕ = 0 is a valid solution, so that the phase field is bounded. This threshold can be derived by inserting the zero solution into Equation 4, and solving for the critical value Wc0 of the crack driving force H.
For the quadratic degradation function g(ϕ) = (1 − ϕ)2, this results into the condition
(6)
from which the critical strain can be identified.
The modified crack driving force is then written as
(7)
This also implies that crack propagation only starts once this threshold is exceeded.
The model geometry and boundary conditions are depicted in Figure 1. Symmetry about the X-axis is used in order to reduce the size of the computational domain.
Plane strain conditions are assumed, and the thickness is arbitrarily set to 1 m.
The pre-existing crack is modeled by prescribing the phase field to a value ϕ = 1 along the lower-left boundary, while being free to deform.
The instantaneous load is applied over the first time step using a smooth step function.
Figure 1: Geometry and boundary conditions for the notched planar tension sample.
The material properties used in this example are inspired by Ref. 1 and listed in Table 1.
Results and Discussion
The damaged stress is shown in Figure 2 at selected time points. The parts of the domain with damage d(ϕ) = 1 − g(ϕ) > 0.95 have been removed from the plot in order to visualize the crack path.
Crack propagation is initiated after approximately 10 μs. First, the crack propagates perpendicular to the applied load, after which the crack branches off at around 33 μs and spreads quickly throughout the sample, as shown in the snapshots at 45 μs and 75 μs.
Figure 2: Damaged stress and crack path at selected times.
The phase field at the end of the simulation is shown in Figure 3. The upper crack branch is overlayed with circular points indicating the location of the crack front as a function of time.
Figure 3: Phase field at the end of the simulation overlayed with colored points indicating the location of the crack front as a function of time.
The total elastic and fracture energy in the sample are plotted as functions of time in Figure 4, and the temporal evolution of the crack length and the crack velocity are shown in Figure 5. The results compare well with those reported in Ref. 1 for the AT2 phase-field damage model.
Figure 4: Total elastic and fracture energy as functions of time.
Figure 5: Crack length and crack tip velocity as functions of time. The crack tip velocity stays below 0.6 times the Rayleigh wave speed vR, as noted in Ref. 1.
Notes About the COMSOL Implementation
To implement the AT1 phase-field damage model, the Solid Mechanics and Phase Field in Solids interfaces are coupled using the Phase-Field Damage multiphysics coupling node.
The crack front is tracked in postprocessing by the use of a heuristic indicator function that is the product between the phase field and its time derivative, two quantities that should be maximal at the crack tip during crack propagation. The crack trajectory is then evaluated using a model method. Note that the method editor is only available in the Windows® version of the COMSOL Desktop.
Reference
1. M.J. Borden, C.V. Verhoosel, M.A. Scott, T.J.R. Hughes, and C.M. Landis, A phase-field description of dynamic brittle fracture, ICES REPORT 11-14, The Institute for Computational Engineering and Sciences, The University of Texas at Austin, 2011.
Application Library path: Nonlinear_Structural_Materials_Module/Damage/dynamic_crack_branching
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D.
2
In the Select Physics tree, select Structural Mechanics > Phase-Field Damage.
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies > Time Dependent.
6
Global Definitions
Parameters 1
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Step 1 (step1)
Define a smooth Step function that will be used to apply the instantaneous tensile load.
1
In the Home toolbar, click  Functions and choose Global > Step.
2
In the Settings window for Step, locate the Parameters section.
3
In the Location text field, type 0[s].
4
Click to expand the Smoothing section. In the Size of transition zone text field, type dtmax.
5
From the Location definition list, choose Beginning of step.
Geometry 1
1
In the Model Builder window, under Component 1 (comp1) click Geometry 1.
2
In the Settings window for Geometry, locate the Units section.
3
From the Length unit list, choose mm.
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type width.
4
In the Height text field, type height/2.
5
Click to expand the Layers section. In the table, enter the following settings:
6
Select the Layers to the left checkbox.
7
Clear the Layers on bottom checkbox.
8
Click  Build All Objects.
Materials
Material 1 (mat1)
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material.
2
In the Settings window for Material, locate the Material Contents section.
3
Solid Mechanics (solid)
1
In the Model Builder window, under Component 1 (comp1) click Solid Mechanics (solid).
2
In the Settings window for Solid Mechanics, locate the Thickness section.
3
In the d text field, type d0.
4
Locate the Structural Transient Behavior section. From the list, choose Include inertial terms.
5
Click to expand the Discretization section. From the Displacement field list, choose Linear.
Symmetry 1
1
In the Physics toolbar, click  Boundaries and choose Symmetry.
2
Boundary Load 1
1
In the Physics toolbar, click  Boundaries and choose Boundary Load.
2
3
In the Settings window for Boundary Load, locate the Force section.
4
Specify the fA vector as
Phase Field in Solids (pfs)
Phase Field Model, AT1
1
In the Model Builder window, under Component 1 (comp1) > Phase Field in Solids (pfs) click Phase Field Model 1.
2
In the Settings window for Phase Field Model, type Phase Field Model, AT1 in the Label text field.
3
Locate the Potential section. From the Potential function list, choose User defined.
4
In the Q text field, type 3*phi/8.
5
Locate the Phase Field Model section. In the lint text field, type lint.
6
In the D text field, type 3/4.
Symmetry 1
1
In the Physics toolbar, click  Boundaries and choose Symmetry.
2
Prescribed Phase Field 1
Model the pre-existing crack by prescribing the phase field to 1, corresponding to a fully damaged material.
1
In the Physics toolbar, click  Boundaries and choose Prescribed Phase Field.
2
3
In the Settings window for Prescribed Phase Field, locate the Prescribed Phase Field section.
4
In the ϕ0 text field, type 1.
Definitions
To set up the AT1 phase-field damage model, start by defining a few help variables as described in the model documentation.
Variables 1
1
In the Model Builder window, under Component 1 (comp1) right-click Definitions and choose Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Integration 1 (intop1)
Define the integration operator intop1 on domain 2 only in order to exclude the pre-existing crack from the evaluation of the total fracture energy. Note that the factor 2 is due to symmetry about the X-axis.
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
In the Settings window for Integration, locate the Source Selection section.
3
From the Selection list, choose Manual.
4
5
Locate the Advanced section. In the Integration order text field, type 2.
6
From the Frame list, choose Material  (X, Y, Z).
Multiphysics
Set up the Phase-Field Damage multiphysics coupling. For the AT1 model, we need to enforce that crack propagation only occurs above the critical strain energy density Wc0. We can do this using the max-operator, the built-in variable pfdmg1.Ws0 for the positive part of the strain energy density, and the help variable Wc0.
Phase-Field Damage 1 (pfdmg1)
1
In the Model Builder window, under Component 1 (comp1) > Multiphysics click Phase-Field Damage 1 (pfdmg1).
2
In the Settings window for Phase-Field Damage, locate the Damage Model section.
3
From the Parent material model list, choose Linear Elastic Material 1.
4
From the Dd list, choose User defined.
5
In the text field, type max(pfdmg1.Ws0,Wc0)*pfdmg1.lint/Gc0.
6
From the Exclude compressive energy list, choose Volumetric only.
7
Click to expand the Advanced section. In the dmax text field, type 1-eta.
Mesh 1
Mapped 1
In the Mesh toolbar, click  Mapped.
Size 1
1
Right-click Mapped 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
5
Locate the Element Size section. From the Predefined list, choose Extremely fine.
6
Click the Custom button.
7
Locate the Element Size Parameters section.
8
Select the Maximum element size checkbox. In the associated text field, type he.
Size 2
1
In the Model Builder window, right-click Mapped 1 and choose Size.
2
In the Settings window for Size, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
5
Locate the Element Size section. From the Predefined list, choose Extremely fine.
6
Click the Custom button.
7
Locate the Element Size Parameters section.
8
Select the Maximum element size checkbox. In the associated text field, type 10*he.
9
In the Model Builder window, right-click Mesh 1 and choose Build All.
Study 1
The crack needs about 75 µs to propagate across the sample. To resolve the dynamics, use the implicit generalized-α method with manual control over the time step. The parameter dtmax is set to be approximately equal to the time required for a Rayleigh wave to propagate across one element. The crack speed is typically slower than the Rayleigh wave.
Step 1: Time Dependent
1
In the Model Builder window, under Study 1 click Step 1: Time Dependent.
2
In the Settings window for Time Dependent, locate the Study Settings section.
3
From the Time unit list, choose µs.
4
In the Output times text field, type range(0,1,75).
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node.
3
In the Model Builder window, expand the Study 1 > Solver Configurations > Solution 1 (sol1) > Dependent Variables 1 node.
4
In the Model Builder window, under Study 1 > Solver Configurations > Solution 1 (sol1) click Time-Dependent Solver 1.
5
In the Settings window for Time-Dependent Solver, click to expand the Time Stepping section.
6
From the Method list, choose Generalized alpha.
7
From the Steps taken by solver list, choose Manual.
8
Clear the Interpolate solution at end time checkbox.
9
In the Time step text field, type dtmax.
10
In the Amplification for high frequency text field, type 0.5.
Use a tolerance-based convergence criterion on the outer segregated solver.
11
In the Model Builder window, under Study 1 > Solver Configurations > Solution 1 (sol1) > Time-Dependent Solver 1 click Segregated 1.
12
In the Settings window for Segregated, locate the General section.
13
From the Termination technique list, choose Tolerance.
Before solving, compute the initial values in order to set up a plot while solving for monitoring the solution process.
14
In the Study toolbar, click  Get Initial Value.
Results
Use a Mirror dataset to visualize the complete sample.
Mirror 2D 1
1
In the Results toolbar, click  More Datasets and choose Mirror 2D.
2
In the Settings window for Mirror 2D, locate the Axis Data section.
3
In row Point 2, set X to 1.
4
In row Point 2, set Y to 0.
5
Stress (solid)
1
In the Model Builder window, under Results click Stress (solid).
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Dataset list, choose Mirror 2D 1.
4
Locate the Plot Settings section. Clear the Plot dataset edges checkbox.
Surface 1
1
In the Model Builder window, expand the Stress (solid) node, then click Surface 1.
2
In the Settings window for Surface, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Solid Mechanics > Damage > solid.misesdGp - von Mises stress, damaged - N/m².
3
Locate the Expression section. From the Unit list, choose MPa.
Deformation
1
In the Model Builder window, expand the Surface 1 node, then click Deformation.
2
In the Settings window for Deformation, locate the Scale section.
3
Select the Scale factor checkbox. In the associated text field, type 10.
Surface 1
Use a Filter to remove the damaged elements from the plot to visualize the crack trajectory.
Filter 1
1
In the Model Builder window, right-click Surface 1 and choose Filter.
2
In the Settings window for Filter, locate the Element Selection section.
3
In the Logical expression for inclusion text field, type solid.dmg<0.95.
Arrow Line 1
1
In the Model Builder window, right-click Stress (solid) and choose Arrow Line.
2
In the Settings window for Arrow Line, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Solid Mechanics > Load > solid.fax,solid.fay - Force per deformed area (spatial frame).
3
Click to expand the Inherit Style section. From the Plot list, choose Surface 1.
4
Clear the Arrow scale factor checkbox.
5
Clear the Color checkbox.
6
Clear the Color and data range checkbox.
Deformation 1
Right-click Arrow Line 1 and choose Deformation.
Line 1
1
In the Model Builder window, right-click Stress (solid) and choose Line.
2
In the Settings window for Line, locate the Expression section.
3
In the Expression text field, type 1.
4
Click to expand the Title section. From the Title type list, choose None.
5
Locate the Coloring and Style section. From the Coloring list, choose Uniform.
6
From the Color list, choose Black.
7
Click to expand the Inherit Style section. From the Plot list, choose Surface 1.
8
Clear the Color checkbox.
9
Clear the Color and data range checkbox.
10
Clear the Height scale factor checkbox.
11
Clear the Tube radius scale factor checkbox.
Deformation 1
Right-click Line 1 and choose Deformation.
Selection 1
1
In the Model Builder window, right-click Line 1 and choose Selection.
2
Filter 1
In the Model Builder window, under Results > Stress (solid) > Surface 1 right-click Filter 1 and choose Copy.
Filter 1
In the Model Builder window, right-click Line 1 and choose Paste Filter.
Study 1
Step 1: Time Dependent
1
In the Settings window for Time Dependent, click to expand the Results While Solving section.
2
Select the Plot checkbox.
3
From the Update at list, choose Time steps taken by solver.
4
In the Study toolbar, click  Compute.
Results
Plot the total elastic strain energy and the regularized fracture energy.
Total Elastic and Fracture Energy
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
2
In the Settings window for 1D Plot Group, type Total Elastic and Fracture Energy in the Label text field.
Global 1
1
Right-click Total Elastic and Fracture Energy and choose Global.
2
In the Settings window for Global, locate the y-Axis Data section.
3
4
Click to expand the Coloring and Style section. Find the Line markers subsection. From the Marker list, choose Cycle.
Total Elastic and Fracture Energy
1
In the Model Builder window, click Total Elastic and Fracture Energy.
2
In the Settings window for 1D Plot Group, click to expand the Title section.
3
From the Title type list, choose Label.
4
Locate the Plot Settings section.
5
Select the y-axis label checkbox. In the associated text field, type Energy (J).
6
Locate the Legend section. From the Position list, choose Upper left.
Crack Front Indicator
Next, evaluate the trajectory, velocity, and length of the crack. To start with, we first need to track the crack front as a function of time. Here, we will identify the crack front using a suitable indicator function. A simple condition that can be used in this case is to look for the maximum of the product between the phase field and its time derivative.
1
In the Results toolbar, click  Evaluation Group.
2
In the Settings window for Evaluation Group, type Crack Front Indicator in the Label text field.
Surface Maximum 1
1
Right-click Crack Front Indicator and choose Maximum > Surface Maximum.
2
3
In the Settings window for Surface Maximum, locate the Expressions section.
4
5
Click to expand the Configuration section. From the Point type list, choose Mesh vertices.
6
Select the Include position checkbox.
7
In the Crack Front Indicator toolbar, click  Evaluate.
Phase Field (pfs)
Visualize the coordinates of the crack front on top of the phase field to verify that the moving front is captured correctly.
1
In the Model Builder window, under Results click Phase Field (pfs).
2
In the Settings window for 2D Plot Group, locate the Data section.
3
From the Dataset list, choose Mirror 2D 1.
4
From the Time (µs) list, choose Last (75).
Table Point 1
1
In the Phase Field (pfs) toolbar, click  More Plots and choose Table Point.
2
In the Settings window for Table Point, locate the Data section.
3
From the Source list, choose Evaluation group.
4
From the x-axis column list, choose x (mm).
5
From the y-axis column list, choose y (mm).
6
From the Data column list, choose Time (µs).
7
In the Phase Field (pfs) toolbar, click  Plot.
Application Builder
With the coordinates of the crack front at hand, the crack length and velocity can be evaluated using a model method. Note that the method editor is only available in the Windows® version of the COMSOL Desktop.
In the Home toolbar, click  Application Builder.
New Method
1
In the Application Builder window, right-click Methods and choose New Method.
2
In the New Method dialog, type computeCrackTrajectory in the Name text field.
3
computeCrackTrajectory
1
In the Application Builder window, under Methods click computeCrackTrajectory.
2
Copy the following code into the computeCrackTrajectory window:
// Extract crack length and crack velocity
EvaluationGroupFeature eg = model.result().evaluationGroup("eg1");
eg.run();
double[][] data = eg.getReal(); // {time, indicator function, x-coordinate, y-coordinate}
int n = data.length;

// Initialize a 2d array for output data
double[][] out = new double[n][5]; // {time, x-coordinate, y-coordinate, crack length, crack velocity}
double xprev = 50.0; // Initial crack front x-coordinate
double yprev = 0.0; // Initial crack front y-coordinate

// Make sure the crack front data are initialized correctly
out[0] = new double[]{data[0][0], xprev, yprev, 0.0, 0.0};
data[0][2] = xprev;
data[0][3] = yprev;

// Compute crack length and velocity
double clength = 0.0;

for (int i = 1; i < n; i++) {
  double x = data[i][2];
  double y = data[i][3];
  if (data[i][1] < 0) {
    // if the indicator function is negative, no crack propagation has been detected
    x = xprev;
    y = yprev;
  }
  // Crack length
  double dx = x-xprev; // x-increment
  double dy = y-yprev; // y-increment
  clength += Math.sqrt(Math.pow(dx, 2)+Math.pow(dy, 2));
  
  // Crack velocity
  double dt = data[i][0]-data[i-1][0]; // Time step
  double vx = dx/dt; // x-velocity
  double vy = dy/dt; // y-velocity
  double cvel = Math.sqrt(Math.pow(vx, 2)+Math.pow(vy, 2));
  
  // Update output data array and previous coordinates
  out[i] = new double[]{data[i][0], x, y, clength, cvel};
  xprev = x;
  yprev = y;
}

String tblTag = "tblcrk";
TableFeature tbl;
if (contains(model.result().table().tags(), tblTag)) {
  tbl = model.result().table(tblTag);
  tbl.clearTableData();
} else {
  tbl = model.result().table().create(tblTag, "Table");
}
tbl.setColumnHeaders(new String[]{"Time (us)", "x-coordinate (mm)", "y-coordinate (mm)", "Crack length (mm)", "Crack velocity (mm/us)"});
tbl.setTableData(out);
tbl.label("Crack Trajectory");
Methods
1
In the Home toolbar, click  Model Builder to switch to the main desktop.
Add a Method Call to computeCrackTrajectory in order to run it.
2
In the Developer toolbar, click  Method Call and choose computeCrackTrajectory.
Global Definitions
Compute Crack Trajectory
1
In the Model Builder window, under Global Definitions click ComputeCrackTrajectory 1.
2
In the Settings window for Method Call, type Compute Crack Trajectory in the Label text field.
3
Results
Crack Propagation
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Crack Propagation in the Label text field.
Crack Length
1
Right-click Crack Propagation and choose Table Graph.
2
In the Settings window for Table Graph, type Crack Length in the Label text field.
3
Locate the Data section. From the Plot columns list, choose Manual.
4
In the Columns list, select Crack length (mm).
5
Locate the Coloring and Style section. Find the Line markers subsection. From the Marker list, choose Cycle.
6
Click to expand the Legends section. Select the Show legends checkbox.
7
Find the Include subsection. Select the Label checkbox.
8
Clear the Headers checkbox.
Crack Velocity
1
Right-click Crack Length and choose Duplicate.
2
In the Settings window for Table Graph, type Crack Velocity in the Label text field.
3
Locate the Data section. In the Columns list, select Crack velocity (mm/us).
Crack Propagation
Compare the crack velocity with the Rayleigh wave speed.
Right-click Crack Velocity and choose Global.
Rayleigh Wave Speed
1
In the Settings window for Global, type Rayleigh Wave Speed in the Label text field.
2
Locate the y-Axis Data section. In the table, enter the following settings:
3
Click to expand the Coloring and Style section. From the Color list, choose From theme.
4
Click to expand the Legends section. Clear the Show legends checkbox.
Crack Propagation
In the Crack Propagation toolbar, click  More Plots and choose Table Annotation.
Table Annotation 1
1
In the Settings window for Table Annotation, locate the Data section.
2
From the Source list, choose Local table.
3
4
Select the LaTeX markup checkbox.
5
Locate the Coloring and Style section. Clear the Show point checkbox.
Crack Propagation
1
In the Model Builder window, click Crack Propagation.
2
In the Settings window for 1D Plot Group, locate the Plot Settings section.
3
Select the x-axis label checkbox. In the associated text field, type Time (µs).
4
Select the Two y-axes checkbox.
5
In the table, select the Plot on secondary y-axis checkboxes for Crack Velocity, Rayleigh Wave Speed, and Table Annotation 1.
6
Locate the Legend section. From the Layout list, choose Outside graph axis area.
7
In the Model Builder window, click Crack Propagation.
8
Click to expand the Title section. From the Title type list, choose Label.
9
In the Crack Propagation toolbar, click  Plot.