Isotropic Phase Field Softening Material

Constitutive Law

This MPM softening material models fracture using a phase field. In phase field fracture modeling, sharp crack surfaces are replaced by diffuse phase fields. The growth of these diffuse cracks is controlled by energy flow through the material and the material's fracture toughness. This material implements phase field fracture methods using methods described in Miehe and extends it in a few areas.

Phase Field Methods

In phase field fracture modeling of small-strain elastic materials, the total strain energy is partitioned into damaging, [math]\displaystyle{ \Psi^{(+)} }[/math], and nondamaging, [math]\displaystyle{ \Psi^{(-)} }[/math], terms by

     [math]\displaystyle{ \Psi = g(\phi)\Psi^{(+)} + \Psi^{(-)} }[/math]

where [math]\displaystyle{ g(\phi) }[/math] is a damage or softening law that depends on damage parameter [math]\displaystyle{ \phi }[/math] given by a phase value tracked on each material point that varies from 0 for undamaged to 1 for complete damage.

Using a viscous regularization (see Miehe ), the phase field evolves by the equation

     [math]\displaystyle{ \eta \frac{d\phi}{dt} = G_c \ell \nabla^2\phi - \frac{G_c}{\ell}\phi - g'(\phi)\mathcal{H} }[/math]

where [math]\displaystyle{ G_c }[/math] is toughness, [math]\displaystyle{ \ell }[/math] is a length scale (describing width of diffuse cracks), and [math]\displaystyle{ \mathcal{H} }[/math] is history variable given by [math]\displaystyle{ \max(\Psi^{(+)}) }[/math]. This equation is a diffusion equation that describes evolution of damage driven by damaging energy. The first two terms lead to diffusive cracks that decay exponentially from complete damage state with exponential decay distance determined by [math]\displaystyle{ \ell }[/math]. The last term is a source term to causes crack propagation. In other words, crack propagation is affected by the methods used to partition energy into [math]\displaystyle{ \Psi^{(+)} }[/math] and [math]\displaystyle{ \Psi^{(-)} }[/math].

Simulations using viscous regularization must couple solution of phase field diffusion equation to the material's mechanics calculations. For that to work, all simulations using this material must include a Diffusion command and choose (style) of "fracture". If needed, you can set grid based phase value boundary conditions using the "fracture" option to define damage conditions. Although you can set damage flux on particle surfaces, such conditions may have few applications in real-world modeling. Note that MPM solutions of transport equations can develop oscillations that prevent output of viable phase field descriptions of cracks. This problem can be solved by using FMPM(k) for transport tasks with the PeriodicXPIC Custom Task

Energy Partitioning

Crack evolution results depend on the method chosen to partition energy into [math]\displaystyle{ \Psi^{(+)} }[/math] and [math]\displaystyle{ \Psi^{(-)} }[/math]. A common method for isotropic materials is to use an eigenstrain analysis and partition energy such that [math]\displaystyle{ \Psi^{(+)} }[/math] depends only on tensile principle strains and tensile trace of the strain tensor. The physical interpretation is that only tensile loading causes fracture while compression does not. A weakness of this approach is that shear loading may not respond as expected.

An alternative method is to partition energy into pressure and deviatoric strains. Now [math]\displaystyle{ \Psi^{(+)} }[/math] depnds on hydrostatic tension (negative pressure) and deviatoric strain while [math]\displaystyle{ \Psi^{(-)} }[/math] depends only on positive pressure and does not promote fracture.

This material supports both these options that can be selected with the partition material property.

Phase Field Damage Law

Crack evolution also depend on the damage (or softening) law used for the phase field material. The vast majority of phase field fracture papers set [math]\displaystyle{ g(\phi) = (1-\phi)^2 }[/math] under the misunderstood concept that [math]\displaystyle{ g'(1) }[/math] needs to be zero to stop dissipating energy or to prevent [math]\displaystyle{ \phi }[/math] from exceeding 1. In dynamic codes, however, any [math]\displaystyle{ g(\phi) }[/math] can be used provided [math]\displaystyle{ g(\phi)=1 }[/math], [math]\displaystyle{ g(1)=0 }[/math], and [math]\displaystyle{ \phi }[/math] is prevented from ever decreasing or ever exceeding 1.

Because [math]\displaystyle{ g(\phi) }[/math] appears in total energy, it implicitly corresponds to stress-strain behavior of the material as a function of damage. The common [math]\displaystyle{ g(\phi) = (1-\phi)^2 }[/math] implies all materials have the same stress-strain law. Reverse engineering show that law further implies that all material are ductile. If phase field materials are expected to apply to a range of materials, the modeling should admit other laws including laws with varying ductility. The available [math]\displaystyle{ g(\phi) }[/math] are selected using the gd material property and using the garg material property to pick an argument [math]\displaystyle{ \chi }[/math]. The current damage laws (by gd setting) are:

0. Quadratic Law: this law generalizes the common quadratic law with an argument:
        [math]\displaystyle{ g(\phi) = (1-\phi)(1-(2\chi-1)\phi) }[/math]
   where [math]\displaystyle{ 0\le\chi\le1 }[/math] varies the material's ductility (higher is more ductile). Choosing [math]\displaystyle{ \chi=1 }[/math] is the standard (and very ductile) law used in the literature. When partitioning energy using an eigenstrain analysis, this law is unstable for [math]\displaystyle{ \chi }[/math] less than about 0.5 (which is when [math]\displaystyle{ g(\phi) }[/math] becomes linear).
1. Exponential Law: a trunctated exponential function to go from 1 to 0 as [math]\displaystyle{ \phi }[/math] goes from 0 to 1:
        [math]\displaystyle{ g(\phi) = \frac{e^{-\chi\phi} - e^{-\chi}}{1-e^{-\chi}} }[/math]
   where [math]\displaystyle{ \chi\gt0 }[/math] increases the ductility. In the limit as [math]\displaystyle{ \chi\to0 }[/math], this law is linear in [math]\displaystyle{ \phi }[/math] but cannot be modeled using this law. To use a linear law, pick quadratic law with [math]\displaystyle{ \chi=0.5 }[/math].
2. Linear Softening Law: if the energy partitioning is pressure and deviatoric stress, one can reverse engineer a linear softening law where stress is linear elastic until failure and then decays linearly to zero as damage increases:
        [math]\displaystyle{ g(\phi) = 1 + \frac{\chi\phi^2}{2} - \phi\sqrt{1+\chi + \frac{\chi^2\phi^2}{4}} }[/math]
   where [math]\displaystyle{ \chi\gt0 }[/math] increases the ductility. The peak stress to initiate damage (if loaded in uniaxial tension) is
        [math]\displaystyle{ \sigma_i = \sqrt{\frac{G_c E}{\ell \chi}} }[/math]
   The area under a linear softening law can be divided into elastic area up to the peak ([math]\displaystyle{ A_1 }[/math]) and area under the post-peak decay to zero stress ([math]\displaystyle{ A_2 }[/math]). For phase field material using this law
        [math]\displaystyle{ A_2 = \chi A_1 = \frac{G_c}{2\ell} }[/math]
   Note that [math]\displaystyle{ \chi\to0 }[/math] corresponds to [math]\displaystyle{ A_2\to0 }[/math] for an elastic-brittle material. This behavior, however, cannot be modeling with this law because the initiation stress becomes infinite. Although this linear softening response was determined by assuming pressure and deviatoric partitioning of energy, this material can use this law even when partitioning by eigenstrain analysis. Simulations show this option provides approximate linear softening response.

Stability Factor

Some phase field literature suggest that letting [math]\displaystyle{ g(\phi) }[/math] reach zero when damage reaches one can cause instabilities. One option to avoid such instability is to revise and damage law to be

     [math]\displaystyle{ g^*(\phi) = (1-k)g(\phi) + k }[/math]

such that [math]\displaystyle{ g^*(1)=k }[/math] where [math]\displaystyle{ k\lt \lt 1 }[/math] is set using the stability material property. Stability factors in dynamic MPM do not seem necessary, but an option is provide if needed, or for comparison to other modeling.

Time-Independent Phase Field Modeling

The viscous regularization methods is a time-dependent solution to the time-independent phase field equation:

     [math]\displaystyle{ 0 = G_c \ell \nabla^2\phi - \frac{G_c}{\ell}\phi - g'(\phi)\mathcal{H} }[/math]

The diffusion equation approaches this solution in steady state (i.e. when [math]\displaystyle{ d\phi/dt\to 0 }[/math]) or in the low viscosity limit (i.e. when [math]\displaystyle{ \eta\to 0 }[/math]). The problem with the time-independent equation is that it does not match well with dynamic MPM modeling. It is possible (and some MPM papers do it) to solve this equation on the background grid using standard linear finite element methods. This approach has two problems. First, one can only use linear finite element analysis when [math]\displaystyle{ g'(\phi) }[/math] is linear, which implies it can only be used for quadratic softening laws (this limitation might be another reason that most of the literature uses quadratic softening laws). Second, this approach in MPM is much slower (more than 30X slower).

Although OSParticulas has a option to solve the time-independent equation, that option is not provided here. Instead, one must use viscous regularization and that approach adds a burden of confirming choice of viscosity has not affected the results. In well behaved problems, results will be independent of viscosity provided viscosity is sufficiently small. If viscosity gets too high, it inhibits crack diffusion, slows down crack propagation, and gives non-physical results. A good rule of thumb is to pick viscosity such that it does not cause the time step needed for diffusion modeling to be smaller than time step needed to solve mechanics equations. A good starting point is to pick:

     [math]\displaystyle{ \eta \sim \frac{2G_c\ell}{v_{max}\Delta x} }[/math]

where [math]\displaystyle{ v_{max} }[/math] is wave speed of the elastic material and [math]\displaystyle{ \Delta x }[/math] cell size in the background grid. Simulations with different viscosites should also be run to determine if this value is small enough or if it needs to be smaller. Reference gives more details on choosing viscosity.

Material Properties

This isotropic phase field fracture modeling using a single energy release rate for toughness that scales evolution of damage. The critical energy release rate is enter using the base material JIc property. The other needed material properties are as follows:

Property	Description	Units	Default
(Isotropic Properties)	Enter all properties needed to define the underlying isotropic material response	varies	varies
ell	Length scale parameter used in variational fracture mechanics. To resolve diffusive cracks, ell should be at least four times material point diameter.	length units	none
viscosity	Viscosity to use when solving coupled phase field evolution in a diffusion tasks. See above for comments on choosing viscosity.	viscosity units	none
partition	Chose the method used to partition energy into energy that causes fracture and energy that does not cause fracture. The options are 0 = use eigenstrain analysis and 1 = divide into pressure and deviatoric strains	none	1
gd	Softening law with options 0 = quadratic, 1 = exponential, 2 = linear softening	none	0
garg	An optional argument for use within the softening law. If not provided, default values depend on gd and are 1, 3, and 4, for gd = 0, 1, or 2, respectively	none	varies
stability	A stability factor thought to stabilize post-failure response	none	0
(other)	Properties common to all materials and must include JIc to provide toughness [math]\displaystyle{ G_c }[/math]. Notice that this material cannot model mixed-mode behavior where results might depend on separate mode I and mode II toughnesses.	varies	varies

The results in Miehe correspond to partition = 0, gd = 0, and garg = 1. These choices give poor results in some problems, especial when shear is important. Other choices are encouraged when these basic ones look questionable.

If needed, you can set set initial phase field value for particles using particle boundary conditions.

History Variables

This material stores several history variables that track the extent of the damage and evolution of the phase field:

1. Maximum energy history term that provides source terms for phase field evolution: [math]\displaystyle{ \mathcal{H}\max(\Psi^{(+)}) }[/math].
2. Particle failure state equal to 0 if not failed and 1 if failed (i.e., phase value has reached 1).
3. Current phase field value
4. Change in phase field since the last time step. It is used in constitutive law modeling and is scaled by 0.5 when using USAVG method.